A semi-implicit unstructured operator-difference scheme for three-dimensional self-gravitating flows

A support operators (operator-difference) method has proven itself well for implicit simulations of different astrophysical fluid flows. Following the operator-difference approach, we construct nodal difference analogues of differential operators in Cartesian coordinates, where the conjugacy properties are the same as for original ones. Using the difference operators, we develop an Eulerian semi-implicit gas-dynamical solver with self-gravity on a three-dimensional collocated unstructured tetrahedral mesh. In the solver, only acoustic waves are treated implicitly, resulting to the only elliptic equation for a pressure on each time-step. The conjugacy properties of derived difference operators allow us to construct symmetric sign-definite matrices for this elliptic equation as well as for Poisson equation for a gravitational potential. The stability condition of the proposed scheme is milder, than the usual Courant-Friedrichs-Lewy condition for explicit solvers, and depends only on the gas velocity. Results of test problems’ simulations of low and high Mach number flows are presented.


Introduction
Numerical simulations of astrophysical fluid flows is one of the key instruments, which allow us to understand better a wide range of astronomical phenomena as well as to compare the theoretical predictions to observational data. An astrophysics is rich for different types of the fluid and plasma flows in the case, when one can use compressible fluid approximation without dissipation, so that the explicit schemes for hyperbolic systems, such as Godunov-type finite volume or SPH solvers, can be successfully applied for such kind of problems [1,2]. But there are the astrophysical problems, such as, for example, star formation from protostellar clouds, collapses of stars' iron cores and supernovae explosions, where the implicit or semi-implicit schemes may possibly outperform the explicit solvers. It is a common situation in astrophysical modelling, when the compact high-density regions of a collapsed matter with high speed of sound (such as protostars, protoneutron stars) appear in a computational domain, where the velocity of the flow is very low compared to the acoustic speed. At the same time, everywhere else in the domain, the flow can be considered as compressible, and strong shocks and other different flow discontinuities may form. The mesh size during such astrophysical simulations is usually refined in the regions of collapsed objects up to several orders of magnitude, compared to the cell size on the outer boundary. It is well known, that the explicit solvers are only conditionally stable, and the Courant-Friedrichs-Lewy (CFL) condition limits the time-step in simulations. In the mentioned types of astrophysical flows, the global CFL condition comes mostly from the 2 acoustic oscillations of these resolved collapsed objects and may be very restrictive, which makes such kind of simulations computationally expensive and challenging. Thus, there is a hope, that the implicit treatment of acoustic waves in hydrodynamical schemes may allow to reduce the computational cost of such astrophysical simulations.
Implicit evaluation of the pressure field is a standard approach in incompressible flow simulations. The idea to extend the classical semi-implicit incompressible numerical scheme by Harlow and Welsh [3] to compressible flows was proposed in mid-80's by Casulli and Greenspan [4], but only quite recently this approach was used to construct the so-called All-Mach number solvers for gas dynamics (see e.g. works [5,6,7,8] on structured and [9, 10] on unstructured grids, and references therein) and for MHD equations [11,12], written in conservative form. They are the pressure-based solvers by its nature, but they can work as high resolution shock-capturing methods as well in presence of discontinuities. The usage of semi-implicit solvers results to the solution of an elliptic equation with a Laplace-type operator for the pressure on each time-step, and it is very important to construct a consistent difference analogue of the solving differential problem, which provides the linearized systems to be cheaply solvable. To obtain a symmetric sign-definite matrices for the linear elliptic problems, we use a Support (basic) Operators Method (operator-difference approach) [14] to construct the difference analogues of differential operators in three dimensions. The support operators method has also many common features with mimetic finite difference [13] approach.
The idea of the support operators method is that the difference operators should be derived in the way to fulfil the same properties as the continuous operators. Thus, the constructed operators should satisfy the difference analogues of integral relations such as Green formula. This approach allows to obtain the so-called completely conservative fully implicit finite-difference schemes for gas dynamics and MHD equations [14] on unstructured grids. For 2D simulations of magnetorotational supernova explosion this method was used e.g. in [15] with a cell-nodal approximation [16,17].
The paper is organized as follows. In the next subsection we derive the nodal difference operators of gradient and divergence in three dimensions. In the 3-rd subsection we provide the description of the solving gravitational gas-dynamics system and derive a semi-implicit numerical scheme for these equations, and in the last subsections we present and discuss the results of numerical simulations of several benchmark problems in low and high Mach number regimes.

