SARAS: A general-purpose PDE solver for fluid dynamics

1 Department of Mechanical Engineering, Indian Institute of Technology Kanpur, Uttar Pradesh 208016, India 2 Department of Physics, Indian Institute of Technology Kanpur, Uttar Pradesh 208016, India 3 Mechanical Engineering, Physical Science and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955, Saudi Arabia 4 Department of Mechanical Engineering, ZHCET, Aligarh Muslim University, Uttar Pradesh 202002, India DOI: 10.21105/joss.02095


Statement of Need
There are a number of open-source solvers for Computational Fluid Dynamics. Some well known solvers include OpenFOAM (Weller et al., 1998), Pencil Code (Brandenburg & Dobler, 2002), Gerris (Popinet, 2003), SU2 (Palacios et al., 2013), etc. to name a few. However, many of them, for instance, SU2, OpenFOAM, etc. use finite-volume method with structured or unstructured meshes. On the other hand, the use of finite-difference method restricts the solver to simpler domains, but offers comparatively higher accuracy with lesser computational effort. In this regard, the Pencil Code also uses finite-difference method, but is used to solve compressible flows, unlike SARAS which is designed for incompressible flows. Moreover, SARAS is written in C++ with an object oriented structure like OpenFOAM.
The design of SARAS is inspired by TARANG (Verma et al., 2013), a pseudo-spectral solver developed in our lab. TARANG has been shown to scale up to 196608 cores (Chatterjee et al., 2018), and SARAS has been designed with the goal of achieving similar scaling performance. We conducted a preliminary scaling analysis of SARAS using up to 1024 cores for the lid-driven cavity (LDC) problem on a 512 3 grid, and observed strong scaling (Verma et al., 2020). However, we need to implement further optimizations to the code before scaling to larger grids and more cores. In solving the Rayleigh Benard Convection problem, we also observed that the computational efficiency of SARAS is at par with that of OpenFOAM. Recently, we also performed a comparative study between spectral and finite difference codes in the context of exascale computing (Verma et al., 2020), using SARAS.
In SARAS, the underlying mathematical constructs like vector and scalar fields are defined as classes. Moreover, vector calculus operations associated with such fields, like gradient and divergence, are defined using these classes. This design makes the code intuitive, allowing users to quickly cast PDEs into readable codes. The initial conditions, boundary conditions, and source/forcing terms are implemented using initial, boundary and force classes. These classes are readily extensible so that users can add custom initial conditions, source terms, and so on.
SARAS includes solvers for hydrodynamic flows, namely the incompressible Navier-Stokes initial value problem (IVP), as well as for scalar convection, like Rayleigh Benard Convection. Presently, we use semi-implicit Crank-Nicholson (Crank & Nicolson, 1947) method for timeadvancing the IVP. The solver uses Marker and Cell method (MAC) (Harlow & Welch, 1965) for discretizing the velocity and pressure fields.

Mathematics
The Navier-Stokes equations, which govern the dynamics of fluid flow, can be written as where u is the velocity field, p is the pressure field, f is the forcing term, and ν is the kinematic viscosity of the fluid. The fluid is assumed to be incompressible. Hence ∇ · u = 0, and density is constant (chosen to be unity).
If the velocity and pressure field at time t = t n are denoted as u n and p n respectively, then the corresponding fields at the next time-step, t = t n+1 , namely u n+1 and p n+1 , can be calculated as described below (Anderson, 1995;Ferziger & Peric, 2001;Patankar & Spalding, 1972). We compute an intermediate velocity field using the known values, u n and p n , as The forcing term has been neglected here for simplicity. Note that the diffusion term (also called the viscous term) is handled semi-implicitly, with equal contribution from u n and u * , as given below: The above equation has to be solved iteratively, and this is achieved through Jacobi iterations. The intermediate velocity field, u * , does not satisfy the continuity equation, and requires appropriate correction. This correction is obtained from the pressure correction term, which is calculated using the pressure Poisson equation, SARAS uses a Geometric Multigrid library to solve the above equation (Briggs et al., 2000;Wesseling, 2004). Presently the library employs the Full Multigrid (FMG) V-Cycle to solve the Poisson equation. Other methods like F-Cycle and W-Cycle are planned updates to the library in future.
Finally, using the above pressure correction, the velocity and pressure fields corresponding to the next time-step are obtained as The numerical implementation of the above procedure will be discussed in the next section.

