Hydrodynamics in full general relativity with conservative adaptive mesh refinement

There is great interest in numerical relativity simulations involving matter due to the likelihood that binary compact objects involving neutron stars will be detected by gravitational wave observatories in the coming years, as well as to the possibility that binary compact object mergers could explain short-duration gamma-ray bursts. We present a code designed for simulations of hydrodynamics coupled to the Einstein field equations targeted toward such applications. This code has recently been used to study eccentric mergers of black hole-neutron star binaries. We evolve the fluid conservatively using high-resolution shock-capturing methods, while the field equations are solved in the generalized-harmonic formulation with finite differences. In order to resolve the various scales that may arise, we use adaptive mesh refinement (AMR) with grid hierarchies based on truncation error estimates. A noteworthy feature of this code is the implementation of the flux correction algorithm of Berger and Colella to ensure that the conservative nature of fluid advection is respected across AMR boundaries. We present various tests to compare the performance of different limiters and flux calculation methods, as well as to demonstrate the utility of AMR flux corrections.

Figure 1

A visualization of a refinement level boundary and its treatment in the flux correction algorithm. The top shows cells in the $x$ direction on refinement level $L$ while the bottom shows equivalent cells for the $L+1$ refinement level (here the refinement ratio is 2). Fluxes are symbolized by arrows. On the bottom level the blue cells (“type B” in the discussion in the text) and those to the left on level $L+1$ are boundary cells and will have their values set by interpolation from level $L$ following an evolution step on level $L$. Because of truncation error, subsequent evolution on level $L+1$ will give a flux differing from that computed on the parent level $L$. Consequently, when the new fine grid solution is injected back to the parent level (in cells to the right of the red/dotted pattern cell), the solution about the boundary will no longer be consistent with the flux previously computed there. To correct this, the fluid quantity in the red/dotted pattern cell is adjusted to exactly compensate for the difference in flux computed between the coarse and fine levels.

Figure 2

Top: The natural log of the $L^2$ norm of the constraint violation, $C_a≔H_a-\square x_a$, for a Brill wave evolution (i.e. natural log of $\sqrt{\int |C_a{|}^2{d}^{2}x/\int d^{2}x}$). The three resolutions shown are scaled assuming fourth-order convergence. Time is shown in units of $M$, the total ADM mass of the spacetime, and the constraints are multiplied by $M$ to make them dimensionless. The lowest resolution has a grid spacing of $h=1.56M$. Middle/Bottom: The value of the metric component $g_{tt}$ evaluated at $(x,y,z)=(0,50M,0)$ (middle) and the difference in this quantity between low and medium resolution, and medium and high-resolution (bottom), the latter scaled so that the two curves should coincide for fourth-order convergence.

Figure 3

The $L^2$-norm of the constraint violation ($C_a≔H_a-\square x_a$) in the equatorial plane for a boosted BH simulation with $v=0.25$. The three resolutions shown are scaled assuming fourth-order convergence. Time is shown in units of $M_{\mathrm{BH}}$, the ADM mass of the BH in its rest frame, and the norm of the constraints is multiplied by $M_{\mathrm{BH}}$ so as to make it dimensionless.

Figure 4

Density at $t=0.4$ for Riemann Test 1 (RT1) with different reconstruction methods and the HLL flux scheme at resolution $N=400$. The inset shows the shock and contact discontinuity. The exact solution was generated using the code of 104.

Figure 5

Density at $t=0.4$ for Riemann Test 2 (RT2) at different resolutions for HLL-WENO5. The average convergence rate in this case is 0.76. The thin shell of material between the shock and the contact is particularly difficult to resolve.

Figure 6

Density at $t=0.4$ for Riemann Test 3 (RT3), a collision problem, for different reconstruction methods with HLL at resolution $N=400$. The post-shock oscillations are largest in the MC case and smallest for PPM.

Figure 7

Density at $t=0.4$ for the transverse test (TVT) at different resolutions for HLL-WENO5. The average convergence rate in this case is 0.83.

Figure 8

Density at $t=2.0$ for the shock-heating test with different reconstruction methods and the HLL flux scheme. The inset focuses on the shock front.

Figure 9

Density contours (30 equally spaced in the logarithm) for the Emery step problem. The upper (lower) two plots show results for resolution $240\times 80$ ($480\times 160$). For each resolution, the upper plot shows results for WENO5 reconstruction, and the lower for PPM. The respective minimum and maximum densities, $({\rho }_{\min },{\rho }_{\max })$, are (1.0, $1.0\times {10}^2$), (0.55, $1.1\times {10}^2$), (0.82, $1.1\times {10}^2$), and ($6.8\times {10}^{-2}$, $1.1\times {10}^2$).

Figure 10

Density contours (30 equally spaced in the logarithm) for the 2D Riemann problem with, from left to right, top to bottom: HLL and MC, HLL and WENO5, Marquina and MC, and HLL and MC with mesh refinement. The respective minimum and maximum densities, $({\rho }_{\min },{\rho }_{\max })$, are ($1.1\times {10}^{-2}$, 7.0), ($8.2\times {10}^{-3}$, 8.1), ($9.1\times {10}^{-3}$, 7.1), and ($7.6\times {10}^{-3}$, 7.0). For the first three simulations a resolution of $400\times 400$ was used. For the final simulation, a refinement region (red box) was placed in the middle with equivalent resolution, while the remaining grid has half the resolution (i.e. this simulation has lower resolution overall). A CFL factor of 0.5 was used throughout.

Figure 11

The $L^2$ norm of the difference between the numerical and exact value of the fluid quantity $D$ (divided by the norm of the exact solution) for Bondi accretion with MC and WENO5. The low resolution was run with a grid spacing of $h=0.078M_{\mathrm{BH}}$ while the medium and high has twice and 4 times the resolution, respectively. The results are scaled for second-order convergence.

Figure 12

The $L^2$-norm of the constraint violation ($C_a≔H_a-\square x_a$) in the equatorial plane for a boosted NS simulation with $v=0.5$ (using HLL flux calculation and WENO5 limiter). The three resolutions shown are scaled assuming second-order convergence. Time is shown in units of $M_0$, the ADM mass of the NS in its rest frame, and the norm of the constraints is multiplied by $M_0$ so as to make it dimensionless.

Figure 13

The relative variation of the maximum central density from its initial value (${\rho }_{\max }/{\rho }_{\max }(t=0)-1$) during a boosted NS simulation with $v=0.5$ for various reconstruction methods and for two different resolutions in the case of WENO5.

Figure 14

The relative variation of the fluid rest-mass ($\overline{M}/\overline{M}(t=0)-1$) for a boosted NS test with (dotted, red line) and without (solid, black line) flux corrections at mesh refinement boundaries. For this test, the mesh hierarchy is fixed so that at $t\approx 20M_0$ the NS has moved from being contained entirely on the finest resolution to being contained entirely on the second finest resolution. At $t\approx 60M_0$ the NS has moved from the second to the third finest resolution.