Implementation of a stable solver for the full (3+1D) vacuum Einstein equations

[JonathanGorard]

This branch contains an implementation of a numerically stable solver for the vacuum ADM evolution equations for the components of the spatial metric tensor $\gamma_{i j}$

$$ \left( \partial_t - \mathcal{L}_{\boldsymbol\beta} \right) \gamma_{i j} = - 2 \alpha K_{i j} $$

and its corresponding canonical momenta, namely the components of the extrinsic curvature tensor $K_{i j}$:

$$ \left( \partial_t - \mathcal{L}_{\boldsymbol\beta} \right) K_{i j} = - {}^{\left( 3 \right)} \nabla_i {}^{\left( 3 \right)} \nabla_j \alpha + \alpha \left[ {}^{\left( 3 \right)} R_{i j} + \mathrm{tr} \left( K \right) K_{i j} - 2 K^{i j} K_{i j} \right] $$

in first-order, flux-conservative form. The form of the equations being solved is based on the first-order Bona-Massó formalism presented in https://arxiv.org/abs/gr-qc/9412071 and https://arxiv.org/abs/gr-qc/9709016, in which symmetric hyperbolicity is maintained by additionally evolving $A_k$ (the derivative of the lapse function), $B_{k}^{i}$ (the derivative of the shift vector), and $V_k$ (used for the indirect enforcement of the momentum constraint). The wv_vacuum_einstein equation struct implements the Bona-Massó equations more-or-less directly as written, while the wv_vacuum_einstein_conformal struct implements a conformal decomposition of the spatial metric tensor:

$$ \gamma_{i j} = \frac{1}{\chi^2} \widetilde{\gamma_{i j}}, $$

with conformal factor:

$$ \chi = \left( \det \left( \gamma_{i j} \right) \right)^{\frac{1}{6}}, $$

similar to what is done within BSSN and CCZ4 codes. We find that evolving the conformally-decomposed forms of the variables increases the numerical robustness of the solver in tests involving singularities (such that both single and binary black hole tests can run for a significantly longer period of time without crashing). We augment the system with the following standard evolution equation for the conformal factor (as used in standard BSSN and CCZ4 codes, e.g. https://arxiv.org/abs/1106.2254), appropriately recast into first-order, flux-conservative form:

$$ \partial_t \chi = \frac{1}{3} \alpha \chi \mathrm{tr} \left( K \right) - \frac{1}{3} \chi \partial_k \beta^k + \beta^k \partial_k \chi. $$

Errors in the ADM Hamiltonian:

$$ \mathcal{H} = R + \left( \mathrm{tr} \left( K \right) \right) ^2 - K_{i j} K^{i j} $$

and momentum:

$$ \mathcal{M}_i = \gamma^{j k} \left( \partial_l K_{i j} - \partial_i K_{j l} - \Gamma_{j l}^{m} K_{m i} + \Gamma_{i j}^{m} K_{l m} \right) $$

constraints are corrected using hyperbolic divergence error correction, analogous to what is done in the Munz approach to electromagnetism. With the inclusion of the conformal factor, we are able to apply gauge conditions in both the alpha-driver:

$$ \partial_t \alpha = - \mu_{\alpha_1} \alpha^{\mu_{\alpha_2}} \mathrm{tr} \left( K \right) + \mu_{\alpha_3} \beta^i \partial_i \alpha, $$

and gamma-driver:

$$ \partial_t \beta^i = \eta_1 B^i, $$

$$ \partial_t B^i = \mu_{\beta_1} \alpha^{\mu_{\beta_2}} \partial_t \hat{\Gamma}^i - \eta_2 B^i, $$

forms, for contracted conformal connection coefficients $\hat{\Gamma}^i$. By default, we select hyperbolic gamma-driver conditions with $\eta_1 = \frac{3}{4}$, $\mu_{\beta_1} = 1$, $\mu_{\beta_2} = 0$, and $\eta_2 = 1$, and through the gkyl_spacetime_slicing struct we provide support for geodesic (GKYL_GEODESIC_SLICING), harmonic (GKYL_HARMONIC_SLICING), and $1 + \log$ (GKYL_1PLUSLOG_SLICING) slicing conditions on the lapse function, with the latter (corresponding to $\mu_{\alpha_1} = 2$, $\mu_{\alpha_2} = 1$, and $\mu_{\alpha_3} = 1$) being particularly well-suited to both single and binary black hole spacetime evolutions.

Validations have been performed against the 1D Apples with Apples tests of https://arxiv.org/abs/0709.3559 (specifically the linearized gravitational wave test and the highly stringent expanding Gowdy wave test, with considerably improved performance in the latter test from the inclusion of the conformal factor evolution), as well as against both static single black hole tests (Schwarzschild and Kerr) and dynamic binary black hole tests (Brill-Lindquist initial data, both with and without gamma driver gauge conditions, and therefore both with and without moving punctures and dynamic excision boundaries). In all cases, we find agreement with analytic or numerical (BSSN) reference solutions wherever the code doesn't crash, and we find generally improved stability in the conformally-decomposed form of the equations.

[ammarhakim]

Ready to merge