Numerical Hydrodynamics: Part 3

Ian Hawke

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

The present

Ideal GRMHD;
exploring parameter space in
- masses
- mass ratios
- compositions
- magnetic field strengths
- magnetic field configurations
- EOSs and parameters
worrying about
- turbulence and small scale features
- finite temperature
- reactions and transport.

Current

Steiner, GitHub.; Page, UNAM.

Goal

Multiscale theory

Let ${\cal L}_{1,2}$ describe different bits of physics. Then $$ \begin{aligned} \partial_t q &= {\cal L}_1 q + \varepsilon {\cal L}_2 q, \\ \partial_t \tilde{q} &= {\cal L}_1 \tilde{q} \end{aligned} $$ are indistinguishable for times $\sim \varepsilon^{-1}$.

As $\varepsilon$ is a timescale, look for physics/problem $\tau$, or $v / L$, or $D / L^2$, or...

Argument fails for

nonlinear operators;
resonances.

There's been a few questions about how important individual bits of physics are. There are two responses to that.
The first is an argument for multiscale theory. This assumes that there's two bits of input physics each described by a differential operator, ${\cal L}_{1,2}$. We've seen the differential operators relevant for hydrodynamics and MHD and can imagine how MHD can be written as the hydro differential operator plus a separate operator. The assumption is then that the physics acts of different scales, with the second physical model being "small", hence multiplied by $\varepsilon$.
In many cases it can be proved that the full models with both sets of input physics is indistinguishable from a model using just the first set of input physics, provided we're allowed to change parameters of the first model a little, and provided we only look on timescales proportional to $1/\varepsilon$.
Now, we want to work on timescales of around a second. So, if the relative difference between two models gives a parameter whose timescale is less than this then we can probably ignore it. To do this we'd look for appropriate dimensionless constants and physical timescales. However, this needs doing with care. For example, bulk viscosity will depend on the Reynolds number and the size of the NS. But that's an average over all viscosity contributions. We know that viscosity typically acts strongly in boundary layers, which would appear - for example - near the crust-core interface. In this case the length scale is much shorter - maybe centimetres or less - which hugely boosts the impact.
In addition to these problems with applying the multiscale theory, it's known that the theory fails in particular circumstances. What the theory is doing is assuming the second operator acts on short, fast scales compared to the first, and then it's averaging (in some sense) over those scales. This only works as long as there's no resonances - that is, no feedback from the large scale to the small. However, this feedback most definitely occurs in fluid models through turbulence. So we would expect the details of this approach to fail in the post merger regime. At that point nonlinear numerical simulation is the only way to tell when the new physical model is quantitatively needed.

Non ideal MHD

Ideal MHD assumes $\sigma \to \infty$ to enforce $$ {\bf E} = {\bf v} \times {\bf B}. $$

Detailed calculations show $\sigma$ drops when

density drops;
temperature increases;
magnetic field increases.

All possible in (post) merger.

Harutyunyan & Sedrakian, 1605.07612

Our current models rely on the ideal MHD approximation. As we've seen in the previous two lectures, this simplifies the full Maxwell equations by imposing an Ohm's law that relates the electric field to the magnetic field and fluid velocity. This is only valid when there is no resistance to the charge current. This means we're requiring the conductivity to be infinite.
Detailed calculations have been done of the conductivity within neutron stars, with one recent example concentrating on the outer layers of the neutron star shown in the figure by Harutyunyan and Sedrakian. We see that the conductivity is large, typically around $10^{20} \text{s}^{-1}$ or bigger. However, the conductivity drops off sharply at low densities and high temperatures. It also drops off sharply as the magnetic field strength increases.
This means that during a neutron star merger we should expect the resistivity to become important in the core (where the temperature and the magnetic field are high) and in the remnant disk (where the density is low and the temperature high). So for multimessenger simulations we need to check if resistivity has an impact.

Plasma models

