Neutron star evolution by combining discontinuous Galerkin and finite volume methods

TL;DR

The paper addresses the challenge of simulating general relativistic hydrodynamics for neutron stars on modern HPCs by marrying high-accuracy discontinuous Galerkin (DG) methods with a compact finite-volume/finite-difference (FV/FD) scheme. The hybrid DG+FV/FD approach uses a Troubled Element Indicator to switch elements between DG in smooth regions and FV/FD near shocks or star surfaces, with a new WENOZ2 reconstruction, atmosphere treatment, positivity limiters, and a robust primitive-conservative conversion, all implemented in the Nmesh code. Key contributions include a compact FV/FD scheme exchanging data only with six nearest neighbors, a post-step a posteriori switching strategy, and extensive 3D GRHD NS tests (including unstable migrating and boosted stars) demonstrating improved stability and scalability over pure DG. The results show that the hybrid method reduces Gibbs oscillations, preserves mass and constraints better, and scales efficiently on large HPCs, enabling more realistic 3D neutron star and potentially binary neutron star simulations relevant to gravitational-wave and multimessenger astrophysics.

Abstract

We present here a new hybrid scheme that combines a discontinuous Galerkin (DG) method with compact finite volume (FV) and finite difference (FD) methods. The computational mesh is divided into smaller elements that touch but do not overlap. Like a pure DG method, our new hybrid scheme requires information exchange only at the surface of neighboring elements. This avoids the need for ghost zones that are usually many points deep in traditional FV implementations. Furthermore, unlike traditional FV implementations, that need information exchange between each element and its 26 surrounding neighbors on noncuboid meshes, our new hybrid method exchanges information only between each element and its six nearest neighbors. With this reduced communication, we aim to retain the high scalability of DG when using large supercomputers. In addition, the information exchange between adjacent elements is much simpler than in a traditional FV implementation, because we always have grid points at the interface, so that only surface interpolation is required. As a result it is much easier to implement adaptive mesh refinement. The goal is to use DG in elements with smooth matter fields and to fall back onto the more robust FV/FD method in elements that contain nonsmooth shocks or star surfaces. For this we devise trouble criteria to decide whether an element should be evolved with DG or FV/FD. We use the Nmesh program to implement and test the new scheme. We successfully evolve various single neutron star cases. These include the challenging cases of a neutron star initially in an unstable equilibrium migrating to a stable configuration and a boosted neutron star. These cases are simulated for the first time here in full 3D with general relativistic hydrodynamics using DG methods. We also describe additional numerical methods, such as the limiters and the atmosphere treatment we need for our simulations.

Figures (20)

• Figure 1: Grid points in 1D. On the left is one DG element, with two FV/FD elements to its right that have uniform grid spacing.
• Figure 2: 1D grid point arrangement for FV with WENO reconstruction. The numbering at the top shows the grid point $q_1$ = 0, 1, … . The dashed lines indicate the midpoints between grid points. The region between two consecutive midpoints is called a cell. For each cell we need to reconstruct fluxes at its left and right end. In the interior of the grid (in cells $q_1$ = 2, …, $N-2$), when we have five consecutive grid points available for interpolation, we use WENOZ reconstruction. In cells $q_1 = 1$ and $q_1 = (N-1)$ we drop the order of interpolation and use WENO3 instead. One stencil for each of WENOZ and WENO3 has been shown in the figure to illustrate their actual usage of the grid points. Even closer to the boundary, e.g. in the cell between grid point 0 and midpoint $\frac{1}{2}$ (marked by $\times$), we perform linear interpolation using fields at $q_1 = 0$ and $q_1 = 1$ to obtain fluxes at midpoint $\frac{1}{2}$. The same is done at midpoint $N-\frac{1}{2}$. The fluxes at the endpoints ($q_1 = 0$ and $q_1 = N$) can be directly computed from the fields there, without any interpolation. From the fluxes at the cell ends we construct a numerical flux as usual from data in this cell and the neighboring cell in the adjacent element. From the numerical fluxes at the cell ends we then obtain flux derivatives at the cell centers. For all interior cells the centers coincide with grid points, so that we immediately obtain flux derivatives at the grid points. However, the cell centers near the left and right end (marked by the $\times$) do not coincide with grid points. To obtain flux derivatives at grid points $q_1 = 0$ and $q_1 = N$ we use a weighted average of the flux derivatives at the two closest cell centers of this element.
• Figure 3: This figure shows the L2 norm of the error in the scalar field for a mesh that contains both DG and FV elements. We show the results for seven resolutions corresponding to h-refinement levels of $l = 0, 1, 2, 3, 4, 5, 6$ and $7$. The error is consistent with 2nd order convergence.
• Figure 4: The plot shows blast wave profiles for pressure $P$, rest-mass density $\rho_0$ and speed $v_x$ (times 10) at $t = 0.4$ after evolving the initial shocks in $P$ and $\rho_0$.
• Figure 5: This figure shows central density $D_\mathrm{C}$ for (S) with pure DG simulations corresponding to h-refinement levels of $l = 4, 5$ and $6$, normalized to the initial value. Large amplitude Gibbs oscillations originating from the star surface spread across it generating large amounts of noise. This is especially true near the beginning, but continues to plague the lowest resolutions throughout. As resolution increases, the star settles down into more systematic oscillations consistent with the star's fundamental oscillation frequency.