Numerical Hydrodynamics: Part 1

Ian Hawke

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

GW170817

The only NS merger seen in GWs:

high SNR;
covers frequencies from 10s-1000s Hz;
observed in multiple detectors, allowing localisation.

For many of us studying numerical relativity and relativistic hydrodynamics the motivation is gravitational waves. The goal was to observe these and use them to tell us about gravity, relativity, and the physics of extreme objects.
With the detection of black hole mergers in 2015, and neutron star mergers from 2017, the goal is now reality. The object now is to squeeze the maximum possible information from these observations, through better modelling. Full, quantitative modelling of neutron star mergers needs numerical simulation, which is going to be the focus of these lectures.
Our first goal in numerical modelling is therefore to predict the gravitational wave signal as seen here: a "chirp" lasting seconds with frequencies in the hundreds to thousand of Hz.

GW170817

NS constraints from GWs:

Mass $\sim 1.4 M_{\odot}$;
Radius $\sim 11 \pm 1 \text{ km}$;
weak constraints on EOS through tidal compressibility.

LIGO, 1710.05832

Summit supercomputer. Photo: Carlos Jones/ORNL

Sun: $\sim 10^{60}$ elementary particles.
Computer resource on Earth: $\sim 10^{25}$ FLOPs.

Averaging

If we can't treat each particle separately, we'll have to average over them in some way. The hydrodynamics approximation averages over them all. Our $10^{60}$ separate particles, each with their own positions, velocities, and properties, are averaged to give a single (vector) field. At each point in spacetime this field will describe the number density of particles, and the direction in which they flow.
This averaging process gives us a small set of differential equations to solve. With some standard approximations, it actually reduces to five PDEs (on top of those describing the evolution of the spacetime). This makes the numerical simulation practical. However, it has also thrown away (or compressed) a lot of physics.
It also means that there is a crucial difference of principle and of scale between vacuum relativity and relativistic hydrodynamics. In the vacuum case we have Einstein's field equations which, in principle, hold on all scales within our simulations and models. In hydrodynamics, however, our model and description breaks down on certain scales. We must keep in mind that as computer power and numerical accuracy increases our results may not get better, as we may hit the limitations of our models.

Modelling: inspiral

For inspiral only bulk properties matter:

mass;
spin;
compactness (EOS);
bulk internal magentic field.

Leads to key numerical limits:

$\lambda_{\text{GW}} \sim 10^6 \text{ m} \implies L_{\text{grid}} \sim 10^7 \text{ m}$;
Phase accuracy $\implies \Delta x \sim 100 \text{ m}$.
Stability: $\Delta x \sim 100 \text{ m} \implies \Delta t \sim 10^{-6} \text{ s}$.

Now we have a starting point for our simulation. We're going to solve a set of PDEs within a domain: a box that we'll put around the neutron stars. Within this box we'll evolve the spacetime and the matter, and extract observable signatures like gravitational waves and emissions in the electromagnetic spectrum. We'll do this by splitting the box into small pieces - grid cells, or points, or elements - and replacing the PDE with discrete approximations that can be solved on a computer. The smaller the small pieces, the more accurate the solution, but the more computer time and memory we need.
So, what's our minimum requirements?
First, let's look at the gravitational waves detected from GW170817. These were detected in the hundreds to thousands of Hz. As the waves travel at the speed of light we can convert the frequencies to a wavelength, suggesting that the gravitational wavelength here is around $10^6$m. So, our computational box needs to be big enough to fit in a gravitational wavelength. In fact, it should be big enough to fit in multiple wavelengths, so we can cleanly extract the wave far from both the source - that is, far from the neutron stars - and also far from the artificial boundaries of our box. So our minimum requirement is that the box should be around $10^7$m to a side.
Next, let's think about the size of the grid cells. We're going to use a hydrodynamic model of the star. In a fluid model like this, information travels back and forth in waves. The numerical accuracy depends both on how accurately we capture the amplitude of the waves, and how accurately we capture their phase. A back of the envelope calculation tells us how many grid cells we need to use to capture a wave of wavelength the size of a neutron star (whose diameter is around $25$km) oscillating say $30$ times (which captures roughly the last $10$ orbits) to around $1$% accuracy, giving a grid cell size of around $10^2$m.
Finally, let's think about how many steps we have to take. When evolving forward in time we can think about the simulation in relativistic terms. Take a grid cell at a fixed time and think about its past domain of dependence: that is, which grid cells (at some previous time) could possibly affect the value of this cell? This suggests that the maximum timestep we can take is roughly the grid cell size divided by the speed of light, which is roughly $10^{-6}$s.
Putting this all together, we need a grid of roughly $10^{15}$ cells to solve a merger with reasonable accuracy. To solve for around $1$s - and remember that the GRB after GW170817 was detected a few seconds after the peak gravitational wave amplitude - we need to take $10^6$ timesteps. Combining this, we're looking at updating a small cell at least $10^{20}$ times. That's a lot of work.
Now, no current simulation evolves for a full second. However, simulations do evolve for up to around $100$ms, and do typically use resolutions around $100$m, and sometimes considerably finer. We need to consider what physics we can capture within a model on these scales, what behaviour we should expect, and how we can simulate that numerically.