Difference analogues of differential operators
In Support operators method the difference analogues of differential operators are built in pairs. Further we are going to consider only the nodal gradient and divergence operators, while the other difference analogues can be obtained in a similar way. To introduce difference GRAD and DIV operators one should define finite dimensional linear grid spaces and scalar products of grid functions in them. The first operator can be obtained using any numerical differencing technique, while the second one has to satisfy the grid analogue of the Green formula in corresponding grid space. Thus, the resulting operators are conjugated to each other. The GRAD operator can be derived in two different ways: either using the finite element approach with linear basis functions inside the tetrahedra (for 2D case see e.g. [18]) or with a finite volume approach using piecewise-constant basis functions on the nodal edge-based (median-dual) cells, as in Godunovtype methods on unstructured grids (see e.g. [19]). The grid gradient operator reads: where p is a nodal scalar function, G j is a set of tetrahedral cells adjoined to sought j-th node, i is a cell index, k is a node index inside the i-th cell, n k and S k are unit normal and area to the tetrahedron's triangular face, which does not contain the node k in its vertexes, and To obtain the consistent approximation of the divergence operator, the Green formula and its difference analogue should be used: where ∇ and ∇· are the gradient and the divergence operators correspondingly, v is a vector function, N n is a total number of the nodes, and V is a volume of the computational domain. The surface integral in the right hand side of the formulae (2) correspond to the so-called Boundary operator Φ in a formalism of the support operators method. Its evaluation can be found in [16,17,20]. The boundary operator is used to formulate the boundary-value problems in the operator-difference form. We do not consider it here, because in the solver we use ordinary finite volume boundary conditions via ghost cells approach.
To derive the DIV operator in the interior nodes of the computational domain, the Green formula in (2) is written without the integral in the right hand side. By substituting the grid gradient operator (1) in the analogue of the Green formula (2) and rearranging the terms, the following approximation for the divergence operator can be obtained We note, that the latter operator has the same form, as the initially derived gradient one (1) due to the usage of purely nodal approximation, while its form and stencil is different from the cell-nodal DIV operator. The divergence operator (3) is conjugated to (1), and thus, the grid Laplacian transforms into the symmetric and sign-definite matrix. The usage of the support operators technique allows to derive the full set of vector calculus operators, which satisfy the conjugacy properties. Finally, the difference analogues of differential operators can be written in a finite-volume form [21] where H j is a set of neighbouring nodes to the sought node with index j, S nj and n nj are the area and the unit normal to the n-th face of the node-centred median-dual control volume. This approximation is exactly the same as (1) [21].
3. Semi-implicit numerical scheme for the equations of gas dynamics with gravitation 3.1. Euler equations with self-gravity Gas dynamical equations with self-gravity can be written in the following conservative formulation: In the system above, the first three equations without source terms in right hand sides correspond to ideal Euler equations, while the last one is the Poisson equation for a gravitational potential Ψ, and the source term corresponds to presence of an acceleration by gravity, G is a gravitational constant. In ideal Euler equations (5) ρ is a mass density, v is a flow velocity; P is a pressure, E = v 2 2 + E th is a total energy density, E th is an internal energy density. A thermodynamic equation of state is needed to close the system of equations (5). In this paper we consider only a perfect gas equation of state with an adiabatic index γ: Extension of a semi-implicit solver (see next subsection) to an arbitrary equation of state is straightforward and can be found in e.g. [6,7]

