Gmunu: Toward multigrid based Einstein field equations solver for general-relativistic hydrodynamics simulations

We present a general relativistic hydrodynamics code Gmunu (General-relativistic multigrid numerical solver) which uses a multigrid method to solve the elliptic metric equations in the conformally flat condition (CFC) approximation. Most of the existing relativistic hydrodynamics codes are based on formulations which rely on a free-evolution approach of numerical relativity, where the metric variables are determined by hyperbolic equations without enforcing the constrained equations in the evolution. On the other hand, although a fully constrained-evolution formulation is theoretical more appealing and should lead to more stable and accurate simulations, such an approach is not widely used because solving the elliptic-type constrained equations during the evolution is in general more computationally expensive than hyperbolic free-evolution schemes. Multigrid methods solve differential equations with a hierarchy of discretizations and its computational cost is generally lower than other methods such as direct methods, relaxation methods, successive over-relaxation. With multigrid acceleration, one can solve the metric equations on a comparable time scale as solving the hydrodynamics equations. This would potentially make a fully constrained-evolution formulation more affordable in numerical relativity simulations. As a first step to assess the performance and robustness of multigrid methods in relativistic simulations, we develop a hydrodynamics code that makes use of standard finite-volume methods coupled with a multigrid metric solver to solve the Einstein equations in the CFC approximation. In this paper, we present the methodology and implementation of our code Gmunu and its properties and performance in some benchmarking relativistic hydrodynamics problems.


