Elliptic Constraint Solvers

Numerical initial data must satisfy the Hamiltonian and momentum constraints before an evolution begins. After freely specified fields and a decomposition are chosen, these constraints usually become coupled nonlinear elliptic equations. Solver convergence is only the first requirement: the discrete residual, continuum constraint, boundary error, and global charges must all be shown to converge for the intended slice. This page develops that computational evidence while leaving the continuum decompositions and their physical choices to Formulations of General Relativity.

Required background. Hypersurfaces and Foliations supplies the induced metric, extrinsic curvature, and constraint equations on a slice.

Helpful background. Formulations of General Relativity develops the continuum constraint systems, Ordinary Differential Equations supplies nonlinear iteration and tolerance concepts, and Relativistic Matter and Compact Objects supplies the matter models whose sources enter the constraints.

The elliptic system and its data

Use $K_{ij}=-\tfrac12\mathcal L_n\gamma_{ij}$ for the future-directed slice normal, and take four-dimensional general relativity with $\Lambda=0$. With a conformal split $\gamma_{ij}=\psi^4\tilde\gamma_{ij}$ and $K_{ij}=A_{ij}+\gamma_{ij}K/3$, define the conformal trace-free curvature by $A^{ij}=\psi^{-10}\tilde A^{ij}$, with indices on $\tilde A^{ij}$ moved using $\tilde\gamma_{ij}$. One representative vacuum Hamiltonian equation is

$$\tilde\nabla^2\psi-\frac18\tilde R\psi -\frac1{12}K^2\psi^5 +\frac18\tilde A_{ij}\tilde A^{ij}\psi^{-7}=0.$$

The momentum constraint supplies coupled equations for a vector potential or shift-like field. Matter adds energy and momentum sources with convention-dependent powers of $\psi$. The computational input record must identify the decomposition, freely specified conformal fields, inner boundaries or punctures, asymptotic conditions, matter data, coordinate domains, and regularity assumptions. Cook 2000, §§2.2.1–2.3 and §§3.1–3.4 lays out the major elliptic constructions and their boundary data.

The powers and signs can be checked before discretization. The vacuum constraint is $R+K^2-K_{ij}K^{ij}=0$, while

$$R=\psi^{-4}\tilde R-8\psi^{-5}\tilde\nabla^2\psi, \qquad K_{ij}K^{ij}=\psi^{-12}\tilde A_{ij}\tilde A^{ij}+\frac13K^2.$$

Substitute these expressions and multiply by $-\psi^5/8$ to recover the displayed elliptic equation. Every term has dimensions of inverse length squared because $\psi$ is dimensionless and $K_{ij}$ has inverse-length units in an orthonormal frame. This short derivation catches a mistaken conformal power that a converged nonlinear solver would faithfully solve. The curvature transformation and constraint reduction are given in Baumgarte and Shapiro 2010, §3.1, Eqs. (3.11)–(3.12) and (3.37).

Finite differences with multigrid, finite elements with adaptive residual control, and multidomain spectral collocation are all viable; Grandclément and Novak 2009, §§4–6 develop elliptic and spectral constructions used in numerical relativity. Newton’s method linearizes $F(u)=0$ as $J(u_k)\delta u=-F(u_k)$, but a small algebraic update does not imply a small continuum constraint. Control the linear solve relative to the current nonlinear residual, then show that tightening both solver tolerances changes the field or observable negligibly compared with grid refinement. A residual becomes a solution-error estimate only after accounting for the inverse operator or Jacobian sensitivity; the conditioning discussion explains why raw residuals and discretization errors cannot simply be compared.

Residuals need an independent operator

The solver naturally reports a discrete residual $F_h(u_h)$. Re-evaluating it with the same matrix mostly tests the linear algebra. Interpolate the solution to a denser or staggered grid and evaluate the original Hamiltonian and momentum constraints with an independently assembled operator. Report volume-weighted $L_2$ and $L_\infty$ norms, excluding no difficult region without an explicit geometric reason.

For a method of order $p$, a smooth solution should satisfy

$$\lVert u_h-u\rVert\sim Ch^p, \qquad p_{\rm obs}=\log_2 \frac{\lVert u_h-u_{h/2}\rVert} {\lVert u_{h/2}-u_{h/4}\rVert}.$$

