Numerical Hydrodynamics: Part 2

Ian Hawke

ianhawke.github.io/slides/icts-2020
Additional material at github.com/IanHawke/icts-2020

What we covered

Balance laws are generic;
lead to shocks;
shocks appear in mergers;
start from Newtonian CFD;
GR gives
- increased cost;
- increased complexity;
Merger problem gives
- surface/atmosphere;
- lots more physics.

In the first lecture I wanted to emphasize the modelling approximations and limitations of what goes in to a numerical simulation. These come from physical and resource limitations. They lead us to consider the most important features of our models, and so what aspects we need to model in what order. The key points from the first lecture were that

conservation and balance laws are generic features of matter simulations;
these lead to shocks, and these shocks do appear in mergers;
using techniques from Newtonian CFD is a good place to start, but we need to consider issues specific to NSs and relativity, like artificial atmospheres, and the more complex equations of motion.

Today's lecture is going to be a deep dive into the numerical methods. We'll look at the standard techniques used in the field. The exercises associated with these lectures ask you to implement and test some of these techniques: that's the best way of getting to grips with the difficulties, limitations, costs, and complexities.

Balance laws

All the terms: $$ \partial_t {\bf q} + \nabla \cdot {\bf f} + A \cdot \nabla {\bf q} = \nabla \left( D \cdot \nabla {\bf q} \right) + {\bf s}. $$

In mergers

ignore diffusive term;
simple models don't have $A$ term;
source terms often local, well-behaved.

The general form of the equations of motion of relativistic matter can be framed as here. Some quantity is evolved using its flux, ${\bf f}$, additional non-conservative terms that push the matter around through the matrix $A$, a second order diffusive derivative operator with dissipation matrix $D$, and a source term ${\bf s}$ that typically represents exchange of "stuff" to other parts of the model. In the most complex cases the functional form of the matrices, sources or fluxes can depend not only on nonlinear algebraic operators of the evolved vector ${\bf q}$, but even on nonlocal operators. An example would be radiation transport where the optical depth is needed.
However, this is far more complicated than is generally needed for NS mergers, particularly through the inspiral stage. In these stages we are dominated by the flux terms. In the standard models (hydrodynamics, maybe including EM fields through the MHD approximation) we can write the equations in conservation law form. And, when radiation transport can be neglected, we can assume that the source terms are given in an algebraic form.
This means we can concentrate on the fluxes. From a mathematical and implementation point of view the most important terms are the principle part with the highest derivatives, so we will (for now) look at the conservation law form and the numerical methods needed for that.

Characteristics

For $\partial_t q + \partial_x f = 0$ get local speed $$ \partial_t q + \color{red}{\partial_q f} \partial_x q = 0. $$

For system $\partial_t {\bf q} + A \partial_x {\bf q} = 0$ diagonalize to get $$ {\bf w} = L {\bf q}, \quad \partial_t {\bf w} + \Lambda \partial_x {\bf w} = 0. $$ Extend to nonlinear case with care.

As a quick recap, remember that for the advection equation (where the flux is a constant velocity $v$ multiplying the quantity $q$) we can define characteristics: lines, of slope $v$, along which the solution is unchanged. We can use this locally for nonlinear problems like Burgers equation: with a general flux $f$, the characteristic speed is $\partial_q f$. This only makes sense away from discontinuities.
When we extend to the system case it's easiest to consider the linear problem first. When the flux is defined by a linear system, ${\bf f} = A {\bf q}$, then we can change variables to get the solution. A crucial aspect of the conservation laws is that they are hyperbolic: the essential consequence of this for now is that $A$ will be diagonalizable. This means we can define and work with the characteristic variables ${\bf w} = L {\bf q}$, where $L$ is the matrix of left eigenvectors. The characteristic variables obey uncoupled advection equations, with the advection velocity being $\lambda_i$, the eigenvalues of the matrix $A$. So we can use the solution from the advection equation: each characteristic variable propagates unchanged with a fixed speed. We can then find the solution for ${\bf q}$ by transforming back using the matrix of right eigenvectors.
Locally, in smooth regions, we can see how this would work in the nonlinear case. The Jacobian matrix $J = \partial_{\bf q} {\bf f}$ plays the role of $A$. The characteristic variables can therefore be defined pointwise, noting that the eigenvalues and eigenvectors are going to vary with space. These are then propagated forward and recombined, in principle. In practice the re-combination is difficult, as the correct eigenvectors to use depend on the new data in ways that can be hard to disentangle.
However, the viewpoint of using characteristics is fundamental to a lot of numerical algorithms. Whilst the full procedure outlined here is just about never done, the information is always useful to interpret how the algorithm is working.

