Fast hyperbolic relaxation elliptic solver for numerical relativity: Conformally flat, binary puncture initial data

We introduce nrpyelliptic, an elliptic solver for numerical relativity (NR) built within the nrpy+ framework. As its first application, nrpyelliptic sets up conformally flat, binary black hole (BBH) puncture initial data (ID) on a single numerical domain, similar to the widely used twopunctures code. Unlike twopunctures, nrpyelliptic employs a hyperbolic relaxation scheme, whereby arbitrary elliptic partial differential equations (PDEs) are trivially transformed into a hyperbolic system of PDEs. As consumers of NR ID generally already possess expertise in solving hyperbolic PDEs, they will generally find nrpyelliptic easier to tweak and extend than other NR elliptic solvers. When evolved forward in (pseudo)time, the hyperbolic system exponentially reaches a steady state that solves the elliptic PDEs. Notably nrpyelliptic accelerates the relaxation waves, which makes it many orders of magnitude faster than the usual constant wave speed approach. While it is still $\sim 12\mathrm{x}$ slower than twopunctures at setting up full-3D BBH ID, nrpyelliptic requires only $\approx 0.3\%$ of the runtime for a full BBH simulation in the einstein toolkit. Future work will focus on improving performance and generating other types of ID, such as binary neutron stars.

Figure 1

Curves of constant ${\mathtt{x}}_{\mathtt{1}}$ (red) and constant ${\mathtt{x}}_2$ (blue) using a cell-centered grid structure with ${\mathtt{x}}_{\mathtt{3}}=0$.

Figure 2

Calibration of twopunctures and nrpyelliptic solutions, axisymmetric ID study. Top: Relative errors between a super-high resolution (${196}^2\times 4$) and lower-resolution results from the twopunctures code. Bottom: Calibration of nrpyelliptic solutions at various resolutions (${\mathtt{N}}_{\mathtt{1}}\times {\mathtt{N}}_{\mathtt{2}}\times 6$), comparing against the trusted solution (i.e., the twopunctures result at ${96}^2\times 4$ resolution). The horizontal axis is logarithmic in the range $|z|>1$ and linear otherwise. Further, the black dots on the horizontal axes denote the puncture BH positions (also the locations of the coordinate foci). twopunctures data were rotated so that the punctures and the foci lie on the $z$ axis.

Figure 3

Axisymmetric ID study. Top: Snapshots of the relative error between the nrpyelliptic solution at ${128}^2\times 4$ resolution and the trusted solution at different times. twopunctures data were rotated so that the punctures and the foci lie on the $z$ axis. Bottom: Black dashed curve: $L^2$-norm of the same relative error as the top panel, computed within a sphere with radius $r=100M$ and as a function of time. Blue (dotted) and red (solid) curves: $L^2$-norms of residuals [second term in Eq. (18)] computed within the same sphere at two different finite-difference orders as a function of time.

Figure 4

Same as Fig. 2, but for GW150914-like, full-3D BBH initial data.

Figure 5

Same as Fig. 3, but for GW150914-like, full-3D BBH initial data, with nrpyelliptic resolution of ${128}^2\times 16$.

Figure 6

Initial data quality of the GW150914-like BBH 3D scenario evolved with the einstein toolkit, showing the relative error of $u$ against the trusted twopunctures solution. Here twopunctures and nrpyelliptic ID used resolutions of ${38}^2\times 16$ and ${128}^2\times 16$, respectively. Note that twopunctures data were rotated so that the punctures are positioned on the $z$ axis.

Figure 7

Hamiltonian constraint violation ${\log }_{10}|\mathcal{H}|$, as defined in Eq. (4), in orbital plane after interpolation to AMR grids at $t=0$ (left) and after one orbit (right). Note that nrpyelliptic data were rotated so that the punctures are positioned on the $x$ axis.

Figure 8

Simulation results with ID by twopunctures (solid, blue) and nrpyelliptic (dashed, orange-red). Left: Dominant mode ($\ell =2$, $m=2$) of the Weyl scalar ${\psi }_4$ at extraction radius $r_{\mathrm{ext}}=500M$. Right: Trajectory of the less-massive BH.

Figure 9

Schematic of seven-grid BiSpheres numerical meshes used for blackholes@home $\sim 6:1$ mass-ratio BBH simulations during the long inspiral phase. BHs are represented as black dots.

Figure 10

Top: Propagation of $u$ along gridpoints closest to the $x$-axis during the first relaxation-wave crossing time for the 3D ID. Bottom: Same as top panel, except for gridpoints closest to the $z$-axis.