When no exact $u$ is available, restrict fine solutions to a common grid before differencing. Spectral domains may show exponential decrease until limited by non-smooth puncture fields, surface discontinuities, or roundoff. A nonlinear residual driven below $10^{-12}$ on every grid while the continuum constraint stalls is evidence of discretization or boundary error, not success.

A Schwarzschild benchmark

Time-symmetric, conformally flat vacuum data have $K_{ij}=0$ and reduce to $\tilde\nabla^2\psi=0$ away from the puncture. Let $M_{\rm phys}$ denote the physical Schwarzschild mass and $M_{\rm geom}=G M_{\rm phys}$ its geometric length. The isotropic solution

$$\psi=1+\frac{M_{\rm geom}}{2r} =1+\frac{G M_{\rm phys}}{2r}$$

tests singular handling, asymptotic falloff, and ADM mass normalization. Solve it on at least three grids or basis orders, move the finite outer boundary, and compare $\psi$, the independent Hamiltonian residual, ADM mass on several coordinate spheres, and—if an inner surface is represented—the horizon mass. The expected physical ADM mass is the input $M_{\rm phys}$; extrapolation in outer radius separates truncation from finite-domain error.

For a regular exterior-domain test, excise the ball $r\lt a$ with $a=M_{\rm geom}/2$. The exact throat condition is $\partial_r\psi+\psi/(2a)=0$ at $r=a$, and the exact monopole outer condition is $\partial_r\psi+(\psi-1)/R_{\rm out}=0$ at $r=R_{\rm out}$. Here $\partial_r$ points toward increasing radius at both surfaces; it is opposite to the computational domain’s outward normal at the inner surface. These Robin conditions admit the analytic solution at every finite $R_{\rm out}$, so this test isolates discretization and nonlinear error. An approximate outer condition for nonspherical data requires a separate radius study. Prescribing $\psi=1$ at a finite outer radius instead changes the problem and introduces a systematic mass error.

For coupled binary or stellar data, no closed form may exist. Then use a second discretization or domain layout, global charge relations, and evolution stationarity as additional evidence. Cook’s discussion of black-hole and neutron-star constructions makes clear that free data affect physical content even when constraints vanish Cook 2000, §§3.1–4.3. Constraint convergence does not prove low eccentricity or low junk radiation.

Environment and rejection criteria

Pin the decomposition, domain maps, mesh or basis, quadrature, boundary conditions, nonlinear and linear solver versions, preconditioner, precision, stopping norms, initial guess, and code revision. A reproducible interface could be solve-constraints --config schwarzschild-puncture.toml --output run/, emitting fields, independent constraint norms, iteration histories, charge surfaces, refinement data, and environment metadata. This specifies a rerun record; it does not claim that an executable is supplied here.

Reject a solution if constraints fail to decrease at the designed order, results depend on nonlinear tolerance after it is well below truncation error, the outer-boundary sequence lacks a stable extrapolation, or the isolated Schwarzschild benchmark has ADM and horizon masses that disagree outside estimated error. For binary or radiative data, binding energy and exterior radiation can produce a genuine mass difference; compare it with the appropriate physical budget. Other failure modes include Newton convergence to an unintended branch, multigrid stagnation on unresolved punctures, spectral ringing at stellar surfaces, and a preconditioned residual that hides a large physical residual.

Diagnostic

Three grids reach a nonlinear residual of $10^{-13}$. Their independently evaluated Hamiltonian $L_2$ norms are all about $10^{-5}$, and the ADM mass changes by two percent when the outer boundary doubles. Has the initial-data calculation converged?

Solution

No. The algebraic solve has converged on each discrete problem, but neither the continuum constraint nor the global charge has entered a convergent regime. Refine the discretization, improve puncture or surface treatment, and perform an outer-radius sequence. Accept only when independent constraint norms decrease with the expected order and the extrapolated ADM mass is stable within the declared tolerance.

References

• Baumgarte, T. W., and Shapiro, S. L. (2010). Numerical Relativity: Solving Einstein’s Equations on the Computer. Cambridge University Press. DOI.
• Cook, G. B. (2000). “Initial data for numerical relativity.” Living Reviews in Relativity 3, 5. DOI.
• Grandclément, P., and Novak, J. (2009). “Spectral methods for numerical relativity.” Living Reviews in Relativity 12, 1. DOI.