Wave types

Rarefaction:
- nonlinear
- continuous;
Shocks:
- nonlinear
- discontinuous;
Contacts:
- linear
- discontinuous.
Compound waves:
- nonlinear, a mess;
- MHD, phase transitions.
Martí & Müller, Liv. Rev.

From the characteristic point of view we can distinguish different generic types of behaviour. If we start from simple initial data, such as the piecewise constant data illustrated on the left, then there's typically three different outcomes. For complex systems we'll see them all, and some more!
The first type of behaviour, on the left, is where the characteristics separate. This leads to a continuous spreading of the data. In hydrodynamics this is referred to as a rarefaction wave: the material is rarefied as it expands.
The second type of behaviour, on the right, is where the characteristics converge. This leads to a discontinuity, even when the initial data is smooth, as discussed in the first lecture. This is a shock, and typically is associated with a jump in all variables. Specifically in the hydrodynamic case we would see entropy and temperature increase through the shock wave, whereas across any other wave they would remain constant.
The final type of behaviour, in the middle, is where the characteristics move parallel to each other, maintaining the jump in the initial data. These linear waves are associated with jumps in material properties, or rotations of vector properties. In hydrodynamics these are the contact waves, and are associated with the advective velocity. The nonlinear waves are associated with the acoustic waves.
In MHD there is a more complex wave structure. There is still a central linear contact wave, but each nonlinear acoustic wave splits into two nonlinear "fast" and "slow" waves separated by a linear wave across which the EM fields can rotate. This additional complexity in the wave structure is harder to resolve numerically.

Questions...

Stretch out;
have a break;
add questions to the chat.

Weih et al, 1912.09340.

Grids and approximations

Finite differences:
store point values $q_i$.
Finite volumes:
store cell averages $\hat{q}_i$.
Finite elements (DG):
store modal coefficients $q^{(m)}_i$.

On to numerical approximation. We want to solve continuum partial differential equations on a computer. The solution, $q(x)$, needs - in principle - an infinite amount of information. That's because the continuum solution could, in principle, take totally different unpredictable values at every separate point $x$. This way of thinking is clearly useless when working with a computer with a finite amount of working memory.
Instead we consider ways of approximating the solution $q(x)$. There are essentially three that are used in the field.
The first is the point value approximation, used in finite differencing. In this case the domain is discretized into points, and we assume that we know (or want to compute) the value of the solution at each of these points. The solution is not known between the points. When we want to link the solution at different points - for example, when approximating a derivative - we have to impose the general behaviour of the solution between those points. For example, if we interpolate between points using a straight line then the derivative is the slope of this line.
The second approach is to split the domain into cells, and within each cell to store the average value. This is the finite volume approach. In this case some information is known at every point in spacetime: the average value near that point. However, the exact value is known nowhere. Again, to get the value of the solution at a specific point we have to impose the general behaviour of the solution. When restricted to a given cell this behaviour has to be consistent with the average value within that cell. The finite volume approach implicitly "smears out" the local behaviour of the solution, which is exactly what we want when dealing with discontinuous solutions like shocks.
The final approach has links to the finite volume approach, but is more often referred to as a finite element method. Instead of storing the average value of the solution within a cell, the moments of the solution with respect to some function basis are stored. Typically the zeroth moment would be the cell average of the solution, the first moment linked to its derivative, the second to its curvature and so on. This is fundamental to the Discontinuous Galerkin method that HP will talk about later. So far this hasn't been used much for relativistic matter simulations, for reasons we'll get to later.

Fluxes and telescoping

Integrate over cell $i \to \hat{q}_i$: $$ \frac{\text{d}}{\text{d} t} \hat{q}_i + \frac{1}{|V_i|} \oint_{\partial V_i} f(q) = 0. $$ Restrict to one dimension: $$ \frac{\text{d}}{\text{d} t} \hat{q}_i = \frac{1}{\Delta x} \left[ f_{i-1/2} - f_{i+1/2} \right]. $$ Gives discrete conservation.

To see the implications of the finite volume approach, where we split the domain into cells, let's use Gauss, or Stokes, to integrate our conservation law over the domain. We see that the integral over all space leads to all the spatial derivatives disappearing. The equation becomes an ODE, where the time derivative of the integral average of the solution over the domain is given by the surface integral of the flux through the domain.
This is the weak form of the conservation law. We'll talk about weak forms more generally later. Crucially, the weak form is used to remove spatial derivatives from the solution (which might be discontinuous). This allows us to talk about discontinuous solutions, like shocks, with confidence that the mathematics will work out.
We see that we're still not at a point where we can solve this equation. The ODE is defined by the surface integral of the flux. Evaluating this flux over the boundary of the domain requires knowing the solution at the boundary of the domain. But we don't know the solution there: we only know the integral average within the domain.
This brings us to the fundamental step in numerical methods for conservation laws. We need a prescription for computing this boundary flux. Once we have that, we can solve the ODE using standard methods.
To be concrete, let's reduce our domain down to one spatial dimension. We split the domain into cells. We label the cells with an integer $i$, and given the cell centre the coordinates $x_i$. The interfaces between the cells and the boundaries of each little subdomain. We label the coordinates of these interfaces with "half integers", so the boundaries of cell $i$ are at $x_{i-1/2}$ and $x_{i+1/2}$. We then write the ODE explicitly, showing that the average solution $\hat{q}_i$ within cell $i$ evolves due to the flux into the cell through the left interface at $i-1/2$ and the flux out of the cell through the right interface at $i+1/2$. Note that "into" and "out of" are conventions without physical meaning here, associated with assuming the material flows from left to right: the flux could be negative, meaning the material is flowing from right to left.
There is a central result from writing the conservation law in this form. If we have a discrete algorithm that updates in this fashion then it is conservative at the discrete level. To see this, imagine summing up all the integral averages in all the cells. This would give the integral average over the entire domain at that time. Now imagine doing this at a later time, after the solution has been updated using the flux differencing form here. Each internal cell interface has a flux that is used twice: the flux through the interface at $i+1/2$ is used to update cell $i$ and cell $i+1$. As each update has different signs, the internal contributions cancel out. This telescoping effect means the total integral average is only changed by what comes in through the boundaries of the full domain. Total conservation is retained. This is essential to ensure the Rankine-Hugoniot conditions are satisfied at the discrete level and any shock propagates at the correct speed.

Computing the intercell flux

Godunov:

$q(x) = \hat{q}$ in cell $i$;
$$ f_{i+1/2} = F(\hat{q}_i, \hat{q}_{i+1}); $$
$$ \partial_q f > 0 \implies f_{i+1/2} = f(\hat{q}_i); $$
Systems: use characteristic variables;
Nonlinear: solve Riemann Problem, usually approximately!

We're now at the point where describing our algorithm boils down to computing the flux through a single cell interface. We'll look at Godunov's method and variants as a starting point.
In Godunov's method we assume the solution is piecewise constant. That is, within each cell the value of the solution is the integral average. This gives us the value of the solution everywhere except at the cell interfaces, where there are two values that, in general, won't agree. We need to take these values and get the intercell flux.
To do this, think about the characteristics. For example, think about the advection equation. If the advection velocity $v$ were positive then the solution at the cell interface would always be given by the value to the left of the interface. More generally, if the characteristic speed is positive, the value of the solution at the interface is given by the solution to the left. From this we can compute the flux.
In the system case this is more complex. For a linear system we can convert to characteristic variables. The eigenvalues tell us which characteristic variables propagate to the right, meaning we should use the solution to the left of the interface, and vice-versa. This gives us the characteristic variables at the interface, from which we can convert back to get the solution, and hence the flux.
In the nonlinear case this gets yet more complex. In this most general case we are looking for a solution to the Riemann problem: what's the solution of the PDEs when the initial data is piecewise constant? This can be hard to compute exactly: in most relativistic cases it's impractical. It can be approximated using some, or all, of the characteristic information as in the linear system case. In relativity people often use simple approximations, as computing the characteristic information needed for more complex approximations is expensive.

Approximate Riemann Solvers

HLLE: to find $$ f_{i+1/2} = F(\hat{q}_i, \hat{q}_{i+1}) = F(q_L, q_R): $$

assume fastest speeds are $\xi_{\pm}$;
impose conservation;
$$ q_* = \frac{\xi_{+} q_R - \xi_{-} q_L - f(q_R) + f(q_R)}{\xi_{+} - \xi_{-}}; $$
get flux from appropriate state.

As an example of a simple approximation we look at the HLLE method. This assumes that there are two characteristic waves and that the solution jumps discontinuously across both. So the approximate solution will have three states: the initial left and right states (either side of the interface), and a central state. If the characteristic speeds have the same sign then the solution along the interface is given by one of the initial states, and everything is straightforward. If they have different signs then the solution along the interface is given by the intermediate state. This can be found by imposing conservation across each wave, leading to a simple formula for the intermediate state: from this the intercell flux follows.
The HLLE method has some nice properties. It's fairly cheap to compute, especially if simple approximations to the characteristic speeds are used. It's not specific to one model. It's linked to positivity preservation: quantities like density that should be positive will remain positive, when HLLE is used in a particular way.
However, HLLE has its problems. In particular, linear waves (like contact discontinuities) tend to be smeared out very badly. For MHD this can be a serious problem, as there's more linear waves.
More complex approximations either work to ensure better behaviour when the waves aren't discontinuities, or to ensure better capturing of the linear waves. For hydrodynamics this helps but isn't essential in many cases: for MHD it's often needed.

Questions...

Stretch out;
take a chance to refocus;
ensure you're ready for more detail;
add questions to the chat.

Bernuzzi, 2004.06419.

Dimensions, costs, and accuracy

Godunov isn't good enough:

Error $\propto (\Delta x)^1$.
Computational cost $\propto (\Delta x)^{-4}$.
Extrema are clipped.
GWs need better phase accuracy.
Need higher order methods. But...
...leads to problems with shocks.

We now have an algorithm that will work. Godunov's method assumes the solution is constant within each cell, and computes the intercell flux using an approximate Riemann Solver. This can be extended to two or three dimensions by solving a Riemann problem for each cell face.
However, the algorithm isn't very accurate. This wouldn't matter if we could throw computing power at it, but as we've seen we're going to be limited in the size of the grid cells we can use. We measure the accuracy by how quickly the error reduces as we reduce the grid size. If the error scales as the cell size $\Delta x$ to some power $p$ then we say the method is $p^{\text{th}}$ order accurate. Godunov's method is first order accurate: if we reduce the cell width $\Delta x$ by a factor of 2 then we reduce the error by a factor of 2.
Unfortunately, reducing the cell width by a factor of 2 does not mean the computational cost goes up by a factor of 2. In three space dimensions we have to reduce the grid size in each direction. The CFL stability limit discussed in the first lecture implies we also have to reduce the time step by a factor of 2. So the computational cost goes up by a factor of 16.
So, for cost and efficiency reasons we want to improve the order of accuracy. If we get up to fourth order then the accuracy will scale linearly with the computational cost, which would be great. Unfortunately, increasing the order of accuracy is hard, and it's hard for a crucial, physical reason: shocks.

Reconstruct-Evolve-Average

Rethink Godunov as three steps:

Reconstruct: $$\hat{q} \to q(x);$$
Evolve: $$q(x);$$
Average: $$q(x) \to \hat{q}.$$

Reconstruction loses accuracy.

A different visual way of thinking about Godunov's method (in one dimensions, at least) is the Reconstruct, Evolve, Average framework. We start with the cell average solution: within each cell the solution is constant. This is the reconstruction step: going from our cell averages to an assumed form for the solution everywhere. We compute the intercell fluxes. Locally, this approximately advects our solution at a certain speed. We now have a (roughly) piecewise constant solution that isn't aligned with the cell boundaries. This is the evolve step. To get back to our cell averages, the average step takes the cell average within each cell of our new solution.
In principle, most of these steps are exact. The averaging step is exact. The evolve step needs doing with some ODE solver, but that can be done to high accuracy. The evolve step also needs the computation of the intercell flux, but it's possible in simple situations to do that exactly. It's the reconstruction step that introduces the first order approximation. By enforcing that the solution is piecewise constant we significantly reduce the accuracy of the method.
To improve this, we need to improve the reconstruction. The next step up would be to use piecewise linear reconstruction: assuming the solution is a straight line, not necessarilly flat, within each cell. Then up to quadratic functions and higher order. This would give the improved accuracy we seek. Unfortunately it also introduces other problems when the true solution is discontinuous.

Monotonicity, Gibbs oscillations, and Godunov's theorem

Fourier Series:
discontinuities $\implies$ oscillations.

Don't converge with $\Delta x$;
Don't converge with more modes;
Don't go away with different function basis.

Monotonicity: scheme doesn't introduce oscillations.

Godunov's theorem: linear monotonic schemes are first order accurate.

Here's a standard discontinuous function: it's piecewise constant with one jump. If we look for a Fourier series representation of this function then the partial sum gives us a representation that oscillates around the true solution. These oscillations are concentrated near the discontinuity.
There are three obvious ways to get around this. First, split the domain into smaller cells. However, we see that there's no scale to the true solution, so these oscillations will still appear. They'll get squeezed closer to the discontinuity, but their magnitude won't reduce.
Second we could include more terms in the partial sum. The result is the same as changing the cell size: the oscillations get squeezed closer to the discontinuity without getting smaller.
Finally, we could change the function basis. Instead of using a Fourier series we could expand in standard polynomials, or Legendre polynomials, or something else. This doesn't help either: these oscillations are a generic feature.
These Gibbs oscillations will destroy our numerical accuracy. They don't converge with resolution: typically they blow up in a few iterations. Most depressing is Godunov's theorem. It asks when the solution might be monotone: that means, when the solution lies between the minimum and maximum of the original solution. Godunov's theorem says that a linear algorithm that is monotone must be first order accurate.
We have one loophole left by Godunov's theorem. It only applies to linear methods: that is, methods where the algorithm essentially does the same thing everywhere, ignoring the values and form of the data. The past 50 years and more of CFD has focused on better ways of locally modifying the algorithm based on the values and form of the solution itself, to avoid Gibbs oscillations.
It's this nonlinearity that makes numerical methods for matter models so much more complex and expensive than numerical methods for smooth solutions (such as the spacetime). It's typical for the computational cost of the matter evolution to be greater than that of the spacetime evolution, even though we're usually evolving far fewer variables, sometimes by an order of magnitude.

Slope limiting

Assume $q(x) = \hat{q}_i + \tfrac{x - x_i}{2} \sigma$.

Slope $\sigma$ could be

Upwind: $\sigma_{\text{u}} = \hat{q}_{i+1} - \hat{q}_i$;
Downwind: $\sigma_{\text{d}} = \hat{q}_{i} - \hat{q}_{i-1}$;
Centered: $\sigma_{\text{c}} = \tfrac{1}{2} ( \hat{q}_{i+1} - \hat{q}_{i-1})$.

All would give oscillations. Limit slope: $$ \bar{\sigma} \equiv \bar{\sigma}(\sigma_{\text{u}}, \sigma_{\text{d}}) \overset{\text{(eg)}}{=} \begin{cases} 0 & \text{if } \sigma_{\text{u}} \cdot \sigma_{\text{d}} \le 0 \\ \sigma_{\text{u}} & \text{if } | \sigma_{\text{u}} | < | \sigma_{\text{d}} | \\ \sigma_{\text{d}} & \text{if } | \sigma_{\text{u}} | > | \sigma_{\text{d}} | \end{cases}. $$

The first nonlinear scheme to look at is slope limiting. We're still doing a reconstruct-evolve-average method. So we'll reconstruct our solution within each cell based on the cell averages within this cell and its neighbours. Then we'll compute the intercell flux, evolve, and average to get new cell averages. The new step is that the reconstruction within each cell will be piecewise linear, not piecewise constant.
We can approximate the slope inside each cell by forward, backward or central differencing. Central differencing uses more information and will be the most accurate in smooth regions. However, if we use two cells across a shock then we'll get oscillations. So what we want to do is compare the three different possible slopes. When it looks like we're in a smooth region we choose the slope from central differencing. When it looks like we're not in a smooth region we use piecewise constant reconstruction: the slope is zero. Then we want some smooth transition between the two extremes.
The easiest check for possible shocks or problems is to look for when forward and backward differencing gives slopes with different signs. The simplest approximate slope that works is to look just at the forward and backward (or upwind and downwind) slopes, and choose the minmod slope. That is, if the two slopes have different signs, the slope is set to zero. If they have the same sign, we choose the slope with smallest magnitude.
Minmod slope limiting combined with a higher-order ODE solver will work to give better than first order accuracy. Unfortunately it's not much better until we get to high resolution: the exercises give examples of this. Better slope limiters are available, but even the best won't do better than second order accuracy, a long way from the fourth order or better that we want.

Finite difference methods

In $N$-d, finite volume needs a surface integral: expensive.

Instead, write a finite difference method as $$ \frac{\text{d}}{\text{d} t} q_i = \frac{1}{\Delta x} \left[ f_{i-1/2} - f_{i+1/2} \right]. $$ Now $f_{i \pm 1/2}$ not intercell fluxes. Directly reconstruct flux.

For stability must split flux: $$ f = f^{(+)} + f^{(-)}, \quad f^{(\pm)} = \tfrac{1}{2} \left( f \pm \max | \lambda | q \right). $$

Trying to move to even higher order methods is complicated by an issue that our focus on one dimension has hidden to now: the computation of the intercell flux. Remember, in our finite volume approach, the cell average is updated by the surface integral of the flux over the cell boundary. Once we move to higher order reconstruction methods the solution, and hence the flux, varies over the surface of the cell. This means we can't evaluate it at a single point, as in one dimension: instead we must evaluate it at multiple points and approximate the surface integral.
The most efficent way of approximating an integral in general is Gauss quadrature, where approximating the integrand at $k$ points in each dimension gives us $(2 k)^{\text{th}}$ order accuracy. For second order accuracy we need only one point - this is the midpoint rule. With two points we get fourth order accuracy, and with three points we get sixth order accuracy. However, this is per dimension: to get fourth order accuracy in three dimensions, where each cell face is two dimensional, we will need to reconstruct to 4 points on each face and solve 4 separate Riemann problems per face. The computational cost and complexity goes up very rapidly.
For this reason, high order finite volume methods are rare (although most methods in the field are finite volume methods, they aim for high absolute accuracy with second order convergence). To get higher accuracy, the standard approach in the field is to use finite differences.
In finite difference methods we still write the update formula as if we were computing intercell fluxes. This form, with its telescoping property, is necessary to get shock speeds correct. However, we're now interpreting these $f_{i \pm 1/2}$ terms differently. We don't need these terms to approximate the intercell fluxes, as in finite differences there are no cells, only the values of the solution at the central point. Instead we allow these half-integer terms to be anything, provided the difference $(f_{i+1/2} - f_{i-1/2})/(\Delta x)$ approximates the derivative of $f$ to the appropriate order. As everything is interpreted at the grid point, we automatically get the accuracy associated with the order of approximation of the derivative.
This appears to be wonderful. By directly dealing with the flux and not going through the solution we are bypassing all the issues with cells, intercell fluxes, and Riemann problems. Unfortunately this simple approach also throws away all the characteristic information and so is horribly unstable. By differencing simply without checking which direction the information is travelling in, we introduce Gibbs-like oscillations that rapidly kill the simulation.
So the standard approach is to use flux splitting. We compute the full flux at the grid point, $f(q_i)$. We then split it into two pieces, $f^{(\pm)}$, one of which travels to the left, and one of which travels to the right. This encodes some of the characteristic information. We then use a suitably upwinded reconstruction method on the split fluxes to get their values at the half-integer locations. Finally, we recombine the split fluxes to get the terms needed in the update formula.
The slide shows a simple flux splitting approach which, like the HLLE method, only uses some of the characteristic information. This sort of finite differencing approach is used by codes like Radice's WhiskyTHC to get the current most accurate simulations.

Questions...

One topic left, so
deep breath;
stretch out;
drink something;
screen break;
add questions to the chat.

Ciolfi, 2001.10241.

MHD

Constraint ${\cal C} = \nabla \cdot {\bf B} = 0$ needed.

Either

Modify equations of motion so $$ \partial_t {\cal C} \sim -\alpha {\cal C} \quad \implies \quad {\cal C} \sim e^{-\alpha t}; $$
Find discrete scheme so, when $\nabla \to D$, $$ D \cdot {\bf B} = 0. $$

In this final part of the lecture we need to look at the extra steps needed to evolve MHD, which is the main model for current simulations. The problem with MHD is not the flux terms, which can be evolved with the exact same techniques as we've discussed. The problem is not with the source terms (unless we include resistive effects). The problem is with the "div B" constraint, $\nabla \cdot {\bf B} = 0$.
Constraints are a fact of life in numerical relativity. They're useful as an independent measure of the accuracy of our simulation. However, some constraints, when violated, can destroy the simulation accuracy. It's been known for decades (even before Brackbill and Barnes did their key investigation) that MHD simulations violating the div B constraint would rapidly fail.
There are essentially two types of approach to solving this problem. We can modify the evolution equations so that the constraint violations are kept "small", by damping and propagating these constraints away. Alternatively, we design our numerical scheme so that the constraints are zero at the discrete level.
Constraint damping is the simplest approach. In essence, it works by adding additional terms, proportional to the constraint itself, so that the evolution equation for the constraint has its time derivative evolved by a term proportional to the constraint itself. Choosing the proportionality constant to be negative means the constraint is exponentially damped. To avoid feedback at specific points it's typical to add additional terms so that the constraint violations propagate as a wave equation.
The main advantage of constraint propagation is its simplicity. Adding one or two equations of motion for the constraints and a couple of extra source terms is easy to implement. The disadvantages can be substantial. It requires hand-tuning a parameter, the damping time. Its mathematical status is dubious - there's some results suggesting the modified system isn't properly well posed. The interaction of the violations with grid boundaries, particularly AMR grid boundaries, can be very messy. And finally, it's not clear how small the constraint violations need to be to not cause problems.

Vector Potentials

Use vector potential ${\bf A}$, so ${\bf B} = \nabla \times {\bf A}$.

Automatically preserves constraint.
EM gauge freedom. Prefer Lorenz gauge, $$ \begin{aligned} \partial_t {\bf A} + \nabla \Phi & \sim {\bf E}, \\ \partial_t \Phi + \nabla \cdot {\bf A} &\sim 0. \end{aligned} $$
Grid structure is conceptually more complex.

The preferred method now is to enforce the constraint at the discrete level. The favourite method at the moment is to not evolve the magnetic field directly, but instead to evolve the vector potential ${\bf A}$. We can always compute the magnetic field from the curl of the vector potential, and because the divergence of a curl automatically vanishes, this enforces the constraint.
The vector potential is not unique: by modifying the EM gauge we can change its value. This means different choices for the evolution equation for the vector potential can be made. As with constraint damping, it turns out that propagating the vector potential is better than trying to fix it for simplicity. This means the preferred method is to use a Lorenz (or generalized Lorenz) gauge.
It's useful to note that the vector potential must be continuous, as the magnetic field must be $C^0$, and the magnetic field goes like a derivative of the vector potential. So standard methods can be applied to the vector potential, or even methods that work for $C^1$ functions. However, in order to cleanly get the magnetic fields from the vector potential, and to get the terms needed for evolving the vector potential itself, its typical to discretely locate the vector potential at different locations to the other variables. Specifically, when thinking of the problem in finite volume form, the vector potential is stored at cell edges.

Summary

We have discussed

the three key approaches,
- finite volumes;
- finite differences;
- finite elements,
discrete flux conservation;
reconstruction, monotonicity and Gibbs oscillations;
MHD and constraints.

Next lecture: some aspects of the future.

Alford et al, 1707.09475.

The flux conservative form that we saw in the first lecture is essential for building numerical methods that deal with shocks correctly. Once that's in place we can improve are method by carefully considering how and where we are representing which variables. I think it's really important to learn the finite volume approach as a way of thinking about the key concepts. However, I also think that the medium to long term future of the field is moving towards finite difference or DG type methods, inspired by finite volume concepts.
For MHD and other models including EM fields it's important that the constraints that follow from Maxwells equations are imposed. When implementing a new code I always first reach for constraint damping methods due to their simplicity. However, it's definitely true that the most accurate and useful codes available now are using methods that enforce the constraints at the discrete level. Combining high order accuracy with constraint enforcement is a challenge.