Introduction
In the past decade, numerical relativity has matured to the state that stable and robust numerical calculations of the Einstein equations with or without matter has become feasible. The spacetimes of many interesting astrophysical systems such as stellar core collapses and binary systems of compact objects have been accurately modeled (see, e.g., [43,14,13] for recent reviews). In the standard 3 + 1 decomposition of spacetime, the Einstein equations are split into a set of evolution equations and constraint equations. Nevertheless, one still has the freedom to choose the basic variables to evolve and reformulate the resulting systems of differential equations in order to improve the stability and accuracy of numerical simulations. This results in different formulations of numerical relativity, such as the so-called BSSN [49,2], CCZ4 [7], and Z4c [6] schemes, which are popular choices for numerical modelings. The practical applications of these different formulations are based on a free-evolution approach where the constraint equations are first solved for preparing the initial data and used subsequently only as an indicator to monitor the numerical accuracy during the evolution (see, e.g., [37]).
Alternatively, one can also formulate the Einstein equations based on a fully constrained-evolution approach where the constraint equations are solved and fulfilled to within the discretization errors during the evolution. Despite the fact that a constrained-evolution approach is theoretical appealing, such an approach is not popular among numerical relativists since solving the elliptic-type constraint equations during the evolution is generally computational expensive. In contrast to the active development and applications of free-evolution formulations, the last proposed constrained formulation of the Einstein equations was already 15 years ago due to Bonazzola et al. [8]. The fully constrained-evolution formulation of Bonazzola et al. [8] has been employed to simulate pure gravitational wave spacetime [8], and also an oscillating neutron star by ignoring the back-reaction of the gravitational waves into the fluid dynamics [16]. However, the application of this constrained scheme and assessment of its performance in modelling more generic dynamical spacetimes without symmetry is still a largely unexplored area.
It is worth to point out that the fully constrained scheme of Bonazzola et al. [8] automatically reduces to the so-called conformally flat condition (CFC) approximation to general relativity [55,32] if a tensor field h ij introduced in their formulation is set to zero (see [15] for a detailed discussion). The CFC approximation results in a simpler set of elliptic equations for the metric sector. Numerical simulations based on the CFC scheme have been successfully carried out for various astrophysical problems [18,40,48,4,3,5,38] and the scheme has also been shown to be a good approximation to full general relativity in rotating iron core collapses [41]. However, the original CFC scheme suffers from mathematical non-uniqueness problems when the system is too compact. In order to overcome the non-uniqueness issue, the scheme was reformulated and extended to the so-called extended CFC (xCFC) scheme so that the modelling of extreme spactimes such as black hole formation becomes possible [15,39,16].
We have in our mind a motivation to experiment and develop our own general relativistic hydrodynamics code based on the fully constrained formulation of Bonazzola et al. [8] (or other similar constrained formulations if available in the future) which maximizes the use of elliptic-type equations for the metric sector of the system in the evolution. In this paper, we take a first step along this direction by developing a relativistic hydrodynamics code based on the xCFC scheme. Although it is not fully general relativistic, the xCFC scheme contains a set of similar, but simpler, elliptic equations as the fully constrained formulation. We can thus use the xCFC scheme to evaluate the performance and robustness of our metric solver.
As already pointed out, it is known in general that solving the elliptic equations frequently in a fully constrained or xCFC scheme during a simulation is computationally expensive. Many numerical methods have been explored to deal with such elliptic systems, including finite-difference methods, different types of iterative solvers, and spectral methods (see [17,19,12] and references therein). A seminal work is due to Dimmelmeier et al. [19] which combines a finite-difference grid and a spectral grid, on which the hyperbolic hydrodynamics and elliptic metric equations are solved, respectively. However, even though the spectral method is known to be extremely fast and accurate, the metric solver is still one of the bottlenecks to slow down a hydrodynamics simulation as the communication between variables defined on the two different grids is time consuming especially in the multidimensional cases [20,19]. In this work, we propose and demonstrate that Multigrid method is an efficient strategy to solve the elliptic metric equations in hydrodynamical simulations.
Multigrid methods solve differential equations with a hierarchy of discretizations and its computational cost is generally lower than other methods such as direct methods, relaxation methods, and successive over-relaxation [28]. The multigrid strategy has been employed in a wide range of problems and it has also been used to generate initial data in numerical relativity [31,22,36]. However, multigrid methods have not been applied in any constrained-evolution schemes for numerical relativity. In order to couple to the matter directly, nonlinear cell-centred multigrid (CCMG) and the corresponding boundary treatments are needed, the latter of which is more complicated than the vertex-centred multigrid and is still being actively studied in the computational physics and applied mathematics.
Our aim is to construct a direct, rapid, and robust multidimensional metric solver which can be easily coupled to the matter based on the non-linear cell-centred multigrid strategy. In this paper, we present the methodology and implementation of our general relativistic hydrodynamics code Gmunu (General-relativistic MUltigrid NUmerical solver), which solves the hydrodynamics equations using standard finitevolume methods and the xCFC metric equations using a multigrid approach. To the best of our knowledge, this is the first relativistic hydrodynamics code that makes use of a multigrid solver in dynamical simulations. We also perform various benchmarking tests in relativistic hydrodynamics to assess the performance and robustness of our code.
The paper is organised as follows. In section 2 we outline the formalism we used in this work. The details of the numerical settings and, the methodology, implementation of our hydrodynamics solver and our multigrid solver are presented in section 3 and section 4 respectively. The code tests and results are presented in section 5. This paper ends with a discussion section in section 6.

