Multi-Domain and Spectral Element Decomposition

Multi-domain methods divide a complicated spacetime region into elements whose coordinates and approximation spaces match local geometry. Spectral elements can achieve rapid convergence for smooth fields, represent spherical boundaries cleanly, and distribute work across processors. Their advantage disappears when a shock, stellar surface, puncture singularity, or map irregularity lies inside an element. Domain design is therefore part of the numerical model, and convergence must be tested under both polynomial enrichment and geometric subdivision.

Helpful background. Formulations of General Relativity supplies the evolved geometric fields, Ordinary Differential Equations supplies time integration and residual concepts, Computational Relativity supplies shared evidence standards, and Probability and Statistical Inference supplies error comparison across adaptive choices.

Elements, bases, and interfaces

Map each reference element $\boldsymbol\xi\in[-1,1]^d$ to a physical block. Expand a smooth field in tensor-product Chebyshev polynomials or another declared basis,

$$u_N(\boldsymbol\xi)=\sum_{n_1=0}^{N_1}\cdots \sum_{n_d=0}^{N_d}\hat u_{n_1\ldots n_d} T_{n_1}(\xi^1)\cdots T_{n_d}(\xi^d).$$

For a smooth analytic field, coefficient magnitudes and truncation error often decrease exponentially until roundoff. Limited differentiability produces algebraic decay. Record collocation points, quadrature, filtering, tensor-component frame, map Jacobians, and how neighboring elements communicate.

Continuous Galerkin methods share degrees of freedom; discontinuous Galerkin methods exchange characteristic numerical fluxes; multidomain collocation may impose penalty or matching conditions. Interface normals, characteristic orientations, penalty strengths, and mortar projections must be consistent. Grandclément and Novak 2009, §§2.6.1–2.6.6 and §§3.1.1–3.1.3 develops multidomain spectral approximation and elliptic applications in numerical relativity.

Choose h or p from regularity

$p$-refinement increases polynomial order on an element. $h$-refinement subdivides it. For smooth waves, $p$-refinement is efficient; for localized features or non-smoothness, subdivision and alignment are usually necessary. A stellar surface should coincide with an element boundary when a sharp EOS creates derivative jumps. A shock in GRMHD needs limiting or a finite-volume subcell treatment, not an unrestricted global polynomial.

Use a truncation indicator based on modal coefficient decay, but check it against an observable error. Coefficients can be small because a filter removed content, while the solution is biased. An hp controller should record thresholds, hysteresis, maximum orders, split directions, and every adaptation event.

A simple example explains the role of interfaces. The function $u(x)=|x|$ is continuous on $[-1,1]$ but its derivative jumps at zero. One polynomial on the whole interval cannot reproduce that corner exactly and its spectral coefficients decay algebraically. Split at zero: $u=-x$ on the left and $u=x$ on the right are each degree-one polynomials. The field is now represented exactly within each element. Whether its derivative jump is permitted by a PDE remains a separate interface-condition question; a good approximation space does not supply missing physical jump conditions.

Smooth and non-smooth benchmarks

Use a periodic interface-crossing wave as the quantitative smooth record:

$$\partial_tq+\partial_xq=0, \qquad q_{\rm ex}(t,x)=e^{2\pi i(x-t)}.$$

Evolve on $[0,1]^3$ with periodic boundaries and four equal hexahedral elements along $x$. Use nodal Legendre–Gauss–Lobatto DG, an upwind numerical flux, RK4, and $\Delta t=0.05/(p+1)^2$ for polynomial orders $p=4,6,8,10$ to $t=1$. The initial data and analytic reference are $q(0,x)=e^{2\pi ix}$ and $q_{\rm ex}$ above. Compare with the quadrature-weighted relative $L_2$ error $E_2$, the fundamental-mode phase error $\Delta\phi$, and the maximum normalized jump $J_\infty$ between independently evaluated left and right traces at every interface.

Accept the $p=10$ result only if $E_2\leq10^{-9}$, $\Delta\phi\leq10^{-9}\,\mathrm{rad}$, and $J_\infty\leq10^{-8}$, with each of the last two increments in $p$ reducing $E_2$ by at least a factor of ten. Independently sample 32 points per element to exclude agreement confined to collocation nodes. These limits specify the comparison; they are not observations.

Repeat each order with half and quarter time steps before interpreting this as spatial spectral convergence. RK4 contributes an algebraic $O(\Delta t^4)$ error floor; once that floor dominates, further $p$-refinement cannot maintain a factor-ten decrease. Apply the spectral-ratio requirement only while an independently reduced temporal error is below the spatial error. A plateau at controlled roundoff or time error is reported as such, rather than misclassified as a failed spatial discretization.

