The new discontinuous Galerkin methods based numerical relativity program Nmesh

In Nmesh, we use a discontinuous Galerkin (DG) method to discretize evolution equations. Often, these evolution equations come from general relativistic conservation laws of the form

$$\nabla_{\mu}J^{\mu}=S,$$

(1)

where $S$ is a possible source term. The covariant divergence on the left hand side can be written in terms of coordinate derivative using $\nabla_{\mu}J^{\mu}=\frac{1}{\sqrt{|g|}}\partial_{\mu}(\sqrt{|g|}J^{\mu})$, where $g$ is the determinant of the 4-metric $g_{\mu\nu}$. In terms of the standard 3+1 decomposition [53], the 4-metric is written as

$$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=-\alpha^{2}dt^{2}+\gamma_{ij}(dx^{i}+\beta^{i}dt)(dx^{j}+\beta^{j}dt),$$

(2)

where $\gamma_{ij}$ is the spatial metric on $t=const$ slices, and $\alpha$ and $\beta^{i}$ are called lapse and shift. We can show that $\sqrt{|g|}=\alpha\sqrt{\gamma}$, where $\gamma$ is the determinant of the 3-metric $\gamma_{ij}$. Thus, Eq. (1) is equivalent to

$$\partial_{t}(\sqrt{\gamma}\alpha J^{t})+\partial_{i}(\sqrt{\gamma}\alpha J^{i})=\sqrt{\gamma}\alpha S.$$

(3)

Usually, we introduce the new variables

$$u=\sqrt{\gamma}\alpha J^{t},\ \ \ f^{i}=\sqrt{\gamma}\alpha J^{i},\ \ \ s=\sqrt{\gamma}\alpha S,$$

(4)

so that Eq. (3) finally yields

$$\partial_{t}u+\partial_{i}f^{i}=s.$$

(5)

Note that the flux vector $f^{i}$ is usually a function $f^{i}(u)$, that depends on $u$. For brevity, we omit this dependence in most equations.

To discretize the spatial derivatives in Eq. (5), we first integrate against a test function $\psi$ with respect to the coordinate volume $d^{3}x$ over a certain region $\Omega$. We obtain

$$\int(\psi\partial_{t}u+\psi\partial_{i}f^{i})d^{3}x=\int\psi sd^{3}x.$$

(6)

The second term is now integrated by parts using

$$\int\psi\partial_{i}f^{i}d^{3}x=\oint\psi f^{i}n_{i}d^{2}\Sigma-\int f^{i}\partial_{i}\psi d^{3}x,$$

(7)

where $d^{2}\Sigma$ is the surface element for integrating over the boundary $\partial\Omega$, and $n_{i}$ is normal to $\partial\Omega$. In the surface integral, the flux $f^{i}n_{i}$ at the boundary appears. To incorporate numerical boundary conditions, $f^{i}n_{i}$ in the surface integral over the boundary is replaced by the so-called numerical flux $(f^{i}n_{i})^{*}$ (See Sec. 2.3 below). It contains any information we need from the other side of the boundary (such as incoming characteristic modes). The replacement of $f^{i}n_{i}$ by $(f^{i}n_{i})^{*}$ in the surface integral yields

	$\displaystyle\int\psi\partial_{i}f^{i}d^{3}x$	$\displaystyle\rightarrow$	$\displaystyle\oint\psi(f^{i}n_{i})^{*}d^{2}\Sigma-\int f^{i}\partial_{i}\psi d^{3}x$		(8)
		$\displaystyle=$	$\displaystyle\oint\psi[(f^{i}n_{i})^{*}-f^{i}n_{i}]d^{2}\Sigma+\int\psi\partial_{i}f^{i}d^{3}x,\ \ \$

where in the last step we have used integration by parts again to eliminate derivatives of $\psi$. With this replacement, Eq. (6) becomes

$$\int(\psi\partial_{t}u+\psi\partial_{i}f^{i})d^{3}x=\int\psi sd^{3}x-\oint\psi[(f^{i}n_{i})^{*}-f^{i}n_{i}]d^{2}\Sigma.$$

(9)

We now introduce a coordinate transformation

$$x^{i}=x^{i}(x^{\bar{i}}),$$

(10)

such that the volume we integrate over extends from $-1$ to $+1$ for all three $x^{\bar{i}}$ coordinates. In these coordinates, Eq. (9) reads

$$\int\psi\left(\partial_{t}u+\frac{\partial x^{\bar{i}}}{\partial x^{i}}\partial_{\bar{i}}f^{i}\right)Jd^{3}\bar{x}=\int\psi sJd^{3}\bar{x}-\oint\psi Fd\bar{A},$$

(11)

where we have defined

$$F:=(f^{i}n_{i})^{*}-f^{i}n_{i},$$

(12)

$$J=\left|\det\left(\frac{\partial x^{i}}{\partial x^{\bar{i}}}\right)\right|,$$

(13)

where $J$ is called the Jacobian, and $d\bar{A}$ is the surface element on one of the six surfaces $x^{\bar{i}}=\pm 1$, but now expressed in $x^{\bar{i}}$ coordinates. For example on $x^{\bar{3}}=-1$,

$$d\bar{A}=\sqrt{{}^{(2)}\bar{\gamma}}dx^{\bar{1}}dx^{\bar{2}},$$

(14)

where ${}^{(2)}\bar{\gamma}$ is the determinant of the 2-metric induced on this coordinate surface by the flat 3-metric $\delta_{ij}$. The flat $\delta_{ij}$ results from our choice to integrate over the coordinate volume $d^{3}x$ in Eq. (6), without including $\sqrt{\gamma}$. Below these $x^{i}$ coordinates will be chosen to be Cartesian-like, so that they can cover all numerical domains.

Next, we expand both $u$ and $f^{i}$ in terms of Lagrange’s characteristic polynomials

$$l_{q}(\bar{x})=\prod_{r=0,r\neq q}^{N}\frac{\bar{x}-\bar{x}_{r}}{\bar{x}_{q}-\bar{x}_{r}}$$

(15)

so that, e.g.,

$$u(\bar{x})=\sum_{r_{1}=0}^{N}\sum_{r_{2}=0}^{N}\sum_{r_{3}=0}^{N}u_{r_{1}r_{2}r_{3}}l_{r_{1}}(x^{\bar{1}})l_{r_{2}}(x^{\bar{2}})l_{r_{3}}(x^{\bar{3}}).$$

(16)

The $\bar{x}_{r}$ are $N+1$ grid points that we choose in the interval $[-1,1]$. For the test function $\psi$, we use these same basis polynomials, i.e.,

$$\psi=l_{q_{1}}(x^{\bar{1}})l_{q_{2}}(x^{\bar{2}})l_{q_{3}}(x^{\bar{3}}).$$

(17)

The final step is to approximate all integrals in Eq. (11) using Gaußian quadrature (specifically Lobatto’s Integration formula on p. 888 of [54]), which, in one dimension, is given by

$$\int_{-1}^{1}d\bar{x}\ g(\bar{x})\approx\sum_{q=0}^{N}w_{q}g(\bar{x}_{q}).$$

(18)

Here we use the $N+1$ Legendre Gauß-Lobatto grid points $\bar{x}_{q}$, that are defined as the extrema of the standard Legendre polynomial $P_{N}(\bar{x})$ in $\bar{x}\in[-1,1]$. The integration weights are then given by [54]

$$w_{q}=\frac{2}{N(N+1)P_{N}(\bar{x}_{q})^{2}}.$$

(19)

The integrals in Eq. (11) then turn into sums over products of $l_{p}(\bar{x})$, or products of $l_{p}(\bar{x})$ and its derivative. If we define

$$(u,v):=\sum_{r=0}^{N}w_{r}u(\bar{x}_{r})v(\bar{x}_{r}),$$

(20)

the products we encounter are $(l_{q},l_{r})$ and $(l_{q},\partial_{\bar{x}}l_{r})$. Since

$$l_{q}(\bar{x}_{r})=\delta_{qr},$$

(21)

we find

$$(l_{q},l_{r})=w_{q}\delta_{qr}$$

(22)

and

$$(l_{q},\partial_{\bar{x}}l_{r})=w_{q}\partial_{\bar{x}}l_{r}(\bar{x}_{q})=:w_{q}D_{qr}^{\bar{x}}.$$

(23)

The differentiation matrix is given by

$$D_{qr}^{\bar{x}}=\frac{c_{r}}{c_{q}}\frac{1}{\bar{x}_{q}-\bar{x}_{r}},$$

(24)

where the Lagrange interpolation weights are

$$c_{q}=\left[\prod_{s=0,s\neq q}^{N}(\bar{x}_{q}-\bar{x}_{s})\right]^{-1}.$$

(25)

Putting all this together Eq. (11) becomes

	$\displaystyle w_{q_{1}}w_{q_{2}}w_{q_{3}}J_{q_{1}q_{2}q_{3}}\Big{(}\partial_{t}u_{q_{1}q_{2}q_{3}}+\frac{\partial x^{\bar{1}}}{\partial x^{i}}\sum_{r=0}^{N}D_{q_{1}r}^{\bar{1}}f^{i}_{rq_{2}q_{3}}$
	$\displaystyle+\frac{\partial x^{\bar{2}}}{\partial x^{i}}\sum_{r=0}^{N}D_{q_{2}r}^{\bar{2}}f^{i}_{q_{1}rq_{3}}+\frac{\partial x^{\bar{3}}}{\partial x^{i}}\sum_{r=0}^{N}D_{q_{3}r}^{\bar{3}}f^{i}_{q_{1}q_{2}r}\Big{)}$
	$\displaystyle=w_{q_{1}}w_{q_{2}}w_{q_{3}}J_{q_{1}q_{2}q_{3}}s_{q_{1}q_{2}q_{3}}$
	$\displaystyle-w_{q_{2}}w_{q_{3}}F_{q_{1}q_{2}q_{3}}\sqrt{{}^{(2)}\bar{\gamma}_{q_{1}q_{2}q_{3}}}(\delta_{q_{1}0}+\delta_{q_{1}N})$
	$\displaystyle-w_{q_{1}}w_{q_{3}}F_{q_{1}q_{2}q_{3}}\sqrt{{}^{(2)}\bar{\gamma}_{q_{1}q_{2}q_{3}}}(\delta_{q_{2}0}+\delta_{q_{2}N})$
	$\displaystyle-w_{q_{1}}w_{q_{2}}F_{q_{1}q_{2}q_{3}}\sqrt{{}^{(2)}\bar{\gamma}_{q_{1}q_{2}q_{3}}}(\delta_{q_{3}0}+\delta_{q_{3}N}),$		(26)

where the Kronecker deltas on the right hand side come from $l_{q}(+1)=\delta_{qN}$ and $l_{q}(-1)=\delta_{q0}$. We now divide Eq. (2.1) by $w_{q_{1}}w_{q_{2}}w_{q_{3}}J_{q_{1}q_{2}q_{3}}$ and use the fact that

$$\sqrt{{}^{(2)}\bar{\gamma}}/J=\sqrt{\bar{\gamma}^{\bar{i}\bar{i}}}$$

(27)

holds on the surface $x^{\bar{i}}=const$. Here $\bar{\gamma}^{\bar{i}\bar{i}}$ is a diagonal component of the inverse 3-metric obtained by transforming the flat 3-metric $\delta_{ij}$ from $x^{i}$-coordinates to $x^{\bar{i}}$-coordinates. This results in

	$\displaystyle\partial_{t}u_{q_{1}q_{2}q_{3}}+\sum_{r=0}^{N}\Big{(}\frac{\partial x^{\bar{1}}}{\partial x^{i}}D_{q_{1}r}^{\bar{1}}f^{i}_{rq_{2}q_{3}}+\frac{\partial x^{\bar{2}}}{\partial x^{i}}D_{q_{2}r}^{\bar{2}}f^{i}_{q_{1}rq_{3}}+\frac{\partial x^{\bar{3}}}{\partial x^{i}}D_{q_{3}r}^{\bar{3}}f^{i}_{q_{1}q_{2}r}\Big{)}$
	$\displaystyle=s_{q_{1}q_{2}q_{3}}-\frac{\sqrt{\bar{\gamma}^{\bar{1}\bar{1}}_{q_{1}q_{2}q_{3}}}}{w_{q_{1}}}F_{q_{1}q_{2}q_{3}}(\delta_{q_{1}0}+\delta_{q_{1}N})$
	$\displaystyle-\frac{\sqrt{\bar{\gamma}^{\bar{2}\bar{2}}_{q_{1}q_{2}q_{3}}}}{w_{q_{2}}}F_{q_{1}q_{2}q_{3}}(\delta_{q_{2}0}+\delta_{q_{2}N})$
	$\displaystyle-\frac{\sqrt{\bar{\gamma}^{\bar{3}\bar{3}}_{q_{1}q_{2}q_{3}}}}{w_{q_{3}}}F_{q_{1}q_{2}q_{3}}(\delta_{q_{3}0}+\delta_{q_{3}N}),$		(28)

which is the version we use in Nmesh’s DG method. Notice that the derivation of this DG method has mostly followed the one introduced by Teukolsky in [48, 55] and tested extensively in [40], except for one important difference. Since we integrate over $d^{3}x$ without including the determinant of the physical metric, we use the flat metric $\delta_{ij}$ when we construct $\bar{\gamma}^{\bar{i}\bar{i}}$ or when we normalize the normal vector $n_{i}$. This same $n_{i}$ also enters the calculation of the fluxes and eigenvalues discussed in subsection 2.3. Recall that $n_{i}$ arose in Eq. (7) after using Gauß’s theorem. As shown in A, Gauß’s theorem can be used with any metric, as long as we normalize $n_{i}$ with this same metric. The advantage of $\delta_{ij}$ is that it is constant and thus cannot have any discontinuities. Our approach simplifies the formalism as one does not have to worry about possible discontinuities in the physical metric or the normal vector across domain boundaries, which is an issue in the other approach [40]. In B we discuss the difference between both approaches for the well understood case of an advection equation. We find that our approach together with the flux of Eq. (33) yields the correct upwind result, while the approach in [40] does not, if the physical metric is discontinuous across the boundary. However, the true physical solution of the Einstein equations is a continuous metric, even in the presence of shocks in the matter. Thus we expect that both approaches will converge to the same result, because any discontinuities in the metric should converge to zero. We have also tried both approaches by evolving a black hole, where the physical metric is far from flat, and where discontinuities in the physical metric arise due to numerical errors. We find no important differences in the numerical solution or its rate of convergence. The discontinuities in the physical metric for this black hole case are described in Sec. 4.2 and shown in Fig. 8.

Let us also consider the field $u$ from Eq. (5) for the case where $s=0$. Then $\int ud^{3}x$ is exactly conserved. In formulations that directly use Eq. (3), and do not include $\sqrt{\gamma}$ in the field definition, the conserved quantity is $\int U\sqrt{\gamma}d^{3}x$, where $U=\alpha J^{t}$. In fact, for the exact solution of Eq. (5) both integrals are identical. However, once $u$ and $U$ are expressed in terms of basis functions (or equivalently in terms of values at grid points), the numerical (Gaußian quadrature) integrals over $u$ and $U$ will yield answers that differ at the level of the numerical truncation error. Nevertheless, a correct DG formulation will preserve these numerical integrals. Furthermore, any differences between the two numerical integrals will converge away with increasing resolution. Thus at any finite resolution the two approaches conserve different quantities, but in the continuum limit both converge to the same result.

So far, we have only considered evolution equations that can be written in conservative form as in Eq. (5), i.e., in terms of a flux vector $f^{i}$. However, the equations describing general relativistic gravity are often not available in this form. Rather they take the form

$$\partial_{t}u+A^{i}(u)\partial_{i}u=s.$$

(29)

Note that here the matrix $A^{i}(u)$ depends on $u$ itself. In this case, we can still integrate against a test function $\psi$, as before. The crucial introduction of a numerical flux in Eq. (8) now takes the form

	$\displaystyle\int\psi A^{i}\partial_{i}ud^{3}x$	$\displaystyle\rightarrow$	$\displaystyle\oint\psi(n_{i}A^{i}u)^{*}d^{2}\Sigma-\int u\partial_{i}(A^{i}\psi)d^{3}x$		(30)
		$\displaystyle=$	$\displaystyle\oint\psi[(n_{i}A^{i}u)^{*}-n_{i}A^{i}u]d^{2}\Sigma+\int\psi A^{i}\partial_{i}ud^{3}x.$

Thus, the surface integral has almost the same form, with $(n_{i}A^{i}u)^{*}$ playing the role of the numerical flux. If we again expand in Lagrange’s characteristic polynomials, and retrace our previous steps, we find the equivalent of Eq. (2.1). We obtain

	$\displaystyle\partial_{t}u_{q_{1}q_{2}q_{3}}+\sum_{r=0}^{N}A^{i}_{q_{1}q_{2}q_{3}}\Big{(}\frac{\partial x^{\bar{1}}}{\partial x^{i}}D_{q_{1}r}^{\bar{1}}u_{rq_{2}q_{3}}+\frac{\partial x^{\bar{2}}}{\partial x^{i}}D_{q_{2}r}^{\bar{2}}u_{q_{1}rq_{3}}+\frac{\partial x^{\bar{3}}}{\partial x^{i}}D_{q_{3}r}^{\bar{3}}u_{q_{1}q_{2}r}\Big{)}$
	$\displaystyle=s_{q_{1}q_{2}q_{3}}-\frac{\sqrt{\bar{\gamma}^{\bar{1}\bar{1}}_{q_{1}q_{2}q_{3}}}}{w_{q_{1}}}G_{q_{1}q_{2}q_{3}}(\delta_{q_{1}0}+\delta_{q_{1}N})$
	$\displaystyle-\frac{\sqrt{\bar{\gamma}^{\bar{2}\bar{2}}_{q_{1}q_{2}q_{3}}}}{w_{q_{2}}}G_{q_{1}q_{2}q_{3}}(\delta_{q_{2}0}+\delta_{q_{2}N})$
	$\displaystyle-\frac{\sqrt{\bar{\gamma}^{\bar{3}\bar{3}}_{q_{1}q_{2}q_{3}}}}{w_{q_{3}}}G_{q_{1}q_{2}q_{3}}(\delta_{q_{3}0}+\delta_{q_{3}N}),$

where $G$ is defined by

$$G:=(n_{i}A^{i}u)^{*}-n_{i}A^{i}u.$$

(32)

When we compare with the surface term $F$ defined in Eq. (12), appearing in the analog Eq. (2.1), we see that $G$ can be obtained from $F$ if we replace $f^{i}$ by $A^{i}u$.

In the interior of the domain, the flux vector $f^{i}(u)$ is simply computed from the field values $u$ in the interior. However, the numerical flux that is used in the surface integral over the boundary is computed from the field values on both sides of the boundary. In many cases, we will use the Rusanov or local Lax-Friedrichs (LLF) flux. It is given by

$$(f^{i}n_{i})^{*}=\frac{1}{2}\left[f^{i}(u_{\mathrm{in}})n_{i}+f^{i}(u_{\mathrm{adj}})n_{i}+|\lambda|_{\mathrm{max}}\left(u_{\mathrm{in}}-u_{\mathrm{adj}}\right)\right].$$

(33)

Here $n_{i}$ is the outward pointing normal to the boundary, $u_{\mathrm{in}}$ is the field value at the boundary using grid points that belong to the domain enclosed by the boundary, $u_{\mathrm{adj}}$ is the field value at the boundary using grid points that belong to the adjacent domain on the other side of the boundary, and $|\lambda|_{\mathrm{max}}$ is the absolute value of the eigenvalue of the characteristic mode with the largest eigenvalue magnitude, considering eigenvalues from both sides.

In the case where our system of equations takes the form of Eq. (29), we often also use another numerical flux, called the upwind flux. It is constructed from the orthonormalized eigenvectors and eigenvalues of the matrix $A^{i}n_{i}$ appearing in the surface term (32). Let the matrix $S$ contain the eigenvectors as its columns. Then we can write

$$A^{i}n_{i}=S\Lambda S^{-1},$$

(34)

where $\Lambda$ is a diagonal matrix that contains the corresponding eigenvalues. The eigenvectors with positive eigenvalues correspond to modes going along the direction on $n_{i}$, while the ones with negative eigenvalues correspond to modes going in the opposite direction. This means positive and negative eigenvalues are associated with modes that are outgoing and incoming through the boundary of the domain. As is usually the case, we wish to impose conditions only on the incoming modes. So, we define the upwind numerical flux that appears in Eq. (32), as

$$(n_{i}A^{i}u)^{*}=(S(\Lambda^{+}+\Lambda^{-})S^{-1}u)^{*}:=S(\Lambda^{+}S^{-1}u_{\mathrm{in}}+\Lambda^{-}S^{-1}u_{\mathrm{adj}}).$$

(35)

Here $\Lambda=\Lambda^{+}+\Lambda^{-}$ with $\Lambda^{+}$ and $\Lambda^{-}$ containing the positive and negative eigenvalues, and $u_{\mathrm{in}}$ and $u_{\mathrm{adj}}$ are the field values from the current domain and the adjacent domain.

To write Eq. (5) or (29), we use particular coordinates that are chosen to be Cartesian-like, and we call them $x^{i}=(x,y,z)$ in Nmesh. As already explained before we map these globally used Cartesian coordinates to local coordinates $x^{\bar{i}}$ via Eq. (10). This mapping is usually carried out in two steps. We first map them into a particular region or patch via

$$x^{i}=x^{i}(X^{j}).$$

(36)

For example, we can use standard spherical coordinates $X^{i}=(r,\theta,\varphi)$ with a range $r\in[r_{\mathrm{min}},r_{\mathrm{max}}]$, $\theta\in[\theta_{\mathrm{min}},\theta_{\mathrm{max}}]$, $\varphi\in[\varphi_{\mathrm{min}},\varphi_{\mathrm{max}}]$, so that we cover a certain section of a shell. Next we use

$$X^{i}=\frac{1}{2}\left[(X^{i}_{\mathrm{max}}-X^{i}_{\mathrm{min}})\bar{X}^{i}+X^{i}_{\mathrm{max}}+X^{i}_{\mathrm{min}}\right]$$

(37)

to map each $X^{i}$ into an $\bar{X}^{i}$ that has the standard range $\bar{X}^{i}\in[-1,1]$. These $\bar{X}^{i}$ are what have been denoted by $x^{\bar{i}}$ in Eq. (10). Each patch is thus described by the particular transformation (36) and range we use for the $X^{i}$ coordinates. In some cases, we only need Cartesian coordinates so that we use the identity transformation in Eq. (36), but we have also implemented the transformation to the cubed sphere coordinates $X^{i}=(\lambda,A,B)$ described in [56]. We then arrange our various patches such that they touch and cover the region of interest.

An example is shown on the left side of Fig. 1. Here we have one central cube that is covered by Cartesian coordinates. This cube is surrounded by six cubed sphere patches. Five of these seven patches intersect the $xy$-plane and are shown on the left of Fig. 1. Note that each of the six cubed sphere patches shares one face with the central cube. The face on the opposite side is curved and arises by deforming one side of a larger cube into a spherical surface via a coordinate transformation of the form in Eq. (36). It thus comprises one sixth of the spherical outer boundary. The remaining faces of each cubed sphere patch touch four other cubed sphere patches, along flat surfaces. All patches are touching each other without any overlap, so that all interior patch faces have their entire face in common with one other patch face. In the next section we explain how information is exchanged via numerical fluxes between adjacent patches.

The patches described before can be directly covered with the Legendre Gauß-Lobatto grid points introduced above. Yet, to gain more flexibility, we can further refine each patch in Nmesh. This is achieved by identifying each patch with a so-called root node that can be further refined. When we refine this root node, we cut the original ranges of all three $X^{i}$-coordinates in half so that we end up with eight touching child nodes that now cover the original root node or patch. Each of these new nodes can then be further refined by again dividing it into eight child nodes. In this way we can refine each patch as often as we want. This can be done in an irregular way, where we further refine only the nodes in certain regions of interest. We end up with a node tree called an octree, where each node has either zero or eight child nodes. The nodes without children are called leaf nodes. Together, these leaf nodes cover the entire patch and are thus the nodes in which we perform any calculations. For this reason, the leaf nodes are often called computational elements or just elements. However, in this paper, we will simply call them leaf nodes or just nodes. We note that, in the context of finite volume methods, the word “node” is also sometimes used in the literature to denote a grid point. Yet, in this paper the word “node” will always refer to a node in our octree.

The $X^{i}$-range of each leaf node is covered by grid points that correspond to the Legendre Gauß-Lobatto points discussed above. This means that we have grid points on each node face. This simplifies any calculations that depend on the values of fields on both sides of a node boundary. The number of grid points in each node can be freely chosen. When we increase it, we obtain higher order accuracy if the fields are smooth within the node. This is called p-refinement because increasing the number of grid points corresponds to an increase in the number of basis polynomials we use to represent a field within a node. Of course, p-refinement is most useful for smooth fields. For non-smooth fields it is often better to refine a node by splitting it into eight child nodes, which is known as h-refinement as it refines the resolution even if each child node has still the same number of grid points as the parent node. In Nmesh both p- and h-refinement can be performed whenever desired. Together we call this Adaptive Mesh Refinement (AMR).

On the right of Fig. 1, we show an example where we h-refine the top node from the left side of Fig. 1. The resulting child nodes now cover the top patch. As we can see, the grid points of the h-refined nodes along e.g. the left patch boundary no longer all coincide with the grid points of the unrefined node covering the left patch. This is a general phenomenon, whenever two neighboring nodes differ in their h- or p-refinement, many of the surface grid points of one node do not coincide with the surface grid points of the touching adjacent node. Furthermore, the surface of one node may be touching several adjacent nodes, as is the case for the node covering the left (orange) patch. This complicates the calculation of numerical fluxes such as Eq. (33) or Eq. (35) in one of our nodes, because we need both the fields $u_{\mathrm{in}}$ and $u_{\mathrm{adj}}$ at every surface grid point of the current node. We already have $u_{\mathrm{in}}$ at every point of the current node. However, $u_{\mathrm{adj}}$ only exists at the grid points of the adjacent nodes, which in general do not coincide with the surface grid points of the current node. To obtain $u_{\mathrm{adj}}$ at one of the surface grid points of the current node, we interpolate the $u_{\mathrm{adj}}$ data from the adjacent node onto this point. For this we currently use Lagrange interpolating polynomials constructed from the 2-dimensional surface data of the adjacent node. To easily find adjacent neighbors each node has an associated data structure that keeps track of all adjacent neighbor nodes. When p-refinement is applied to a node this data structure can remain unchanged since the size of the nodes and therefore the number of neighbors does not change. However, the data structure has to be updated whenever a node is h-refined. In this case the structure gets updated on this one node and on all of its neighbors. In this way Nmesh is able to accommodate arbitrary levels of h-refinement. For example, it is possible to h-refine a node into eight child nodes, and to then repeat this as often as desired with any of the child nodes, without at the same time refining any of the original neighbor nodes. In each case the final result is a number of touching leaf nodes that cover each patch. Since the patches themselves are also touching, the collection of all leaf nodes from all patches forms the mesh on which we perform our calculations.

The word Adaptive in AMR usually also implies an algorithm that automatically chooses p- or h-refinement. We currently have implemented only one such algorithm. Within each node, it expands a chosen quantity, such as the matter density $\rho_{0}$, in terms of Legendre polynomials. The coefficients in front of each of the Legendre polynomials can then be used to judge the smoothness of $\rho_{0}$. If $\rho_{0}$ is perfectly smooth, we expect the coefficients to fall off exponentially with increasing polynomial order. In our algorithm we then fit the logarithm of the coefficient magnitudes to a linear function. If the slope values $b_{i}$ in all three directions ($i=1,2,3$) of this linear function are not negative enough, we consider $\rho_{0}$ to be not smooth. We then h-refine the node. This algorithm is in principle geared toward dealing with the non-smooth behavior of matter across a neutron star surface. We have, however, not had any real success with this algorithm yet. We can turn it on and evolve neutron stars with it, but we have not been able to tune the parameters, that decide when the $b_{i}$ are not negative enough, to values that work well after some matter leaves the star surface. We mention this particular algorithm here only to show that Nmesh has AMR capabilities in principle. We note, however, that these capabilities were not used in the simulations discussed below, where we use uniform h-refinement.

Figure. 2 shows the mesh for a single neutron star in a plane through its center. The different levels of h-refinement shown here are obtained from the coefficient drop off based algorithm mentioned above. It also refines neighbor nodes if their level of refinement is more than one below the just refined node. The purpose of the latter is to avoid abrupt changes in resolution, but is technically not required by Nmesh.

Note that Eqs. (2.1) or (2.2) still contain time derivatives, since up to this point we have only discretized spatial derivatives. This means that these equations represent a set of coupled ordinary differential equations for the fields $u_{q_{1}q_{2}q_{3}}$ at the grid points. In this paper, we use standard Runge-Kutta time integrators to find the solution of these ODEs from the $u_{q_{1}q_{2}q_{3}}$ at the initial time. Such Runge-Kutta methods are only stable when the Courant-Friedrichs-Lewy (CFL) condition is satisfied, i.e., if the time step $\Delta t$ is small enough. In this work we use

$$\Delta t=\Delta x_{\mathrm{min}}/v,$$

(38)

where $\Delta x_{\mathrm{min}}$ is the distance between the two closest grid points, and $v$ is a number that needs to be greater than the largest characteristic speed.

Even when the time step satisfies the CFL condition, instabilities can still occur in some cases, e.g., on non-Cartesian patches. To combat such instabilities, we filter out high frequency modes in the evolved fields. This is achieved by first computing the coefficients $c_{l_{1}l_{2}l_{3}}$ in the expansion

$$u(\bar{x})J=\sum_{l_{1}=0}^{N}\sum_{l_{2}=0}^{N}\sum_{l_{3}=0}^{N}c_{l_{1}l_{2}l_{3}}P_{l_{1}}(x^{\bar{1}})P_{l_{2}}(x^{\bar{2}})P_{l_{3}}(x^{\bar{3}})$$

(39)

of each field $u$, where the $P_{l}(\bar{x})$ are Legendre polynomials. The coefficients are then replaced by

$$c_{l_{1}l_{2}l_{3}}\to c_{l_{1}l_{2}l_{3}}e^{-\alpha_{\mathrm{f}}(l_{1}/N)^{s}}e^{-\alpha_{\mathrm{f}}(l_{2}/N)^{s}}e^{-\alpha_{\mathrm{f}}(l_{3}/N)^{s}},$$

(40)

and $u$ is recomputed using these new coefficients. Note that we typically use $\alpha_{\mathrm{f}}=36$ and $s=32$, so that the coefficients with the highest $l=N$ are practically set to zero, while all others are mostly unchanged.

One of the simplest systems one can evolve is a scalar wave. Here we consider a single scalar field obeying the wave equation

$$\partial_{t}^{2}\phi=\delta^{ij}\partial_{i}\partial_{j}\phi.$$

(41)

The DG method described earlier cannot be applied directly to systems with second order derivatives. We therefore introduce the extra variables

$$\Pi:=\partial_{t}\phi$$

(42)

and

$$\chi_{i}:=\partial_{i}\phi.$$

(43)

This results in the following system of first order equations

	$\displaystyle\partial_{t}\Pi+\partial_{j}f_{\Pi}^{j}$	$\displaystyle=$	$\displaystyle 0,$
	$\displaystyle\partial_{t}\chi_{i}+\partial_{j}f_{\chi_{i}}^{j}$	$\displaystyle=$	$\displaystyle 0,$
	$\displaystyle\partial_{t}\phi$	$\displaystyle=$	$\displaystyle\Pi.$		(44)

Here we have defined the flux vectors

	$\displaystyle f_{\Pi}^{j}$	$\displaystyle:=$	$\displaystyle-\chi_{j},$
	$\displaystyle f_{\chi_{i}}^{j}$	$\displaystyle:=$	$\displaystyle-\Pi\delta_{i}^{j},$
	$\displaystyle f_{\phi}^{j}$	$\displaystyle:=$	$\displaystyle 0.$		(45)

As we can see, the system in Eq. (4.1), consists of two coupled partial differential equations and one ordinary differential equation.

To evolve this system with our DG method, we also need to provide initial values and boundary conditions. Since

$$\phi=\sin(k_{i}x^{i}-\omega t)$$

(46)

is a solution for any $k_{i}$ and $\omega=\sqrt{\delta^{ij}k_{i}k_{j}}$, we initialize the system according to this equation at $t=0$. We also use Eq. (46) at the outer boundary so that this sine wave is continuously entering through the outer boundary.

The DG method requires numerical fluxes. We have successfully used both the LLF flux of Eq. (33) as well as the upwind flux of Eq. (35). To impose Eq. (46) at the outer boundary, we use the same numerical flux as in the interior, but we set $u_{\mathrm{adj}}$ to the value coming from Eq. (46).

The wave vector $k_{i}$ in Eq. (46) is arbitrary. For the test results presented here, we choose $k_{i}=(0.7,-2,4.3)$, so that it represents the general case where $k_{i}$ is not aligned with any coordinate direction. As we can see in Fig. 4, we get a sinusoidal wave that propagates through our numerical domain, with this $k_{i}$ vector. In this case we have used seven patches and the figure shows the scalar field $\phi$ in the $xy$-plane after evolving up to time $t=4.5$. For the test case depicted in Fig. 4, we have used an equal number of grid points ($19\times 19\times 19$) in all directions, without any h-refinement applied to the root nodes. However, we have also performed tests with an unequal number of grid points and with the root nodes h-refined. We have obtained stable evolution for the system for all these cases with both the LLF and upwind numerical fluxes. The choice of numerical flux did not have any significant effect in any of the scalar wave evolution tests, as both cases yield results that have errors of the same order.

To demonstrate the convergence of our new Nmesh program in this scalar wave evolution test, we have applied p-refinement and h-refinement separately. In Fig. 5 we plot the L2-norm of the error in $\phi$ for both. For the case of p-refinement (top) we increase the number of points $n$ in all directions, setting them to $n=10,12,14,16,18,20$, with no h-refinement applied to the root nodes. We observe an exponential drop in the L2-norm error, as expected in this case. For the case of h-refinement (bottom), we apply $l=0,1,2,3,4,5$ levels of refinement to the root nodes, and have $n=10$ points in each direction in each node. Again, we observe convergent behavior for the L2-norm of the error, when increasing the number of h-refinement levels. Figure 5 only shows the results for the LLF numerical flux, since the results for the upwind flux are very similar. For all the runs shown in Figs. 4 and 5, the time step was set to $\Delta t=\Delta x_{\mathrm{min}}/5$ in accordance with Eq. 38.

Even though all the convergence results stated above have been obtained for a mesh using six cubed sphere patches that surround one central cube, we also obtain convergence for a mesh covered by a single cubic Cartesian patch. The key difference between these is that we need the filters of Eq. 40 to stabilize the evolution on meshes that contain both Cartesian and cubed sphere patches, while this is not necessary on a purely Cartesian patch, even if it is h-refined. For the runs of Figs. 4 and 5 we set the filter parameters of Eq. 40 to the values $\alpha_{\mathrm{f}}=36$ and $s=32$.

For the gravitational part, we have implemented the first-order reduction of Generalized Harmonic Gauge (GHG) formulation [58, 31]

		$\displaystyle\partial_{t}g_{ab}$	$\displaystyle-(1+\gamma_{1})\beta^{k}\partial_{k}g_{ab}=-\alpha\Pi_{ab}-\gamma_{1}\beta^{k}\Phi_{kab},$		(49)
		$\displaystyle\partial_{t}\Pi_{ab}$	$\displaystyle-\beta^{k}\partial_{k}\Pi_{ab}+\alpha\gamma^{ki}\partial_{k}\Phi_{iab}-\gamma_{1}\gamma_{2}\beta^{k}\partial_{k}g_{ab}=$
			$\displaystyle 2\alpha g^{cd}\left(\gamma^{ij}\Phi_{ica}\Phi_{jdb}-\Pi_{ca}\Pi_{db}-g^{ef}\Gamma_{ace}\Gamma_{bdf}\right)$
			$\displaystyle-2\alpha\nabla_{(a}H_{b)}-\frac{1}{2}\alpha n^{c}n^{d}\Pi_{cd}\Pi_{ab}-\alpha n^{c}\Pi_{ci}\gamma^{ij}\Phi_{jab}$
			$\displaystyle+\alpha\gamma_{0}\left[2\delta^{c}{}_{(a}n_{b)}-g_{ab}n^{c}\right]\left(H_{c}+\Gamma_{c}\right)-\gamma_{1}\gamma_{2}\beta^{k}\Phi_{kab}$
			$\displaystyle-16\pi\alpha\left(T_{ab}-\frac{1}{2}g_{ab}g^{cd}T_{cd}\right),$
		$\displaystyle\partial_{t}\Phi_{iab}$	$\displaystyle-\beta^{k}\partial_{k}\Phi_{iab}+\alpha\partial_{i}\Pi_{ab}-\alpha\gamma_{2}\partial_{i}g_{ab}=$
			$\displaystyle\frac{1}{2}\alpha n^{c}n^{d}\Phi_{icd}\Pi_{ab}+\alpha\gamma^{jk}n^{c}\Phi_{ijc}\Phi_{kab}-\alpha\gamma_{2}\Phi_{iab}.$

Here $g_{ab}$ is the spacetime metric, $n_{a}$ is the unit normal to the hypersurface of constant coordinate time $t$, and $\Gamma_{a}=g^{bc}\Gamma_{abc}$ is the contracted Christoffel symbol. The equations are written in terms of the extra variables $\Pi_{ab}:=-n^{c}{\partial}_{c}g_{ab}$ and $\Phi_{iab}:={\partial}_{i}g_{ab}$, that have been introduced to make the original second-order GHG system first-order in both time and space. Gauge conditions in the GHG system are specified by prescribing the gauge source function $H_{a}$. The lapse $\alpha$, shift $\beta^{i}$ and spatial metric $\gamma_{ij}$ come from the $3+1$ decomposition in Eq. (2). The GHG evolution equations also contain extra terms that are multiplied with the parameters $\gamma_{0}$, $\gamma_{1}$, and $\gamma_{2}$. In this paper we set $\gamma_{1}=-1$, and choose $\gamma_{2}=\gamma_{0}=1$ for the constraint damping parameters.

To test the gravitational part of Nmesh, we evolve a black hole spacetime. As initial data, we choose the metric of a single Schwarzschild black hole in Kerr-Schild coordinates [59],

$$g_{ab}=\eta_{ab}+\frac{2M}{r}l_{a}l_{b},$$

(50)

where $\eta_{ab}$ is the Minkowski metric, and $M$ is the mass of the black hole. In the Cartesian coordinates, $r=(x^{2}+y^{2}+z^{2})^{1/2}$, and $l_{a}=(1,x/r,y/r,z/r)$. The gauge source function is initialized based on the above metric (50), and is left constant during the simulation [60],

$$H_{a}(t=0)=-\Gamma_{a}(t=0),\quad\partial_{t}H_{a}=0.$$

(51)

With this initial condition, the analytic solution of the evolution equations is simply the static Schwarzschild metric, so that all evolved fields should be constant. Of course evolution will lead to some amount of numerical errors. Thus we test here if Nmesh can stably evolve this setup and whether the numerical evolution will settle down to a stable state.

As computational domain, we choose a spherical shell that extends from $r=1.8M$ to $11.8M$, and is covered by six cubed sphere patches. The inner boundary is thus inside the black hole horizon of $r=2M$. The speeds of all characteristic modes at the inner boundary are such that every mode is moving towards the black hole center and thus leaving the computational domain. We therefore do not impose any boundary conditions at the inner boundary. The situation at the outer boundary is different as we have both incoming and outgoing modes. We impose no condition on the outgoing modes, but we keep the incoming modes constant at their initial values (consistent with the static analytic solution). This is done using Eq. (35), where $u_{\mathrm{adj}}$ is set to the analytic Schwarzschild solution, $u_{\mathrm{in}}$ are the evolved fields at the boundary, and $S$, $\Lambda^{+}$, and $\Lambda^{-}$ come from the characteristic modes and their eigenvalues [58], calculated from normals that are normalized with respect to the flat metric. To improve the accuracy we use either 2 or 3 levels of h-refinement in each patch. We choose the time step according to Eq. (38), with $v=4$. We find that, with this setup, no filters are necessary to stabilize our runs. As in the scalar field test cases discussed above, filters only become necessary when both Cartesian and cubed sphere patches are present. In the latter case, the filter of Eq. (40) is again sufficient for obtaining stable runs. We have evolved this setup using both the LLF and upwind fluxes of Eqs. 33 and 35 at inter domain boundaries. As described below, both fluxes work about equally well.

In order to demonstrate stability of our runs at high resolution, we have evolved the black hole with three levels of h-refinement for three different numbers of grid points. We find that the system of equations reaches a state where the time derivatives $\partial_{t}g_{ab}$, $\partial_{t}\Pi_{ab}$, and $\partial_{t}\Phi_{iab}$ all approach zero (up to machine precision), as expected for a static black hole. As an example, we show $\partial_{t}\Pi_{xx}$ in Fig. 6 when evolved with $4\times 4\times 4$, $5\times 5\times 5$, and $6\times 6\times 6$ grid points in each node. As we can see, this time derivative falls exponentially until it settles down to below $10^{-12}$. Beyond this point, the terms that determine the time derivatives in the GHG evolution equations (49), (49), (49) add up to almost zero, and deviate from zero only because of roundoff errors due to the use of floating point numbers. As expected, at higher resolution this steady state is reached earlier, because our numerical method then has smaller discretization errors.

The Hamiltonian and momentum constraints of general relativity read as

	$\displaystyle H$	$\displaystyle=$	$\displaystyle R-K_{ij}K^{ij}+K^{2}-16\pi\rho_{\mathrm{ADM}},$		(52)
	$\displaystyle M^{i}$	$\displaystyle=$	$\displaystyle D_{j}(K^{ij}-\gamma^{ij}K)-8\pi j^{i},$		(53)

where $R$ and $D_{j}$ are the Ricci scalar and derivative operator associated with the 3-metric $\gamma_{ij}$, $K_{ij}=-\frac{1}{2}\pounds_{n}\gamma_{ij}$, $\rho_{\mathrm{ADM}}=T_{ab}n^{a}n^{b}$, and $j^{i}=-T_{ab}n^{a}\gamma^{bi}$ (see e.g. [61]). General relativity dictates $H=M^{i}=0$ for all time. In Fig. 7, we show the infinity norm of $H$ over the grid when we evolve the black hole with 2 levels of h-refinement for various numbers of grid points per node. As we can see, $H$ stabilizes after a short time and then stays practically constant, which again indicates stability. As expected $H$ converges to zero exponentially as we increase the number of grid points. We also see that the results for the LLF flux (on the left) and the upwind flux (on the right) are almost the same. In our simulations the momentum constraint $M^{i}$ behaves just like $H$ and also converges to zero exponentially as we increase the number of grid points.

Since the initial data is given by the Kerr-Schild metric, the analytic solution is a static black hole. The numerical solution, however, evolves for a while until it settles down into a stable configuration (see Fig. 6). This happens because the analytic solution does not exactly satisfy the discretized GHG equations. As already mentioned in Sec. 2 we use domain normals that are normalized with respect to the flat metric for our black hole evolutions.

In Fig. 8 we plot the metric component $g_{xx}$ at different times (for the $n=4$ case of Fig. 6). We have zoomed in onto a region close to a domain boundary that is in the strong field region close to the horizon at $x=2M$. We find that the metric rapidly evolves away from the continuous Kerr-Schild initial data. During this rapid evolution discontinuities develop across domain boundaries due to numerical errors. These discontinuities persist throughout the evolution, but do not negatively affect its stability or convergence, even though dynamic evolution takes place. This indicates that such discontinuities in the physical metric are not a problem for a DG method that uses numerical fluxes and eigenvalues, which are computed from normals that are normalized with respect to the flat metric.

To treat neutron star matter we use the Valencia formulation [62]. Matter is thus described as a perfect fluid, where the stress-energy tensor is given by

$$T^{\mu\nu}=\left(\rho+P\right)u^{\mu}u^{\nu}+Pg^{\mu\nu}.$$

(54)

Here $\rho$ is the energy density, $P$ is the pressure, and $u^{\mu}$ is the four-velocity. The total energy density is written as

$$\rho=\rho_{0}(1+\epsilon),$$

(55)

where $\rho_{0}$ is the rest-mass energy density, and $\epsilon$ is the specific internal energy. We express the four-velocity $u^{\mu}$ in terms of the three-velocity given by

$$v^{\mu}=u^{\mu}/W-n^{\mu}$$

(56)

and also introduce the Lorentz factor

$$W=-n_{\mu}u^{\mu}=\alpha u^{0}.$$

(57)

Together $(\rho_{0},\epsilon,Wv^{i},P)$ are known as the primitive variables.

The matter equations follow from the conservation law for the energy-momentum tensor and the conservation law for the baryon number. In order to obtain flux conservative evolution equations of the form (5), one introduces the conserved variables

	$\displaystyle D$	$\displaystyle=$	$\displaystyle\rho_{0}W,$		(58)
	$\displaystyle\tau$	$\displaystyle=$	$\displaystyle\rho_{0}hW^{2}-P-\rho_{0}W,$		(59)
	$\displaystyle S_{i}$	$\displaystyle=$	$\displaystyle\rho_{0}hW^{2}v_{i}.$		(60)

Here $D$ is rest-mass density, $\tau$ the internal energy density, $S_{i}$ the momentum density as seen by Eulerian observers. The last two equations also contain the specific enthalpy given by

$$h=1+\epsilon+P/\rho_{0}.$$

(61)

The conserved variables are then

$$u=\sqrt{\gamma}\left(\begin{array}[]{ccc}D\\ \tau\\ S_{l}\\ \end{array}\right).$$

(62)

They satisfy Eq. (5), with the flux vectors and sources given by

$$f^{i}=\sqrt{\gamma}\left(\begin{array}[]{ccc}(\alpha v^{i}-\beta^{i})D\\ (\alpha v^{i}-\beta^{i})\tau+\alpha Pv^{i}\\ (\alpha v^{i}-\beta^{i})S_{l}+\alpha P\delta^{i}_{l}\\ \end{array}\right)$$

(63)

and

$$\frac{s}{\sqrt{\gamma}}=\left(\begin{array}[]{ccc}0\\ T^{00}(\beta^{i}\beta^{j}K_{ij}-\beta^{i}\partial_{i}\alpha)+T^{0i}(2\beta^{j}K_{ij}-\partial_{i}\alpha)+T^{ij}K_{ij}\\ T^{00}(\frac{\beta^{i}\beta^{j}}{2}\partial_{l}\gamma_{ij}-\alpha\partial_{l}\alpha)+T^{0i}\beta^{j}\partial_{l}\gamma_{ij}+T^{0}_{i}\partial_{l}\beta^{i}+\frac{T^{ij}}{2}\partial_{l}\gamma_{ij}\\ \end{array}\right).$$

(64)

The components of the stress-energy tensor appearing here can be expressed in terms of the primitive variables as

	$\displaystyle T^{00}$	$\displaystyle=$	$\displaystyle(W^{2}h\rho_{0}-P)/\alpha^{2},$		(65)
	$\displaystyle T^{0i}$	$\displaystyle=$	$\displaystyle Wh\rho_{0}u^{i}/\alpha+P\beta^{i}/\alpha^{2},$		(66)
	$\displaystyle T^{ij}$	$\displaystyle=$	$\displaystyle h\rho_{0}u^{i}u^{j}+P(\gamma^{ij}-\beta^{i}\beta^{j}/\alpha^{2}),$		(67)
	$\displaystyle T^{0}_{i}$	$\displaystyle=$	$\displaystyle h\rho_{0}W^{2}v_{i}/\alpha,$		(68)

where

$$u^{i}=Wv^{i}-W\frac{\beta^{i}}{\alpha}.$$

(69)