Introducing Flow-er: a Hydrodynamics Code for Relativistic and Newtonian Flows

1. Introducing Flow-er: a Hydrodynamics Code for Relativistic and Newtonian Flows Patrick M. Motl (motl@phys.lsu.edu), Joel E. Tohline, & Luis Lehner (Louisiana State University) Features of Flow-er Abstract • Implements the Kurganov-Tadmor scheme for systems of hyperbolic conservation laws • Flow-er can calculate flows in one, two or three spatial dimensions • User specifiable spatial accuracy of advection operators with choice of six reconstruction schemes (1) linear reconstruction with minmod function (2) linear reconstruction with VanLeer slope calculation (see Anninos & Fragile 2003) (3) linear reconstruction with UNO limiter (Harten & Osher 1987) (4) quadratic reconstruction with Liu’s limiter (Liu & Tadmor 1998) (5) quadratic reconstruction with Kurganov’s limiter (Kurganov & Petrova 2001) (6) PPM reconstruction with optional discontinuity detection (Colella & Woodward 1984) • User specifiable spatial accuracy for finite difference operators used in evaluating source terms • Utilize either second or third order formulations of Runge-Kutta time integrators (Shu & Osher 1988) • Currently implement three equation sets (1) Euler’s equations for ideal Newtonian Flows (2) Relativistic hydrodynamic equations following from T = 0 for an ideal fluid where g = det(g), h = 1 +  + P/, P = (-1), Vi = ui/ut, D = W, W = -gut (3) Relativistic hydrodynamic equations following from T = 0 • The relativistic formulations can evolve systems with a general, static metric We present a new numerical code (Flow-er) for calculating astrophysical flows in 1, 2 or 3 dimensions. We have implemented equations appropriate for the treatment of Newtonian gravity as well as the general relativistic formalism to treat flows with either a static or dynamic metric. The heart of the code is the recent non-oscillatory central difference scheme by Kurganov and Tadmor (2000; hereafter KT). With this technique, we do not require a characteristic decomposition or the solution of Riemann problems that are required by most other high resolution, shock capturing techniques. Furthermore, the KT scheme naturally incorporates the Method of Lines, allowing considerable flexibility in the choice of time integrators. We have implemented several interpolation kernels that allow us to choose the spatial accuracy of an evolution. Through the Cactus framework or independent code, Flow-er serves as a driver for the hydrodynamical portion of a simulation utilizing adaptive mesh refinement or on a unigrid. In addition to describing Flow-er, we present results from several test problems. We are pleased to acknowledge support for this work from the National Science Foundation through grants PHY-0326311 and AST-0407070. ∂t + ∂i(vi) = 0 ∂t(vj) + ∂i(vjvi + ijP) = 0 ∂tE + ∂i(vi(E + P)) = 0 Non-Oscillatory Central Difference Schemes, Briefly We want to solve a conservation law (or set of conservation laws) of the form, ∂tq +∂i (q) =  (eq1) where q is the conserved quantity,  is the corresponding flux function and  is the source function which for the present we will assume is zero. In the Nessyahu-Tadmor (1990; NT) scheme, this equation is solved numerically by constructing the updated value of q at a staggered grid location. As illustrated in Figure 1, the cell averages at timestep n, qnj are used to build qn+1j+1/2 as follows: qn+1j+1/2 = qnj+1/2 - 1/x [tt+t j+1d - tt+t j d ] (eq2) where qnj+1/2 is reconstructed to the desired order of accuracy and the flux integrals are evaluated via quadrature or using Runge-Kutta integrators. For second order accuracy in time, the integrals can be evaluated using midpoint values predicted from the conservation law to yield the original NT scheme qn+1j+1/2 = qn j+1/2 - t /x [n+1/2j+1 - n+1/2j] (eq3) ∂tD + ∂i(DVi) = 0 ∂tSj + ∂i(SjVi + -g(Pgij - PgtjVi)) = --gjT ∂tE + ∂i(EVi + -g(Pgit - PgttVi)) = --g tT Sj = hWuj + -gPgtj, E = hWut + -gPgtt ∂tD + ∂i(DVi) = 0 ∂tSj + ∂i(SjVi + -gPij) = -gjT ∂tE + ∂i((E - -gP)Vi) = -g tT Sj = hWuj, E = hWut + -gP qn+1j+1/2 q nj-1 Example Calculation q nj q nj+1 x j x j+1 x j-3/2 x j-1/2 x j+1/2 x j+3/2 minmod,  = 1.5 VanLeer Figure 1 The NT scheme in one spatial dimension. The updated value qn+1j+1/2 is calculated at a staggered grid location from the reconstructed value of qnj+1/2 and time-centered fluxes through the faces at xj+1 and xj. Due to the Courant condition, q will remain smooth during the timestep at xj+1 and xj and a Riemann solver is not required to properly upwind the solution to compute fluxes. The updated, staggered, solution can then be averaged back onto the original grid or one may alternate between staggered and non-staggered grids. Kurganov & Tadmor (2000) refined this scheme by noting that the finite signal speed limits the required size of the staggered grid cells. They conducted an analogous calculation to (eq2) and took the limit of t  0 to yield a system of equations that were discrete in space and were ordinary differential equations with respect to time. The staggered grid for constructing the KT scheme is illustrated in Figure 2. The left and right limits of the non-uniform grid that the solution is realized on have width given by the largest signal speed at each cell interface. UNO PPM q nj-1 q nj q nj+1 Liu Kurganov xLj-1/2 xRj-1/2 xLj+1/2 xRj+1/2 x j-3/2 x j-1/2 x j+1/2 x j+3/2 Figure 2 The KT scheme grid in one spatial dimension. Because of the Courant condition, a wave from the possibly discontinuous region at xj-1/2 can travel only as far as xRj-1/2 into cell j in a timestep whereas it was effectively assumed to be able to travel to xj in the NT scheme. A solution can be calculated across these variable width cells using (eq2), forming the basis for the KT scheme. Figure 3 The mass density from a simple shock tube calculation utilizing the six reconstruction methods incorporated in Flow-er. All calculations use the same number of points and all use the R-K3 integrator.The minmod, VanLeer and UNO calculations are linear reconstructions that limit the slopes near local extrema. The PPM, Liu and Kurganov calculations use quadratic reconstructions and limit the interpolation coefficients to avoid introducing spurious extrema in the reconstructed distributions. Discontinuity detection was not used in the PPM calculation. After taking the limit of t  0, the resulting KT scheme takes the very simple form, dqj(t)/dt = - [Hj+1/2(t) - Hj-1/2(t)]/x (eq4) where the flux, H, is given by Hj+1/2(t) = 1/2[(q+j+1/2(t)) + (q-j+1/2(t))] - 1/2 aj+1/2(t) [q+j+1/2(t) - q-j+1/2(t)] (eq5) aj+1/2 is the signal speed at the zone interfaces and q+/-j+1/2 are the edge values of the reconstruction polynomials. Future Work In future work we will implement an additional set of evolution equations for Newtonian systems in general curvilinear coordinates. We are also developing a general Poisson solver to allow us to simulate self-gravitating, Newtonian flows. Finally, we are incorporating Flow-er as one of the physics applications in the adaptive mesh refinement code HAD. When this is finished, we will also be able to evolve systems with a dynamic metric by coupling Flow-er to a finite difference code for solving Einstein’s equations. Advantages of the Kurganov-Tadmor Scheme • No staggered grid as required by NT scheme, all variables remain naturally cell centered • Algorithm does not require solution of Riemann problems or characteristic tracing • Spatial accuracy of advection dictated by choice of reconstruction polynomial • Temporal accuracy given by ODE integrator • Simple, conservative flux formula • Numerical viscosity within the scheme is independent of the timestep size, allowing small timesteps with minimal numerical diffusion • Enables simple, modular programming with significant code reuse for different conservation law systems • Solution quality for 2nd and higher order accurate KT formulations can be comparable to more expensive techniques such as PPM with Riemann solvers References Anninos, P. & Fragile P.C. 2003, ApJS, 144, p243 Colella, P. & Woodward, P.R. 1984, JCP, 54, p174 Harten, A. & Osher, S. 1987, SINUM, 24, p279 Kurganov, A. & Tadmor, E. 2000, JCP, 160, p241 Kurganov, A. & Petrova, G. 2001, Numer. Math., 88, p683 Liu, X.-D. & Tadmor, E. 1998, Numer. Math., 79, p397 Nessyahu, H. & Tadmor, E. 1990, JCP, 87, p408 Shu, C.-W. & Osher, S. 1988, JCP, 77, p439 An online archive of papers developing and using Non-Oscillatory Central Difference schemes is available from Eitan Tadmor’s website at http://www.cscamm.umd.edu/people/faculty/tadmor/centralstation