Tie strong scaling to the same exact wave by holding a $16\times8\times8$ element mesh at $p=10$ fixed for $P/P_0=1,2,4,8$. For weak scaling, tile $1,2,4,8$ unit-periodic cubes, each containing $4^3$ elements at $p=10$, and assign one cube per $P_0$ resources. Before collecting timings, qualify every layout against the exact wave using the same quadrature, interface flux, and timestep rule: its $E_2$, $\Delta\phi$, and $J_\infty$ must satisfy the limits above, including the independent off-node and timestep checks. The lower-order $p=8$ tile need not meet the stated $10^{-9}$ field-error target even when time error is reduced. After qualification, use one warm-up and seven measured repetitions without changing the numerical parameters.

Observed discrepancy: not available—specification only.

A complementary elliptic test should solve an analytic smooth scalar field on spherical shells and distorted blocks, compare the field and derivatives with the exact solution, and evaluate residuals between collocation points. Vary element boundaries and map distortion. Exponential convergence should persist until a precision plateau.

Next place a function with a derivative discontinuity inside an element, then repeat with the discontinuity aligned to an interface. The unaligned case should show slow convergence or Gibbs oscillations; the aligned case should recover rapid convergence within each smooth subdomain. This negative control demonstrates the method’s validity regime rather than advertising only its best case.

For hyperbolic evolution, propagate a wave across multiple interfaces and measure phase, reflected amplitude, constraint flux, and conservation. Kidder et al. 2017, §4, §§5.1–5.3, Figures 4–6 and §6 describes task-based DG infrastructure, representative tests, and scaling; it does not supply the hp controller or a shock-subcell algorithm. Fully adaptive spectral-element hp refinement and its accuracy tests are documented by Renkhoff et al. 2023, §§II–VI and Figures 1–14. An a posteriori DG–finite-difference subcell method for GRMHD is developed and tested by Deppe et al. 2022, §§4–6 and Figures 2–14.

Cross-check a selected solution with a uniformly refined finite-difference or finite-volume discretization. Agreement at matched error is more informative than agreement at matched point count because spectral degrees of freedom have different resolution power. For elliptic problems, also compare surface charges and an independently assembled continuum residual. For evolutions, vary time-step error separately from spatial order; an overly loose integrator tolerance can hide the expected spectral improvement.

Environment and failure criteria

Pin element topology, maps, bases, quadrature, polynomial orders, numerical fluxes or penalties, filters, hp indicators, subcell trigger, time integrator, CFL policy, precision, and code revision. A clean interface could be multidomain-test --config shell-wave.toml --output run/, emitting modal spectra, interelement jumps, residuals, adaptation history, reflection measures, refinement sweeps, and environment metadata. This specifies an implementation record without claiming that software was executed for this page.

Reject a result when convergence is inferred from filtered coefficients alone, interface jumps fail to decrease, a map loses regularity, or a non-smooth feature remains inside an unlimited spectral element. Other failure modes include wrong normal orientation, underintegrated nonlinear products and aliasing, penalty parameters outside their stable range, and adaptation that changes the physical boundary. Higher polynomial order can worsen an aliased or ill-conditioned problem.

Diagnostic

A stellar density profile shows exponential modal decay on every element except one crossed by the surface, where coefficients decay algebraically and oscillations persist as order increases. Should the polynomial order be raised further?

Solution

No. The surface limits differentiability inside that element, so pure $p$-refinement cannot restore exponential convergence. Move or split the element boundary to align with the surface, or use a robust subcell representation. Then repeat field, derivative, and global-charge convergence tests and verify interface conditions.

References

• Deppe, N., Hébert, F., Kidder, L. E., and Teukolsky, S. A. (2022). “A high-order shock capturing discontinuous Galerkin–finite difference hybrid method for GRMHD.” Classical and Quantum Gravity 39, 195001. DOI.
• Grandclément, P., and Novak, J. (2009). “Spectral methods for numerical relativity.” Living Reviews in Relativity 12, 1. arXiv.
• Kidder, L. E., Field, S. E., Foucart, F., et al. (2017). “SpECTRE: A task-based discontinuous Galerkin code for relativistic astrophysics.” Journal of Computational Physics 335, 84–114. arXiv.
• Renkhoff, S., Cors, D., Hilditch, D., and Brügmann, B. (2023). “Adaptive hp-refinement for spectral elements in numerical relativity.” Physical Review D 107, 104043. DOI.