A more general Ohm's law: $$ \def\n{{\rm n}} \def\p{{\rm p}} \def\e{{\rm e}} \def\l{\Lambda} \def\c{{\rm c}} \def\h{\Sigma} \def\s{\mathrm{s}} \newcommand{\A}{{\cal A}} {\scriptsize \begin{aligned} e n_\e \mathcal E_b &- \left( 1 - {n_\e \mu_\e \over p+\varepsilon}\right) \epsilon_{bac} J^a b^c - {1\over n_\e e} \left( \hat {\mathcal R} - \Gamma_\e s \A^{\e\s} \right) J_b \\ &= - en_\e \epsilon_{bac} v_\p^a b^c + \mathcal R_{\e\n} w^{\n\p}_b + \left( \mathcal R_{\e\s} - \Gamma_\e s \A^{\e\s} \right) \left( {q_b \over sT} - v^\p_b\right) \\ &-n_\e \mu_\e \left[ \left( v_\p^a - {J^a\over e n_\e} \right) \nabla_a u_b + \perp^c_b u^a \nabla_a \left( v^\p_c - {J_c\over e n_\e} \right) + \left( v^\p_b - {J_b \over e n_\e} \right) u^a {1\over \mu_\e} \nabla_a \mu_\e \right] \\ &+ 2 n_\e u^a \nabla_{[a} s\A^{\e\s} w^{\s\e}_{b]} - e \Gamma_\e \left[ \perp^a_b+ u^a \left( v^\p_b -{J_b\over en_\e}\right) \right] A_a . \end{aligned} } $$

From Andersson et al, 1610.00449. Only considers $e, p, n$ plus heat.

Resistive MHD

Instead impose $j^a = q n^a + J^a$, $$ J^i = q v^i + W \sigma \left[ E^i + \epsilon^{ijk} v_j B_k - (v_k E^k) v^i \right]. $$

Can now solve full Einstein-Euler-Maxwell (11 PDEs), given constituitive relations (EOS, $\sigma \equiv \sigma(\rho, T)$).

However, near ideal MHD write $\eta = \sigma^{-1} \ll 1$, to see $$ \partial_t {\bf E} \sim \frac{1}{\eta} \left[ {\bf E} + \dots \right]: \quad \text{stiff}. $$

Dionysopoulou et al, 1502.02021

In current simulations using resistivity people don't start from the plasma model. Insted the Ohm's law given introduces the conductivity tensor explicitly, performs a $3+1$ split of the charge current little $j^a$, and imposes that the spatial part of the current big $J^i$ is linked to the velocity and the EM fields. All current simulations, starting from the work of Palenzuela, and through to work by Lehner, Dionysopoulou and Alic, restrict to scalar conductivities.
We now have a prescription for the charge current meaning we can solve the full Einstein-Euler-Maxwell equations, which are 11 PDEs: one from continuity, four from generic energy-momentum conservation, and six from Maxwell's equations. In addition we have to provide the constituitive relations for the equation of state (fixing, for example, the pressure) and for the conductivity itself. For now we can think of the conductivity as a constant, but as we've seen above, in general it will be a function of density, temperature, composition and magnetic field strength.
We now run into a problem, not at the analytic level, but at the numerical. Consider how the electric field evolves. We know its time derivative is given by the curl of the magnetic field and is sourced by the charge current. The conductivity in the charge current is large, particularly near the ideal MHD limit. To match up with how numerical analysts talk about this problem we rewrite in terms of the scalar resistivity, $\eta$, which tends to zero in the ideal limit. So, in this limit, the problem is stiff. That is, it's trying to react on a timescale proportional to the resistivity $\eta$, but this is very small. This has consequences at the level of ODEs.

Stiffness and ODE solvers

Solution of $$ \frac{\text{d} q}{\text{d} t} = F(q) = -\alpha q, \qquad q(0) = 1 $$ is $\exp[-\alpha t]$. Can use $$ \begin{aligned} \text{Forward Euler:} && q^{n+1} &= q^n + \Delta t \, F(q^n), \\ \text{Backward Euler:} && q^{n+1} &= q^n + \Delta t \, F(q^{n+1}). \end{aligned} $$

Explicit: fast, unstable if $\Delta t \gtrsim \alpha$;
Implicit: slow (nonlinear algebraic equation at every step), stable.

The problem of stiffness is directly linked to the solution of the ODE. So far we haven't been too explicit about how to solve the ODE - HP has gone into more detail. The essential point is that we want to be sufficiently accurate whilst keeping it computationally cheap. This means we're usually using third or fourth order explicit Runge-Kutta methods.
However, when dealing with stiff problems we have the source term reacting on timescales that are comparable to or smaller than the timestep used by the ODE solver. This happens in resistive MHD, where near the ideal MHD limit the source term is trying to force the electric field to behave as in the traditional Ohm's law, and it's trying to force it to react very quickly.
On the exercise sheet there's an analytical calculation for the linear problem shown here. The exact solution is a rapid exponential decay. The explicit forward Euler scheme only works when the timestep is small enough: that is, when the timestep is roughly smaller than the decay timescale $\alpha$. This is a qualitative feature of all explicit methods, including those used in our standard simulations.
The alternative shown here is the implicit backward Euler scheme. For a general ODE this is called implicit as, in order to compute the time update $q^{n+1}$, we need to solve the implicit nonlinear algebraic equation specified by the function $F$. In our full PDE problem we're going to have to do this at every grid cell, at every update. These implicit schemes are stable, but much more expensive, often by orders of magnitude.
At first glance this makes resistive MHD impractical. However, there are two tricks we can play.

IMEX

Solves stiff system $ \partial_t q + \partial_x f = \varepsilon^{-1} R(q)$ by writing a Runge-Kutta step as $$ \begin{aligned} q^{(i)} &= \color{green}{q^n} - \Delta t \, \color{green}{\sum_{j=1}^{i-1} \left[ \tilde{a}_{ij} \partial_x f(q^{(j)}) - \frac{a_{ij}}{\varepsilon} R(q^{(j)}) \right]} + \\ &\qquad \Delta t \, \color{red}{\sum_{k=i} \frac{a_{ik}}{\varepsilon} R(q^{(k)})}. \end{aligned} $$

Choose coefficients $a_{ik}$ to minimize implicit cost.

Ponce et al, 1404.0692

The first workaround for the computational cost of implicit schemes is to note that for many problems, including ours, the stiff terms only appear in a few equations. In our case we're solving 11 PDEs and the stiff source term only affects three: those linked to the electric field. We could therefore try a mixed scheme, where the bulk of the system is updated explicitly, and only the electric field is treated implicitly.
This is the fundamental idea behind the implicit-explict (or IMEX) Runge-Kutta methods of Pareschi and Russo. They've been used for resistive MHD by Palenzuela and others. This allows for practical simulations with resistivities up to $\sim 10^{11} \text{s}^{-1}$. There's been considerable work done to accelerate the implicit step: see, for example, the work of Cordeiro-Carillon who's used IMEX schemes to deal with issues in spherical coordinates.
In the end, however, the implicit step is still limited by the size of the resistivity. It's still impractical to run full merger simulations with physical values of the resistivity.

Relaxation systems

Look at $$ \begin{aligned} \partial_t q + \partial_x v &= 0 \\ \partial_t v + a \partial_x q &= \varepsilon^{-1} \left[ f(q) - v \right]. \end{aligned} $$

Equilibrium is $v = f(q)$. Remove $v$ from system.
Chapman-Enskog expansion $v = f(q) + \varepsilon v_1$ gives $$ \partial_t q + \partial_x f(q) = \varepsilon \partial_x \left\{ \left[ a - \left( f'(q) \right)^2 \right] \partial_x q \right\}. $$
Smaller system with diffusive correction.

As the purely numerical IMEX approach doesn't yet allow for practical simulations of the full resistive simulations in mergers, we instead look for an analytic trick. This comes from work on relaxation systems, where the stiff source term is meant to force the large system towards an equilibrium. This is the situation we've got with resistive MHD.
The problem we're looking at here is a toy problem. The variable $q$ is the one we really care about. It evolves slowly, but is coupled to another variable $v$. $v$ has a stiff source term that tries to force $v$ to its equilibrium solution, and because it is stiff this causes $v$ to evolve quickly. At equilibrium $v$ is given by a function of the slow variable $q$: out of equilibrium the source tries to rapidly push it back.
In the limit where $\varepsilon \to 0$ the equilibrium constraint must be enforced. At that point we can drop the evolution equation for the fast variable $v$ and only consider the slow variable $q$. Away from this limiting case we can expand our fast variable in the small timescale $\varepsilon$. This is a Chapman-Enskog expansion. The leading order term gives us the equilibrium solution. The next term can be determined explicitly by the full evolution equations.
Plugging in this linear term we get an equation of motion for the slow variable that decouples from the fast variable. It looks like the equilibrium conservation law with a correction term. This correction term is small - it's multiplied by the small timescale $\varepsilon$. Importantly this correction term is a second derivative: it's a diffusive correction to the original system.
This toy problem leads to a scalar equation. We can extend this to systems to deal with resistive MHD. In resistive MHD the fast variables are the electric fields. The equilibrium solution will be ideal MHD. The short timescale is given by the resistivity. The diffusive correction will capture the long-lengthscale corrections coming from EM resistivity.

REGIME

Chapman-Enskog expansion on resistive MHD: $$ {\small \newcommand{\pdv}[2]{\frac{\partial {#1}}{\partial {#2}}} \newcommand{\bm}[1]{{\bf {#1}}} \begin{aligned} \partial_t (\sqrt{\gamma} {\bf q}) &+ \partial_i (\sqrt{-g} {\bf f}^i_0 + {\bf F}^i ) = \sqrt{-g} {\bf s}_0 + {\bf S} + \partial_i {\bf D}^i, \\ {\bf D}^i &= - \pdv{\bm{f}^i_0}{\overline{\bm{q}}} \bigg( \pdv{\overline{\bm{s}}_0}{\overline{\bm{q}}} \bigg)^{-1} \bigg[ \partial_j(\sqrt{-g} \overline{\bm{f}}^j_0) + \dots \bigg]. \end{aligned} } $$

Corrections to fluxes and sources;
Not stiff as $\varepsilon \to 0$;
Works for NS merger regime;
No free lunch: new timestep restriction!

Wright & Hawke, 1906.03150

Alex Wright has done the hard work of extending the Chapman-Enskog approach to relativistic resistive MHD. The result is conceptually similar to the toy problem, but practically much more complex. The full expressions are found using computer algebra, but their general form is as here: there are corrections to the fluxes, the geometric sources, and there's the expected diffusive correction. The leading order expansion also gives us the standard ideal MHD Ohm's law.
With this approach we can now use standard explicit solvers. We're back to only solving 8 PDEs, whilst still getting the resistive corrections to leading order. It's more expensive that ideal MHD, but only by constant factors, not by orders of magnitude, because it's not stiff.
However, there's no free lunch, so we should expect some problems. In particular, numerical algorithms for solving diffusive terms typically have stricter timestep restrictions for stability, with the timestep scaling proportional to $(\Delta x)^2$. That also happens with the correction term here. That would significantly restrict our simulations at very fine resolutions. However, for current NS merger simulations we're not hit by that stability limit.

Questions...

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

Radice et al, 1809.11161

Discontinuous Galerkin

Revisit the weak form. Multiply by $\phi$, integrate: $$ \int_V \phi \partial_t {\bf q} + \oint_{\partial V} \phi {\bf f} - \int_V {\bf f} \cdot \nabla \phi = \int_V \phi {\bf s}. $$

Moves all spatial derivatives from ${\bf q}$ ($C^0$) to $\phi$ (piecewise $C^\infty$).
Information in $q$ is lost by reconstruction.
Surface integral only coupling to neighbours. Finite volumes/differences couple more.

There are two key computational problems with the higher order methods we've described to now.
The first is waste. The reconstruction used in both finite difference and finite volume methods takes a limited amount of information about the solution - the cell averages, or the point values - and produces a form of the solution everywhere, usually as a piecewise polynomial. This form is then used to compute some terms, such as the intercell fluxes, and then all the high order information is thrown away. We only store the values of the solution or its cell average.
The second problem is communication. To get enough information for a high-order approximation to the solution we need to look not just at the cell and its immediate neighbours, but to ever further neighbouring cells. In general, we need to look at $k$ neighbours on either side to get $(2 k - 1)^{\text{th}}$ order accuracy. This is a real issue for big simulations on modern supercomputers.
As a quick digression. Simulations have essentially four things that can slow them down. The first is how fast it does operations: the FLOP count. The second is how much memory it has. The third is how fast results can be saved to disk. As the fourth, when we split the calculation across multiple processes, is how fast the different bits of the simulation communicate with each other.
On modern and near-future machines the main problem is communication speed. High order finite volume and finite difference schemes communicate too much information with cells too far away from themselves. This stops our simulations using the computational power available.
One answer to both these issues is to use a Discontinuous Galerkin method. HP has covered this in his talk on vacuum numerics, so here I'll just touch on the features important to hydrodynamics.
First we have to revisit the weak form. In DG and other finite element methods we take the equations of motion, multiply by a test function $\phi$, and integrate by parts. This moves the spatial derivatives from the solution (which might be discontinuous) to the test function (which we can choose to be sufficiently differentiable).

Function basis expansion

Expand both ${\bf q}, \phi$ as eg $q(x, t) = \hat{q}_m(t) P_m(x)$. The weak form gives $$ \int_V P_m P_n \hat{\phi}_m \partial_t \hat{q}_n + \oint_{\partial V} P_m P_n \hat{\phi}_m \hat{f}_n - \int_V P_m \hat{f}_m \hat{\phi}_n \nabla P_n = \int_V P_m P_n \hat{\phi}_m \hat{s}_n. $$

Simplifies to $$ M \partial_t \hat{{\bf q}} + S^T f(\hat{{\bf q}}) = - [\phi {\bf F}]_{x_{i-1/2}}^{x_{i+1/2}}. $$

More accurate.
Less communication - exascale!
Issues with shocks.

We then expand the solution and the test function in terms of some function basis: think of Fourier Series, or Legendre polynomials. We're going to store the modes: the coefficients of the solution with respect to the function basis expansion. We end up with evolution equations for the modes, but these are coupled through the mass matrix and stiffness vector. These matrices and vectors can be pre-computed. We see there is a term that looks like the standard boundary flux, for which we'll need something like a Riemann solver.
DG methods don't have the communication or waste problems seen in finite volume or finite difference methods. They only couple a cell to its neighbours through the intercell flux. The high order reconstruction is automatic thanks to the mode information stored in each cell, which is evolved, not discarded between steps.

Limiting...

Simple DG methods have Gibbs oscillations;
$h$-refinement localises, but doesn't remove, them;
the DG moments can be limited to stop oscillations;
limiting reduces DG advantages.
Better DG implementations needed for matter.

HP showed the wonderful convergence properties of DG methods. They don't converge as rapidly as spectral methods, but they are more flexible.
However, I personally am sceptical about the performance of DG methods when it comes to discontinuous solutions. Here I'm showing the problems that arise when evolving discontinuous problems. I'm evolving the "top hat" profile using the advection equation.
As expected the DG method gives Gibbs oscillations. HP noted that you can use $h$-refinement to resolve the shock: essentially dropping to a low order scheme and using small cells around the discontinuity. I am sceptical that this can be done dynamically in a nonlinear problem where the shocks develop without forcing the cell size, and hence the timestep, to be impractically small.
An alternative is to extend the techniques in slope limiting to limit the moments of the solution: that is, directly limiting the mode coefficients. This is sufficient to eliminate oscillations, but increases the communication by coupling more cells, in the same way that slope limiting does. However, it's still better in terms of communication that say a high order WENO scheme.
However, when comparing the efficiency and convergence rate of DG methods with tradiational high order methods, like the WENO based finite difference methods, we should compare to the limited DG method (or some other DG method that is avoiding oscillations). The second figure does this by comparing computational cost to error when advecting a smooth solution, but still using limiting.
We see that, when comparing like with like, the efficiency of the two methods is comparable. At this point I wouldn't be willing to commit to say which method is better. There are two points I would make.
The high order finite difference methods have been shown to work with NS mergers. They're mature technology.
The DG methods are new and so may improve substantially, once robust ways to avoid Gibbs oscillations are known. Right now they're not obviously better, but this may change quickly.

The future

We have discussed

relaxation systems and methods;
Discontinuous Galerkin methods;

Further possibilities for the future include

Path conservative schemes;
Well balancing;
Adaptive model refinement;
Uncertainty quantification.

Kiuchi et al, 1509.09205