Numerical Method and Implementation
SARAS uses finite-difference method (Ferziger & Peric, 2001) to calculate derivatives of the field variables. Presently the solver uses second-order central difference stencils to compute first and second derivatives. We use a highly optimized and fast array manipulation library, Blitz++ (Veldhuizen, 1998), for all array operations. Blitz++ also offers finite-difference stencils on uniformly spaced grids. SARAS augments this feature with grid-transformation terms (Anderson, 1995) to perform simulations on grids with non-uniform spacing.
SARAS offers a set of extensible boundary and initial conditions, as well as source terms.
Presently the solver can switch between periodic, Neumann, and Dirichlet boundary conditions. There also exists a mixed boundary condition class that can simulate many practical applications, such as a conducting heating plate on an adiabatic wall (Teimurazov et al., 2017). Currently the solver is restricted to Cartesian grids, but it will be extended to cylindrical, toroidal, and spherical grids in the future. The solver also supports adaptive time-stepping, where the Courant-Friedrichs-Lewy (CFL) condition (Courant et al., 1928) is used to dynamically compute the appropriate time-step.

Results
We validate our code using two very well-known test problems: lid-driven cavity and decaying turbulence. We simulate these problems using SARAS and compare the results with standard and validated solutions.

Problem 1
We solve the two-dimensional lid-driven cavity (LDC) problem using SARAS, and compare the results with those of Ghia et al. (1982). LDC is an important fluid system and it serves as a benchmark for testing numerical methods. This system consists of a square cavity of dimension 1 × 1 with no-slip boundary conditions on all the four walls. However, the top wall moves laterally to the right with a constant velocity of U = 1.0 that serves as the reference velocity for non-dimensionalization of the problem. The length of the side of the cavity, L = 1.0, is the reference length. The Reynolds number of the flow is approximately Re ≈ 1000.
At the start of the simulation, the fluid, which is at rest, is driven impulsively by the top lid. This results in the formation of a vortex at the upper-right corner of the cavity, and this vortex rapidly grows in size and occupies the entire region of the cavity. For the simulation we employ a 129 × 129 grid, and carry it out till t = 30. The solver computes this solution using 4 MPI processes in approximately 12 minutes on an Intel workstation. The output by SARAS at the final time is used for comparison with the results of Ghia et al. (1982) In the website we supply the Bash and Python scripts to compile and run this test case. After automatic execution of SARAS, a Python script reads the output from SARAS and compares the horizontal and vertical velocity profiles across the geometric center of the square cavity. The corresponding results from Ghia et al. (1982) are also available with the installation. We compare the two results and exhibit them in Figure 1. We observe that the profiles computed by SARAS match very well with those by Ghia et al. (1982), thus providing a strong validation for our code. This quick validation of the solver, which can be done as a part of its installation, is one of the strengths of this package.

Problem 2
We simulate decaying turbulence using SARAS with Taylor-Green vortex (Taylor & Green, 1937) as the initial condition. That is, where u 0 = 1 and k 0 = 1. We perform our simulation in a periodic box of size 1 × 1 × 1 (L = 1) with a grid resolution of 257 3 up to t = 3.0. Here, we nondimensionalize time using L/u 0 . The initial Reynolds number of the flow is Re = 1000. We choose a constant dt = 0.001 for time-integration. Besides, we use a uniform mesh along all the three directions.
To validate the accuracy of the finite difference scheme of SARAS, we compare the results of SARAS with those of a pseudo-spectral code TARANG, which has been benchmarked and scaled up to 196608 cores of Cray XC40, Shaheen II of KAUST (Chatterjee et al., 2018;Verma et al., 2013). Note that pseudo-spectral method yields very accurate derivatives (Canuto et al., 1988). We perform our spectral simulation using the same initial condition and grid resolution as above, up to t = 3.0. The comparison of the results from the two codes is discussed below. Comparison of SARAS results with those of TARANG provides another validation of SARAS.
The results of TARANG and SARAS are quite similar, as exhibited in Figures 2 and 3. In Figure  2(a) we exhibit the plots of the total energy ( ∫ dru 2 /2) vs. time for both the runs. As is evident from the plots, both the codes yield very similar evolution profiles for the total energy. In addition, we also compute the energy spectra for the results at t = 1. As shown in Figure  2(b), both the energy spectra are very similar. Interestingly, in both the plots, the energy spectra in the inertial range are close to the k −5/3 power law (Kolmogorov, 1941;Verma, 2019).   : For the simulation of decaying turbulence on a 257 3 grid, vector plots of the velocity field, and the density plots of the magnitude of velocity (|u|) computed at the horizontal mid-plane: for the data from TARANG(a, c), and SARAS(b, d) at t = 1 (top row) and t = 3 (bottom row).

Conclusions
This paper provides a brief description and validations tests of a new parallel finite-difference code SARAS. SARAS has been designed as a general PDE solver using the object-oriented features of C++. We validate SARAS using two test cases: lid-driven cavity and decaying turbulence.