have matured into standard numerical methods to simulate a wide variety of partial differential equations. In the context of numerical relativityBaumgarte and Shapiro (2010), discontinuous Galerkin methods have shown advantages for relativistic hyperbolic problems over traditional discretization methods such as finite difference, finite volume and spectral finite elements Field et al. (2010); Brown et al. (2012); Field et al. (2009); Zumbusch (2009); Radice and Rezzolla (2011); Mocz et al. (2014); Zanotti et al. (2014); Endeve et al. (2015); Teukolsky (2016); Bugner et al. (2016); Schaal et al. (2015); Zanotti et al. (2015); Miller and Schnetter (2017); Kidder et al. (2017); Fambri et al. (2018) by combining the best aspects of all three methods. As computers reach exa-scale power, new methods like DG are needed to tackle the biggest problems in numerical relativity such as realistic supernovae and binary neutron-star merger simulations on these very large machines Kidder et al. (2017).
DG efforts in numerical relativity have so far targeted evolutionary problems Field et al. (2010); Brown et al. (2012); Field et al. (2009); Zumbusch (2009); Radice and Rezzolla (2011); Teukolsky (2016); Bugner et al. (2016); Schaal et al. (2015); Zanotti et al. (2015); Miller and Schnetter (2017); Fambri et al. (2018). Radice and Rezzolla Radice and Rezzolla (2011) showed for spherically symmetric problems that the discontinuous Galerkin method could handle strong relativistic shock waves while maintaining exponential convergence rates in smooth regions of the flow. Teukolsky Teukolsky (2016) showed how to develop discontinuous Galerkin methods for applications in relativistic astrophysics. Bugner et al. Bugner et al. (2016) presented the first three-dimensional simulations of general relativistic hydrodynamics with a fixed spacetime background using a discontinuous Galerkin method coupled with a WENO algorithm. Kidder et al. Kidder et al. (2017) showcased the first discontinuous Galerkin code to use a task-based parallelism framework for applications in relativisitic astrophysics. Kidder et al. tested the scalability and convergence of the code on relativistic magnetohydrodynamics problems. Miller et al. Miller and Schnetter (2017) developed a discontinuous Galerkin operator method for use in finite difference codes and used it to solve gauge wave problems involving the BSSN formulation of the Einstein field equations. Hebert et al. Hébert et al. (2018) presented the first discontinuous Galerkin method for evolving neutron stars in full General Relativity. Finally, Fambri et al. Fambri et al. (2018) used an ADER discontinuous Galerkin scheme to solve general relativistic ideal magnetohydrodynamics problems in fixed spacetimes. They compared their DG method to a finite volume scheme and showed that DG is much more efficient.
This paper explores DG for elliptic problems in numerical relativity. We develop an elliptic solver with the following primary features: (i) It operates on curved meshes, with non-conforming elements. (ii) It supports adaptive h and p refinement, driven by a posteriori error estimators. (iii) It employs multi-grid for efficient solution of the resulting linear systems. (iv) It uses compactified domains to treat boundary conditions at infinity. While each of these features has appeared in the literature beforeArnold et al. (2002); Stiller (2017); Hesthaven and Warburton (2008); Kronbichler and Wall (2018); Fick (2014); Kozdon and Wilcox (2018); Kozdon et al. (2019), to our knowledge, our solver for the first time combines all these elements simultaneously and demonstrate their effectiveness on difficult numerical problems.
Specifically, the present article is a step toward a solver for the Einstein constraint equations, which must be solved to construct initial data for the evolution of compact binary systems Pfeiffer (2005); Cook (2000); Baumgarte and Shapiro (2010). The constraint equations are generally rewritten as elliptic equations, and depending on detailed assumptions, this results in one or more coupled non-linear elliptic partial differential equations. Construction of initial data is arguably the most important elliptic problem in numerical relativity, but not the only one: Elliptic equations also occur in certain gauge conditions Baumgarte and Shapiro (2010) or for implicit time-stepping to alleviate the computational cost of high-mass-ratio binaries Lau et al. (2009, 2011) .
The motivation for developing a new solver is multi-fold. First, current spectral methods have difficulty obtaining certain initial data sets, such as binaries at short separation containing a neutron star, where the neutron star has high compactness and a realistic equation of state Henriksson et al. (2016). Furthermore, there is a need for a solver which can routinely obtain high-accuracy. Errors from inaccurate initial data sets creep into the evolutions with sizeable effect: Ref. Tsokaros et al. (2016) shows that despite only global (local) differences of () in the initial data of the two codes COCAL and LORENE for irrotational neutron-star binaries, the gravitational wave phase at the merger time differed by 0.5 radians after 3 orbits. The design of a more accurate code requires adaptive mesh refinement, load-balancing and scalability which a DG code potentially can provide. The present work also complements the DG evolution code presented in Ref. Kidder et al. (2017), leading to a complete framework for solving both elliptic and hyperbolic PDEs in numerical relativity.
The organization of the paper is as follows: Section II presents the components of our discontinuous Galerkin code. Section III showcases our hp-adaptive multigrid solver on increasingly challenging test-problems, each illustrating the power of the discontinuous Galerkin method. This section includes a solution for the Einstein constraint equations in the case of a constant density star. This problem has a surface discontinuity which mimicks Neutron stars with phase transitions. Lastly, this section also presents solutions for initial-data of two orbiting non-spinning black-holes to showcase and compare it with the elliptic solver Spells Pfeiffer et al. (2003a) as well as solutions for the initial data of three black-holes with random locations, spins and momenta. We close with a discussion in Sec. IV.
Ii Numerical Algorithm
ii.1 DGFEM discretization
We base our nomenclature on Arnold et al. (2002); Sherwin et al. (2006); Di Pietro and Ern (2011); Fick (2014); Kozdon and Wilcox (2018); Kozdon et al. (2019), which we summarize here for completeness.
Our purpose is to solve elliptic equations in a computational domain in dimensions. To allow for non-trivial topology and shape of , we introduce a macro-mesh (also known as a multi-block mesh) consisting of macro-elements (also known as blocks) such that (1) the macro-elements cover the entire domain, i.e., , (2) the macro-elements touch each other at complete faces and do not overlap; and (3) each is the image of the reference cube under a mapping . As an example, Fig. 1 shows a macro-mesh of five elements covering a two-dimensional disk. The macro-mesh represents the coarsest level of subsequent mesh-refinement.
The macro-mesh is refined by subdividing macro-elements into smaller elements along faces of constant reference cube coordinates. Refinement can be multiple levels deep (i.e. refined elements can be further subdivided), and we do not assume uniform refinement. The refined mesh is referred to as the fine mesh, which is exemplified in panel (b) of Fig. 1. In contrast to the macro-mesh, the fine mesh will generally be non-conforming at element boundaries, both within one macro-element and at boundaries between macro-elements. We assume that there is at most a 2:1 refinement difference at any element boundary, i.e. the boundary of a coarse element faces at most two smaller elements (per dimension). Each face of each element is endowed with an outward-pointing unit-normal . The map in macro-element induces a map on each fine element within , denoted
The reference coordinates of are linearly related to the reference coordinates of .
Turning to boundaries, we define the set of all element boundaries (internal and external), , where is the boundary of element e of the fine mesh. is called the mortar, and it decomposes into a finite set of -dimensional mortar elements , arising from the faces of each mesh element , s.t. and each mortar element intersects at most at the boundary of two elements. splits into interior mortar elements and exterior mortar elements , with . Similarly, splits into interior mortar and exterior mortar , cf. Fig. 1, which intersect at -dimensional edges where the interior mortar touches the exterior boundary. We partition the exterior mortar elements further into elements where we apply Dirichlet boundary conditions, , Neumann boundary conditions, , and Robin boundary conditions, , respectively. This induces sets of boundary points via , and . We assume , and are disjoint, i.e. one type of boundary condition is employed on each connected part of the boundary. Finally, we introduce two definitions to help us simplify equations later in the text. Firstly, because internal boundaries and external Dirichlet boundaries are often treated similarly, it is convenient to define and . Secondly, we refer to the set of mortar elements surrounding an element as .
ii.1.2 Weak Equations
On the fine mesh , we wish to discretize the following model elliptic problem with the symmetric interior penalty discontinuous Galerkin method Arnold et al. (2002):
In Eqs. (2) and in the following, we employ the Einstein sum convention, so that the first term in Eq. (2a) represents , and similar for the left hand side of Eqs. (2c) and (2d). In accordance with the underlying ideas of DG, we seek to approximate the solution of Eqs. (2) with polynomials on each element , without strictly enforcing continuity across elements. Let denote the set of polynomials on the reference cube up to degree mapped to element . We assume the same maximum polynomial order along each dimension; it is straightforward to extend to different polynomial order along different dimensions. The functions in are understood to be extended by outside of (i.e. on other elements). The global function space is the direct sum of the per-element polynomial spaces,
where the polynomial order may vary between elements. Because neighboring element touch, on internal boundaries the discretized solution will be represented twice on touching elements with generally different values on either element (this is origin of the term ’discontinuous’ in DG).
The discretized solution is determined, such that the residual of Eqs. (2a) is orthogonal to the function space . Within the symmetric interior penalty discontinuous Galerkin discretization Arnold et al. (2002); Di Pietro and Ern (2011), this yields
For internal boundaries, the operators and are defined by
Here is a scalar function, ’+’ and ’-’ is an arbitrary labeling of the two elements and touching at the interface, and and are the function values and the outward pointing unit normal on the two elements that share the interface. These operators are extended to external boundaries by
For later reference, we will refer to the first four integrals in Eq. (II.1.2) as the stiffness, consistency, symmetry and penalty integrals respectively and denote them , , and so that we may write
ii.1.3 Basis Functions
So far, we have not yet specified a concrete basis for the polynomial spaces introduced in Sec. II.1.2. We do so now. Recall that curvilinear elements are mapped to reference cubic elements . That is, each point corresponds to a reference cube coordinate , cf. Eq. (1).
Along each dimension , where the subscript denotes dimension, we first choose Legendre-Gauss-Lobatto collocation points . The index identifies the point along dimension . We take so that the points with and , that lie on the faces of the cube, are always collocation points. The collocation points in all
dimensions form a tensor product grid of-dimensional collocation points that we index by
that identifies a point regardless of dimension. With respect to these collocation points we can now construct the one-dimensional interpolating Lagrange polynomials
and employ their tensor product to define the -dimensional basis functions
When evaluating in the physical coordinates , one uses to map the function argument. For instance, the set of test-functions on element becomes
Furthermore, because the form a nodal basis, the expansion coefficients are the function values at the nodal points , and each test-function can be written as
ii.1.4 Semi-Discrete Global Matrix Equations
The global solution over the entire mesh is the direct sum of the solutions on each element, that is
Thus, with a chosen global ordering of elements and local ordering of basis functions , we may prescribe a global index to each of the local expansion coefficients . With this, we may now turn Eqs. (II.1.2) and (II.1.2) into a global linear system of equations. Let , and , then the global linear system is
Instead of forming the full global matrix
and performing the matrix-vector operatorover the global space, we can restrict the integrals in (II.1.2) and (II.1.2) to those pertaining to element and the mortar elements of , and perform elemental matrix-vector operations.
with an inverse given by
The Jacobian determinant is
Similarly, the surface Jacobian for a mortar element is
where is the index of the reference coordinate which is constant on the surface under consideration (no sum-convention). For a mortar on a macro-element boundary, may change depending on which of the two macro-element mappings are used to compute Eqn. (23). This ambiguity is not problematic as long as all quantities in a mortar integrand are transformed to the macro-element coordinate system in which is being computed. The normal is computed by
where sgn is the signum function, where labels the dimension that is constant on the face under consideration (no sum-convention).
Let us now consider the stiffness integral in Eq. (II.1.2). Because basis-functions are local to each element, we can consider each element individually, and we will omit the superscript to lighten the notational load. Substituting the expansion of the solution , as well as , the stiffness integral becomes
We recall that we employ the Einstein sum convention, i.e. Eq. (25) has implicit sums over and . The derivatives are just polynomials in , therefore, they can be re-expanded in our nodal basis as with . The physical derivative is then . The tensor-grid in Eq. (15) implies that the matrices are sparse for . Using these expressions, and converting to an integral in the reference coordinates, we obtain
We evaluate the integral in Eq. (26) with Gauss-Legendre (GL) quadrature. This choice follows Mengaldo et al. (2015) in the use of a stronger set of quadrature points (higher order GL grids) to decrease the error in geometric aliasing. We denote the GL abscissas and weights by and , respectively, where . In multiple dimensions these will be a tensor product of the 1-d GL abscissas and weights denotes by and , respectively. All non-polynomial functions, including geometric quantities such as , are evaluated directly at , whereas polynomial functions like the trial function and the test function —which naturally are represented on the Legendre-Gauss-Lobatto grid— are interpolated to the GL–quadrature points. Denoting the interpolation matrix by , we find
Here is an interpolation matrix from the GLL points to the GL points. The other volume integrals in Eq. (II.1.2) are computed in a similar manner.
The mortar integrals are a bit more involved owing to the extra book-keeping arising from the two elements (named ’+’ and ’-’) that touch at the boundary , each with their own local basis-functions, denoted by and . Taking the penalty-integral as an example, for test-functions that are a basis-function of the ’-’ element, the definition Eq. (7) yields
Equation (31) contains both and ; therefore, when substituting in their respective local expansions, , and pulling the coefficients outside the integral, we see that the penalty term will result in entries of the global matrix equation (19) that couple the two elements . Once the coefficients are moved outside the integral, the remaining integral only depends on basis-functions and geometric quantities. These integrals are evaluated with GL-quadrature, which is expressed in terms of interpolation matrices , where labels the GL collocation points on the boundary and labels the local basis-functions in , which are to be evaluated in the suitably mapped locations.
ii.2 Penalty Function
where and represent a typical polynomial degree and a typical length-scale of the elements touching at the boundary, respectively, and where the parameter is large enough such that is coercive. We choose and we set as this has yielded the best results empirically. Here is the volume Jacobian on mapping to , and is the surface Jacobian mapping to the boundary face of .
Throughout the following sections we will use the energy norm,
and the norm,
Here, is a scalar function on .
ii.4 Multigrid-Preconditioned Newton-Krylov
In general, a set of elliptic partial differential equations can be written in the form
where is the solution and R is a non-linear elliptic operator. Defining the Jacobian of the system by , Newton-Raphson iteratively refines an initial guess for the solution by solving the following linear system for the correction :
Upon discretization as described in Sec. II.1, Eq. (36) results in a linear system. Once Eq. (36) is solved, the improved solution is given by . The Newton-Raphson iterations are continued until the residual is sufficiently small. At each step, the linear system is solved by a linear solver we abbreviate as . The full algorithm is described in Algorithm 1.
For three-dimensional problems with a large number
of degrees of freedom, the Jacobianis too large to form fully and partition across processes. For building a scalable solver, we are therefore limited to iterative solvers, the most popular class of which are the Krylov solvers Saad (2003) such as the Conjugate Gradient method or GMRES, which only involve matrix-vector operations. In this paper we use the flexible conjugate gradient Krylov Notay (2000) solver at each Newton-Raphson step, as summarized in Algorithm 2, where denotes the initial guess of the linear solver, and is the number of iterations.
A typical Krylov solver will take iterations to reach a desired accuracy, where the condition number
is defined as the ratio of the maximum eigenvalue and minimum eigenvalue of. For the discontinuous galerkin method, the discretized Laplacian matrix has a condition number that grows with p-refinement and h-refinement (see Antonietti and Houston (2011); Hesthaven and Warburton (2008) for example) and therefore the number of iterations will grow with each AMR step unless Eq. (36) is preconditioned.
ii.4.2 Multigrid preconditioner
We use a multigrid V-cycle Briggs et al. (2000) as a preconditioner for each FCG solve (called on line 4 in Algorithm 2). Multigrid is an multi-level iterative method aimed at achieving mesh-independent error reduction rates through a clever method of solving for the error corrections on coarser grids and prolonging the corrections to the finer grids. For a detailed overview of multigrid methods see Briggs et al. (2000). The main drawback of multigrid is its complexity. In this section we will briefly describe the components of the multigrid algorithm employed in solving problems in this paper with hp-grids in parallel with the interior penalty discontinuous Galerkin method.
ii.4.3 Multigrid Meshes
Multigrid uses a hierarchy of coarsened meshes, labeled with , where represents the fine mesh on which the solution is desired, and represents the coarsest mesh. One V-cycle proceeds from the fine grid to the coarsest grid, and back to the fine grid, as indicated in Fig. 3. We construct the coarse meshes by successively coarsening. Coarsening by combining elements into one element can lead to interfaces with a 4:1 balance, even if the original mesh is 2:1 balanced. We avoid such 4:1 balance with surrogate meshes, which have been used before in large-scale FEM multigrid solvers, see for example Sampath et al. (2008); Sundar et al. (2012). A surrogate mesh is the naively h-coarsened mesh, as indicated in blue in Fig. 3. If the surrogate mesh has indeed interfaces with 4:1 balance, the desired 2:1 balance is enforced by refining at the unbalanced interfaces. This results in the coarsened mesh on Level , and iteratively down to . It is easy to show that at each coarsening the coarse level mesh is strictly coarser than the level fine mesh. Following this, we must make a decision as to what polynomial degree we choose for the coarsened elements. We take the minimum polynomial degree of the children elements for each parent element of the coarsened mesh, ensuring that functions on the coarse grid can be represented exactly on the fine mesh, a property we will utilize below in Eq. (42). Lastly, it is important to note that we never purely p-coarsen on any of the multigrid levels. We do not store the multigrid meshes because we have empirically shown that the coarsening, refining and balancing is not the bottleneck for the multigrid algorithm.
ii.4.4 Multigrid algorithm
Having constructed the coarse meshes, we can now turn to the multigrid algorithm, summarized in Algorithm 3. It consists of several important components most importantly smoothing and inter-mesh operators like coarsening, balancing of the mesh, restriction and prolongation. We now look at each of these in turn.
ii.4.5 Multigrid Smoother
We explore two different smoothers, Chebyshev smoothing and a Schwarz smoother. Both of these avoid explicit storage of the matrix, as required by smoothers like the Gauss-Seidel method Shang (2009).
Chebyshev smoothing Li et al. (2011), a type of polynomial smoother Saad (2003) requires only matrix-vector operations and has been shown to work satisfactorily in scalable multigrid solvers Ghysels et al. (2012). Our implementation of Chebyshev smoothing Li et al. (2011) is presented in Algorithm 4. For this algorithm there are two user-defined parameters: and . Typical values we use for are in the range and is usually set in the range .
Chebyshev polynomial smoothers require a spectral bound on the eigenvalues of the linear operator. We use conjugate gradients to estimate the eigenvalues of a general symmetric linear operator, as implemented in Algorithm 5. It can be shown Scales (1989) that each iteration of conjugate gradients obtains one row of the underlying linear operator in tri-diagonal form:
The values and in this expression are obtained from the conjugate gradients algorithm 5 on lines 6 and 9. Furthermore, the eigenvalues of each of the sub tri-diagonal matrices is a subset of the eigenvalues of the full matrix. So at each iteration we can get an estimate of the bound by using the Gershgorin circle theorem Bell (1965). Our Alg. 5 combines the CG steps with the estimation of the bound of the largest eigenvalue.
Besides estimating the spectral radius, we also use the CG iterations to further smooth the solution. This adds robustness to multigrid algorithms Elman et al. (2001).
While the Chebyshev smoother is easy to implement, it is reliant on a robust estimate of the largest eigenvalue and this may not be always true in our case. Thus we also implement an additive Schwarz smoother in the manner of Stiller (2017), which is much more robust. The Schwarz smoother works by performing local solves on element-centered subdomains and then adding a weighted sum of these local solves to obtain the smoothed solution. A simple example of one of these subdomains (where the grid is both p and h-uniform) is shown in Figure 4. More generally, the Schwarz subdomain centered on element of the mesh is constructed as follows: Starting with all collocation points on , one adds all boundary points of neighboring elements which coincide with faces, vertices or corners of . Around these collocation points, one then adds layers of additional collocation points (in Fig. 4, ). If the mesh has non-uniform or refinement, the resulting set of collocation points will have ragged boundaries.
The solutions on the individual Schwarz subdomains are combined as a weighted sum. The weights differ with each collocation point of each subdomain. In 1-D, for a Schwarz subdomain centered on element with left and right neighbours and we define an extended LGL coordinate for a collocation point as follows: