Framework for large-scale relativistic simulations in the characteristic approach

We present a new computational framework (LEO) that enables us to carry out the very first large-scale, high-resolution computations in the context of the characteristic approach in numerical relativity. At the analytic level, our approach is based on a new implementation of the eth formalism, using a nonstandard representation of the spin-raising and spin-lowering angular operators in terms of nonconformal coordinates on the sphere; we couple this formalism to a partially first-order reduction (in the angular variables) of the Einstein equations. The numerical implementation of our approach supplies the basic building blocks for a highly parallel, easily extensible numerical code. We demonstrate the adaptability and excellent scaling of our numerical code by solving, within our numerical framework, for a scalar field minimally coupled to gravity (the Einstein-Klein-Gordon problem) in 3 dimensions. The nonlinear code is globally second-order convergent, and has been extensively tested using as a reference a calibrated code with the same initial and boundary data and radial marching algorithm. In this context, we show how accurately we can follow quasinormal mode ringing. In the linear regime, we show energy conservation for a number of initial data sets with varying angular structure. A striking result that arises in this context is the saturation of the flow of energy through the Schwarzschild radius. As a final calibration check, we perform a large simulation with resolution never achieved before.

Figure 1

The cubed sphere: covering of the sphere with six nonoverlapping gnomic patches.

Figure 2

Convergence rate of the $ð$ operator, built upon angular derivatives of order 2, 4, 6, and 8 (indicated in the graph as circles, squares, diamonds, and triangles, respectively), acting on ${}_{2}Y_{43}$, and with grid sizes ranging from $N_{\zeta }=16$ to $N_{\zeta }=128$.

Figure 3

Convergence rate of the various interpolation schemes used. Shown in the graph are the third-order (circles), fifth-order (squares), seventh-order (diamonds), and ninth-order (triangles) interpolators. For the highest order interpolator used (ninth order), the error goes down to double-precision round-off level ($\sim {10}^{-16}$) when more than $N_{\zeta }=100$ angular points per patch are used.

Figure 4

Convergence rate of the integral of the area element over the sphere, for grid sizes ranging from $N_{\zeta }=8$ to $N_{\zeta }=512$. The markers indicate the error of the area element at the corresponding resolution; the line is the least-squares fit, yielding a convergence rate of 2.0.

Figure 5

Convergence of the orthonormality condition, illustrated here by computing the convergence rate of ${\int }_S{}_{0}Y_{lm0}{\overline{Y}}_{lm}\equiv 1$, for the case $l=6$, $m=0\dots 6$, on grid sizes ranging from $N_{\zeta }=32$ to $N_{\zeta }=64$.

Figure 6

Error in the coefficients $c_{lm}$ when the function being projected is the spherical harmonic $Y_{63}$, on a grid of $N_{\zeta }=64$ points. Coefficients whose errors are at round-off level are not shown.

Figure 7

The function $\chi (u)$ at $\mathcal{I}$ as a function of Bondi time, showing the quasinormal mode regime oscillations for $\ell =1$, $m=0$. The solid line is the output from LEO for $N_{\zeta }=11$ and $N_x=1001$, when the initial data and boundary conditions are given as for the convergence test; the dashed line is the quasinormal mode extracted from the data.

Figure 8

Log of the absolute value of the function $\chi (u)$ at $\mathcal{I}$ as a function of Bondi time. Parameters and conditions are the same as in Fig. 7. The solid line is the output from LEO; the dashed line is the quasinormal mode extracted from the data.

Figure 9

The function $\chi$ at $\mathcal{I}$ as a function of Bondi time, showing the quasinormal mode regime oscillations for $\ell =2$, $m=0$. The solid line is the output from LEO for $N_{\zeta }=11$ and $N_x=1001$, when the initial data and boundary conditions are given as for the convergence test, except that $r_{\mathrm{in}}=2.13M$; the dashed line is the quasinormal mode extracted from the data.

Figure 10

Log of the absolute value of the function $\chi (u)$ at $\mathcal{I}$ as a function of Bondi time. Parameters and conditions are the same as in Fig. 9. The solid line is the output from LEO; the dashed line is the quasinormal mode extracted from the data.

