Numerical Methods in GR

Ian Hawke

github.com/IanHawke
orcid.org/0000-0003-4805-0309
STAG, University of Southampton

ianhawke.github.io/slides/et-2022

The aim

Numerics:

given model,
find observable
to given accuracy
using minimum resources.

ETK:

within GR,
find GW signal from binary merger
to 1 radian phase accuracy
using minimum CPU time.

Hammond+, 2205.11377.

There is no generic best model. We want to find a specific model prediction as efficiently as possible, often to estimate parameters of the model. It's clear that as we change model, or the observable of interest, or the level of accuracy required by observations, or the resources available, then the best numerical method could change.
The ETK is founded on design decisions taken 20-25 years ago, built around assumptions from LIGO. It is designed to work within GR or similar theores, for finding gravitational wave signals, with accuracy assumptions from the early LIGO era, and working with HPC resources built around CPUs.
We have to remember that "resource" here includes human time. We have to balance the time needed to improve the methods with "just" throwing more compute time at a problem.

The models

EFEs $\implies$ wave equation $$ \partial_{tt} \phi = c^2 \nabla^2 \phi; $$
Charge conservation $\implies$ advection $$ \partial_t \phi + \nabla_k \left( \phi v^k \right) = 0; $$
$\nabla_a T^{ab} \; \implies$ balance laws $$ \partial_t q + \nabla_k f^{(k)}(q) = s(q); $$
Exchange $\implies$ relaxation $$ \dot{q} = \tau^{-1} s(q). $$

When looking at the numerical methods literature, and developing and testing methods, it's really useful to look at toy models and know why they're relevant. A number are clear: information is pushed around by waves and by matter, so the wave equation is obviously important for gravitational waves, and the advection equation needed for matter models. Most matter models will also give nonlinear balance laws, for which the standard toy model is Burgers' equation where the flux term is quadratic in terms of the conserved quantity.
In addition, information is exchanged between spacetime and matter, or between different components of the matter. This often leads to relaxation towards an equilibrium on a specific timescale. This timescale can be really important for the numerics, as we'll talk about later.

Split out time

Method of Lines (thorn MoL):

start from spacetime;
discretize space, giving ODEs;
solve ODEs.

(Dis)advantages:

MoL simplifies coupling;
causality restricts $\Delta t$;
width of stencil affects cost;
accuracy, efficiency issues.

Evolving just the spacetime is complex enough. However, there are many weakly coupled systems and models that we care about, including fluids, magnetic fields, radiation, gauge, and more. Monolithic codes that deal with everything at once are hard to write, hard to maintain, and hard to analyse. Instead it's best to split the code into separate sections - thorns, in the Cactus language - which are easier to work with.
The central idea is to split space and time at the discrete level. This is semi-discretization, or the Method of Lines. First we introduce some finite representation of the spatial dependence of the fields, whilst keeping time continuous. Think about introducing a grid of points on each slice of the manifold. The result of this is to convert the PDEs describing the model to ODEs. This system is large, and couples (for example) all the neighbouring points on the grid, but is a lot easier to solve numerically in a consistent and accurate fashion.
Solving the ODEs introduces the discretization in time, and hence the timestep. Causality means that the representative lengthscale introduced by the spatial discretization will be linked to the maximum stable timestep that can be taken: a given point on the future slice needs to be updated using all points from the previous slice that could potentially be causally connected to it. The more timesteps you have to take, the higher the cost and the lower the accuracy.
The ETK implements this idea through the MoL thorn. It's central to the ideas I'll present here. When absolute efficiency and accuracy are essential, particular when working with less accurate methods, it can be better to avoid using this approach.

Fields as data

To represent $\phi(x)$ as finite data:

Finite differences. Points $x_i$, $$ \phi_i = \phi(x_i). $$
Finite volumes. Cells $[x_{i-1/2}, x_{i+1/2}]$, $$ \hat{\phi}_i = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} \text{d} x \, \phi(x). $$
Finite elements. Basis functions $P_m(x)$, $$ \phi(x) = \sum_m \hat{\phi}_m P_m(x). $$

ETK works best with FD and FV.

With the time advance dealt with by converting everything to ODEs, specifying the rest of the numerical method depends on how the spatial derivatives are represented discretely. This is directly linked to how the continuous fields describing the model are represented discretely.
The first, and usually simplest, option is to represent fields by introducing a set of gridpoints covering the slice and storing discretely the values of the fields at these points. Regularly spaced grids in some coordinate system are simplest and often lead to the fastest expressions to numerically evaluate; however, they may not be the most efficient as we'll touch on later.
Another approach that's often used for matter models is to think of the slice as split into sub-volumes called cells, and to store the average value of the fields within each cell. This has specific advantages when the fields lose differentiability, for a slight increase in complexity and cost.
The third approach is to express the fields in terms of a function basis. Standard examples are Fourier series, particularly truncated Fourier series, where a function is expressed as a finite sum of trigonometric functions. We then store the mode coefficients. We can do this over the whole domain at once, or apply it again on subdomains. This approach needs much more set-up and can be very costly, but can also give huge gains in efficiency.
There are also particle methods. Despite some recent successes of SPH within numerical relativity - see Sphincs-BSSN - from the point of view of this classification SPH is closer to a function basis method than a traditional particle method. In my opinion, particle methods have very limited applicability to our problems.
The ETK was designed with finite difference methods in mind, using the grid point representation. Adapting it to finite volume methods is not difficult. Function basis methods are used in places, particularly for initial data, but are not a natural choice within the ETK for evolution.

Finite differencing

Approximate:

fields $\phi$ as point values $\phi_i = \phi(x_i)$;
polynomial interpolation gives $\hat{\phi}(x)$ from $\{\phi_i\}$;
derivative from $\partial_x \phi \simeq \partial_x \hat{\phi}$.

See Kranc, NRPy.

Example:

$$ \begin{aligned} && \partial_t \phi + \partial_x (v \phi) & = 0 \\ \rightarrow && \partial_t \phi_i & = \frac{-v \left( \phi_{i+1} - \phi_{i-1} \right)}{2 \, \Delta x}. \end{aligned} $$

Start from the "field values at grid points" viewpoint. Finite differencing is the standard method used in many codes, particularly for evolving the spacetime within the ETK. Given the field values at the points, we can fit an interpolating polynomial through the values and then approximate the derivative of the field as the derivative of the interpolating polynomial. For most cases we can explicitly calculate the derivative in terms of the neighbouring points, so this is very cheap.
The more points we use the higher the order of the interpolating polynomial and hence the better the accuracy of the approximation, in general, but the higher the cost.
If developing a new code you can find the finite differencing coefficients on Wikipedia, for example. However, it's simplest - within the ETK - to use one of the code generation frameworks such as Kranc or NRPy to generate efficient finite differencing code.

Order of accuracy

Finite differencing errors $\propto \Delta x^p = N^{-p}$:

$$ \partial_t \phi_i + \frac{v \left( \phi_{i+1} - \phi_{i-1} \right)}{2 \, \Delta x} = \mathcal{O} \left( \Delta x^2 \right). $$

$\Delta x$ "small": higher order usually better. However, costs increase:

$$ \partial_x \phi \simeq \frac{-\phi_{i+2} + \phi_{i-2} + 8 (\phi_{i+1} - \phi_{i-1})}{12 \, \Delta x}. $$

Usual to plot convergence with $N$ or $\Delta x$. Really resource (eg CPU time) that matters.

The order of the interpolating polynomial is directly linked to how the error of the method scales as we change the grid resolution. It's essential to check the correctness of the method by measuring how the error changes on a known solution. On a log-log plot, the power-law scaling of the error shows up as a straight line. Once we have confidence in the implementation, this allows us to estimate the error in simulations where the exact solution is not know.
However, higher accuracy using higher order methods means higher costs. We need to check that the extra cost is worth it. By plotting error against CPU time we see that - for this model and example, and ignoring implementation overhead - for any moderate level of error, the higher cost is worth the effort.

Phase errors

Finite differencing errors:

dissipation, reducing amplitude;
dispersion, changing phase.

Phase error dependence on $\Delta x$ can be estimated.

For 1kHz wave for 1s, 1% phase error:

$2^{\textrm{nd}}$ order	$4^{\textrm{th}}$ order	$6^{\textrm{th}}$ order
$10$m	$265$m	$784$m

Phase errors dependence on $\Delta x$ can be estimated.

For 10kHz wave for 1s, 0.01% phase error:

$2^{\textrm{nd}}$ order	$4^{\textrm{th}}$ order	$6^{\textrm{th}}$ order
$0.3$m	$47$m	$247$m

Looking at how the error scales with the grid resolution is not the only way of approaching matters. If we look at a single Fourier mode - essentially a sine wave - then we can look at how the amplitude and phase of the wave changes with resolution and order of the method. For gravitational waves in particular the phase accuracy is crucial for parameter estimation.
We can see from evolutions how the waves are damped by the numerical error, and that this is the most obvious effect. However, as the mode number increases the phase error from dispersive effects also increases, and this can be the most signficant effect.
Kreiss and Oliger showed back in the 70s that, for linear finite difference methods, the phase error can be approximated as a function of the grid spacing and the order. Using this we can bound the required grid spacing for a given wave. As a rough estimate for a LIGO wave, we're looking for a 1kHz wave over 1 seconds, which leads to a required grid spacing of around 10 metres for the simple second order method. This needs to hold over all the grid between its generation and extraction, which leads to an absurdly large grid and is totally impractical. We see the bound is much more relaxed with higher order. Unsurprisingly, production codes using finite differences use sixth or eighth order schemes as a matter of course.
While this is bad enough, the task for numerical relativity in the next decade will switch to producing waveforms for the Einstein Telescope and Cosmic Explorer. In this case we may need to resolve waves to both higher frequencies and better phase errors, making higher order even more important.

Spectral

Higher order is better: use all information. Mode expand:

$$ \phi = \sum_i \hat{\phi}_i \sin \left(\frac{2 \pi i x}{L} \right). $$

Then $\partial_x \phi \simeq D \hat{\phi}$.

Spectral error $\sim \exp(-p N)$.

Accuracy of full evolution often constrained.

ETK has spectral initial data (TwoPunctures, LORENE), no evolution.

If higher order is better, why not use all the information available to get the highest order possible? This is the essential idea of a spectral method, where the approximate derivative computed anywhere depends on every single data point.
In order to practically use spectral methods for nonlinear problems it's usual to use the function basis approach rather than grid points, also known as collocation points for spectral methods. In either case, the result of using all the information means the discrete spectral method can be written in terms of a differentiation matrix, which is typically dense.
The key advantage over finite differencing methods is the speed with which the error converges with the amount of computational steps. The comparison is between N grid points for finite differencing compared to N mode coefficients for spectral methods. But the convergence rate is exponential - the second plot is linear in the degrees of freedom rather than logarithmic - which is a huge advantage.
However, eventually the error from the spatial derivatives is not the main source of error. Because I'm using a toy model here this comes in earlier than it would with a more complex system, but it's always going to be the case that eventually the nonlinear couplings and error from the time integrator will stop the extremely rapid convergence. From the efficiency plot, however, we should expect spectral methods to be better than finite differencing. However, they are much more complex to set up and get working.
As noted earlier, the ETK isn't best set up to do evolutions with spectral methods, although they appear in initial data and analysis codes which are used. The standard spectral code in the field for evolutions is SpEC.

Discontinuities

Any high order method shows Gibbs' oscillations at jumps.

No convergence with resolution.

Use piecewise polynomial reconstruction (eg WENO) to avoid oscillations.

Converges with resolution.
Formally first order at jumps.
Much more expensive.
Data dependent - hard to analyse.
Extension to higher dimensions an issue.

See GRHydro, IllinoisGRMHD, Spritz, WhiskyTHC.

So, why don't we all use spectral methods? Well, they have major disadvantages. Here we'll mention two. The first is when the fields lose differentiability. This occurs naturally with matter models, particularly fluids.
This is why we think about the methods in terms of interpolating functions or polynomials. The standard Gibbs' effect is usually introduced in terms of Fourier series but can be seen in other function basis representations. That is, when a finite amount of data is used to represent a discontinuous jump in a field, the truncated series approximation will oscillate near that jump. When we differentiate the approximation in order to approximate the derivative of the original field, those oscillations are amplified. As Gibbs' oscillations do not converge with amplitude as the order of the truncation is increased, this means the oscillations in the approximate derivative get worse as more computational effort is used.
For finite difference methods these oscillations exist but are confined to a region near the jumps. As spectral methods rely on a dense differentiation matrix, the oscillations induced contaminate the entire domain.
The alternative that avoids, or at least minimises, the oscillations, is to change the data used in the interpolating polynomial to avoid the jump. This ensures that Gibbs' oscillations are never introduced, so the derivative is fine. However, it means that the reconstruction step now depends on the data itself. There are many ways of doing this, from second order slope limited methods to high order WENO methods. Whilst they are essential for discontinuous data they are much more costly - in some cases an order of magnitude or more - than finite differencing.
These sort of reconstruction methods can be found in a number of thorns in the ETK. If developing your own code then for initial tests I would recommend taking the slope limited methods from, for example, GRHydro, but higher order methods will be essential eventually. For those I would recommend starting with the WENO methods which are most generically implemented in WhiskyTHC.

Parallel compute

Classic CPU (Carpet etc): split domain.

Each CPU fast, high memory.
Need ghostzones, linked to stencil size.
Communication is bottleneck.
Problem for higher order schemes.

GPUs (CarpetX etc): create tasks.

Each GPU slow, low memory.
$N_{\textrm{GPU}} \gg N_{\textrm{CPU}}$.
Overload GPUs via pipeline.
Task based approach hides communication.

The second problem with spectral methods is an implementation question. The number of fields that need evolving in a modern GR model, and the tight accuracy constraints we need, mean there's too much computation and memory to do a simulation on a single core of any machine. This means we have to use parallel computing.
The standard model for parallel computing, which is built into the ETK, is domain decomposition. Take the grid covering the spatial slice. Split the grid into chunks, or patches, and put each patch onto a different core. The question is what happens at the boundaries of the patch. In order to compute the spatial derivative we need the data from the neighbouring points. The standard approach is to introduce ghostzones: points which will live on both cores, used on both, but evolved only on one. After the points have been evolved, the ghostzone data needed to be communicated from one core to the other.
Compared to computation, communication is really expensive. It's often the main bottleneck in an evolution, and as more cores are used it becomes worse. The higher the order of the numerical method, the wider the computational stencil, and hence the more the communication costs. This is significant problem for spectral methods in principle, although there are ways around it in real codes.
It's important to note that the next generation of HPC machines are likely to be built around GPUs rather than CPUs. GPUs are slower, cheaper, and with much less memory than CPUs. However, there will be orders of magnitude more available. Here, a pure domain decomposition technique will fail because of the communication. Instead it's standard to split the problem both into patches and into tasks on each patch. By ensuring the number of tasks is greater than the number of cores, the communication between patches can be hidden behind the computation taking place in the pipeline of tasks.
For now the ETK is a pure CPU framework, although (with a lot of effort) a user could put tasks onto GPUs within their own code. The CarpetX project should build GPU capabilities into the ETK much more systematically. The key takeaway is the importance of communication in the cost of HPC simulations.

Nonlinearity

Formally need to work with weak form

$$ \begin{aligned} && \partial_t \phi + \nabla_k f^{(k)}(\phi) & = 0 \\ \implies && \frac{\text{d}}{\text{d}t} \int_V \phi + \oint_{\partial V} \hat{n}_k f^{(k)}(\phi) & = 0. \end{aligned} $$

Discrete version, 1d:

$$ \frac{\text{d}}{\text{d}t} \hat{\phi}_i + \frac{1}{\Delta x} \left[ f_{i+1/2} - f_{i-1/2} \right] = 0. $$

Finding the flux $f_{k \pm 1/2}$ depends on the model (for fluids see GRHydro etc).

All the previous discussion has focused on linear models such as the wave and advection equations. However, matter models such as fluids naturally lead to nonlinear balance laws. The nonlinearity means that information is propagated at a speed that depends on the data, meaning that generic initial data will lead to discontinuities forming. At this point the strong form of the PDE does not make sense. Instead we need to use the weak form which follows from the finite volume representation shown earlier, where we integrate the PDE over the cell volume.
We saw earlier that it's essential to reconstruct the data using a data-dependent method such as a WENO scheme. When applied to the weak form we use this to compute the flux through the cell boundaries. This can be done in multiple ways. The first standard one is to compute a flux at the cell centre, reconstruct to the boundaries, and average appropriately. The other standard approach is to reconstruct the fields to the boundaries and then solve the Riemann Problem to say how the fields evolve when jumping at the boundary, from which the flux is computed. The precise approach does not matter - the key is that the same cell boundary flux is used when updating neighbouring cells. This ensures conservation and is essential for a discontinuity to propagate at the right speed.

Consistency

Using the right weak form is essential.

$$ \begin{aligned} && \partial_t q^n + \frac{n}{n+1} \partial_x q^{n+1} & = 0 \\ \implies && \partial_t q + q \partial_x q & = 0. \end{aligned} $$

Strong solutions agree when continuous; inconsistent at shocks.

Can use entropy pairs, path-consistent methods for complex cases. Links to SummationByParts.

Getting shock propagation speeds right is essential and depends both on how the model is written and how the numerical method implements it. A standard example is shown on the slide, where - if we were only dealing with differentiable solutions - the results should be identical for all values of n. However, as soon as a discontinuity forms, its speed depends strongly on the value of n. This means that manipulating a balance law for numerical advantage has to be done with care, as both the model and the numerics need to be consistent.
This links to other quantities or concepts that should be constrained or bounded. As we've seen, in MHD the divergence constraints are also crucial in order to get good numerical results. Other examples would be the second law of thermodynamics, or the constraints resulting from the 3+1 split. However, conservation constraints are probably the most extreme, and must be enforced at the discrete level.
For more complex matter models it is possible to have shocks even when the model has equations not in conservation law form. In this case we need to look for entropy functions which are non-increasing and change only at shocks, and which encode the underlying dissipation of the model. It is then possible to construct path consistent schemes. This is often linked to methods that ensure that energy functions are conserved discretely, which are implemented within the SummationByParts thorn of the ETK.

Spectral element

Want high order, low communication. Apply weak form to basis expansion $q(x, t) \sim \hat{q}_m(t) P_m(x)$, getting

$$ M \partial_t \hat{q} + S^T f(\hat{q}) = - \left[ F \right]_{i-1/2}^{i+1/2}. $$

Spectral accuracy with modes.
Only communication through RHS.
Expensive, complex.
(Solvable?) issues with shocks.

No code in ETK (?) - see Spectre.

For linear, smooth, models we saw the huge advantages of the spectral methods. There's been a lot of work tackling their main disadvantages - the communication for parallel computing, and dealing with shocks.
This leads to spectral element methods, of which the currently focus is on Discontinuous Galerkin. As in finite volume methods the domain is split into subdomains. The PDEs are multiplied by a test function - a smooth function whose precise form does not matter. The fields and the test function are then expressed in terms of basis functions. Integration by parts is then used over the subdomain to shift spatial derivatives from the fields onto the test functions. The resulting discrete scheme still depends on dense matrices, but all the coefficients are local to the element - the subdomain. The only coupling to neighbouring elements is through the flux term, so communication is minimized.
In a sense, DG methods are better than, for example, WENO methods, as the "reconstruction" step is not needed. By keeping all the modal coefficients around the fields are known everywhere all the time. However, this comes with significant extra cost. The number of degrees of freedom is now the number of elements times the number of modes. Care is needed to interpret DG results - never trust comparisons of gridpoints to elements.
Dealing with shocks remains an issue, although limiters and hybrid methods have been used with a lot of success to get around the oscillations. The downside of most of these fixes is they break the locality advantages by coupling more elements, or otherwise increase communication. However, these problems should be confined to near shocks. It seems likely that, in the long run, DG methods will be better than the alternatives. However, their complexity and cost currently work against them.
As with spectral methods, the ETK isn't best suited to DG or spectral element approaches. The Spectre code is the field-specific framework to look at.

Mesh refinement

Errors $\sim \Delta x^p$, but also local.

Refine mesh where errors large (Carpet(X)).

Boundaries need interpolation - source of error.

Physics not Cartesian: overlapping patches complex, but help (Llama, BH@Home).

So far I've emphasized how the errors scale with the grid size, but haven't specified how the grid itself is set up. The finer the grid the higher the computational cost. Due to the timestep restriction, halving the spatial grid step means increasing the computational cost by a factor of 16.
By concentrating the computational grid where the error is largest we can minimise the cost. This typically focuses on the compact objects, near the neutron stars or black holes. This is done via mesh refinement, where multiple grids are used with different grid spacings. In the ETK, this is implemented by Carpet. The data from the finer grids is the most accurate, so used by overwriting the data on coarser grids. The boundary data on the finer grids needs to be interpolated from the coarser grids.
The boundaries can be a significant source of error. In addition to the interpolation issues, consider the phase errors discussed earlier. Waves propagate at slightly different speeds depending on the numerical resolution. This leads to reflections and diffractions from the boundaries. Whilst this can be mitigated by high order methods, its impact on lower order schemes, particularly when matter is involved, can be severe.
Some of these issues can be significantly reduced by changing the coordinate system the grid is in - the wave will be roughly spherical, so a spherical grid is better. However, getting a useful grid set up can be technically painful. Look carefully at the Llama project within the ETK, but also BlackHoles@Home.

Stiff sources

For $\dot{q} = -\tau^{-1} f(q)$,

$$ q^{\color{green}{n+1}} = q^{\color{green}{n}} - \frac{\Delta t}{\tau} f(q^{\color{green}{n}}). $$ Explicit: fast, unstable for $\tau \lesssim \Delta t$.
$$ q^{\color{red}{n+1}} = q^{\color{red}{n}} - \frac{\Delta t}{\tau} f(q^{\color{red}{n+1}}). $$ Implicit: slow, stable $ \forall \tau$.

Stable does not mean accurate!

IMEX methods can be added to MoL.

When there is a big difference in the fastest and slowest timescales in the problem it is said to be stiff. This is a problem, particularly when the physics on the slowest timescale is interesting or essential, as the simplest and fastest numerical methods need to resolve the fast timescale for stability. This can make the source terms, which we have so far ignored, crucially important for the stability and accuracy of the problem. Here we're thinking of matter problems in particular, where reactions or radiation can add a relaxation term to the problem, and the field is expected to relax exponentially quickly on a particular timescale to equilibrium.
An example would be the Euler methods. The explicit Euler method solves the ODE simply and cheaply. However, it is only stable if the timestep is comparable to or smaller than the relaxation timescale. If the relaxation scale is faster then the numerical method will overshoot and the oscillations will explode with time.
The alternative is to use an implicit method. As the unknown data appears nonlinearly within the algorithm we need to use a nonlinear root-finder to solve each point at each step, which is expensive - usually an order of magnitude worse than an implicit scheme. However, implicit methods can be stable for arbitrary timesteps. For really stiff problems they can be much more efficient.
However, an implicit method being stable does not mean it is accurate. It's usual for the relative error of the implicit scheme being large even when the result is stable. If the impact of the fast behaviour on the observable of interest is small then implicit schemes are a good choice. However, this can be hard to judge or predict - see the current disagreements about the impact of bulk viscosity.
There are methods that split the system of equations into slow parts that can be solved explicitly and fast parts that need an implicit solver. The implicit-explicit (IMEX) methods can then be used to get the best efficiency for the whole system. Within the ETK and particularly MoL the focus is on explicit methods. It is possible to include IMEX methods within this framework - it has been done - but they're not included by default.

Toy heat

The Cattaneo model for heat is stiff:

$$ \begin{aligned} \partial_t T & = \partial_x q, \\ \partial_t q & = - \tau^{-1} \left( q - \kappa \partial_x T \right). \end{aligned} $$

For $\tau \ll 1$ write $q = \kappa \partial_x T + \mathcal{O}(\tau)$, finding

$$ \partial_t T = \kappa \partial_{xx} T + \tau \kappa^3 \partial_x^{(4)} T. $$

Solve heat equation explicitly whilst $\Delta t \gtrsim \kappa \Delta x^2 / 2$.

Tricks exist to adapt the model to the numerics. I mentioned bulk viscosity. Resistivity in non-ideal MHD models is another stiff source that can be approached this way. Here I'll use the example of heat conduction in non-ideal hydro, for which there's a simple toy model.
The Cattaneo equation is a causal model of heat propagation. Temperature evolves through the gradient of the flux; the flux evolves in relaxation form so that Fourier's law is approached on a relaxation timescale. When the timescale is small we can write the flux in the Fourier law form plus a small correction. Plugging this back in gives us the standard heat equation plus a correction term proportional to the relaxation scale.
The point is that standard explicit numerical methods work well for the heat equation when the grid spacing is large. So when the relaxation scale is much smaller than the timestep, the corrected heat equation is physically acceptable and numerically preferrable. However, the timestep restriction for the heat equation, due to its parabolic nature, is much more restrictive. Therefore, when the relaxation timescale is resolvable it will be more efficient to solve the original Cattaneo equation.

Adaptive Model

Changing model can reduce costs by

reducing size of state vector;
increasing allowed timestep.

Choose model adaptively to improve numerics while capturing physics.

Examples:

Low Mach approximation (Harpole);
Fast reactions as bulk viscosity (Hammond, Celora, Most+);
Resistive MHD as ideal plus corrections (Wright).

Asymptotic preserving models and schemes are needed.

This leads into the idea of Adaptive Model Refinement. In different regions of spacetime it may be possible to use simpler, less accurate, models to capture the physics to sufficient accuracy, purely for numerical advantage.
To some extent adaptivity is a core part of every large scale simulation we do. It's also standard to use the results of a simulation with a complex model as initial data for a simpler, cheaper model for long term simulations. This is known as co-simulation. Dynamics adaptation of the models and numerical methods used within the simulation is a known, but technically difficult, approach, that has been discussed in the past but never fully implemented within the ETK.
There's a range of examples in the literature where it's worked - the use of adaptivity for resistive relativistic MHD is the case that I think is the most successful. However, I want to highlight the Low Mach approach, which is a well-known application of the AMReX code. Here, acoustic waves are eliminated by solving elliptic equations. This only works when the Mach number is low. However, this holds true for neutron stars all the way up until merger, if the coordinate system is chosen appropriately. We would also expect it to be true throughout much of the disk. Adaptive model could give us order of magnitude increases in the timestep, and hence the associated accuracy, but the technical difficulty will be high.

Summary

The ETK is optimized for

finite difference/volume methods
using block-structured mesh refinement
on CPUs
for computing gravitational waves.

CarpetX should add GPU support.

Adaptive and hybrid methods becoming common, but complex.

Uncertainty quantification a key buzzword - lots of formal (non-GR) work to understand.

Hammond+; Most+.