Metric equations and Conformal flatness approximation
We use the standard ADM 3+1 formalism [27,1]. The metric can be written as where α is the lapse function, β i is the spacelike shift vector and γ ij is the spatial metric.
In the 3+1 formalism, the Einstein equations are split into a set of constraint equations which must be satisfied on every hypersurface and a set of the evolution equations for γ ij and the extrinsic curvature K ij where ∇ i is the covariant derivative with respect to the three-metric γ ij , R ij is the corresponding Ricci tensor, R is the scalar curvature and K is the trace of the extrinsic curvature K ij . For the matter sources, E := n µ n ν T µν , S i := −n µ γ i ν T µν and S ij := γ i µ γ j ν T µν , where T µν is the energy-momentum tensor and n µ is the unit normal vector of a spacelike hypersurface.
It is difficult to maintain the constraint equations in the numerical evolution of the evolution equations above because these ADM equations are numerically unstable. There are serveral different re-formulations of 3 + 1 numerical relativity that can lead to stable evolutions [49,2,7,6]. However, these schemes are based on a free-evolution approach where the Einstein equations are evolved with hyperbolic-type equations. The constraint equations are only used for solving the initial data and serve as a monitor for numerical errors during the simulations. On the other hand, a fullyconstrained evolution approach where the constraints are enforced at each time step is generally not favored as solving the elliptic-type constraint equations is computational expensive comparing to hyperbolic equations [45]. We have a motivation to develop and experiment efficient multigrid solvers for elliptic-type metric equations that one needs in order to carry out fully-constrained evolutions for numerical relativity. As a first step towards this goal, our relativistic hydrodynamics code employs the xCFC scheme which is an improved version of the CFC approximation to general relativity.
In a CFC approximation [17,12], the three metric γ ij is assumed to be decomposed according to where f ij is a time independent flat background metric and ψ is the conformal factor which is a function of space and time. Another assumption is the maximum slicing condition of foliations K = 0.
With these conditions, one can derive the time derivative of the conformal factor ψ and also the extrinsic curvature K ij The CFC approximation of the ADM equations can be reduced into five coupled nonlinear elliptic equations where∇ i and∆ are the covariant derivative and the Laplacian with respect to the flat three metric f ij , respectively. The original CFC scheme suffers from mathematical non-uniqueness problems. The CFC scheme was later reformulated so that the elliptic equations are fully decoupled and the local uniqueness of the solution is guaranteed. The reformulated CFC scheme is the so-called xCFC scheme [15], which is the scheme that we implemented in Gmunu. In the xCFC scheme, one introduces a vector potential X i , and the metric can be solved by the following equations: whereẼ := ψ 6 E,S i := ψ 6 S i andS := ψ 6 S are the rescaled fluid source terms. The tensor fieldÃ ij can be approximated on the CFC approximation level by (see the Appendix of [15]

General relativistic hydrodynamics equations
The evolution equations for the matter are derived from the local conservations of the rest-mass and energy-momentum: where ρ is the rest-mass density of the fluid and u µ is the fluid four-velocity. For a perfect fluid, the energy-momentum tensor is given by T µν = ρhu µ u ν + P g µν , where P is the pressure, h = 1 + + P/ρ is the enthalpy, and is the specific internal energy. The three-velocity v i of the fluid as measured by the Eulerian observers of four-velocity n µ is given by In the flux-conservative Valencia formulation (e.g., [25]), the set of hydrodynamics equations are given by where As shown in Eq. (22), the source terms Q contain the time derivatives of the metric quantities. In order to reduce the accumulated error due to the time update, it is good to avoid the time derivatives in the code. We can rewrite the Q j terms into compact form [19] It is also possible to bypass the time derivatives in the Q τ term (i.e., the last element of the vector in Eq. (22)): In order to adapt to the extended CFC scheme, we have to evolve the conformal transformed conserved quantities, the (non-conformal transformed) conserved variables and thus the primitive variables will be updated once the conformal factor ψ is solved.
In particular, we defineŨ ≡ ψ 6 U ,F i ≡ ψ 6 F i andQ ≡ ψ 6 Q. We can then reformulate the hydrodynamics equations as In the case that we need to solve the metric, we pass the conformal conserved quantities into the metric solver. Both metric quantities and primitive hydrodynamic variables will then be updated.

Numerical methods and implementation
We use spherical polar coordinates {r, θ, φ} and adopt the axial symmetry with cellcentered discretization. In particular, the coordinate grid covers 0 < r < r max and 0 < θ < π/2 and we discretize it into n r × n θ cells with uniform coordinate grid spacing, i.e. ∆r = r max /n r and ∆θ = π/2n θ . Gmunu solves the general relativistic hydrodynamic equations by using standard high-resolution shock-capturing (HRSC) schemes [33]. In particular, various different cell-interface reconstruction methods and Riemann solvers are implemented and tested. For the cell-interface reconstruction, we have implemented piecewise constant scheme (PC), monotonized-central limiter (MC), 5 th order weighted-essentially nonoscillatory scheme (WENO5) [50] and 5 th order monotonicity preserving scheme (MP5) [52]. For the Riemann solver, the Rusanov flux (also known as Total variation diminishing Lax-Friedrichs scheme (TVDLF)) [53], Harten-Lax-van Leer (HLL) [30], Harten-Lax-van Leer-Einfeldt (HLLE) [23,29], Marquina flux formula [21] have been implemented. For the recovery of the primitive variables (ρ, v i , P ) from the conservative variables (D, S i , τ ), we follow the formulation presented in the Appendix C in [26] and use Regula-Falsi method to find the root. This formulation was shown to be not only robust, accurate and efficient, but also suitable for different kinds of relativistic hydrodynamical simulations [47,26]. For the region outside the star, we fill with an "atmosphere" of density ρ atmo = 10 −6 ρ max (t = 0), where ρ max (t = 0) is the maximum density over the whole computational domain initially at time t = 0, and use the standard atmosphere treatment during the simulation. In particular, if the density at a particular grid drops below ρ atmo , the density at that grid is reset to ρ atmo . The velocity of that grid is also set to be zero and other primitive variables such as the specific internal energy is updated accordingly.
The simulations reported in this paper were preformed with HLLE Riemann solver and the MC reconstruction method. We also use a 3 rd order Runge-Kutta integrator for the time integration. It is not necessary to solve the metric at each time step [17,12]. To speed up the simulations, we solve the metric at every 50 time steps, and extrapolate the metric in between. A practical difficulty related to spherical coordinates is that the convergence of grid points near the pole axis puts a severe constraint on the time step imposed by the Courant-Friedrichs-Levy stability condition in numerical simulations. In order to increase the size of the time steps in our simulations, we treat the first 10 grid points, which cover about 5% of the stellar radius, as a spherically symmetric core (i.e., only radial motions are allowed). This should be a good approximation as non-radial fluid motions in the core is negligibly small. We used the open-source code XNS [12,42,43,44], which is also based on the CFC approximation, to generate initial neutron-star models for our dynamical simulations.

Metric solver with xCFC scheme
By following [15], the metric equations Eqs. (13)-(16) can be decoupled, thus they can be solved in a hierarchical way and the local uniqueness is guaranteed. Once the time integration at each time step for the hydrodynamics equations is completed, the conformally rescaled hydrodynamical conserved variables (D,S i ,τ ) are updated, and they will be used to solve the metric equations. The steps for solving the metric equations are summarized in the following: (i) Solve Eq. (13) for the vector potential X i from the conserved variablesS i .
(ii) Calculate the tensorÃ ij in Eq. (17) from the vector potential X i .
(iv) With the updated conformal factor ψ, calculate the conserved variables (D, S i , τ ) and thus convert the conserved variables to the primitive variables (ρ, v i , P ). Theñ S can be worked out consistently.
To solve the metric equations, appropriate boundary conditions at the origin and outer computational boundary are required. In order to let the solution fall off as the Schwarzschild solution at large distance, we require ψ = C r + 1, α = C r + 1 and X i = β i = 0 at the outer boundary (r = r max ). In Gmunu, we impose the boundary conditions More detailed implementation of the metric equations and boundary conditions at the axis and the origin can be found in Appendix A.

An overview of multigrid
Multigrid is an efficient method for solving elliptic partial differential equations with low memory and work complexity. Different modes are filtered out with different rates at different resolutions. For some low-frequency modes, it is computationally expensive to compute directly on a high-resolution discretization. However, it can be done efficiently at a low-resolution discretization. The main concept of the multigrid method is to solve the problem recursively with a series of coarse grids. It is noted that the multigrid method is not a single method but a strategy with many possible implementations. However, the elements needed to construct a multigrid solver is more or less the same for different implementations. For instance, as shown in figure 1, the key ingredients of multigrid solver includes (1) a cycle framework as the backbone of the whole solver,  Figure 1: Three types of (4-grid) cycles can be used in multigrid methods. "S" denotes smoothing while "Sol" denotes solving the equation directly. Each descending line \ represents restriction and each ascending line / represents prolongation. The key ingredients of multigrid solver includes a cycle framework, restriction and prolongation operators, smoothers and solvers.

Nonlinear multigrid
Due to the nonlinearity of the Einstein equations needed to be solved (i.e. Eqs. (13)-(16)), nonlinear multigrid method is required. The implementations for the multigrid for non-linear elliptic equations are different from the linear cases. Two well-known methods for solving non-linear partial differential equations with multigrid techniques are the Full Approximation Scheme (FAS) [10] and Newton-multigrid (Newton-MG) [11,54]. The two methods are widely used and obtained successes in various problems. We refer the interested reader to [9] for a detailed comparison of the two methods. Due to the fact that the memory used is low in FAS, the current version of Gmunu adopts the FAS algorithm (see [10,9,46] and references therein for more details). Algorithm 1 shows the pseudocode of a single cycle of the non-linear multigrid elliptic partial differential equation solver implemented in Gmunu.

Cell-centered discretization and intergrid transfer operators
As mentioned in section 2.1, the source terms of the metric equations consist of the hydrodynamical variables. Meanwhile, since Gmunu solves the hydrodynamics with finite volume approach, the grids are discretized with cell-centred discretization. On the other hand, the discretization for the metric solver is in general different from the hydrodynamics sector. For instance, for the pseudospectral method, the choice of grid points has to be consistent with the basis functions which can not be chosen arbitrarily. The corresponding interpolation or extrapolation are needed so that the hydrodynamical variables can be passed correctly into the metric solver (e.g. [19,12]), which might be another bottleneck of the computational time of the simulations. In order to adapt the grid of the hydrodynamics sector so that the hydrodynamical variables can be passed into the metric solver without any interpolation or extrapolation, we implemented the cell-centred multigrid (CCMG) [35,28], which is one of the novelties of Gmunu.
Constructing a cell-centred multigrid solver is non-trivial. Unlike the vertex-centred case, in which a node of the coarse grid is also a node of the fine grid, the nodes on coarser grids do not form a subset of the fine grid nodes in the case of cell-centred discretization. The choices of inter-grid transfer operators and the boundary condition implementation are different from the vertex-centred cases. There are also many possible approaches for different situations. Indeed, constructing problem-independent efficient cell-centred transfer operators is still under an active research area in computational physics and applied mathematics (see [35] and references therein).
There are many possible choices of restriction and prolongation operators and they cannot be chosen arbitrarily [35,28]. In Gmunu, we adopt the most standard combination, i.e., piecewise constant restriction (figure 2a) and bi-linear prolongation (figure 2b).

Key features of the nonlinear cell-centered multigrid metric solver in Gmunu
This section shortly summarizes the key features of the metric solver implemented in Gmunu. In our multigrid metric solver, we adopt the Full Approximation Storage (FAS) to deal with the nonlinear metric equations with V-, W-and F-cycle implemented. For the smoother and solvers, we use the standard red-black Gauss-Seidel relaxation. In particular, the smoother consists of 15-times relaxations and the direct solver consists of 200-times relaxations. For the inter-grid transfer operators, we adopt piecewise constant restriction (figure 2a) and bi-linear prolongation (figure 2b). Note that multigrid solvers are iterative solvers. Practically, this function has to be called until the solution converges, i.e., when the L ∞ or L 1 norm of the residual is below some chosen threshold value. The "*" denote the location of the coarse grid note. The notation shows the weighting of the value which are the neighbor of the coarse grid node "*".

Convergence properties
Before applying Gmunu to any astrophysical problems, we first demonstrate the convergence properties and performance of our multigrid metric solver. We solve the metric of model BU8 (see table 1), which represents a rapidly rotating neutron star and far from spherically symmetric, with our metric solver with the resolution n r × n θ = 640 × 64. The computational domain covers r = [0, 30] and θ = [0, π/2]. To compare the convergence rate, we focus on solving the lapse function α with the flat space initial guess α = 1. Figure 3 shows the L 1 norm of the residual of Eq. (15) as a function of the number of iterations with different level of V-cycle. The number represents how deep the solver goes when solving for the lapse function. The horizontal black dashed line represents the threshold tolerance. As we can see in figure 3, the solver convergence faster when it goes to a deeper level. Also, it takes O(10 5 ) interactions (not shown in the plot) to reach the threshold tolerance for the V1 (or Gauss-Seidel) case, while it takes only 37 steps for the V6 case.
In practice, at the beginning of the simulation, we use the initial data provided by XNS as initial guess. During the evolution, we use the previous solution as initial guess for the next iteration. This makes the solver converge much faster as the solutions on previous time step are usually good approximation to the solution.

Code tests
We present a set of numerical tests of Gmunu. Some standard hydrodynamics tests in a static background metric such as shocktube [34] and Cowling approximation [51] can be found in the Appendix B. Here we will mainly focus on hydrodynamic evolution with dynamical spacetime. Table 1 lists the models we used in various tests. The name of the models are originally defined in the literature [20,15]. All the models are constructed with the polytropic equation of state with Γ = 2 and K = 100 (in units of  (c = G = M = 1)). Table 1: The equilibrium neutron star models we used in this paper. The name of the models are originally defined in the literature [20,15]. "BU" represents a sequence of fixed central rest-mass density ρ c = 1.28 × 10 −3 uniformly rotating models, and "SU" is a nonrotating unstable model. All the models are constructed with the polytropic equation of state with Γ = 2 and K = 100. Ω is the angular velocity; M is the gravitational mass; r e and r p are the equatorial and polar radii, respectively. Unless otherwise noted, we use units where c = G = M = 1.

Stability of a stable TOV star
The first test we perform is the stability of a stable spherically symmetric neutron star BU0. Even though BU0 is a 1D model, we run this test with a 2-dimensional setup. The resolution of the simulation is n r × n θ = 640 × 64, where r = [0, 30] and θ = [0, π/2]. While BU0 is a static and stable configuration, the discretization errors and the diffusion at the contact discontinuity of the neutron star surface trigger stellar oscillation modes. The upper panel of figure 4 shows the relative variation of the central density ρ c as a function of time. The relative variation of the density is of the order 10 −4 . The code is able to evolve BU0 stably for more than 10 ms as expected since BU0 is a stable configuration. The Fourier transform of the radial velocity v r (t) and non-radial velocity v θ (t) at r = 5, θ = π/4 (inside the neutron star). The vertical lines represent the known eigenmodes frequency. Our results agree with the known eigenmodes.
One way to test if the code handles a dynamical spacetime correctly, at least in the linear regime for small perturbations, is to extract the eigenmode frequencies from the simulations and compare them with the known values from perturbative calculations. Note that the eigenmode frequencies of an oscillating neutron star in a dynamical spacetime are quite different from the Cowling approximation where the metric is kept fixed in time (see [20] and also compare the results in Appendix B.2). In order to compare the eigenmodes more clearly, instead of applying the Fourier transform on the central density directly, we analyse the radial component of the velocity v r and the θ-component v θ at r = 5, θ = π/4 (inside the neutron star). The lower panel of figure 4 shows the Fourier transform of v r and v θ . The vertical dashed lines represent the known and well-tested eigenmode frequencies. Our results agree with the known eigenmode frequencies [24,20,12].

Stability of rotating neutron stars
The tests we perform in this subsection are the stability of stable rotating neutron stars. The resolution of these simulation is again n r × n θ = 640 × 64, with r = [0, 30] and θ = [0, π/2].
The upper panel of figure 5 shows the evolution of a rotating stable neutron star BU2. This model is stably evolved for more than 20 ms, where the relative variation of ρ c is about the order of O(10 −4 ). The Fourier transforms of v r and v θ are shown in the lower panel of the figure 5. The extracted eigenmodes again agree with the known results [20].  [20]. Our results agree with the known eigenmodes.
Comparing the density profiles and the rotational velocity profiles at later times with the initial profiles serves as another indicator for code performance. Figure 6 shows the comparison between the initial density profiles and the rotational velocity (solid lines) with the same quantities (dashed lines) at T max = 20 ms of BU2. It can be seen that the density profiles along the equatorial and polar radii are maintained very well during the evolution. However, the rotational profile is slightly suppressed at the surface of the star. This is due to the fact that the Riemann solver HLLE we used in these tests is known to be too diffusive to deal with the star surface or any discontinuous surface [12]. The same test with the same numerical setup has been done for the model BU8. Unlike the moderately rotating case BU2, model BU8 is a rapidly rotating neutron star which is close to the mass shedding limit and therefore this test is more demanding. Even for this demanding case with the diffusive Riemann solver HLLE, Gmunu is able to maintain this model for more than 20 ms. The relative variation in ρ c is about O(10 −4 ) as shown in the upper panel of figure 7. The oscillation modes extracted from the simulation also agree with the results reported by other groups [20], as shown in the lower panel in figure 7. Moreover, as shown in the left panel in figure 8, the density profiles are well-preserved even up to T max = 20 ms. However, unlike the previous cases, the average value of the central density increases slowly during the evolution, and the angular velocity profile is slightly distorted at the surface of the star.

Migration of an unstable TOV star
Previous tests are all based on the evolution of stable neutron stars, where the configurations are almost stationary with some perturbations. In this subsection, we report the performance of the code in the fully non-linear regime with significant  Figure 7: The evolution of a stable rapidly rotating neutron star (BU8) with the resolution n r × n θ = 640 × 64. Long term evolution of a rapidly rotating neutron star in a dynamical spacetime is known as a challenging test. Even though the HLLE solver might be too diffusive, Gmunu is still able to maintain this model for more than 20 ms with small variations and the fluid behaves correctly and the extracted eigenmodes agree with the results from the other groups [20]. Upper panel : The relative variation of the central density ρ c is about O(10 −4 ) and the average value ρ c increases slowly during the evolution. Lower panel : The Fourier transform of the radial velocity v r (t) and non-radial velocity v θ (t) at r = 5, θ = π/4 (inside the neutron star). The vertical lines represent the known and well-tested eigenmode frequencies [20].
changes and coupling in the metric and fluid variables. One of the standard tests for hydrodynamical evolution coupled with dynamical spacetime in the fully non-linear regime is the migration of an unstable neutron star [24,6,15,12]. Following [15], we use model SU for this migration test, which lies on the unstable branch of the massradius curve. As the star evolves and migrates to the corresponding stable configuration ρ c = 1.346 × 10 −3 with the same mass, the radius of the star expands to a large value. In this test, we setup a 1D run with a simulation box r = [0, 30] with 1024 grid points. Unlike the previous cases, we adopt the ideal gas (gamma-law) equation of state P = (Γ − 1)ρ with K = 100 and Γ = 2 for the fluid so that we can also capture the shock heating effect. Figure 9 shows the evolution of the central density ρ c as a function of time. The oscillations are damped due to the fact that shock waves are formed at every pulsation and some kinetic energy is dissipated into thermal energy. Our results agree with the T max = 20 ms equatorial initial equatorial Figure 8: Comparison between the initial density profiles and the normalized rotational velocity (solid lines) and the some quantities (dashed lines) at T max = 20 ms of a stable rotating neutron star BU8. The left panel shows a comparison between the initial normalized density profiles along the polar and equatorial radii and the same quantities at T max . The right panel compares the rotational velocity profiles at times t = 0 and t = T max . Due to the fact that the HLLE solver is too diffusive to deal with the surface of a rapidly rotating star, the angular velocity profile is slightly distorted at the star surface.

Conclusion
We present the methodology and implementation of Gmunu, a new general-relativistic hydrodynamics code which makes use of cell-centered nonlinear multigrid methods to solve the elliptic-type metric equations in the extended conformally flat condition (xCFC) approximation to general relativity. The set of hydrodynamics equations are solved with standard high-resolution shock-capturing schemes. Four different Riemann solvers have been implemented in the code: TVDLF [53], HLL [30], HLLE [23,29], and Marquina flux formula [21]. For the cell-interface reconstruction, various options are also available: PC, MC, WENO5, and MP5.
We have tested Gmunu with some benchmarking tests for relativistic hydrodynamics codes such as the relativistic shocktube problem, the evolution of rapidly rotating neutron stars, and the migration from an unstable TOV star to a corresponding stable solution with the same mass.
The main novelties of Gmunu are the following: • Although the system is highly nonlinear and fully coupled, our multigrid solver is robust and converges rapidly. In practical use, the computational time needed for solving the (elliptic-type) metric equations is comparable to the hydrodynamics step.
• In contrast to the code presented in [19] where the hydrodynamic and metric variables are defined in two different grids, our multigrid metric solver uses the same grid as the hydrodynamics sector. This avoids the need to perform interpolations for the variables defined in two different grids.
• Even using spherical polar coordinates, besides standard boundary conditions, our code does not require special treatment or regularization near the origin and pole axis.
We have demonstrated that multigrid method is an efficient strategy to solve nonlinear elliptic metric equations in hydrodynamical simulations. It seems to be promising that the method would make fully-constrained evolution scheme in numerical relativity become more affordable computationally. In the future, we shall extend Gmunu to a fully-constrained scheme in exact general relativity such as the formulation of Bonazzola et al. [8]. While numerical-relativity codes based on a free-evolution approach are standard choices for modelling dynamical spacetimes, it is still challenging for these codes to preform stable and accurate long-term evolutions of compact objects. It would be interesting to see whether a fully-constrained evolution code could improve the situation without requiring much more computational resources in the future.

Acknowledgements
PCKC thanks Elias Most for the useful discussion about the details of conserved to primitive variables, and Ninoy Rahman for the detailed discussion about the special treatment at the centre of the simulation box. This work was partially supported by grants from the Research Grants Council of the Hong Kong (Project No. CUHK 14310816 and CUHK 24304317), the Croucher Innovation Award from the Croucher Fundation Hong Kong and by the Direct Grant for Research from the Research Committee of the Chinese University of Hong Kong. Table B1: Initial conditions for the relativistic shocktube problem.
x < 0.5 x > 0.5 ρ = 10 ρ = 1 P = 13. Einstein equations, the diffusion at the contact discontinuity of the neutron star surface triggers the natural oscillation modes which were well studied for this kind of polytropic star [51]. The oscillation modes of a neutron star can be obtained approximately by perturbation calculations in the Cowling approximation [56,51,57]. They can also be extracted by nonlinear hydrodynamical simulation with a fixed background metric. By comparing the simulated results with the oscillation modes obtained in the Cowling approximation, we can thus focus and check the correctness of our hydrodynamics solver on a fixed background spacetime. In this test, we simulate the stable spherically symmetric neutron star BU0 model as mentioned in the table 1. We setup a 1D run with the resolution n r × n θ = 640 × 1 and keep the background metirc fixed during the simulation.
The upper panel of figure B2 shows the relative variation of the central density ρ c as a function of time. The relative variation of the density is of the order 10 −4 . The oscillation modes extracted from the simulation also agree very well with those obtained in the Cowling approximation.