Figure 11

Energy conservation as a function of Bondi time for $\ell =0$ (solid line), $\ell =1$ (dotted line), and $\ell =2$ (long dashed line). This calculation was done using the same grid parameters as for Fig. 9 except for $r_{\mathrm{in}}=2.3$. For each specific $\ell$, the descending curve corresponds to energy given by Eq. (54). The ascending curve corresponds to the algebraic sum of $E_{\mathrm{in}}=-\int P_{\mathrm{in}}du$ and $E_{\mathrm{out}}=\int P_{\mathrm{out}}du$. Thus, in accordance with Eq. (56), the horizontal curve represents the global conservation of energy.

Figure 12

Percentage variation in $\Sigma (u)$ with respect to $\Sigma (0)$ as a function of Bondi time for $\ell =0$ (circles), $\ell =1$ (squares), and $\ell =2$ (triangles). The graph shows that energy is conserved to within less than 0.2% of the energy content of the initial surface.

Figure 13

Energy content $E(u)$ as a function of Bondi time for $\ell =0$ (circles), $\ell =1$ (squares), $\ell =2$ (diamonds), $\ell =3$ (triangles), and $\ell =4$ (stars). This calculation was done using the grid parameters $N_{\zeta }=11$ and $N_x=1001$. The initial data and boundary conditions are the same as in the convergence test.

Figure 14

Energy flow to infinity $E_{\mathrm{out}}=\int P_{\mathrm{out}}(u)du$ as a function of Bondi time for $\ell =0$ (circles), $\ell =1$ (squares), $\ell =2$ (diamonds), $\ell =3$ (triangles), and $\ell =4$ (stars). This calculation was done using the same conditions as in Fig. 13.

Figure 15

Energy flow into the black hole, $E_{\mathrm{in}}=-\int P_{\mathrm{in}}(u)du$, as a function of Bondi time for $\ell =0$ (circles), $\ell =1$ (squares), $\ell =2$ (diamonds), $\ell =3$ (triangles), and $\ell =4$ (stars). Both curves for $\ell =3$ and $\ell =4$ eventually saturate, without crossing, for for $u>30$. This calculation was done using the same conditions as in Fig. 13.

Figure 16

Surface plots of $\chi$ at $\mathcal{I}$ for $u=2.5M$ (top graph) and $u=30M$ (bottom graph). The parameters of the initial data are $\lambda ={10}^{-4}$, $r_0=3M$, $\sigma =0.5M$, $\ell =8$, $m=6$. The grid size is $N_{\zeta }=93$, $N_x=1501$.

Figure 17

Surface plots of $J\overline{J}$ at $\mathcal{I}$ for $u=2.5M$ (top graph) and $u=30M$ (bottom graph). The parameters of the initial data and grid size are the same as in Fig. 16.

Figure 18

Surface plots of $\beta$ at $\mathcal{I}$ for $u=2.5M$ (top graph) and $u=30M$ (bottom graph). The parameters of the initial data and grid size are the same as in Fig. 16.

Figure 19

Surface plots of $U\overline{U}$ at $\mathcal{I}$ for $u=2.5M$ (top graph) and $u=30M$ (bottom graph). The parameters of the initial data and grid size are the same as in Fig. 16.

Figure 20

Surface plots of $W$ at $\mathcal{I}$ for $u=2.5M$ (top graph) and $u=30M$ (bottom graph). The parameters of the initial data and grid size are the same as in Fig. 16.

Figure 21

Energy conservation for the simulation which generated the results shown in Figs. 16, 17, 18, 19, 20. The solid line corresponds to the energy content $E(u)$ at successive times, the dashed line to the sum of the energy radiated through the inner [$E_{\mathrm{in}}(u)$] and outer [$E_{\mathrm{out}}(u)$] boundaries, and the dot-dashed line to the sum, $\Sigma (u)=E(u)+E_{\mathrm{in}}(u)+E_{\mathrm{out}}(u)$, respectively. The insert graph shows the percentage variation in $\Sigma (u)$, relative to its final value at $u=30$.