Multi messenger

Want to complement GW observations with

electromagnetic signals (radio, X-ray, $\gamma$, ...);
neutrino observations;
anything else we can get.

Done with GW170817: inferred jet properties suggest more constraints on progenitors.

LIGO, Fermi, INTEGRAL et al, 1710.05833

Let's look at a simulation using a single perfect fluid with magnetic field. This simulation is from Bruno Giacomazzo's group. It ignores dissipation, viscosity, heat transport, magnetic resistivity, superfluidity, superconductivity, the elastic crust of the neutron star, and any radiation transport. It includes full general relativity, models the star as a hot fluid, and adds an ideal electromagnetic field.
Initially the stars inspiral due to the emission of gravitational waves. Around $15$ms after the start of the simulation the stars merge. At this point things get messy. When the stars smash into each other a shock propagates through the merged object, massively increasing the temperature. At the same time the magnetized material "slips" past its companion, winding up the total magnetic field in a dynamo effect. The remnant tries to settle down, with some material expelled from the system. But the object is too massive to support itself, and at around $30$ms the core collapses to a black hole. Most of the dense matter is rapidly swallowed into the black hole, leaving a small fraction at larger radius to settle into a disk. It's believed that the kilonova occurs when this remaining matter collapses back on the black hole.
We see that we need to model the fluid interior and the vacuum exterior. We need to model the shock as the stars smash together, and the resulting jump in temperature. We need to model the small scale structure in the interior as the magnetic fields wind up, forming the precursor to the jet. And we need to do this on top of evolving the spacetime, tracking any black holes that might appear.

Questions...

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

Shibata & Uryū 1999

Conservation

EFEs imply $\nabla_a T^{ab} = 0$.

Pick a tetrad, $e_b^{(j)}$ to get $$ \begin{aligned} && \nabla_a \left[ e_b^{(j)} T^{ab} \right] &= \tfrac{1}{\sqrt{-g}} \partial_a \left( \sqrt{-g} e_b^{(j)} T^{ab} \right) = -T^{ab} \nabla_a e_b^{(j)} \\ \implies && \partial_t {\bf q} + \partial_i {\bf f}^{(i)}({\bf q}) &= {\bf s}. \end{aligned} $$ Balance law form.

Only four equations: need other constituitive equations for, eg, EM, particle number, etc.

We now know the outlines of our problem. We want a model of matter, coupled to GR, that we can write as a set of PDEs. This model will encode the detailed physics we're interested in, at the level we can simulate with our given resources. But before we try and be specific, let's look at what we can say in general.
We know the field equations require total stress energy conservation: the divergence of $T^{ab}$ is zero. We can use this to get four equations of motion. The most useful form for our purposes is to use the standard identity for the 4-divergence of a vector, which can be written as the partial 4-divergence of that vector, weighted by the metric determinant. To get from the stress-energy, which is a two-tensor, to a vector, we contract the free index with one member of a tetrad, which is a set of orthonormal (in our case) vectors that we can link to the coordinates. This then allows us to write stress energy conservation as four PDEs, one for each member of the tetrad.
The key point is that the only derivatives of the matter appear in total divergence form. The term on the right-hand-side which isn't in this form only involves derivatives of the tetrad, which can be written as derivatives of spacetime quantities. We write this abstractly in balance law form: all derivatives of the matter are total derivatives, but these fluxes must balance off against the source terms which come from the geometry.
If we simplify to Minkowski space in standard Cartesian coordinates then the source terms vanish and these balance laws reduce to the conservation laws for energy-momentum. However, we need to be careful of this interpretation (locally): if we write Minkowski spacetime in spherical coordinates then there will be source terms, but that certainly does not mean that energy-momentum is being transferred between the matter and the spacetime.

Shock formation

Advection equation $$ \partial_t q + \partial_x (v q) = 0. $$ Information moves right, speed $v$.

Burgers equation $$ \partial_t q + \tfrac{1}{2} \partial_x q^{2} = 0. $$ Information moves right, speed $q$. Shocks form.

We see that generically the matter will obey conservation laws. We expect the fluxes ${\bf f}$ to depend nonlinearly on the conserved variables ${\bf q}$. This immediately tells us something very important.
The simplest conservation law is the advection equation, where the solution is propagated to the right with constant speed $v$. For a nonlinear conservation law we can locally approximate it as an advection-like equation, where the solution propagates with speed $\partial_q f$. This means the propagation speed depends on the data. Depending on the form of the flux function, some gradients will steadily get steeper, and some will get shallower. In the absence of dissipation (and crucially the generic total energy-momentum equations have no dissipative terms) these steepening effects will lead to discontinuities: a generic shock will form.

Shocks and uniqueness

Shock speed $V_s$ from Rankine-Hugoniot, $V_s \left[ q \right] = \left[ f \right]$. Burgers equation $$ \partial_t q + q \partial_x q = 0 $$ is equivalent to any of $$ \partial_t q^n + \tfrac{n}{n+1} \partial_x q^{n+1} = 0, $$ but they all have different shock speeds.

Total derivative form crucial: otherwise need more information to fix solution.

From the conservation law form we can integrate over a small region of spacetime (see the exercises) to find the Rankine-Hugoniot conditions, telling us how fast the shock travels. It's important to note that this depends on the form of the flux function.
However, this highlights how important the form of the equations of motion is. Take the example of Burgers equation where the flux function is quadratic. This is equivalent, when solutions are differentiable, to a PDE where the local advection velocity is given by the data $q$ itself. However, this isn't unique. The conservation laws at the bottom are all equivalent to this non-conservative form of Burgers equation for positive $n$. Crucially, they all have different shock speeds, so different solutions in the discontinuous case.
For our purposes this means that when the shock forms as the neutron stars merge, and as it propagates around in the post-merger state, it will be essential that we have the correct form of the conserved variables and the fluxes. If we use the wrong conservative form then we will get the wrong solution. If we use the non-conservative form then we need to explicitly include a dissipative term and take appropriate limits, which is hard. If the model can't be phrased in a balance law form, getting correct shock forms will be complex.

Euler equations in relativity

Need: $\rho, u^a, \epsilon$.
Particles conserved: $$\nabla_a ( \rho u^a ) = 0.$$
Stress-energy: $$T_{ab} = \rho (1 + \epsilon) u_a u_b + p \perp_{ab}.$$
EOS: implies $p = p(\rho, \epsilon, \dots)$, and more.
Mixture models: $Y_{\text{x}}$ from $$\nabla_a ( \rho Y_{\text{x}} u^a ) = 0.$$

Total stress-energy conservation gives us four equations. When we talk about a standard ideal fluid it depends on more information. We'll need the number density: the number of particles within a local volume. We'll need its velocity, which is three bits of information making up a unit 4-vector. And we'll need a variable describing the internal energy, or temperature: essentially, how much variation there is at the micro-scale that we've averaged over to get our hydrodynamics description.
This is five independent quantities. So, in addition to the four PDEs we get from stress-energy conservation, we'll need one more. We get this by assuming that particles are neither created nor destroyed, so the particle number flux is conserved. This gives us $\nabla_a (\rho u^a) = 0$. Using the same identity for total divergence as in the stress-energy case, we can use this to write the continuity equation as a PDE in conservation law form.
We can then write down the stress-energy tensor for the fluid. There is the total energy piece which moves with the flow. Then there is the isotropic pressure piece which acts normal to the flow. Plugging this into the conservation laws and adding the continuity equation gives us five PDEs for six quantities. To close the system we need to connect the thermodynamic variables: the density, internal energy, and pressure. We do this through the equation of state, which we can take to be a function giving pressure (and other quantities) in terms of density and internal energy, or density and temperature.
Before we move on we should talk about more complex mixture models. The equation of state describes the microscopic interactions between the particles. In the model we've talked about so far we haven't made any allowance for the differences between particles: we've only cared about how many there are, through the number density. This is too simplistic in many cases. Many equations of state really care about the relative proportions of the different particles involved. For example, many care about the lepton fraction $Y_e$, which is (roughly) the fraction of particles that are electrons.
The simplest mixture models are formed by imposing that each individual type of particle is also conserved, and that they all move with the same velocity. We then get the equivalent of the continuity equation for each species, which can be seen as an equation for the species fraction, say $Y_e$. The equation of state can then be extended to depend on these species fractions. In more complex mixture models we can allow particle species creation and destruction by adding source terms, provided the total number of particles stays fixed. This can require some care, as the particle interactions can impact on the energy-momentum.

Slow motion limit

Newtonian: $$ \partial_t \begin{pmatrix} \rho \\ \rho v \\ E \end{pmatrix} + \partial_x \begin{pmatrix} \rho v \\ \rho v^2 + p \\ (E + p) v \end{pmatrix} = {\bf 0}. $$

Propagation speeds $$ v, v \pm c_s: \quad c_s^2 \sim \frac{p}{\rho}. $$

It's useful to compare the resulting equations of motion to various limits. This is particularly true when developing numerical methods, algorithms and codes, and it allows us to test in simpler situations where it's easier to debug.
The particular limit we're most interested in is the Newtonian limit, where the speed of light is very large compared to the fluid velocities, and pressures (for example) are small. The Newtonian limit is particularly useful to us as there's a long and successful history of building simulation codes in academia and industry, and we can (and still do) use the experience and methods of this huge community.
The Newtonian equations in $1+1$ Cartesian coordinates have the following conservation law form. The continuity equation is first, essentially advecting the density along with the flow. The momentum equation is next. The momentum is advected, but there's also a contribution from the gradient of the pressure. The total energy equation is last, which again shows advective-like behaviour. If we were to investigate the speeds of propagation, which is one of the exercises, then we wouldn't be surprised that one is precisely the advective velocity $v$. The other speeds of propagation are the acoustic velocities which are $v \pm c_s$, where $c_s$ is the speed of sound. The sound speed depends only on the thermodynamic quantities and is given by the equation of state. In simple cases it's roughly the square root of the ratio of the pressure to the density.

Back up to relativity

Newtonian	$\partial_t \begin{pmatrix} \rho \\ \rho v \end{pmatrix} + \partial_x \begin{pmatrix} \rho v \\ \rho v^2 + p \end{pmatrix} = {\bf 0}.$
SR	$\partial_t \begin{pmatrix} \rho \color{red}{W} \\ \rho \color{red}{h W^2} v \end{pmatrix} + \partial_x \begin{pmatrix} \rho \color{red}{W} v \\ \rho \color{red}{h W^2} v^2 + p \end{pmatrix} = {\bf 0}.$
GR	$ \begin{aligned} \partial_t \color{gold}{\sqrt{\gamma}} \begin{pmatrix} \rho W \\ \rho h W^2 v_j \end{pmatrix} &+ \partial_i \color{gold}{\sqrt{-g}} \begin{pmatrix} \rho W \left( v^i - \color{gold}{\frac{\beta^i}{\alpha}} \right) \\ \rho h W^2 v_j \left( v^i - \color{gold}{\frac{\beta^i}{\alpha}} \right) + p \delta^i_j \end{pmatrix} = \color{gold}{{\bf s}}, \\ {\bf s} &= \sqrt{-g} \begin{pmatrix} 0 \\ T^{\mu\nu} \left( \partial_\mu g_{\nu j} - \Gamma^\mu_{\nu j} \right) \end{pmatrix}. \end{aligned} $

To see how useful a model the Newtonian model will be to us, first look at the equations in special relativity. That is, take the GR equations in Minkowski spacetime in standard Cartesian coordinates. We see that they're directly comparable to the Newtonian limit.
We can see the additional relativistic effects in red. In the continuity equation we see the appearance of the Lorentz factor, which is linked to the velocity, and appears here because of length contraction or time dilation effects. The Lorentz factor is one when velocities are small but diverges as the velocity approaches the speed of light. This is a purely kinematic effect, exactly mirroring standard SR results on length contraction.
In the momentum equation we also see the Lorentz factor coming in. However, we also see the density term is modified by the appearance of the specific enthalpy $h$. This comes in because of mass-energy equivalence in relativity: the internal energy of the fluid has its own contribution to the mass.
The kinematic modifications here leave the qualitative form of the conservation laws the same, but modify the quantitative details. Possibly the most important for numerical work is the speeds of propagation. As you would expect in relativity the speeds have to be below the speed of light. One speed remains the advective velocity $v$. The acoustic velocities now use the relativistic addition formula, dividing the Newtonian formula by $1 \pm v c_s$. The appearance of the enthalpy in contributions to the mass also modifies the speed of sound.
Finally we look at the equations in full general relativity. These equations look much more complex, but should contain few surprises based on our previous discussions. First, from the form of the general energy-momentum conservation laws, we expect to get terms relating to the determinant of the metric. The difference between the spatial metric determinant which appears in the time derivative, and the total metric determinant appearing in the spatial divergence, is a factor of the lapse. This can be thought of as making the dimensions match up. The geometric source terms we have already discussed; they're very messy but, along with the metric determinants, result from the spacetime curvature. The other difference to note is the way the velocity is modified by the shift vector. This results from the $3+1$ decomposition, so the fluid velocity now has a piece related to the way coordinates move relative to our spatial slice.
The conceptual issues from the transition from SR to GR are (at least in terms of the matter) small. The practical implementation issues are, however, substantial. The speeds of propagation become much more complex to write down. The relation between conserved and primitive variables becomes steadily more complex as we move from Newtonian to relativity. More care is required in dealing with the variables: see that both covariant and contravariant components of the velocity appear in the fluxes, for example. However, many of the key conceptual points can be illustrated and tested in a Newtonian setting and Newtonian codes.

Electromagnetism

Charge current $j^a \sim e (n^a_{\text{e}} - n^a_{\text{p}})$. Source for EM fields described by Faraday tensor $F_{ab} = \nabla_{[a} A_{b]}$, so $$ \begin{aligned} \nabla_a F^{ab} &= 4 \pi j^a \\ \nabla_a {}^* F^{ab} &= 0. \end{aligned} $$ Gives Maxwell's equations: use $E^a, B^a \implies$ balance law form.

$$ T_{ab} = T_{ab}^{\text{hydro}} + T_{ab}^{\text{EM}}. $$

Current diverges as resistivity $\to 0$: MHD limit $\implies {\bf E} = {\bf v} \times {\bf B}$.
Constraint $\nabla_a B^a = 0$ must hold.

Before we move on to talking about how we simulate these PDEs, I want to note how the electromagnetic fields can be included in our models.
When we talked about mixture models we mentioned how we could distinguish different types of particles moving together. Going beyond just the fluid properties, we can also consider the charge of the particles. At any point, the total charge current is the sum over all species of their number current times their charge. This charge current acts as a source for the Maxwell equations. The Maxwell equations result from total charge conservation.
All of this can be expressed in terms of the Faraday tensor. This gives a term which adds on to the stress-energy tensor, modifying the conserved variables and fluxes in our generic conservation equations. In addition, the Faraday tensor can be written in terms of electric and magnetic fields, and in these terms the Maxwell equations can be written in balance law form with the current as a source.
The problem is that a lot of dynamics is acting on very short length and time scales. This makes the evolution of the full Einstein-Euler-Maxwell system very expensive, as we need - for stability - to take tiny timesteps. It is typical to make additional assumptions to restrict the impact of the small scale behaviour, and to enforce ideal conditions where there is no resistance to the charge current. In this case a simplified Ohm's law allows us to replace the electric field with the cross product of the magnetic and velocity fields, allowing us to drop the evolution equation for the electric field completely.
We now have eight PDEs: the continuity equation for total particle number, the generic energy-momentum balance laws, and the three equations for the magnetic field. We have an equation of state closing the system, and constraints on the magnetic field following from Maxwell's equations. The evolution equations are practically more complex than in the hydrodynamic case, but conceptually similar. The constraint equation, however, causes conceptually new problems, as we'll see later.

Questions...

Deep breath;
have a drink;
stretch out;
add questions to the chat.

Kiuchi et al, 1509.09205

Outside the NSs

Now we have equations of motion that we can, in principle, solve. Before doing so we need to check we having all the conditions required to solve a set of PDEs. The crucial things we need to think about are