Poisson solver
The Poisson equation should be solved on each time-step to calculate the value of the gravitational potential, and hence, a gravitational force. The support operators method allows to construct a discrete Laplace operator in the form of a symmetric and sign-definite matrix.
Purely nodal approximation requires the stencil, which consists of the second-order neighbours for each node. The boundary conditions are projected to the boundary nodes and included in the operator-difference form of the solving problem. For example, the Dirichlet boundary-value problem reads [16,17] ( where Ψ γ is a boundary value of the gravitational potential, I is a unit operator, δ is an operator, which is equal to unity in the boundary nodes, and zero -in the interior ones [17]. For other types of boundary conditions (e.g. an absence of the force on the boundary) we solve the Poisson using the ghost cells approach. The resulting linear system can be solved with a matrix-free conjugate gradient method. The gravitational acceleration then can be calculated as a nodal gradient of the potential. In Fig.1  The computational domain is cube (0 < x, y, z < 5) with no force conditions on the symmetry boundaries and Dirichlet conditions on outer ones.

Hydro-solver
The difference gas dynamical system of equations can be written in the following form on a collocated mesh This system describes the flow evolution from time t = t n to t = t n+1 = t n + τ . The last terms in the left parts of the second and third equations in (8) should be evaluated implicitly according to [4,6], which allows to treat implicitly only the acoustic waves in the gas dynamical system (5). The source terms are given by the formulae S v = −ρ n+1 ∇ × Ψ, S e = −ρ n+1 v n · ∇ × Ψ. The gravitational potential is solved for the density distribution on the (n + 1)-th time-step.
The difference operator with tilde ∇ × · in explicitly evaluated advective fluxes corresponds to an upwind monotonic scheme for advection. It is written together with a Rusanov-type numerical dissipation: where dissipation is written in a finite-volume form (4), in which u and F are the vector of conservative variables and the flux vector correspondingly: We use | ∂F ∂u |= max(|v| i , |v| j ), thus, dissipation coefficient is independent on the speed of sound [6]. The explicit advective scheme is implemented in the same manner, as the edge-based Godunov-type finite volume solvers on median-dual grids [22], so that the advective fluxes F in (9) are calculated at the centres of the edges between the nodes. The numerical scheme is conditionally stable with a milder CFL condition, which depend only on the flow velocity and not on the speed of sound: where ∆X is a characteristic local mesh size. By substituting the momentum from the momentum equation in (8) to the energy equation, we can obtain the following nonlinear equation: where h = E th + P ρ = γP (γ−1)ρ is a thermodynamic enthalpy. The latter system is nonlinear, and the usual approach to solve it is a usage of nested Newton-type iterations (r-iterations, Picard iterations) [6,7], which converge rapidly (usually 2-3 iterations are enough). The resulting linearised system reads: To obtain the solution from n-th time level to n + 1-th, the following procedures should be done: (1) solving the mass conservation equation in (8) with the explicit scheme and obtaining ρ n+1 (2) solving the Poisson equation (7) with the density ρ n+1 for the gravitational potential and calculating the forces in the source terms of gas-dynamical system (8) (3) defining P n+1,1 = P n , h n+1,1 = h n and v n+1,1 = v n (4) r-iterations, r = 1, r p (a) Solving (13) and obtaining P n+1,r+1 (b) recalculating h n+1,r+1 and calculating v n+1,r+1 from the momentum equation (8) using P n+1,r+1 (5) defining P n+1 = P n+1,rp , h n+1 = h n+1,rp and v n+1 = v n+1,rp (6) calculating the energy using a semi-implicit equation from the system (8).
This procedures allow to obtain a numerical solution of the gas-dynamical system (5) with a first order in space and time. In this study we use r p = 3 Picard iterations.
To extend this scheme to higher order in space, we use the same approach as in explicit finite volume Godunov-type schemes on unstructured meshes, as it is done in e.g. [22] (see also [23] for a different approach). We apply a piecewise-linear reconstruction of conservative gas-dynamical variables in advection upwind fluxes (9) with Barth-Jespersen slope limiters [19], thus, the explicit advective part of the solver is calculated with the second order in space. For this reason, the resulting order of spatial approximation of the scheme is higher, than one. Note, that the difference operators like (4) provide the second order in space on the mesh, consisting on equilateral tetrahedra (or equilateral triangles in 2D), because these operators are analogous to the usual central differences on the uniform grid, so that on such kind of meshes the solver has the second approximation order in space.
The linearised problem (13) has a matrix of a Laplace-type second order difference operator and, in addition, a diagonal, filled with the positive elements. The conjugated difference operators allow to obtain an easily reversible symmetric and positive definite matrix. We use a matrix-free conjugate gradient method to solve the system (13), which is very efficient.

Test computations
The code was developed in two variants: 2D and 3D Cartesian versions. We have tested our code in several benchmark problems in high and low Mach number regimes as well as simulated a test problem of a dust collapse to check the code's ability to simulate the flows with self-gravity.  Figure 2. Velocity, pressure and density fields (the first three figures, blue lines) of the gas in the Sod shock tube at time T = 0.45, compared to the reference solution given by a second-order explicit 2D local Lax-Friedrichs solver (black lines) on a finer mesh with N x = 500 nodes. On the last picture the density plot is shown for 2D explicit (red line) and 2D semi-implicit (blue line) solvers with the same mesh resolution.
In Fig.2 the usual Sod shock tube calculation is shown for the 3D (the first three pictures) and 2D versions of our code at the time T = 0.45 for the mesh with N x = 200 nodes in X-direction and N y = N z = 10 ones in Y-and Z-directions. The initial conditions are as usual: ρ(x < 1) = 0.125, P (x < 1) = 0.1 and ρ(x ≥ 1) = 1, P (x ≥ 1) = 1. The velocity is   Figure 3. Velocity, pressure and density fields of the gas in a slowly moving contact discontinuity problem at time T = 50. The blue lines correspond to our numerical solution, while the red dash-dotted line and dashed purple line are initial and final location of the contact discontinuity correspondingly.
In Fig.3 the density profile of a slowly moving contact discontinuity is presented. The initial conditions are the following: v = (0.01, 0, 0), P = 3 everywhere, while the density ρ(x < 0.25) = 50, ρ(x ≥ 0.25) = 1. The problem is solved until the final time T = 50, which corresponds to the contact discontinuity located at x = 0.75. The profile is smeared over ∼ 10 grid points. On this problem, the semi-implicit scheme greatly outperforms the explicit solver by a factor of ∼ 40. The time step τ adv (11) is equal to ∼ 200τ CF L .
The next test is the Sedov-Taylor planar blast wave solution, which was performed in a 2D cylindrical geometry. The shock front location is R ≈ √ T . Fig.4 shows the density and pressure profiles at T = 0.2. The 2D mesh, used in the simulations, is refined at the centre of explosion by a factor of 10. It has 8756 nodes. The numerical solution shows a good agreement with analytics. In the refined central part (r = x 2 + y 2 1) of the computational domain we have a region with a high sound speed and low velocity, while at r 1 the strong shock wave is pronounced. At the end of the simulation, the advective time-step was about 50 times higher, than the CFL time-step τ CF L . In real astrophysical conditions this ratio of time-steps could be even higher.   We have tested our code including the self-gravity module and simulated a well-known problem on the spherical collapse of a pressureless star (dust collapse), which was originally considered by Colgate and White [24,25]. The octant-symmetric 3D results for the density at dimensionless time T = 1.791 are presented in Fig.5.We use the system of units, which removes the multiplier 4πG in the Poisson equation (7). In these units the free-fall time t f f = π √ 3 2 √ 2 ≈ 1.9238, the initial cloud radius and density are equal to unity. We use the mesh, which is refined to the centre of the computational domain. The snapshot of our simulation shows a physically correct evolution of collapsing matter, although at the developed stage of the collapse the density field has a discrepancy from the analytical solution (blue line in Fig.5, right) due to the high spatial and temporal steps.

Conclusion
In this paper we have constructed the semi-implicit nodal (edge-based) pressure-based solver for gravitational gas dynamics on unstructured tetrahedral meshes. The usage of the Support operators technique results to linear systems with symmetric and sign-definite matrices for the pressure system in semi-implicit gas-dynamical solver and for the Poisson equation for the gravitational potential. The code works efficiently both in low and high Mach number regimes, and the spatial accuracy is close to the second order.