General-relativistic neutron star evolutions with the discontinuous Galerkin method
GGeneral-relativistic neutron star evolutions with the discontinuous Galerkin method
François Hébert,
1, 2, ∗ Lawrence E. Kidder, and Saul A. Teukolsky
1, 2 Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York 14853, USA Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125, USA (Dated: August 29, 2018)Simulations of relativistic hydrodynamics often need both high accuracy and robust shock-handlingproperties. The discontinuous Galerkin method combines these features — a high order of convergencein regions where the solution is smooth and shock-capturing properties for regions where it is not— with geometric flexibility and is therefore well suited to solve the partial differential equationsdescribing astrophysical scenarios. We present here evolutions of a general-relativistic neutron starwith the discontinuous Galerkin method. In these simulations, we simultaneously evolve the spacetimegeometry and the matter on the same computational grid, which we conform to the spherical geometryof the problem. To verify the correctness of our implementation, we perform standard convergenceand shock tests. We then show results for evolving, in three dimensions, a Kerr black hole; a neutronstar in the Cowling approximation (holding the spacetime metric fixed); and, finally, a neutron starwhere the spacetime and matter are both dynamical. The evolutions show long-term stability, goodaccuracy, and an improved rate of convergence versus a comparable-resolution finite-volume method.
I. INTRODUCTION
Numerical simulations are a crucial tool in the studyof core-collapse supernovae, compact binary mergers, ac-cretion disks with relativistic jets, and other energeticastrophysical sources. In these events, the dynamics aregoverned by the high-density matter and its coupling tothe strong gravitational field. Nuclear reactions, neu-trino physics, and magnetic fields can also play significantroles. Because of the highly nonlinear nature of the un-derlying general-relativistic hydrodynamics (GR-hydro),simulations are necessary to obtain observable predic-tions from physics models. Achieving sufficient accuracyin the simulation outputs (e.g., gravitational waveforms,ejected masses, and nucleosynthesis products) remains achallenge, however. High resolution is needed to resolvemultiscale fluid flows, and the presence of shocks in thematter reduces the accuracy of the numerical schemes.The standard approach taken in present-day GR-hydrocodes is to cast the partial differential equations into con-servative form and discretize them using a finite-volume(FV) method (see reviews [1–3] for an overview and his-tory). FV methods are favored for their robustness andthe various “shock-capturing” schemes that enable themto handle fluid shocks and stellar surfaces. The Einsteinequations for the spacetime geometry are typically solvedwith a finite-difference method on the same grid but caninstead be solved with a spectral method on a differentcomputational grid [4]. Over the past decade, the applica-tion of improved high-resolution shock-capturing schemes(e.g., the piecewise parabolic method (PPM) [5, 6] and theweighted essentially non-oscillatory (WENO) scheme [7])and higher-order difference schemes has led to significantadvances in the accuracy and stability of the numericalresults (e.g., for core-collapse supernovae [8, 9], binary ∗ [email protected] mergers [10–12], and accretion flows [13–15]). In spite ofthese successes, FV methods have inherent limitationswhen used as high-order methods: the large stencils re-quired for the differencing and shock-capturing schemesmake it difficult to adapt the grid to the problem geometryand can also lead to challenges in efficiently parallelizingthe algorithm.In the pursuit of improved accuracy and efficiency, dis-continuous Galerkin (DG) methods have recently emergedas a promising contender for astrophysical problems. DGmethods share properties with both spectral methods andFV methods — they inherit the high-order accuracy ofthe former for smooth solutions while maintaining therobust shock-handling properties of the latter. They aregeometrically flexible, enabling the use of grids adaptedto the problem geometry. They are well suited to hp -adaptivity, where the grid resolution can be set eitherby adjusting the order of the polynomial approximationwithin an element ( p -refinement) or by adjusting the sizeof the element ( h -refinement). Finally, DG methods arelocally formulated, enabling efficient parallelization andgood scaling.The application of DG methods to problems in rela-tivistic astrophysics is recent and remains exploratoryin nature. With several of these explorations focusingon the evolution of the spacetime geometry, different for-mulations of Einstein’s equations have been investigated.In an early study, Zumbush [16] obtained a space-timeDG scheme for the linearized vacuum Einstein equationsin harmonic gauge. For the commonly used Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of theEinstein equations, Field et al. [17] and Brown et al. [18] developed DG schemes in spherical symmetry. Morerecently, Miller and Schnetter [19] developed an operator-based (vs the typical differential equation-based) DGdiscretization of the BSSN equations and showed successin evolving three-dimensional (3D) test problems. Usinga new first-order form of the constraint-damping Z4 for-mulation (FO-CCZ4), Dumbser et al. [20] evolved a single a r X i v : . [ g r- q c ] A ug black hole (BH) spacetime using a puncture and showeda short-timescale “proof of concept” evolution of a binaryBH system.Efforts on the hydrodynamics side began with Radiceand Rezzolla [21], who presented a formulation of DG forthe evolution of fluids in curved spacetimes and evolveda neutron star (NS) in spherical symmetry. In their work,the spacetime is treated self-consistently by satisfying aradial constraint equation. In Ref. [22], Zhao and Tangimplemented DG with a WENO shock-capturing schemefor special-relativistic hydrodynamics in one and two di-mensions. Bugner et al. [23] were the first to apply DGto a 3D astrophysical fluid problem, evolving a NS in theCowling approximation (i.e., fixed background metric).In a DG code using a task-based parallelism paradigm( SpECTRE ), Kidder et al . [24] showed special-relativisticmagnetohydrodynamics tests in two and three dimensions.Anninos et al . [25] and Fambri et al . [26] (see also Ref. [27])implemented DG schemes with adaptive mesh refinement(AMR) for applications to special- and (fixed-background)general-relativistic magnetohydrodynamics, and showedresults in two and three dimensions.In this paper, we use a DG method to evolve a NSin coupled GR-hydro in three dimensions (prior effortsin this direction are the subject of theses by Hébert [28]and Bugner [29]). As tests of our implementation, wealso evolve a NS in the Cowling approximation and aKerr BH. In these simulations, we investigate the use ofcubed-sphere grids conforming to the spherical geometryof the BH and NS problems. We adopt the DG formula-tion described by Teukolsky [30], using the generalizedharmonic formulation of Einstein’s equations [31–33] andthe València formulation [1] of the general-relativistichydrodynamics.We implement our DG code in the framework of theSpectral Einstein Code [34] (
SpEC ). SpEC combines amultidomain penalty spectral method to evolve binaryBH spacetimes [35–37] with a FV method to evolve thematter in BH-NS [4] and NS-NS [12, 38] systems. Our DGGR-hydro code is independent from
SpEC ’s FV compo-nent and is instead built on the algorithms from
SpEC ’svacuum spectral code: spectral bases and differentiations,domain mappings, communication, etc.There are two main goals of this work:1. Explore the DG method as a means of solving theGR and hydrodynamics equations simultaneously.As we will see below, the equations of the two theo-ries take fundamentally different forms (conservativevs nonconservative), so it is not a priori obviousthat solving them on the same grid with the sametechnique will work.2. Explore the use of conforming grids for BH andNS applications. In these grids, cubical elementsare mapped to match the spherical geometry of anexcision boundary inside the BH or the sphericalboundary at large distances from the BH or NS. This paper is organized as follows. We first summarizethe formulation of our DG method in Sec. II. We discussour use of geometrically adapted grids, “manual” mesh re-finement, and limiters in Sec. III. We detail the GR-hydroequations and associated algorithms of our numerical im-plementation in Sec. IV. To validate our code, we performstandard test cases; we show these in Sec. V. We presentour results — NS evolutions using the DG method — inSec. VI, before concluding in Sec. VII.
II. DISCONTINUOUS GALERKINFORMULATION
Our code uses a DG method to solve conservation lawsin curved spacetimes and also to evolve the spacetimeitself. We express the spacetime metric g µν using thestandard 3+1 form ds = g µν dx µ dx ν = − α dt + γ ab ( dx a + β a dt )( dx b + β b dt ) , (1)where α is the lapse function, β a is the shift vector, and γ ab is the spatial metric (with determinant γ ) on hyper-surfaces of constant time t . Our index convention is asfollows. Greek indices ( µ , ν ,...) refer to spacetime com-ponents and range from 0 to d in d spatial dimensions.Latin indices ( a , b ,...) refer to spatial components andrange from 1 to d . Repeated indices are summed over.We denote by x the spatial point with coordinates x a .We use units where G, c = 1 . We additionally set M (cid:12) = 1 for the NS simulations in Sec. VI.A conservation law in this curved spacetime can bewritten as a 4-divergence ∇ µ F µ = s , where ∇ µ is thecovariant derivative, F µ encodes the conserved quantity u = F and its corresponding spatial flux vector F a ( u ) ,and s is the source term for u . Separating the time andspatial components gives the more common form √ γ ∂ t ( √ γu ) + 1 √ γ ∂ a ( √ γF a ) = s, (2)which we aim to solve for √ γu ( x , t ) given initial condi-tions √ γu ( x , and suitable boundary conditions. Whensolving a system of conservation laws (e.g., for mass, en-ergy, and momentum in hydrodynamics), u is a vector ofseveral conserved quantities, and F a is a vector of fluxvectors.We numerically solve the conservation law using astrong-form, nodal DG method on square/cube elements.In this section, we summarize the method and give thespecifics of our implementation. We follow the formulationgiven by Teukolsky in Ref. [30], in which greater detailmay be found. The conservation law is discretized (see Sec. II B) and solved fora numerical approximation to the true solution u . We do notmake the distinction between the approximate and true solutions. A. Representing the solution
We divide the spatial domain into K elements. On eachelement, we expand the quantities u , F a , s , etc., over aset of polynomial basis functions φ i , e.g., u ( x , t ) = (cid:88) i u i ( t ) φ i ( x ) . (3)We adopt a nodal representation: we evolve the values u i ( t ) = u ( x i , t ) at the nodes x i of the computational grid,and the φ i interpolate between these grid nodes. Below,we define these quantities; more detailed discussion canbe found in textbooks [39, 40].The partition into elements is chosen so that each el-ement is a mapping of a topologically simple referenceelement: a cube (in three dimensions), square (in twodimensions), or interval (in one dimension). The mappingfrom the reference element coordinates ¯ x to the com-putational coordinates x = x (¯ x ) of each element has aJacobian matrix J = ∂x a ∂x ¯ a (4)and Jacobian J = det J .In each direction, the x ¯ a coordinate spans the inter-val [ − , , and on this interval, we place the nodes x ¯ ai of a Gauss-Legendre-Lobatto quadrature. The one-dimensional (1D) Lagrange interpolation polynomials (cid:96) j ( x ¯ a ) are defined on these nodes and satisfy (cid:96) j ( x ¯ ai ) = δ ij .In the full d dimensions, we construct a tensor-productgrid — we obtain the grid nodes ¯ x i from the direct prod-uct of the x ¯ ai and the basis functions φ i from the productof the (cid:96) i ( x ¯ a ) , e.g., (with some abuse of indices to indicatethe tensor product) φ i (¯ x ) → φ ijk (¯ x ) = (cid:96) i ( x ¯1 ) (cid:96) j ( x ¯2 ) (cid:96) k ( x ¯3 ) . (5)With N p nodes in the x ¯ a coordinate, (cid:96) i ( x ¯ a ) is a polynomialof degree N = N p − . When N is the same in alldirections, we say we have an N th -order DG element.We will occasionally use a modal representation inwhich the solution is expanded over a basis of orthonormalpolynomials, e.g., u (¯ x , t ) = (cid:88) i ˆ u i ( t ) ψ i (¯ x ) . (6)The ˆ u i are the expansion weights, and the ψ i are obtainedfrom the tensor product of orthonormal 1D basis functions,the Legendre polynomials P l . The Vandermonde matrix V ij = P j ( x i ) gives the transformation between the nodaland modal representations, u i = (cid:88) j V ij ˆ u j . (7) B. DG for conservation laws
We impose the conservation law (2) in a Galerkin sense,by integrating the equation against each basis function φ i on each element. We integrate over proper volume √ γd x , giving (cid:90) [ ∂ t ( √ γu ) + ∂ a ( √ γF a ) − √ γs ] φ i ( x ) d x = 0 . (8)To establish the flow of information between neighboringelements, we integrate the flux divergence term by parts,and apply Gauss’s law to the resulting boundary term(see Ref. [30]), (cid:90) ∂ a ( √ γF a ) φ i ( x ) d x = − (cid:90) √ γF a ∂ a φ i ( x ) d x + (cid:73) F a n a φ i ( x ) d Σ . (9)Here, d Σ is the proper surface element on the element’sboundary, and n a is the outward-directed unit normal.The flux vector F a is double valued on the boundarybecause of the local (i.e., discontinuous) nature of thesolution. However, for the scheme to be conservative, aunique flux must cross the boundary between two adjacentelements — this is the so-called numerical flux F a ∗ . Thenumerical flux is computed from the data on both sidesof the boundary and so requires the communication ofboundary data between nearest-neighbor elements. Wesubstitute F a → F a ∗ in the last term of (9).We now undo the integration by parts, using (9) toeliminate the second (i.e., ∂ a φ i ) term (this time, however,we do not substitute in the numerical flux) and obtain (cid:90) ∂ a ( √ γF a ) φ i ( x ) d x → (cid:90) ∂ a ( √ γF a ) φ i ( x ) d x + (cid:73) ( F a ∗ − F a ) n a φ i ( x ) d Σ . (10)The surface integral term provides a boundary conditionon the element and serves to connect the solution betweenneighboring elements of the domain. Defining F = ( F a ∗ − F a ) n a and putting the terms back together, we get theDG equation in integral form, (cid:90) [ ∂ t ( √ γu ) + ∂ a ( √ γF a ) − √ γs ] φ i ( x ) d x = − (cid:73) F φ i ( x ) d Σ . (11)To obtain a form more suitable for computation, we firstexpand each term of (11) using the nodal expansion (3).We rewrite the integrals in the reference coordinates ¯ x ,where d x → Jd ¯ x and d Σ → (cid:112) (2) γd ¯ x , with (2) γ thedeterminant of the two-dimensional (2D) metric inducedby γ ab on the surface. Finally, we evaluate the integralswith a Gauss-Lobatto quadrature rule. By using the gridnodes ¯ x i as the quadrature nodes, we can use the identity (cid:96) i ( x ¯1 j ) = δ ij to greatly simplify the scheme. The tradeoffis that the quadrature rule will not be exact — especiallywhen a nontrivial Jacobian J multiplies the integrand— and this can lead to aliasing and introduce numericalinstabilities that require filtering. Finally, after simplifying the geometric factors on theboundary terms (see Ref. [30], Appendix A) and dividingthrough by common factors, we arrive at the evolutionequation, d ( √ γu ) ijk dt + (cid:104) ∂x ¯1 ∂x a (cid:12)(cid:12)(cid:12) ijk (cid:88) l D ¯1 il ( √ γF a ) ljk + ∂x ¯2 ∂x a (cid:12)(cid:12)(cid:12) ijk (cid:88) m D ¯2 jm ( √ γF a ) imk + ∂x ¯3 ∂x a (cid:12)(cid:12)(cid:12) ijk (cid:88) n D ¯3 kn ( √ γF a ) ijn (cid:105) − ( √ γs ) ijk = − (cid:113) γ ¯1¯1 Njk w N ( √ γF ) Njk δ iN − (cid:113) γ ¯2¯2 iNk w N ( √ γF ) iNk δ jN − (cid:113) γ ¯3¯3 ijN w N ( √ γF ) ijN δ kN + (cid:113) γ ¯1¯10 jk w ( √ γF ) jk δ i + (cid:113) γ ¯2¯2 i k w ( √ γF ) i k δ j + (cid:113) γ ¯3¯3 ij w ( √ γF ) ij δ k . (12)Here, D ¯1 il is the differentiation matrix along the x ¯1 direc-tion, given by D ¯1 il = ∂ ¯1 (cid:96) l (cid:0) x ¯1 (cid:1)(cid:12)(cid:12) i . (13)Although our derivation and resulting evolution equa-tion (12) are given for the 3D case, restricting to a lower-dimensional problem is straightforward. For instance, ina 2D problem, the third tensor-product index on eachterm is dropped (e.g., u ijk → u ij ), as are the ¯3 terms ofthe flux derivative and flux boundary terms. C. DG for the Einstein equations
We use a formulation of the Einstein equations, detailedin the next section, that cannot be written in conservativeform. These equations are instead in hyperbolic form, ∂ t u + A a ∂ a u = s, (14)where the matrices A a and the vector s may be functionsof u , but not of derivatives of u . To obtain the corre-sponding DG algorithm, we again multiply by a basisfunction φ i and integrate over the proper volume element.We integrate by parts twice, substituting the numericalflux after the first integration, to obtain the integral formakin to (11), (cid:90) [ ∂ t u + A a ∂ a u − s ] φ i ( x ) √ γd x = − (cid:73) [( A a u ) ∗ − ( A a u )] n a φ i ( x ) d Σ . (15) Evaluating the integrals as before, we find du ijk dt + A aijk (cid:104) ∂x ¯1 ∂x a (cid:12)(cid:12)(cid:12) ijk (cid:88) l D ¯1 il u ljk + ... (cid:105) − s ijk = − (cid:113) γ ¯1¯1 Njk w N (cid:0) [( A a u ) ∗ − ( A a u )] n a (cid:1) Njk δ iN + ... . (16)This result is analogous to (12), so we have reproducedhere only one term of each type.
1. Comparison with
SpEC ’s penalty spectral algorithm
SpEC solves the Einstein equations using a multido-main penalty pseudospectral method (see, e.g., Ref. [41]).This method is closely related to our nodal DG method:the DG boundary term represents a particular type ofpenalty term, one chosen to enforce conservation via thenumerical flux. Indeed, the spectral method in
SpEC takes the form of (16) with an upwind flux, differing onlyin the numerical prefactor multiplying the boundary fluxterm. Where our DG method has a prefactor of /w N arising from the Legendre Gauss-Lobatto quadrature rule,the SpEC penalty method instead uses the prefactor de-rived for stability of a Chebyshev penalty method [42]. Innumerical experiments (not reported in this paper), weobserve a higher order of convergence under h -refinementfrom the DG method (order N + 1 ) than from SpEC ’sspectral method (order N ). III. APPROACH TO GRID STRUCTURE,MESH REFINEMENT, AND LIMITING
Early applications of the DG method to problems in as-trophysics have used uniform grids. We adopt a differentphilosophy and take advantage of the DG method’s geo-metric flexibility to tailor our grid to the problem beingsolved. This approach was also taken by Refs. [25, 27],which use a 2D wedge-shaped domain when evolving gasflows in a BH spacetime. We discuss here our choice ofgrid structures, mesh refinement, and limiting.
A. Grid structure and mesh refinement
It is well known that constructing the computationalgrid to mirror the underlying symmetries of the problemcan greatly increase the accuracy of a numerical method.In astrophysical problems, the symmetry is often spher-ical, reflecting the gravitational potential of a star orBH. The use of a conforming spherical grid comes with aloss of generality: the grid must remain centered on theastrophysical body. This is especially important whentaking advantage of the spherical grid to excise the sin-gularity inside a BH. With the use of moving grids [35]and control systems [37], however, conforming grids canbe successfully used in simulations of binary mergers.The evolutions shown in this paper make use of twobasic types of grid structures:1. Cartesian grids, obtained by a straightforward affinemapping (a translation and a scaling) of the ref-erence element. These grids are used in severalstandard test problems.2. Cubed-sphere grids, obtained by conforming severalcube-like elements to the surface of a sphere, usingmappings detailed in Appendix A and illustratedin, e.g., Fig. 3. These grids are used for problemswith spherical geometry such as single BH or NSevolutions. The cubed-sphere grid may cover ahollow spherical shell, allowing for excision of thespacetime region inside the BH’s event horizon, ora filled ball, for evolution of the full NS. As weconsider isolated systems at rest, moving grids arenot needed.To further take advantage of the geometric flexibilityof the DG method, we use hp -adaptivity to vary the spa-tial resolution across the simulation domain. The AMRinfrastructure of SpEC is designed to operate under arestricted set of conditions and is not general enough tohandle the shocks and surfaces encountered in the hydro-dynamics evolutions. We instead manually set up fixedmesh refinement, where we initially assign the size andorder of the DG elements based on a priori knowledge ofthe solution. When constructing the grid for a NS evo-lution, for instance, we use larger, higher-order elementsinside the star and smaller, lower-order elements at thesurface. We use “higher-order” (“lower-order”) as a quali-tative description of a DG element, typically referring toelements with N (cid:38) ( N (cid:46) ).The SpEC framework, designed and optimized for evo-lutions on O (10 – spectral elements, scales poorly to the large number of elements often used in DG simulations.In spite of several improvements to the data structures,we find that the code’s memory usage and parallelismbecome inefficient when the domain approaches O (10 ) elements. We therefore stay below this threshold in mostof the tests presented. This restriction on the maximumnumber of elements would be problematic for a typicalDG implementation, in which the domain is split intoa regular grid of many small cubical elements. As weinstead conform our grids to the problem geometry, weobtain satisfactory accuracy using many fewer elements. B. Limiting
In DG elements containing a shock or surface in thefluid, the solution is susceptible to spurious oscillations(Gibbs phenomenon) and overshoots. If unaddressed,these overshoots can lead to unphysical fluid states (e.g.,negative densities) in which the fluid equations are nolonger solvable. A limiter controls these oscillations andovershoots by modifying the solution in a way that isconservative and — ideally — does not overly degradethe accuracy of the method.Typical DG implementations apply the limiter agnosti-cally across the uniform grid. A “troubled-cell” detectoridentifies cells containing spurious oscillations and appliesthe limiter to those cells. While this is the most generalway to set up the problem, finding a general troubled-celldetector that does not misidentify smooth extrema in thesolution can be challenging. This can lead to problems,such as a smearing out of the density maximum at thecenter of a star.In the context of an hp -adaptive DG method, however,the AMR criteria can also be used to inform the troubled-cell detector. When the solution is not smooth (i.e., themodal coefficients do not fall off rapidly enough), theAMR algorithm will reduce the order N of the elementand trigger h -refinement. High-order elements, then, havesmooth solutions and do not require limiting. In our man-ually refined grid, we apply the limiter only to elementswith N ≤ in any spatial direction.While our choices of grid setup and limiter applicationare not fully general, they are representative of the out-come from a more general AMR DG code. Our resultsare an exploration and will serve to inform the choicesmade in a future AMR update to SpECTRE (the newDG code mentioned in Sec. I).
IV. EVOLUTION OF GR-HYDROA. Spacetime geometry
1. Generalized harmonic equations
We evolve the spacetime geometry using the generalizedharmonic formulation of Einstein’s equations [31–33]. Weuse a first-order representation of the system [43] in whichthe evolved variables are the spacetime metric g µν , itsspatial first derivatives Φ iµν = ∂ i g µν , and its first deriva-tive Π µν = − t σ ∂ σ g µν along the (timelike, future-directed)normal t σ to the constant- t hypersurface. The completeequations for ∂ t g µν , ∂ t Φ iµν , and ∂ t Π µν in a vacuumspacetime can be found in Ref. [43]; when coupling thespacetime to matter, we add the source term ∂ t Π µν = (cid:18) vacuumterms (cid:19) − α (cid:18) T µν − g µν T ρσ g ρσ (cid:19) . (17)The DG method for this system of equations takes theform (16). The characteristic variables and speeds for thesystem, used in the upwind numerical flux shown below,are also given in Ref. [43].For the cases we present in this paper, the natural co-ordinates of the initial data are well suited to prolongedtime evolution. The generalized harmonic gauge function H σ , which specifies the coordinates, is therefore indepen-dent of time. Its precise form will depend on the databeing evolved. The constraint-damping parameters γ and γ , which constrain the evolution of the coordinatesand the growth of short-wavelength perturbations, respec-tively, are also problem dependent. Following Ref. [43],we fix the parameter γ to − because this makes thegeneralized harmonic system linearly degenerate.
2. Upwind flux
As the solutions to the Einstein equations are smooth,we use an upwind numerical flux, which sets the fluxthrough the boundary according to the propagation di-rection of each characteristic variable. The characteristicdecomposition of the system is given by A a n a u = S Λ S − u, (18)where S diagonalizes the product A a n a ; i.e., the i th col-umn of S is the right eigenvector of A a n a , with eigenvalue λ i . Physically, the S − u are the characteristic variables ofthe system, and λ i are the associated propagation speedswith respect to the normal n a . The diagonal matrix Λ = diag ( λ , ..., λ n ) holds these eigenvalues and can beseparated by the sign of the eigenvalues, Λ = Λ + +Λ − . Ata boundary with two edge states u L and u R and a normal n a directed toward the R state, the upwind numericalflux takes the form ( A a n a u ) upwind = S (cid:0) Λ + S − u L + Λ − S − u R (cid:1) , (19) Where we use g µν , the cited papers use ψ µν to denote the space-time metric. At each point, we treat the background spacetime (i.e., A a n a )as constant and compute the wave decomposition of the statevector u by treating it as a perturbation. so that characteristic variables propagating left to right(in the direction of n a , with λ i > ) are set from the u L state, whereas variables propagating right to left (with λ i < ) are set from u R . B. Hydrodynamics
1. Relativistic fluid equations
We treat the matter as a perfect fluid. Its stress-energytensor takes the form T µν = ρhu µ u ν + pg µν , (20)where ρ is the fluid’s rest-frame mass density, p is thepressure, and h = 1 + (cid:15) + p/ρ is the relativistic specificenthalpy, with (cid:15) the specific internal energy density. Fromthe fluid’s 4-velocity u µ = W (1 , v i ) , we define the lower3-velocity components v i = γ ij v j and the Lorentz factor W = αu = 1 / √ − v i v i . An equation of state (EOS)relates p , ρ , and (cid:15) ; we use an ideal-gas EOS p = (Γ − ρ(cid:15) ,with Γ the adiabatic index. In the absence of shocks, thisis equivalent to a polytropic EOS where p = κρ Γ , with κ some constant.The dynamics of the fluid are governed by the rela-tivistic Euler equations. We use the València form ofthese equations [1], with conserved quantities { D, S i , τ } :the mass-energy density, momentum density, and internalenergy, as measured by a generalized Eulerian observer.These are given by √ γu = ˜ D ˜ S i ˜ τ = √ γW ρ √ γW ρhv i √ γ (cid:0) W ρh − p − W ρ (cid:1) . (21)We follow the convention of using tildes to indicate “den-sitized” variables, ˜ X ≡ √ γX . The corresponding fluxvector and source term are √ γF a = ˜ Dv a tr ˜ S i v a tr + √ γαpδ ai ˜ τ v a tr + √ γαpv a (22) √ γs = α/
2) ˜ S lm ∂ i γ lm + ˜ S k ∂ i β k − ˜ E∂ i αα ˜ S lm K lm − ˜ S l ∂ l α . (23)Here, v a tr = αv a − β a = u a /u is the transport velocityrelative to the coordinates; S lm and E are componentsof the stress energy, ˜ S lm = √ γT lm = √ γρhW v l v m + √ γpγ lm (24) ˜ E = √ γn µ n ν T µν = √ γρhW − √ γp ; (25)and K lm is the usual extrinsic curvature of the constant- t hypersurface. The system of equations is evolved ac-cording to the discretized form (12), with the densitizedconserved variables { ˜ D, ˜ S i , ˜ τ } serving as the primary vari-ables in the code. The characteristic speeds, used in thenumerical fluxes shown below, are given in Ref. [44].Solving for the primitive variables { ρ, v i , (cid:15) } from { D, S i , τ } requires root finding and may additionally re-quire “atmosphere fixing” in regions of low density wherethe inversion may be numerically poorly behaved. Wefollow the inversion and fixing procedure of Ref. [45], Ap-pendix C. This fixing procedure takes grid points wherethe low-density state { D, S i , τ } does not correspond toa physical state { ρ, v i , (cid:15) } and alters the conserved vari-ables to recover a physical state. Additionally, a small(i.e., dynamically negligible) floor ρ atmo is set on the fluiddensity, ensuring that round-off level errors are controlled.In the test problems of Sec. V, fixing is not needed; weset ρ atmo to 0. For the NS evolutions of Sec. VI, fixing isnecessary outside the star; we give the parameters of thefixing within that section.
2. Numerical fluxes
For the fluid, we use a numerical flux chosen to ap-proximately solve the Riemann problem corresponding tothe discontinuity between elements. As before, we labelthe two states at the boundary as u L and u R , and thenormal n a points toward the R state. A popular choice ofnumerical flux, because of its robustness and simplicity, isthe local Lax-Friedrichs (LLF) flux. This flux is computedaccording to ( F a ∗ n a ) LLF = F a ( u L ) n a + F a ( u R ) n a − C (cid:0) u R − u L (cid:1) , (26)where C = max( | λ i ( u L ) | , | λ i ( u R ) | ) is the largest speedacross the interface. The speeds λ i are again the eigenval-ues of the flux Jacobian (see the upwind flux discussionabove, with A a → ∂F a /∂u ). We maximize over the λ i on both sides of the interface, but independently at eachinterface grid point.A more sophisticated numerical flux, which includes anapproximate treatment of the system’s underlying wavestructure, is given by Harten, Lax, and van Leer (HLL)[46, 47], ( F a ∗ n a ) HLL = c max F a ( u L ) n a + c min F a ( u R ) n a c max − c min − c max c min c max − c min (cid:0) u R − u L (cid:1) . (27)Here, c min and c max are estimates for the fastest left-and right-moving signal speeds, respectively. We use thesimple estimates [48], computed pointwise, c min = min (cid:0) λ i ( u L ) , λ i ( u R ) , (cid:1) c max = max (cid:0) λ i ( u L ) , λ i ( u R ) , (cid:1) . (28)Note that the HLL flux reduces to upwinding when all λ i share the same sign, i.e., all characteristic variables arepropagating in the same direction. We find that the LLF and HLL fluxes give very similarresults in most of the cases we tested (for an exception,see the supersonic accretion flow test in Sec. V B) andconclude that the use of an approximate solution to theRiemann problem does not introduce a significant errorin these problems. The results presented in this paperare computed using the HLL flux.
3. Limiters
In this work, we use and compare two limiters. Thefirst is the simple, but also low-order, ΛΠ slope limiter[39, 49], which we will refer to simply as minmod. Thislimiter computes several estimates for the slope of thesolution on each element, and then, in elements wherethese estimates indicate the presence of oscillations in thesolution, it acts to reduce the slope. Taking the 1D caseas example, we write the solution u k on the k th elementas a series expansion, u k = ¯ u k + u ( x − x ) + O ( x − x ) , (29)where ¯ u k is the element-averaged mean of u k , u is themean slope, and x is the center of the element. Theminmod limiter’s slope estimates are a = u , a = ¯ u k +1 − ¯ u k h/ , a = ¯ u k − ¯ u k − h/ , (30)where h is the width of the element. The limiter selectsthe estimate with the smallest absolute value (or 0, if thethree estimates differ in sign). If the selected estimate isnot the original slope u , the limiter activates by reducingthe slope u to the selected estimate (or 0) and discardingany higher-order terms in the approximation. On elementswith order N > , we use the ΛΠ N generalization of thelimiter described in Ref. [39]. We do not use the “totalvariation bound” generalization, which sets a scale belowwhich oscillations are tolerated, since we find that it isnot robust at star surfaces.In higher dimensions, the 1D limiter is applied to eachdirection in turn. After this process, the limited solutionmay occasionally correspond to a nonphysical state. Whenthis occurs, we further reduce the slope until the followingare satisfied throughout the element: min ( D ) > ρ atmo ,min ( τ ) > , and S < τ ( τ + 2 D ) .For evolutions on deformed grids, we apply the 1Dlimiter along each direction of the reference ¯ x coordinates.This choice of coordinates leads to a straightforward com-putation of the minmod slope estimates, because theseries representation of u k takes a simple form, and theelement’s upper and lower neighbors in each direction arewell defined. However, the choice introduces a (relativelysmall) violation of conservation: the limiter will conservethe means ¯ u k with respect to the reference ¯ x coordinates,but the means with respect to the “global” x coordinateswill in general be modified after the limiter activates. Thiscan be understood by noting that the means in the twocoordinate systems are differently sensitive to the shapeof the function u k : ¯ u k | ¯ x = (cid:82) u k d ¯ x (cid:82) d ¯ x vs ¯ u k | x = (cid:82) u k d x (cid:82) d x = (cid:82) u k Jd ¯ x (cid:82) Jd ¯ x . (31)We explored two simple corrections to the minmod lim-iter that restore conservation in the x coordinates. Thefirst correction limits the Jacobian-weighted solution Ju k instead of u k ; the second shifts the postlimiting solution u k → u k + δu k , with δu k a constant computed to re-store the prelimiting mean ¯ u k . Both of these correctionssuccessfully restore the limiter’s conservative properties,but we found that they also introduced long-timescaleinstabilities at the surface of the star — we note thatRadice and Rezzolla [21] also found poor behavior whenusing similar corrections with the simple minmod limiter.Consequently, we do not use these corrections in our simu-lations. Instead, we will quantify the error in maintainingconservation when presenting our results.The second limiter we consider is that of Moe et al .[50], henceforth MRS. This limiter acts by scaling theconserved variables u about their means ¯ u , u → ¯ u + θ ( u − ¯ u ) , (32)with θ ∈ [0 , determined from analysis of the minima andmaxima of the solution in the immediate neighborhoodof the element. A tolerance function α ( h ) sets the scalebelow which oscillations are tolerated; we use the function α ( h ) = 100 h / for the cases presented in this paper, asit performs well on many different test problems.We obtain best results when computing θ from theprimitive variables, as MRS recommend. However, caremust be taken when computing the primitive variables,as the fluid state may be unphysical until limited. We“prelimit” by applying an additional scaling of the form(32) to the conserved variables. The steps below restore aphysical state and ensure the inversion procedure is wellposed:1. If min( D ) < ρ atmo or min( τ ) < , scale to fix theseviolations.2. If S i S i > τ ( τ + 2 D ) at any grid point, scale to fixthis violation. This requires solving a quadraticequation for θ .3. If the inversion to primitive variables encountersany of the errors outlined in Ref. [45], Appendix C(this is rare), scale again with θ = 1 / .This procedure is conservative by construction, and wefind it to be robust. After this prelimiting step, we com-pute the primitive variables and limit according to theMRS prescription. We handle deformed grids as for theminmod limiter, by computing the means in the reference ¯ x coordinates and incurring some error due to loss ofconservation. As with minmod, attempts to reformulatethe limiter to restore conservation (we tried the simple approach of computing the MRS means directly withrespect to the x coordinates as well as the same two re-formulations described for minmod) were not stable atthe star surface.We apply the limiter to the fluid variables at the endof each time-stepper substep. As described in Sec. III B,we may not apply the limiter to every element, choosinginstead to mimic an AMR scheme in which high-orderelements are known to be smooth. The use of morecomplex, higher-order, limiters, e.g., subcell methods [21,23] or the compact-stencil WENO [51] and HWENO [52]limiters, will be the subject of future investigation. C. Combined GR-hydro system
For self-consistent NS evolutions, the equations of thegeneralized harmonic and relativistic Euler systems areeach treated as described above and are evolved in par-allel. The two systems couple via their respective sourceterms and the geometry terms in the hydrodynamics flux F a ( u ) . We compute the characteristic speeds indepen-dently for each system, leaving out the cross-couplingarising from the off-diagonal ∂F a hydro /∂u GR flux Jacobianterms. When the fluid variables require limiting, the lim-iter is applied to the fluid variables only, and the spacetimevariables are left unmodified. D. Filtering
The use of inexact quadratures to obtain an efficientDG scheme may result in numerical instabilities causedby aliasing. Where these numerical instabilities exist,we address them by filtering the higher modes in thesolution’s modal representation. We use an exponentialfilter, e.g., in one dimension, ˆ u i → F ( i )ˆ u i = exp { ( − α ( i/N ) s ) } ˆ u i , (33)where α controls the strength of the filter’s effect and s isan even integer controlling how many modes are affected.In d > dimensions, we take advantage of the tensor-product basis to apply the filter dimension by dimension;this gives d exponentials. On deformed grids, we filterthe Jacobian-weighted solution Ju and then divide by J ,so that the operation remains conservative. We applythe filter at the end of each complete time step to thecomponents of u and on the elements that show numericalinstability. E. Time stepping
We use the third-order strong stability-preservingRunge-Kutta scheme of Ref. [53] for the time integra-tion. Given the solution u n at time t n , the solution u n +1 at time t n +1 = t n + ∆ t is computed as u (1) = u n + ∆ tL ( u n ) u (2) = 34 u n + 14 (cid:104) u (1) + ∆ tL ( u (1) ) (cid:105) u n +1 = 13 u n + 23 (cid:104) u (2) + ∆ tL ( u (2) ) (cid:105) . (34)Here, L ( u ) = du/dt is computed from expressions (12)for the fluid variables or (16) for the spacetime variables.In all cases presented, the initial t = 0 data are com-puted by pointwise evaluation of a known state. Thelimiter is applied to the initial data and at the end ofevery subsequent substep. Filtering is done at the end offull time steps. V. CODE TESTS
In this section, we present a selection of benchmarktests that we use to validate our implementation of theDG method within
SpEC .We first show tests of vacuum spacetime evolution.From a family of gauge wave evolutions at varying resolu-tions, we verify that the method converges to the exactsolution at the expected rate. Next, by evolving a Kerr(i.e., isolated and spinning) BH over long timescales, weshow the stability of the algorithm.We then show our tests of the hydrodynamics imple-mentation. We again verify the convergence rate of theerrors, now with a generalized Bondi problem in whichthe fluid undergoes spherically symmetric accretion ontoa Schwarzschild BH. This test verifies the fluid equationsas well as the sourcing of the fluid by the spacetimecurvature. We then show standard shock tests in oneand two dimensions, comparing the effectiveness of theimplemented limiters.In these tests, whenever possible, we compare the nu-merical solution to an exact solution, and we use theirdifference as an error measure. We report a normalizederror err [ X ] in a quantity X , defined aserr [ X ] = (cid:107) X − X exact (cid:107) (cid:14) (cid:107) X exact (cid:107) . (35)Here, (cid:107) X (cid:107) is the L -norm, evaluated pointwise by directsummation over every node of the computational grid, (cid:107) X (cid:107) = 1 N nodes N nodes (cid:88) i =0 X i . (36)When X is a vector or tensor quantity, we compute acomponent-wise norm (cid:107) X (cid:107) = (cid:107) X (cid:107) + (cid:107) X (cid:107) + ... , ratherthan the physical norm X a X a . When X exact = 0 so thatwe cannot define the normalized error, we instead use (cid:107) X (cid:107) as our error measure. A. Spacetime tests
1. Gauge-wave test
The spacetime of the “apples to apples” gauge-wave test[54], obtained via a nonlinear, plane-wave transformationof Minkowski space, takes the form ds = − (1 + a ) dt + (1 + a ) dx + dy + dz , (37)with a = A sin[2 π ( x − t )] . (38)We show results for a wave of amplitude A = 0 . on aunit-cube domain with extents [0 , . As the gauge waveis harmonic, the generalized harmonic gauge function H σ is zero. We set the constraint-damping parameters ( γ , γ , γ ) to (1 , − , , values that give stable evolutionsover long timescales (up to at least t fin = 1000 , or 1000crossing times). For the convergence study, however, wemeasure the error in the spacetime metric g µν at a finaltime t fin = 10 , after evolution with time steps of size ∆ t =10 − . This time step corresponds to ∆ t/ ∆ x min (cid:39) . for the highest-resolution case in the convergence study( K = 128 , N = 4 ).We show in Fig. 1 the convergence under h -refinement,measured for elements of order N = 2 , , . For a baseresolution, we partition the unit-cube domain into elements along the x direction; we h refine by furthersplitting each element along x , reducing the element’swidth h in half each time. We do not split in y or z —the anisotropic refinement is chosen to match the x -onlydependence in the problem. For each order N of the DGmethod, we compare our measurements to the theoreticalscaling of the error (see, for instance, Ref. [39]),err [ g µν ] ≤ Ch N +1 ∝ /K N +1 , (39)for some constant C . We find excellent agreement be-tween the measured and expected convergence rates. Thehighest-resolution case ( K = 128 and N = 4 ) has aslightly larger error, having reached the round-off levelerror in the derivatives of the spacetime.In Fig. 2, we show the convergence under p -refinement,obtained by increasing the order N of the DG methodwhile maintaining the base resolution of 16 elements. Weexpect the errors to decrease exponentially with the order N and recover this behavior in our measurements. Thisresult demonstrates the spectral convergence of the DGmethod for smooth solutions.
2. Kerr black hole
We next evolve the spacetime of a Kerr BH, describedby the Kerr metric in Kerr-Schild coordinates [55]. TheBH has spin (cid:126)a = (0 . , . , . M BH with magnitude a ≈ . M BH , not aligned with any grid symmetries. We useunits where M BH = 1 .0
16 32 64 128
Number of elements K − − − − − − − − − − − e rr [ g µ ν ] − − − N = 2 N = 3 N = 4 FIG. 1. The error in g µν as a function of the number ofelements ( h -refinement) for the gauge wave test of Einstein’sequations. The symbols indicate the measured error norms formethods of order N = 2, 3, 4. The dashed lines, normalizedto the K = 16 data, indicate the expected error scaling forthird-, fourth-, and fifth-order convergence. Polynomial order N − − − − − − − − e rr [ g µ ν ] K = 16 FIG. 2. The error in g µν as a function of the order of approx-imation ( p -refinement) for the gauge wave test of Einstein’sequations. The number of elements is fixed at K = 16 . Thedots indicate the measured errors; the dashed line is a fitdemonstrating the exponential decrease in the error with N . The domain is a hollow spherical shell that excises thesingularity within the BH. In terms of the coordinateradius r , the domain extends from r in = 1 . (just insidethe event horizon) to r out = 32 . At the inner boundary, allthe characteristics of the system are outgoing (i.e., leavingthe domain, toward the singularity), so no boundarycondition needs to be imposed. Physically, no informationenters the simulation from the interior of the BH. At theouter boundary, we impose the analytic solution as a (a) Projection (b) Equatorial cut FIG. 3. The grid structure for the Kerr BH evolution test.Shown are (a) a projected view and (b) an equatorial cut.The black lines show the element boundaries, and the lightgrey lines show the Gauss-Legendre-Lobatto grid within eachelement for order N = 5 . Dirichlet boundary condition. We choose constraintdamping parameters γ = 3 exp (cid:2) − ( r/ / (cid:3) + 0 . (40) γ = − (41) γ = exp (cid:2) − ( r/ / (cid:3) + 0 . . (42)The generalized harmonic gauge function H σ = Γ σ ≡ g µν Γ σµν is the trace of the Christoffel symbols of theKerr-Schild metric; it is constant in time.We set up a cubed-sphere grid on this domain, usingthe mappings from Appendix A. The wedges of the cubedsphere are split radially into five concentric shells locatedbetween the surfaces r = 1.8, 3.2, 5.7, 10, 18, 32, and thentangentially into × angular portions, for a total of 120elements. The tangential coordinates of each wedge aremapped to obtain an equiangular grid, as this is a moreoptimal distribution for the grid points on the sphericalsurface. We show in Fig. 3 two views of this grid: on theleft a projected view showing the angular structure on aconstant-radius surface, and on the right an equatorialcut showing the radial structure. The increasing densityof grid points toward the center of the domain helps toresolve the stronger spacetime curvature near the BH.In Fig. 4, we show the stability of the Kerr BH evolutionby monitoring the simulation errors over a duration of M BH . We carry out the simulation using elements oforder N = 5, 6, and 7; the time-step size is ∆ t = 10 − ,giving ∆ t/ ∆ x min (cid:39) . for the N = 7 case. The figure’stop panel shows the error err [ g µν ] in the spacetime metric,a measure of the solution’s drift from the exact value. Thebottom panel shows the dimensionless norm (cid:107)C(cid:107) of thegeneralized harmonic energy constraint [43], a measure of We do not use the constraint-preserving boundary conditionstypically used in
SpEC simulations, because these are bound-ary conditions on ∂ t u , rather than u , and so would require amodification of the DG formulation. − − − − e rr [ g µ ν ] N = 5 N = 6 N = 70 2 4 6 8 10 t [10 M BH ] − − − k C k FIG. 4. The errors during the Kerr BH evolution test. Thetop panel shows the error in the spacetime metric g µν for threedifferent orders of the DG method. The lower panel showsthe dimensionless norm of the generalized harmonic energyconstraint at the same three orders. how well the numerical solution at each constant- t slicesatisfies Einstein’s equations. After a rapid settling ofthe solution to its numerical equilibrium, we see clearconvergence in the error quantities. We conclude thatthe method is convergent and stable up to at least t =10 M BH and, we presume, forever. B. Relativistic hydrodynamics tests
1. Spherical accretion onto black hole
In the relativistic Bondi problem, an ideal gas accretesradially onto a nonrotating BH. The feedback from thefluid onto the spacetime is ignored: the BH mass is con-stant, and the spacetime is Schwarzschild. We use Kerr-Schild coordinates, and again we set M BH = 1 . The ana-lytic profile for the fluid flow is presented by Michel [56];following Ref. [44], we pick a solution for a Γ = 5 / idealgas with the sonic point and mass accretion rate given by r crit = 200 and ˙ M = 10 − . We measure the error in theconserved relativistic density ˜ D at a final time t fin = 100 ,after evolution with time steps of size ∆ t = 5 × − . Thistime step corresponds to ∆ t/ ∆ x min (cid:39) . for the highest-resolution case in the convergence study ( K = 120 × ,
120 120 × × Number of elements K − − − − − − e rr [ ˜ D ] − − − N = 2 N = 3 N = 4 FIG. 5. The error in the conserved density ˜ D as a functionof the number of elements ( h -refinement) for the sphericalaccretion test. The symbols indicate the measured error normsfor methods of order N = 2, 3, 4. The dashed lines, normalizedto the K = 120 data, indicate the expected error scaling forthird-, fourth-, and fifth-order convergence. N = 4 ).We evolve the fluid in a hollow spherical shell extend-ing from r in = 1 . (just inside the event horizon), to r out = 12 . The sonic point in the accretion flow is locatedoutside this region, so the flow is smooth and supersonicthroughout the simulation domain. In this test problem,we obtain significantly more accurate results when usingthe HLL numerical flux (vs LLF), as the supersonic flowis best represented by the HLL upwinding limit. At theinner boundary, the characteristics of the fluid systemare outgoing (i.e., leaving the domain into the BH), sono boundary condition needs to be applied. At the outerboundary, we impose the analytic solution as a boundarycondition.We use a cubed-sphere grid similar to that of the KerrBH test above. At the base resolution, we divide thedomain into five spherical shells between the surfaceslocated at radii r = 1.8, 2.7, 4, 6, 9, 12, and we spliteach wedge into × angular portions. The tangentialcoordinates are again mapped to obtain an equiangulargrid.We show in Fig. 5 the convergence under h -refinementof this grid, for elements of order N = 2, 3, 4. We h refineby splitting each element into smaller elements; we splitgeometrically in radius according to r split = √ r lower r upper and linearly in the tangential directions. As the elementsare not uniform, this choice of radial split is not unique,but we find it gives reduced error compared to a linearsplit r split = ( r lower + r upper ) / . We again see the errorsconverging at the expected rate.In Fig. 6, we show the convergence under p -refinement.Again, we use the base configuration of elements and2 Polynomial order N − − − − − − e rr [ ˜ D ] K = 120 FIG. 6. The error in the conserved density ˜ D as a functionof the order of approximation ( p -refinement) for the sphericalaccretion test. The number of elements is fixed at K = 120 .The dots indicate the measured errors; the dashed line is a fitdemonstrating the exponential decrease in the error with N . increase the order N of the method from 2 to 7. Weconfirm that for this smooth fluid evolution problem,the errors decrease exponentially with the order of themethod.
2. 1D shock tube test
We perform a standard 1D relativistic shock test prob-lem, in which a high-density and -pressure fluid expandsinto a low-density and -pressure fluid. Following Ref. [4],we take a
Γ = 5 / ideal gas initially split at x = 0 . intoleft and right states characterized by ( ρ, v x , p ) = (cid:40) (10 , , / , x < . , , , x > . . (43)The simulation domain is an interval x ∈ [0 , , which wedivide into K = 160 elements of order N = 2 . We evolvethe shock until a final time t fin = 0 . , with time steps ∆ t = 4 × − ( ∆ t/ ∆ x min = 0 . ).In Fig. 7, we show the profiles of ρ , v x , and p at thefinal state, comparing the minmod and MRS limiters.Both limiters capture the features of the shock profile.The minmod limiter produces a larger overshoot at themain shock front and increased oscillation at the frontend of the rarefaction fan, a known behavior when apply-ing this limiter to the conserved variables (rather thancharacteristic variables [57]).
3. 2D Riemann shock interaction test
We next study a standard 2D Riemann problem inwhich two shocks and two contact discontinuities interact. . . . . . . x pρ × v x minmodMRS FIG. 7. Snapshot of the fluid variables in the shock tubetest. The fluid pressure p , the rest-mass density ρ , and thevelocity v x (scaled up × ) are plotted after evolutions usingthe minmod and MRS limiters. The mean value on eachelement is shown. The exact solution to the problem is givenby Centrella and Wilson [58] and is plotted here in the solidline. As in the 1D shock test, the fluid is a
Γ = 5 / idealgas. The initial conditions for the problem were firstgeneralized from Newtonian to relativistic hydrodynamicsby Del Zanna and Bucciantini [59] and later modified byMignone and Bodo [6] to give a cleaner wave structure.The initial condition divides the computational domain [ − , into four quadrants, each of which holds a constantfluid state, ( ρ, v x , v y , p ) = (0 . , , , , x < , y < . , , . , , x > , y < . , . , , , x < , y > ρ , , , p ) , x > , y > , (44)where the low-density state in the upper-right quadrant isdefined by ρ = 5 . × − and p = 2 . × − .We partition the domain into × elements of order N = 2 , and we evolve until a final time t fin = 0 . withtime steps ∆ t = 10 − ( ∆ t/ ∆ x min = 0 . ).In Fig. 8, we show contour plots of the density ρ atthe final state. We interpolate the evolved ρ onto ahigh-resolution uniform grid on which the contours arecomputed. The results in the top panel are computedwith a minmod limiter, and those in the bottom panelare computed with MRS. We find, qualitatively, excellentagreement between the results from the two limiters; onlythe jet feature (in the lower-left quadrant) shows a cleardifference in resolution, with the MRS limiter producinga cleaner structure.3 − . − . . . . − . − . . . . − . − . . . . − . − . . . . FIG. 8. The density ρ in the 2D Riemann problem. The toppanel is computed with the minmod limiter, and the bottompanel is computed with the MRS limiter. The plots each show30 contour lines, equally spaced in log ρ . VI. NEUTRON STAR EVOLUTIONS
Having verified the convergence and shock-capturingproperties of our code, we now present our main results:evolutions of an isolated, spherical NS using the DGmethod. We first evolve the NS under the Cowling approx-imation, i.e., keeping the background spacetime fixed tothe Tolman-Oppenheimer-Volkoff (TOV) solution. Thisremains a challenging test of the hydrodynamics code’sability to handle the discontinuity at the stellar surface.We then evolve the NS self-consistently using the coupledGR-hydro system.The initial data for the NS fluid and spacetime arefound by integrating the TOV equations [55, 60, 61] forthe mass-energy density ρ E ( R ) ≡ ρ ( R )(1+ (cid:15) ( R )) , enclosedADM mass m ( R ) , and metric potential φ ( R ) in terms of the areal radius R . The spacetime metric is given by ds = − e φ dt + (cid:18) − mR (cid:19) − dR + R d Ω . (45)In computing the TOV solution, we describe the NSmatter by a polytropic EOS. When time evolving thesolution, we return to the corresponding ideal-gas EOS.The results presented throughout this section are for astar with κ = 100 and Γ = 2 . The star has central massdensity ρ c = 1 . × − , giving a stable, nonrotatingTOV solution with ADM mass M NS (cid:39) . M (cid:12) and arealradius R NS (cid:39) . M (cid:12) (cid:39) km. Its radius in the isotropiccoordinates used during evolution is r NS (cid:39) . M (cid:12) . Inthis section, we use units where M (cid:12) = 1 .For the NS evolutions, we use the atmosphere fixingfrom Ref. [45], Appendix C. We set the density cut-off ρ cutatmo = 10 − so as to resolve 12 orders of mag-nitude in density. Where the density falls below thiscutoff, we set the fluid to the “atmosphere” state where ρ = ρ atmo = 10 − , v i = 0 , and (cid:15) = 0 . Elsewhere,we constrain the specific internal energy to the range κρ ≤ (cid:15) ≤ κρ , with κ from the polytrope describing theinitial conditions. These bounds serve to control the fluidentropy in the region around the star surface, by prevent-ing numerical errors from causing an entropy decreaseand allowing heating only within a reasonable range. Tocheck that our results are not influenced by the choice ofthese atmosphere fixing thresholds, we evolved a few com-parison cases in which we increased the densities ρ cutatmo and ρ atmo by a factor of 10. These comparison evolutionsdeviated only slightly from the primary evolutions, con-firming that our atmosphere treatment does not stronglyimpact the neutron star simulations. A. Cowling neutron star in spherical symmetry
We begin with 1D evolutions in spherical symmetry.For these simulations, we rewrite the conservation law (2)and the relativistic Euler equations in terms of sphericalcoordinates { r, θ, φ } . The DG formulation takes a formsimilar to (12) in one dimension, but with a spherical di-vergence ∂ r ( r u r ) /r replacing the Cartesian divergence ∂ x u x . The fluid equations pick up an additional momen-tum source term: s ( S r ) = s ( S x ) + αp ( g rr ∂ r g rr + 2 /r ) .To avoid the coordinate singularity at r = 0 , we set upa symmetric domain on the interval [ − , and use astaggered grid so that no nodes are located at the origin.On this domain, we consider three grids with differentresolutions. The first two, which we name I1 and I2 , havecomparable resolutions to the grids of our 3D simulations.These two grids differ in the order of the DG elementsnear the surface of the star; linear elements are used in I1 vs quadratic elements in I2 . The third grid, I1R , The grid names are structured as follows: the first letter encodes TABLE I. The structure of the spherically symmetric NSgrids I1 , I2 , and I1R . For each grid, the parameters definingthe elements in the interior, right-side surface, and right-sideexterior regions are given; the elements in the left-side surfaceand left-side exterior regions are obtained by symmetry. Theinterior and exterior regions of I2 are identical to those of I1 .Extents K region N region I1 Interior [ − . , .
25 3Surface (right side) [7 . ,
10 1Exterior (right side) [10 , I2 Surface (right side) [7 . , I1R
Interior [ − ,
101 3Surface (right side) [8 ,
20 1Exterior (right side) [9 ,
30 3 has a higher resolution and is more aggressively refinedaround the surface of the star. In all three grids, we dividethe domain into five regions: the interior of the star, thesurface on the left/right, and the exterior on the left/right.We use larger, higher-order elements in the interior andexterior regions and smaller, lower-order elements in theneighborhood of the star’s surface. The number andorder of the elements within each region are listed inTable I. We evolve the system until t = 10 (cid:39) ms. Onthe lower-resolution grids I1 and I2 , we use time steps ∆ t = 0 . corresponding to ∆ t/ ∆ x min = 0 . . On thehigher-resolution grid I1R , we use time steps ∆ t = 0 . with ∆ t/ ∆ x min = 0 . .We now compare evolutions of the spherically symmet-ric NS for different choices of the grid and the limiter —specifically the I1 or I2 grid, and the minmod or MRSlimiter. We plot in Fig. 9 the normalized density errorerr [ ˜ D ] for each case over the duration of the simulation.We first examine the two minmod cases. Here, the datareveal two components in the dynamics: a short-periodoscillatory behavior and a gradual drift as the star set-tles to its numerical equilibrium configuration on muchlonger timescales. The I2 case has a much higher initialerror and increased dissipation, as indicated by the morerapid decay of the oscillatory component. The increasederror and dissipation occur because the minmod limiterlinearizes the solution on the quadratic elements at thesurface, resulting in the loss of information. In the twoMRS cases, we find very different behavior: the densityerror is roughly an order of magnitude larger (vs minmod)and grows over time, and the high-frequency oscillationsare not damped on the timescale of the simulation. Al-though the MRS evolutions are stable (on this timescale),the star does not settle to an equilibrium. the domain’s topology (“ I ”: interval; “ B ”: ball), the integer givesthe (radial) order of approximation of the elements near thesurface of the star, and a final “ R ” indicates a refined, higher-resolution grid. t [10 M (cid:12) ] − − − e rr [ ˜ D ] I1 minmod I2 minmod I1 MRS I2 MRS
FIG. 9. The density errors in the spherically symmetricCowling NS evolution. The four cases correspond to differentchoices of grid ( I1 or I2 ) and limiter (minmod or MRS) atthe star surface. As some of the curves are highly oscillatory,we plot for each case the mean error in a solid line, andthe envelope as a light-colored shaded region. The mean iscomputed by applying a Gaussian smoothing to the data,with half-width σ (cid:39) ; the envelope minimum/maximum arecomputed in bins of width ∆ (cid:39) . To better understand this difference in behavior be-tween the minmod and MRS limiters, we now look at thedistribution of the errors across the star. In Fig. 10, weshow the error in the fluid density vs the stellar radius,at time t = 2000 . The solid lines in this plot show theangle-averaged errors (i.e., the average of the left- andright-side data), and the lighter filled region shows thespread in error values at fixed radius. From this plot, wemake two observations. First, while the minmod limitermaintains excellent symmetry across the star, the MRScase shows a large spread in the error values, indicatinga loss of spherical (i.e., reflection) symmetry. Second,while the density and velocity errors in the minmod caseare largest at the surface of the star, denoted by a verti-cal dotted line in the figure, the fluid remains confinedwithin the true surface of the star. When using MRS, thestar instead extends significantly beyond the true surface;matter with non-negligible densities and large ( v > . )velocities is present out to r (cid:39) . Our interpretationis that the MRS limiter provides insufficient dampingof small-scale fluctuations in the atmosphere near thestar. These slowly grow, leading to the expansion of thestar beyond its true surface and the contamination of the Moe et al. [50] discuss interpolating the solution to a finer grid(to better sample its shape) before computing the maxima andminima used by the limiter. We found no significant improvementin behavior when using this interpolation. r [ M (cid:12) ] − − − − − | ρ − ρ e x a c t | I1 minmod I2 minmod I1 MRS I2 MRS
FIG. 10. The rest-frame density error vs the radius, at time t = 2000 , in the spherically symmetric Cowling NS evolution.The shaded region shows the minimum/maximum errors, ateach radius, from the left and right sides of the symmetricdomain. We plot the mean error in the solid line. The verticaldotted line indicates the location of the TOV star surface at r NS (cid:39) . . solution inside the star.We conclude from these comparisons that the MRSlimiter — although effective at handling shocks — ispoorly suited to the task of controlling a stellar surface onthe conforming grids that we are using. In the remainderof this paper we therefore only show results obtainedwith the minmod limiter. We will also restrict to resultsfrom grids with linear-order elements at the star surface,as these are much less dissipative, given our use of thislow-order limiter.In Sec. V, we showed that our DG implementation hasthe expected convergence properties for test problemswith smooth solutions. For the NS problem, the expectedconvergence behavior is less clear; the convergence ratewill be degraded by the discontinuity at the stellar surface,and, furthermore, our use of a geometrically adapted gridwith elements of different sizes and different orders makesit more difficult to quantify the resolution. In spite ofthis, we examine the convergence behavior with a seriesof short evolutions in which we successively h refine the I1 grid, and the time step, by factors of 2 and 4 (in theinterior region of the grid, we refine from 25 to 49, then97, elements; this maintains a grid that straddles theorigin). From these evolutions, we measure the error in ˜ D at time t = 1000 ; we expect the error decrease to liebetween a fourth-order convergence (corresponding to thecase where the error is dominated by the N = 3 interior)and first-order convergence (corresponding to the casewhere the error is dominated by the discontinuity at thesurface). We find (but do not show) that as the grid isrefined, the error decrease is consistent with a third-order t [ M (cid:12) ] r [ M (cid:12) ] C o n v . o r d e r FIG. 11. The local convergence order in the sphericallysymmetric Cowling NS evolution. The convergence order ismeasured from the decrease in the element-average of err[ ˜ D ]between simulations on two computational grids: the I1 gridand a denser grid obtained by refining each I1 element into4. Lighter colors correspond to a higher order of convergence,i.e., a more rapid decrease in error. convergence.To investigate the degradation in accuracy caused bythe stellar surface, we plot in Fig. 11 the spatial andtemporal variation of a locally defined convergence or-der. At each point on the plot, the convergence orderis computed from Eq. (39) by measuring the decrease inthe element-averaged error in ˜ D between simulations ontwo computational grids: the I1 grid and the denser gridobtained by a four times refinement of I1 . The interiorof the star initially shows the expected fourth-order con-vergence, but during the first 50 M (cid:12) the order quicklydecreases to roughly second order as lower-accuracy datafrom the stellar surface propagate into the interior. Wefind, in agreement with Refs. [23, 62], that as the star set-tles to its numerical equilibrium, the order of convergenceincreases again to roughly third-order convergence. Thisshows that the degradation in convergence caused by thestellar surface is limited, and so the high-order qualitiesof the DG method, to a large degree, continue to apply.We now take a second, closer look at the spherically sym-metric NS evolution. In Fig. 12, we compare evolutions ofthe spherically symmetric NS on the I1 and I1R grids. Weshow, in the top two panels, the errors in the conservedquantities ˜ D and ˜ S r over the first M (cid:12) (cid:39) ms ofevolution time. The errors are lower by one or two ordersof magnitude in the I1R case. In the bottom panel, weplot the time dependence of the central density ρ c as afractional error with respect to the initial central density ρ c, . We see in ρ c a qualitative difference between thetwo evolutions, with the I1R case showing a clearly pe-riodic structure corresponding to the crossing time forperturbations seeded at the surface of the star. In thefull evolution to t = 10 , not shown in the figure, thehigh-resolution evolution maintains its equilibrium, withthe remaining oscillations in ˜ S r and ρ c slowly decaying.In the low-resolution evolution, the density error err [ ˜ D ] asymptotes to roughly − , the oscillations in (cid:107) ˜ S (cid:107) slowlydecay, and the central density continues to slowly drop,6 − − − − e rr [ ˜ D ] I1I1R − − − − − − k ˜ S k t [ M (cid:12) ] − . − . − . − . − . . . . . ( ρ c − ρ c , ) / ρ c , / − − . . . FIG. 12. The errors in the spherically symmetric CowlingNS evolution. The top (middle) panel shows the error in theconserved density ˜ D (conserved momentum ˜ S i ) for evolutionsusing the minmod limiter on the grids I1 and I1R . The bot-tom panel shows the evolution of the central density ρ c as afractional error with respect to its initial value ρ c, . The insetin the bottom panel zooms in to better show the initial ρ c evolution in the I1R case; the I1 curve is omitted from theinset for visual clarity. reaching a . deficit at t = 10 .While the DG method is fundamentally conservative,we have discussed (in Sec. IV B 3) how our use of limiterson deformed elements can violate this property. The spher-ically symmetric simulations are also affected, even thoughthe 1D elements are themselves undeformed, because thelimiter does not account for the spherical Jacobian πr that takes the 1D volume element to the spherical volumeelement. The corrections to restore conservation exploredin Sec. IV B 3 are no more effective in the 1D case; con-servation is restored at the expense of stability near thestar surface. We quantify the conservation error by track-ing the NS’s baryon mass M b ≡ (cid:82) ˜ D πr dr during theevolution, as this should be a conserved quantity. Wefind (but do not show) that M b slowly grows. On the I1 grid, the relative error in M b (with respect to its initialvalue) reaches roughly − at t = 10 , with over 90% ofthis growth occurring over the initial M (cid:12) as the starsettles toward equilibrium. On the I1R grid, the errorgrows to about × − , again mostly over the initialportion of the evolution.We now reconsider Fig. 12, and focus on the oscillatorybehavior seen in the different quantities. These oscilla-tions are triggered by errors from two sources: truncationerrors from the evaluation of the exact TOV solutionon the finite-resolution numerical grid and the action ofthe limiter which modifies the initial solution near thestar’s surface. These errors seed perturbations of thevarious radial eigenmodes of the star, each of which sub-sequently resonates with its corresponding eigenfrequency.A common test of NS evolution codes is to compare thefrequency spectrum of the simulated star against theeigenfrequencies obtained from linear theory.To make this comparison, we compute the frequencyspectrum from the central rest-mass density during thefirst M (cid:12) of evolution time. After subtracting theinitial density offset ρ c, , we apply a Hanning window tothe time interval and compute the discrete Fourier trans-form. We plot in Fig. 13 the absolute value of the Fouriercoefficients against frequency. The dotted vertical linesindicate the (Cowling) NS’s radial eigenmode frequencies,as listed in Table I of Ref. [63]. The evolution on the I1 grid resolves few of the star’s eigenmodes; the spec-trum has sharp peaks corresponding to the fundamentalmode and the first harmonic only. Modes with higherfrequencies (i.e., shorter wavelengths) are not spatiallyresolved by this computational grid, and so the powerthey contain aliases into the lower-frequency modes. Theevolution on I1R , on the other hand, reproduces veryclearly the fundamental mode frequency and the first fourharmonic frequencies. At higher frequencies, the peaksare still identifiable, though they become broader and lessprecisely centered. We note the presence of intermediatepeaks in the spectrum, shaped in a manner suggestiveof sidebands; however, these features are too noisy forunambiguous identification.The evolution on
I1R , with roughly 210 points acrossthe NS’s radius, has a similar resolution to the 75-elementcase presented by Radice and Rezzolla [21]. While thesetwo simulations are not directly comparable (the 1D starin the cited work self-consistently treats the gravity anduses a uniform grid), we see a qualitative agreement in thenumber of resolved modes and the precision with whichthey are resolved.We conclude the discussion of 1D evolutions by notingthat, while the
I1R grid has significantly reduced error,the lower-resolution I1 grid, representative of the 3Dresolution, is sufficient to resolve the important featuresin the evolution. The lower-resolution case remains stableon long timescales, and the oscillations as the star settlesto its numerical equilibrium correctly reflect the low-frequency eigenmodes from linearized theory.7 f [kHz] − − − − − − − − − I1I1R /10
FIG. 13. The Fourier transform of the central rest-mass den-sity ρ c from the spherically symmetric Cowling NS evolutions.The data from evolutions on the I1 and I1R grids are shown— the
I1R curve is shifted downward on the plot, by a factorof 10, for visual clarity. The vertical dotted lines indicate thefrequencies of the fundamental normal mode and the first sixharmonics. The units of the vertical axis are arbitrary.
B. Cowling neutron star in three dimensions
The simulation domain for the 3D star is a filled ballextending to r max = 24 . We consider two different cubed-sphere grids on this domain, constructed using the map-pings detailed in Appendix A. As in the spherically sym-metric case, we adapt these grids to the geometry byusing larger, higher-order elements in the star’s interioras well as outside the star. In the region near the surface,the grids are composed of thin cubed-sphere shells with alinear basis in the radial direction. The first grid, B1 , hasa similar resolution across the radius of the star to the I1 grid used in the 1D evolutions; an equatorial cut throughthis grid is shown in Fig. 14. The second grid, B1R , isobtained by adaptively refining B1 : the large elementsin the interior and exterior regions are p -refined and thethin elements near the surface of the star are h -refined,i.e., are radially split into thinner shells. The complete de-scription of these two grid structures is given in AppendixB. In these 3D evolutions, we apply as before the minmodlimiter to the surface region of the grid only. We evolvethe hydrodynamics system until t = 10 (cid:39) ms, withtime steps ∆ t = 0 . on the B1 grid ( ∆ t/ ∆ x min (cid:39) . )and ∆ t = 0 . on the B1R grid ( ∆ t/ ∆ x min (cid:39) . ).When evolving the star on these cubed-sphere grids, wefind a numerical instability in the conserved momentum ˜ S i that leads to an exponential growth of this quantityon O (100 M (cid:12) ) timescales. This numerical instability iscaused by the aliasing of the spectral modes as a resultof an insufficiently resolved quadrature rule in the DGmethod [see the paragraph below (11)]; we therefore filter FIG. 14. The B1 grid used in the 3D NS evolutions. Thethick dotted circle indicates the location of the star’s surface.The black lines show the element boundaries, and the lightgrey lines show the Gauss-Legendre-Lobatto grid within eachelement. The details of the grid mappings and structure aregiven in Appendixes A and B„ respectively. the ˜ S i variable. In the central (as described in AppendixB) portion of the grid, the filter takes the form (33) with α = 36 and s = 6 . In the cubed-sphere shells that makeup the bulk of the interior, the numerical instability isweaker and is controlled by a milder filter with α = 36 and s = 12 . With these filters, the NS evolutions remainstable until at least t = 10 . As the stars presented inthis paper have a rest state with no velocity, i.e., ˜ S i = 0 ,with dynamics that consist primarily of short-timescaleoscillations while the system settles to the rest state, thefilters have only a minor effect on the long-term evolution.For stars undergoing rotation or pronounced dynamics,the filtering should not qualitatively affect the results butwould reduce the method’s order of convergence.We plot, as before, the errors in ˜ D , ˜ S i , and ρ c duringthe first M (cid:12) of evolution time in Fig. 15. The er-rors in the B1 simulation closely match those seen in thespherically symmetric case (cf. the I1 curves in Fig. 12) —this is expected, given the comparable resolution and thespherically symmetric nature of the problem. We notethat the gradual decrease in central density seen in the1D simulation is not observed in three dimensions, andthe central density instead approaches a constant as thestar settles. The errors in the B1R case are reduced byabout an order of magnitude. Note that this case cannotbe directly compared to the high-resolution 1D case
I1R ,which has a substantially higher resolution across the NS With these values, the filter reduces power in approximately theupper half of the modes. The power in the highest mode isreduced to round off. − − − − e rr [ ˜ D ] B1B1R − − − − − k ˜ S k t [ M (cid:12) ] − . − . − . − . . . . . . ( ρ c − ρ c , ) / ρ c , / − FIG. 15. The errors in the 3D Cowling NS evolution. Thetop (middle) panel shows the error in the conserved density ˜ D (conserved momentum ˜ S i ) for evolutions using the minmodlimiter on the grids B1 and B1R . The bottom panel shows theevolution of the central density ρ c as a fractional error withrespect to its initial value ρ c, . radius. In the full evolution to t = 10 , for both grids,the oscillations damp away, and the errors tend towarda constant equilibrium value. We again check the con-servation error by tracking the baryon mass, which wefind to slowly grow during the evolution. The relativeerror in M b in the B1 evolution is comparable to the I1 case, reaching about − at t = 10 . For the B1R casethe relative error is about × − . Again, most of thesedrifts are accumulated during the initial settling of thestar.We compute the frequency spectrum of the stellar os-cillations from ρ c , using the procedure described for thespherically symmetric case. The results are shown inFig. 16. Comparing the B1 spectrum to the I1 spectrumfrom Fig. 13, we see good agreement: the first two reso-nant frequencies are clearly resolved, and additional peaksat higher frequencies are suggestive but not conclusive. f [kHz] − − − − − − − − − B1B1R /100
FIG. 16. The Fourier transform of the central rest-massdensity ρ c from the 3D Cowling NS evolution. The data fromevolutions on the B1 and B1R grids are shown — the
B1R curveis shifted downward on the plot, by a factor of 100, for visualclarity. The vertical dotted lines indicate the frequencies ofthe fundamental normal mode and the first six harmonics.The units of the vertical axis are arbitrary.
Going from B1 to B1R we see improvement in the moderesolution, with a third and fourth peak appearing in thefrequency spectrum. These new peaks are increasinglyshifted toward higher frequencies, which indicates thatthe corresponding modes are not yet fully resolved. Theintermediate peaks seen in one dimension are still visiblein three dimensions but remain close to the noise level.We also performed (but do not show) simulations of the3D Cowling NS using the same grid and limiter configura-tions that, in the 1D study, were found to be problematic.These configurations use a grid where the surface elementshave a quadratic radial basis, and/or use the MRS limiternear the NS surface. For the grid check, we employed athird cubed-sphere grid on the domain, B2 , that is simi-lar to B1 but uses thicker shells with a quadratic radialbasis (comparable to the I2 grid in one dimension; seeAppendix B for details). We found, as in the 1D case,that evolutions on this grid using the low-order minmodlimiter are stable over long timescales, but with highdissipation and increased error. For the MRS limitercheck, we found, in contrast to the 1D case, that the 3Dsimulations are unstable on O (1000 M (cid:12) ) timescales; thehigh-frequency oscillations seen in the 1D case grow in 3D,presumably due to the limiter’s inability to control theadditional tangential basis modes, and the star rapidlybecomes unstable.9 C. GR-hydro neutron star
For the coupled GR-hydro evolutions, we again use thetwo grids B1 and B1R from above. The hydrodynamics aretreated as for the Cowling star, with a minmod limiterat the star surface. We additionally evolve the spacetimegeometry, with the constraint damping parameters set to γ = 0 . (cid:2) − ( r/ / (cid:3) + 0 . (46) γ = − (47) γ = 3 exp (cid:2) − ( r/ / (cid:3) + 0 . . (48)The gauge function H σ is computed, as for the Kerr BHevolution, from the contraction of the Christoffel symbolsof the exact metric; it is constant in time. We evolvethe combined system until t = 10 (cid:39) ms, with timesteps ∆ t = 0 . on the B1 grid ( ∆ t/ ∆ x min (cid:39) . ) and ∆ t = 0 . on the B1R grid ( ∆ t/ ∆ x min (cid:39) . ).We show in Fig. 17 the errors in ˜ D , ˜ S i , and ρ c forthe self-consistent NS evolution. Comparing the resultsfrom the B1 grid to the Cowling results of Fig. 15, wesee that the self-consistent NS is more dissipative thanthe Cowling one — the oscillations decay quickly andbecome negligible by t ∼ . Additionally, we see thatthe star settles to a different equilibrium, because thegravity responds to the fluid rather than providing a fixedpotential well. The equilibrium central density is higherthan the TOV value, indicating that in its numericalequilibrium, the star has compressed slightly. The errorsusing the higher-resolution grid B1R are, as in the Cowlingcase, significantly reduced as compared to the B1 grid. Inthe full evolution to t = 10 , the B1R case exhibits a slowlygrowing error component at late times: from t (cid:39) onward, the errors increase by order 10%. This growingerror is consistent with a weak numerical instability over O (10 M (cid:12) ) timescales and could presumably be addressedby improved filtering. We again check the conservationerror by tracking the baryon mass during the evolution.The relative error in M b for the B1 case is about − at t = 10 , as in the Cowling case, and as before is mostlyaccumulated during the initial settling of the star. Inthe B1R case, however, the relative error reaches about × − , or three times the value from the Cowling case.Here, the error is accumulated in two phases: first duringthe initial settling of the star and again at the end of theevolution when the slowly growing errors begin to affectthe computation of M b .We compute once more the frequency spectrum of thestellar oscillations from ρ c , and we show in Fig. 18 theresults from the evolutions on both grids. We also in-dicate the first seven radial eigenmode frequencies fromlinear theory by the vertical dotted lines. In the lower-resolution B1 case, we see clear peaks in the spectrum from These eigenfrequencies were kindly provided to us by DavidRadice; see the discussion of Fig. 11 of Ref. [21] for details. − − − e rr [ ˜ D ] B1B1R − − − − − k ˜ S k t [ M (cid:12) ] − ( ρ c − ρ c , ) / ρ c , / − FIG. 17. The errors in the coupled GR-hydro NS evolution.The top (middle) panel shows the error in the conserved density ˜ D (conserved momentum ˜ S i ) for evolutions using the minmodlimiter on the grids B1 and B1R . The bottom panel shows theevolution of the central density ρ c as a fractional error withrespect to its initial value ρ c, . the fundamental mode up through the sixth harmonic.The first three of these peaks are sharpest, indicatingwell-resolved modes; the subsequent peaks become gradu-ally less prominent and increasingly shifted toward higherfrequencies. The B1R case is qualitatively similar — wesee the same seven peaks in the spectrum, and althoughthey are more prominent than in the B1 case becausethe noise floor is lower, the shift toward high frequenciespersists. Compared to the Cowling case, more modes areresolved. We also note that the intermediate peaks seenin the Cowling case are no longer prominent in the fullGR-hydro results.We conclude our analysis by comparing the accuracy ofthe DG and FV methods for the NS problem. We use the SpEC
GR-hydro code — a FV code that takes a dual-gridapproach for coupled GR-hydro problems — to performadditional evolutions of the NS. The spacetime is evolved0 f [kHz] − − − − − − − − − B1B1R /100
FIG. 18. The Fourier transform of the central rest-massdensity ρ c from the coupled GR-hydro NS evolutions. Resultsare shown for evolutions on the B1 and B1R grids — the
B1R curve is shifted downward on the plot, by a factor of 100, forvisual clarity. The vertical dotted lines indicate the frequenciesof the fundamental normal mode and the first six harmonics.The units of the vertical axis are arbitrary. on a high-resolution grid of nested spherical shells using apseudospectral penalty method, closely related to the DGmethod presented in this paper. The matter is evolvedon a Cartesian grid covering the interval [0 , in eachdirection (octant symmetry is imposed), using a fourth-order finite-difference scheme with a WENO reconstructor.We consider two resolutions for the FV grid. For the baseresolution, we require that the FV grid has the samenumber of grid points within the volume of the star asthe B1 grid of the DG evolution. This corresponds to agrid of points on [0 , . The high-resolution griduses points — far more than B1R . These two casesare labeled FV and FVR respectively.In Fig. 19, we compare the central density errors inevolutions with the DG and FV methods. The DG resultsmake use of the grids B1 and B1R (a two times increase inthe number of grid points), and the FV results make useof FV and FVR (an eight times increase) described above.Comparing the results, we find a few differences betweenthe DG and FV evolutions. First, the DG method is moredissipative than the FV method used, with the star’soscillations damping away by t ∼ . A contributingfactor to the increased dissipation is the use of a low-ordershock-capturing scheme in the surface regions for the DGevolutions vs the high-order reconstruction scheme of theFV method. Second, the error in the central density isgreatly reduced in the DG evolution, primarily becauseof the negligible drift rate after the star has settled toits numerical equilibrium. Finally, in going to the refined B1R and
FVR grids, we find that the error decreases morerapidly in the DG case, even though the resolution change t [ M (cid:12) ] − − − − ( ρ c − ρ c , ) / ρ c , / − B1B1RFVFVR
FIG. 19. The central-density error in the coupled GR-hydroNS, for evolutions using the DG and FV methods. For eachmethod, two resolutions are shown: the DG method uses thegrids B1 (base) and B1R (refined, with two times as many gridpoints), and the FV method uses the grids FV (base) and FVR (refined, with eight times as many grid points). is smaller. This is because the DG method has higherorder in the bulk of the star’s interior, so that p -refinementleads to rapid convergence. Precise statements about theorder of convergence for the DG results are difficult tomake, however, because we use geometrically adaptedgrids with elements of different order. VII. CONCLUSIONS
In this paper, we have presented 3D evolutions using aDG method of (a) a Kerr BH and (b) a general-relativisticNS treated self-consistently. We adopted the DG formu-lation of Teukolsky [30] to solve the generalized harmonicformulation of Einstein’s equations and the València for-mulation of general-relativistic hydrodynamics. We usedconforming grids to take advantage of the problem ge-ometries, and we evolved the spacetime and matter si-multaneously on these grids. We implemented the DGmethod in the
SpEC framework and showed convergenceand shock-capturing tests for our code. We also evolveda NS under the Cowling approximation (fixed spacetimemetric) in spherical symmetry and in three dimensions.With the 3D Kerr BH evolution, we showed that the DGmethod is accurate and stable for long-timescale spacetimeevolutions. By adapting the grid to the (nearly) sphericalgeometry of the BH spacetime, we were able to excise thesingularity from the domain — a promising result for thefuture use of the DG method in compact-object binarysimulations. The success of the DG method here draws onprevious successes of the (closely related) spectral penaltymethod for the BH problem.1For the NS, we again showed long and stable evolu-tions, and we additionally recovered the eigenfrequenciesfrom linearized theory. By using domains conforming tothe star’s spherical geometry and adapted to resolve thesurface, we were able to obtain good accuracy with com-paratively few elements and a low-order shock-capturingscheme. We compared the DG evolution to a FV evolu-tion and found significantly lower errors and an improvedrate of convergence from the DG case.One of the advantages of the DG method over the FVmethod is that it is easier to scale the algorithm on largemachines. However, we were not able to show scalingresults from our implementation in
SpEC . As discussed,the
SpEC framework scales poorly to large numbers ofelements. For the NS results shown, the domains are com-posed of over 5000 elements, enough for
SpEC ’s scalingto break down and for timing measurements to lose theirsignificance. We do note that our DG method, which usesthe same grid for the spacetime geometry and the matter,solves the Einstein equations on a denser grid of pointsthan the dual-grid
SpEC
GR-hydro code. This adds a sig-nificant computational cost for the runs presented in thispaper. The added cost would be reduced in the contextof a science-producing simulation with a spacetime gridextending to large radii, as the addition of some extragrid points in the central portion of the domain would beless significant.Improvements to our work will include the adoption ofhigher-order shock-capturing schemes (e.g., WENO) tolower the errors in the treatment of the star surface. Thedevelopment of an adaptive mesh-refinement scheme willallow geometrically adapted grids to be used in systemswith reduced symmetry and/or dynamics. These improve-ments are planned for implementation in the
SpECTRE code, where they will enable evolutions with the DGmethod of dynamical systems such as rotating or dynam-ically unstable stars.
ACKNOWLEDGMENTS
We thank Andy Bohn, Mike Boyle, Nils Deppe,Matt Duez, Francois Foucart, Jan Hesthaven, CurranMuhlberger, and Will Throwe for many helpful conversa-tions through the course of this work. We gratefullyacknowledge support for this research from the Sher-man Fairchild Foundation; from NSF Grants No. PHY-1606654 and No. AST-1333129 at Cornell; and from NSFGrants No. PHY-1404569, No. PHY-1708212, and No.PHY-1708213 at Caltech. F.H. acknowledges supportby the NSF Graduate Research Fellowship under GrantNo. DGE-1144153. Computations were performed atCaltech on the Zwicky cluster, which is supported bythe Sherman Fairchild Foundation and by NSF GrantNo. PHY-0960291, and on the Wheeler cluster, which issupported by the Sherman Fairchild Foundation. (a) Undeformed cube (b) Rounded cube
FIG. 20. Two grids on the filled ball, constructed from acubed sphere with (a) an undeformed cube as central elementand (b) a rounded cube as central element. Both panels showan equatorial cut through the grid. The grids are obtainedfrom the mappings described in Appendix A, with parameters c min = 0 in panel (a) and c min = 0 . in panel (b); in bothpanels, c max = 1 , x min = 0 . , and x max = 2 . The black linesshow the element boundaries, and the light grey lines show theGauss-Legendre-Lobatto grid within each element for order N = 5 . Appendix A: Cubed-sphere mappings
In simulations of systems with spherical geometries,we use grids based on the cubed sphere [64]. The cubedsphere is obtained by projecting the faces of a cube ontoits circumscribed sphere, thereby defining a grid on thesphere composed of six deformed Cartesian grid patches.The radial direction is introduced by a tensor product,giving a grid on a hollow spherical shell composed ofsix mapped cubes; we call each of these mapped cubesa wedge of the spherical shell. For our NS simulations,however, we require a filled ball topology, rather than ahollow spherical shell.To obtain a grid on the filled ball, one possibility isto insert a cube-shaped element at the center of the gridand deform the inner surface of the spherical shell sothat it conforms to this cube. This is shown in panel (a)of Fig. 20. In numerical experiments, we find that thisgrid configuration often suffers from large errors alongthe diagonal edges where three of the wedges meet (e.g.,the line x = y = z ) because of the large grid distortionsat these locations. This source of error can be reducedby inserting a “rounded” cube, which reduces the griddistortion in the wedges, as shown in panel (b) of Fig.20. As we are not aware of previous uses of such a gridconfiguration, we show here the mappings used.
1. Wedges
The geometry of a cubed-sphere wedge is specified byits inner and outer surfaces. Each of these surfaces isdescribed by two parameters — its curvature c and itsposition x . The surface’s curvature c ∈ [0 , controlsthe shape: when c = 0 , the surface is flat (i.e., the sixwedges together form a cube), and when c = 1 , the surface2is spherical. The surface’s position is the “radius” fromthe origin to the center of the surface, i.e., the pointwhere the surface intersects the x/y/z axis. The positions x min and x max of the inner and outer surfaces satisfy < x min < x max .The mapping from the reference element to each cubed-sphere wedge is a radial interpolation between the wedge’sinner and outer surfaces, and is computed by composingfour transformations, x (¯ x ) = ( x rot ◦ x cs ◦ x tan ◦ x affine )(¯ x ) . (A1)The actions of these transformations are:1. x affine shifts and scales the reference cube alongthe + x axis to obtain a parallelepiped spanning < x min ≤ x ≤ x max . The y and z coordinates areunaffected.2. x tan maps the tangential coordinates y and z ac-cording to y → tan( πy/ , and likewise for z . Thisshifts the grid point distribution tangentially inwardto produce a more uniform, equiangular grid whenthe destination surface is spherical. This transfor-mation is optional; we use it for the spherical-shellgrids of the Kerr BH and spherical accretion tests,but elsewhere, we omit it.3. x cs deforms the parallelepiped into one wedge of thecubed sphere, intersecting the + x axis at x min and x max . It is computed with the intermediate steps a = 1 / (cid:112) y + ¯ z (A2) b min = x min [1 + c min ( a − (A3) b max = x max [1 + c max ( a − (A4) ξ = b min + ( b max − b min ) ¯ x − x min x max − x min (A5) x cs (¯ x ) = ( ξ, ξ ¯ y, ξ ¯ z ) . (A6)4. x rot rotates the wedge to its location on the sphere,corresponding to one of the axes + x , − x , + y , − y , + z , or − z .Figure 20 shows two (filled) cubed-sphere grids wherethe outer surfaces are spherical and the inner surfaces have c = 0 and 0.66. Figure 3 shows a cubed-sphere grid whereboth surfaces are spherical, and each wedge is dividedradially and tangentially into several elements. This isachieved by dividing the unit cube into the correspondingelements before applying the chain of maps in (A1).
2. Rounded central cube
The mapping from the reference element to the roundedcentral cube is chosen to conform to the inner bound-ary of the cubed-sphere wedges. The cube is thereforeparametrized by the x min and c min that give the innerboundary of the wedges and by whether the equiangular transformation is applied. The mapping is again obtainedby composition, x (¯ x ) = ( x rc ◦ x tan )(¯ x ) , (A7)with x rc , the transformation that deforms the cube, givenby a = 1 / (cid:112) x ¯ y + ¯ x ¯ z + ¯ y ¯ z − ¯ x ¯ y ¯ z (A8) b min = x min [1 + c min ( a − (A9) x rc (¯ x ) = ( b min ¯ x, b min ¯ y, b min ¯ z ) . (A10)Inverting this mapping for ¯ x = x − rc ( x ) requires rootfinding and is done numerically.The right panel of Fig. 20 shows a cubed-sphere gridwith a rounded central cube. Figure 14 shows a roundedcube as used in the NS simulation grids; just as for thewedges, the division of the central cube into several el-ements is achieved by dividing the unit cube prior toapplying the chain of maps in (A7). Appendix B: Neutron star simulation grids
The simulation domain for the 3D NS evolutions is afilled ball extending to r max = 24 . We use three differentcubed-sphere grids on this domain: B1 , B2 , and B1R . Here,we define each of these grids.The B1 grid is shown in Fig. 14. For the bulk of thestellar interior, the grid is composed of nested, sphericalcubed-sphere shells containing higher-order elements. Inthe center of the domain, the grid transitions to a roundedcube using the mappings described in Appendix A. In theregion near the surface, we use thinner elements with fewerpoints; for B1 there are ten shells of thickness ∆ r = 0 . ,each of which contains elements with a linear basis in theradial direction (we denote this as N r = 1 ). Outside thestar, the grid is again made up of larger, higher-orderelements. The details of this radial structure are given inTable II, which lists the parameters of the cubed-sphereshells that make up the grid. The angular structure of B1 is obtained by tangentially splitting each wedge into × elements, each of which has a basis of order N tan = 3 in the two tangential directions. The resolution of thecentral rounded cube is set by conforming to the angulargrid of the shell. The equiangular tangent mapping is not applied — omitting this mapping gives a more optimalresolution of the cube in the center of the star. The B1 grid has a total of 5184 elements, with ∆ x min (cid:39) . .The B2 grid differs from B1 in the radial resolution ofthe surface region. Where B1 uses ten shells of linear order, B2 instead uses five spherical shells, of thickness ∆ r = 0 . ,each of which contains elements with a quadratic basisin the radial direction (i.e., N r = 2 ). The B2 grid has atotal of 4104 elements, and the same ∆ x min as B1 .The B1R grid is obtained from B1 by selectively refiningto further take advantage of the hp -adaptivity of the DGmethod: h -refinement is used in the neighborhood of the3 TABLE II. The radial structure of the 3D NS grids, B1 , B2 , and B1R . For each region of each grid, the location and curvature ofthe surfaces that bound the cubed-sphere elements are given.Duplicated information is omitted: the unspecified regions of B2 are identical to those of B1 . x i c i N r B1 Center 1.3, 1.9, 2.5 0.55, 0.85, 1 4Int. 2.5, 3.0, 3.6, 4.33, 5.2, 6.24, 7.5 1 4Surf. 7.5, 7.75, 8, ..., 9.5, 9.75, 10 1 1Ext. 10, 12, 15, 18, 21, 24 1 3 B2 Surf. 7.5, 8, 8.5, 9, 9.5, 10 1 2
B1R
Center (see B1 ) (see B1 ) 5Int. (see B1 ) (see B1 ) 5Surf. 7.5, 7.625, 7.75, ..., 9.875, 10 1 1Ext. (see B1 ) (see B1 ) 4 surface where the solution is not smooth, and p -refinementis used in the smooth interior and exterior regions. Theradial parameters are again given in Table II; the angularparameters are as for B1 but with N tan = 4 . This gridhas 7344 elements, with ∆ x min (cid:39) . , and has roughlytwice as many grid points inside the NS ( r (cid:46) . ) asthe B1 grid. REFERENCES [1] J. A. Font, Living Rev. Rel. , 7 (2008).[2] J. M. Martí and E. Müller, Living Rev. Comput. Astro-phys. , 3 (2015).[3] D. S. Balsara, Living Rev. Comput. Astrophys. , 2 (2017),arXiv:1703.01241 [astro-ph.IM].[4] M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer, M. A.Scheel, and S. A. Teukolsky, Phys. Rev. D , 104015(2008), arXiv:0809.0002 [gr-qc].[5] P. Colella and P. R. Woodward, J. Comput. Phys. ,174 (1984).[6] A. Mignone, T. Plewa, and G. Bodo, Astrophys. J. Suppl.Ser. , 199 (2005), astro-ph/0505200.[7] G.-S. Jiang and C.-W. Shu, J. Comput. Phys. , 202(1996).[8] S. W. Bruenn, A. Mezzacappa, W. R. Hix, E. J. Lentz,O. E. B. Messer, E. J. Lingerfelt, J. M. Blondin, E. Endeve,P. Marronetti, and K. N. Yakunin, Astrophys. J. Lett. , L6 (2013), arXiv:1212.1747 [astro-ph.SR].[9] P. Mösta, C. D. Ott, D. Radice, L. F. Roberts,R. Haas, and E. Schnetter, Nature , 376 (2015),arXiv:1512.00838 [astro-ph.HE].[10] K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata,and T. Wada, Phys. Rev. D , 041502 (2014),arXiv:1407.2660 [astro-ph.HE].[11] S. Bernuzzi and T. Dietrich, Phys. Rev. D , 064062(2016), arXiv:1604.07999 [gr-qc].[12] F. Foucart, R. Haas, M. D. Duez, E. O’Connor, C. D.Ott, L. Roberts, L. E. Kidder, J. Lippuner, H. P. Pfeif-fer, and M. A. Scheel, Phys. Rev. D , 044019 (2016),arXiv:1510.06398 [astro-ph.HE].[13] J. C. McKinney, A. Tchekhovskoy, A. Sadowski, andR. Narayan, Mon. Not. Roy. Astr. Soc. , 3177 (2014),arXiv:1312.6127 [astro-ph.CO].[14] O. Porth, H. Olivares, Y. Mizuno, Z. Younsi, L. Rez-zolla, M. Moscibrodzka, H. Falcke, and M. Kramer,Computational Astrophysics and Cosmology , 1 (2017),arXiv:1611.09720 [gr-qc].[15] B. J. Kelly, J. G. Baker, Z. B. Etienne, B. Giacomazzo,and J. Schnittman, Phys. Rev. D , 123003 (2017),arXiv:1710.02132 [astro-ph.HE]. [16] G. Zumbusch, Class. Quantum Grav. , 175011 (2009),arXiv:0901.0851 [gr-qc].[17] S. E. Field, J. S. Hesthaven, S. R. Lau, and A. H. Mroue,Phys. Rev. D , 104051 (2010), arXiv:1008.1820 [gr-qc].[18] J. D. Brown, P. Diener, S. E. Field, J. S. Hesthaven,F. Herrmann, A. H. Mroué, O. Sarbach, E. Schnetter,M. Tiglio, and M. Wagman, Phys. Rev. D , 084004(2012), arXiv:1202.1038 [gr-qc].[19] J. M. Miller and E. Schnetter, Class. Quantum Grav. ,015003 (2017), arXiv:1604.00075 [gr-qc].[20] M. Dumbser, F. Guercilena, S. Köppel, L. Rezzolla,and O. Zanotti, Phys. Rev. D97 , 084053 (2018),arXiv:1707.09910 [gr-qc].[21] D. Radice and L. Rezzolla, Phys. Rev. D , 024010(2011), arXiv:1103.2426 [gr-qc].[22] J. Zhao and H. Tang, J. Comput. Phys. , 138 (2013).[23] M. Bugner, T. Dietrich, S. Bernuzzi, A. Weyhausen,and B. Brügmann, Phys. Rev. D , 084004 (2016),arXiv:1508.07147 [gr-qc].[24] L. E. Kidder, S. E. Field, F. Foucart, E. Schnet-ter, S. A. Teukolsky, A. Bohn, N. Deppe, P. Diener,F. Hébert, J. Lippuner, J. Miller, C. D. Ott, M. A.Scheel, and T. Vincent, J. Comput. Phys. , 84 (2017),arXiv:1609.00098 [astro-ph.HE].[25] P. Anninos, C. Bryant, P. C. Fragile, A. M. Holgado,C. Lau, and D. Nemergut, Astrophys. J. Suppl. Ser. ,17 (2017), arXiv:1706.09939 [astro-ph.IM].[26] F. Fambri, M. Dumbser, S. Köppel, L. Rezzolla, andO. Zanotti, Mon. Not. Roy. Astr. Soc. (2018), 10.1093/mn-ras/sty734, arXiv:1801.02839 [physics.comp-ph].[27] S. Köppel, J. Phys. Conf. Ser. , 012017 (2018),arXiv:1711.08221 [gr-qc].[28] F. Hébert, Exploring Relativistic Gravity with NumericalSimulations , Ph.D. thesis, Cornell University (2017).[29] M. Bugner,
Discontinuous Galerkin Methods for Gen-eral Relativistic Hydrodynamics , Ph.D. thesis, Friedrich-Schiller-Universität Jena (2018).[30] S. A. Teukolsky, J. Comput. Phys. , 333 (2016),arXiv:1510.01190 [gr-qc].[31] H. Friedrich, Commun. Math. Phys. , 525 (1985). [32] F. Pretorius, Class. Quantum Grav. , 425 (2005), gr-qc/0407110.[33] C. Gundlach, G. Calabrese, I. Hinder, and J. M. Martin-Garcia, Class. Quantum Grav. , 3767 (2005), gr-qc/0504114.[34] .[35] M. A. Scheel, H. P. Pfeiffer, L. Lindblom, L. E. Kidder,O. Rinne, and S. A. Teukolsky, Phys. Rev. D , 104006(2006), gr-qc/0607056.[36] B. Szilágyi, L. Lindblom, and M. A. Scheel, Phys. Rev.D , 124010 (2009), arXiv:0909.3557 [gr-qc].[37] D. A. Hemberger, M. A. Scheel, L. E. Kidder, B. Szilágyi,G. Lovelace, N. W. Taylor, and S. A. Teukolsky, Class.Quantum Grav. , 115001 (2013), arXiv:1211.6079 [gr-qc].[38] R. Haas, C. D. Ott, B. Szilágyi, J. D. Kaplan, J. Lippuner,M. A. Scheel, K. Barkett, C. D. Muhlberger, T. Dietrich,M. D. Duez, F. Foucart, H. P. Pfeiffer, L. E. Kidder,and S. A. Teukolsky, Phys. Rev. D D93 , 124062 (2016),arXiv:1604.00782 [gr-qc].[39] J. Hesthaven and T. Warburton,
Nodal DiscontinuousGalerkin Methods: Algorithms, Analysis, and Applications (Springer, Berlin, New York, 2008).[40] D. A. Kopriva,
Implementing spectral methods for par-tial differential equations: Algorithms for scientists andengineers (Springer, Berlin, New York, 2009).[41] J. Hesthaven, S. Gottlieb, and D. Gottlieb,
Spectral Meth-ods for Time-Dependent Problems (Cambridge UniversityPress, Cambridge, UK, 2007).[42] D. Gottlieb and J. S. Hesthaven, J. Comput. Appl. Math. , 83 (2001).[43] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, andO. Rinne, Class. Quantum Grav. , S447 (2006), gr-qc/0512093.[44] F. Banyuls, J. A. Font, J. M. Ibáñez, J. M. Martí, andJ. A. Miralles, Astrophys. J. , 221 (1997).[45] F. Galeazzi, W. Kastaun, L. Rezzolla, and J. A. Font,Phys. Rev. D , 064009 (2013), arXiv:1306.4953 [gr-qc]. [46] A. Harten, P. D. Lax, and B. van Leer, SIAM Review , 35 (1983).[47] E. F. Toro, Riemann solvers and Numerical Methodsfor Fluid Dynamics: A Practical Introduction (Springer,Berlin, New York, 2013).[48] S. F. Davis, SIAM Journal on Scientific and StatisticalComputing , 445 (1988).[49] B. Cockburn, J. Comput. Appl. Math. , 187 (2001).[50] S. A. Moe, J. A. Rossmanith, and D. C. Seal, ArXive-prints (2015), arXiv:1507.03024 [math.NA].[51] X. Zhong and C.-W. Shu, J. Comput. Phys. , 397(2013).[52] J. Zhu, X. Zhong, C.-W. Shu, and J. Qiu, Communica-tions in Computational Physics , 944 (2016).[53] S. Gottlieb, C.-W. Shu, and E. Tadmor, SIAM Review , 89 (2001).[54] M. Alcubierre et al. , Class. Quant. Grav. , 589 (2004),arXiv:gr-qc/0305023 [gr-qc].[55] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravi-tation (Freeman, New York, New York, 1973).[56] F. C. Michel, Astrophysics and Space Science , 153(1972).[57] B. Cockburn, S.-Y. Lin, and C.-W. Shu, J. Comput. Phys. , 90 (1989).[58] J. Centrella and J. R. Wilson, Astrophys. J. Suppl. Ser. , 229 (1984).[59] L. Del Zanna and N. Bucciantini, Astron. Astrophys. ,1177 (2002), astro-ph/0205290.[60] R. C. Tolman, Phys. Rev. , 364 (1939).[61] J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. , 374(1939).[62] A. Schoepe, D. Hilditch, and M. Bugner, Phys. Rev. D , 123009 (2018), arXiv:1712.09837 [gr-qc].[63] J. A. Font, T. Goodale, S. Iyer, M. Miller, L. Rezzolla,E. Seidel, N. Stergioulas, W.-M. Suen, and M. Tobias,Phys. Rev. D , 084024 (2002), gr-qc/0110047.[64] C. Ronchi, R. Iacono, and P. S. Paolucci, J. Comput.Phys.124