do we know the domain on which we're going to solve;
do we know the boundary conditions we're going to impose;
do we know the relationships between all the variables involved in the problem?

We discussed the domain when talking about the practical requirements at the start: a box big enough to hold multiple gravitational wavelengths. That means it will be many orders of magnitude bigger than our neutron star. That makes the boundary conditions we're going to impose deceptively easy: at the edge of the box there's no matter, and never will be any matter, so we can "just" set it to be vacuum.
This is where we should realize that the vast majority of the box will be vacuum, maybe with some propagating EM waves. This is problematic, not for the boundary conditions themselves, but for the equations of motion. All our fluid conserved variables contain the density, which vanishes outside of the star. This means the equations degenerate in this vacuum region. This makes the surface of the neutron star in particular tricky to handle, and the boundary between a region that's mathematically well described by our equations, and one that's not.
A range of analytical and numerical tricks, modifications and fixes have been suggested. The simplest, and still widely used, trick, is to introduce an atmosphere. This purely numerical trick replaces vacuum with a low density region that is left to evolve freely. Any point or region where the density, pressure or energy gets too low is reset to this atmosphere value.
Whilst simple and (reasonably) robust, the physical implications of adding an atmosphere remain worrying. Experiments have shown that the impact of accreting this atmosphere onto the neutron star are negligible, and its interaction with the outer boundaries is also not a problem. However, its interaction with low density matter expelled in the merger process, particularly with magnetic fields, or in regions of high temperature, could have a significant impact on the observables. This is somewhere we must improve.

C2P

Newtonian: $v = S / \rho$, etc.

Relativity: $\{ \rho_0, v^j, \epsilon \} \leftarrow f(D, S_j, \tau)$.

Implicit and nonlinear: guess $\bar{p}$, solve $\bar{p} - p(\bar{\rho}_0, \bar{\epsilon}) = 0$;
Expensive and not always robust.

GRMHD: strong coupling through ${\bf E} = {\bf v} \times {\bf B}$.

Multi-dimensional nonlinear algebraic equation;
Very expensive and sensitive.

The final topic to explore before we work out how to solve our PDEs is to relate the different variables that appear in our equations of motion. We have a number of different classes of variables. The primitive variables are any minimal set of variables that could, in principle, be observed. Typically in hydrodynamics we would think of density, velocity, and (say) internal energy. The conserved variables are those actually evolved. Here that's something like the density, momentum, and energy. The auxilliary variables are other quantities that are useful to compute for analysis or simplicity, such as the enthalpy or Lorentz factor. The fluxes and sources finish the description of the equations of motion, and are usually a combination of primitive, conserved, and auxilliary variables.
In the Newtonian case there are closed-form analytic expressions relating all variables, except those that need the equation of state. There is minimal cost in converting between any sets of variables, but minimising the number of times the equation of state is called is useful for computational performance.
In the relativistic cases things change drastically. There are closed-form analytic expressions to go from the primitive to the auxilliary and conserved variables. However, it is very rare to have closed-form expressions that compute the primitive variables from the conserved variables. So, given initial values for all variables we can cheaply update the conserved variables using the equations of motion. It is then a complex, expensive process to recover the primitive and auxilliary variables that we need in order to compute the fluxes and sources to update to the next step.
The typical approach is to write the relation between the conserved and primitive variables as a nonlinear implicit algebraic equation. Standard algorithms for solving these are known: when it can be written in terms of a single variable, as in the hydrodynamic case, then algorithms like bisection and Newton-Raphson can be used.
Interestingly, when more physics is added in a minimally coupled fashion, the conversion does not get more complicated. You can see this from the stress-energy tensor. Minimal coupling essentially means the new physics is added in, meaning the new terms are added in to the conserved variables. Separating these terms that are linearly combined is possible in, for example, Einstein-Euler-Maxwell.
However, this is not true in MHD. In using Ohms law to eliminate the electric field we strongly couple the Euler and Maxwell parts of the system, through "E equals v cross B". The MHD system becomes a multi-dimensional root-find, which is much more expensive and much less robust. There are available codes based on lengthy papers for doing this inversion: I strongly recommend relying on them.

Summary

We have discussed

what's practical;
the hydrodynamic description;
balance laws;
the remaining EOM;
implementation issues.

Next lecture: how we actually evolve.