Module equations 

The equations

Merlin solves systems of conservation laws

\partial_t \mathbf{u} + \nabla\!\cdot\!\mathbf{F}(\mathbf{u}) = 0,
\qquad \mathbf{u}(\mathbf{x},0)=\mathbf{u}_0(\mathbf{x}),

where $\mathbf{u}$ holds the conserved quantities and $\mathbf{F}$ is the flux. Every flux equation in this book (all but the level set) is of this form, possibly with the geometric terms of curvilinear coordinates. They are hyperbolic: the Jacobian $\partial\mathbf{F}\cdot\mathbf{n}/\partial\mathbf{u}$ has real eigenvalues (the wave speeds) and a full set of eigenvectors for every direction $\mathbf{n}$ (LeVeque 2002, Chapter 2). In the code an equation is a type that supplies the physical flux and a bound on the wave speed: [merlin_systems::Equation].

Solutions with shocks

Solutions of nonlinear conservation laws develop discontinuities even from smooth data, so the differential equation is understood in its integral (weak) form: for every cell $[x_{i-1/2},x_{i+1/2}]$,

\frac{d}{dt}\int_{x_{i-1/2}}^{x_{i+1/2}} u\,dx = F(u_{i-1/2}) - F(u_{i+1/2}).

Two consequences matter for validation.

Rankine-Hugoniot. A discontinuity separating $\mathbf{u}_L$ from $\mathbf{u}_R$ moves at a speed $s$ fixed by conservation,

s\,(\mathbf{u}_R-\mathbf{u}_L) = \mathbf{F}(\mathbf{u}_R)-\mathbf{F}(\mathbf{u}_L).

A scheme that is not written in conservation form puts shocks in the wrong place; the shock-position tests measure this directly.

Entropy. Weak solutions are not unique. Only the one that is the limit of vanishing viscosity is physical; it admits compressive shocks and rejects expansion shocks (LeVeque 2002, Section 11.14). A scheme can satisfy Rankine-Hugoniot and still converge to the wrong weak solution, which is what the deliberately unfixed Roe flux does at a transonic rarefaction.

Linear advection

u_t + \mathbf{a}\cdot\nabla u = 0

with a constant velocity $\mathbf{a}$. The exact solution is the initial data translated, $u(\mathbf{x},t)=u_0(\mathbf{x}-\mathbf{a}t)$. On a periodic domain and after a whole number of periods the exact answer is the initial condition, which makes the error trivial to measure and is how the schemes’ orders are established in 1D.

Burgers’ equation

u_t + \left(\tfrac12 u^2\right)_x = 0 .

The simplest nonlinear law, with wave speed $f'(u)=u$, so that larger values catch up with smaller ones. It has exact solutions for exactly the cases that test a shock capturing scheme:

• Smooth data before the shock forms. Along the characteristic leaving $\xi$ the solution is constant, so $u(x,t)=u_0(\xi)$ with $x=\xi+t\,u_0(\xi)$. For $u_0=\sin 2\pi x$ the characteristics first cross at $t_b=-1/\min u_0' = 1/(2\pi)\approx0.159$. After that the shock position follows from Rankine-Hugoniot, and for this symmetric profile it sits at $x=1/2$ with values $\pm u_s$, where $\xi + t\sin 2\pi\xi = 1/2$ and $u_s=\sin 2\pi\xi$.
• A shock, $u_L>u_R$: a jump moving at $s=(u_L+u_R)/2$.
• A rarefaction fan, $u_L<u_R$: $u=(x-x_0)/t$ between $u_L$ and $u_R$. If $u_L<0<u_R$ the fan is transonic: it contains the sonic point $u=0$ where the characteristic speed changes sign.

The exact solvers are in src/pymerlin/exact_solutions/burgers.py (riemann, sine) and in crates/merlin-finite-volume/tests/common/burgers_exact.rs.

Euler equations of an ideal gas

For density $\rho$, velocity $\mathbf{v}$, pressure $p$ and total energy density $E$, in $D$ dimensions,

\mathbf{u}=\begin{pmatrix}\rho\\ \rho\mathbf{v}\\ E\end{pmatrix},\qquad
\mathbf{F}_d=\begin{pmatrix}\rho v_d\\ \rho v_d\mathbf{v}+p\,\mathbf{e}_d\\ (E+p)v_d\end{pmatrix},\qquad
E=\frac{p}{\gamma-1}+\tfrac12\rho|\mathbf{v}|^2 .

The sound speed is $c=\sqrt{\gamma p/\rho}$ and the waves along direction $d$ travel at $v_d-c$, $v_d$ (the contact, carrying the density jump and the transverse velocity) and $v_d+c$. The gas law is that of an ideal gas with one $\gamma$ per material (set in the problem file; the shock-tube problems use 1.4). The implementation is [merlin_systems::ideal_gas::IdealGas].

Other fluids close the same equations with another law $p=p(\rho,e)$ (the fluid model of a material in the problem file; [merlin_materials::FluidModel]). The stiffened gas $p=(\gamma-1)\rho e-\gamma p_\infty$ stands for a liquid (water: $\gamma=4.4$, $p_\infty=6\times10^8$ Pa); its sound speed is $c=\sqrt{\gamma(p+p_\infty)/\rho}$, and with the shifted pressure $p+p_\infty$ and the energy $E-p_\infty$ its equations are those of the ideal gas, which gives the exact solution used in the water tube. The ideal gas is $p_\infty=0$. The Mie-Gruneisen law describes a metal under shock, a pressure linear in the energy about a reference curve, $p=p_{\text{ref}}(\rho)+\Gamma(\rho)\rho\,(e-e_{\text{ref}}(\rho))$ (Menikoff and Plohr 1989, Section VI), here with the shock Hugoniot of the material at rest as the reference curve, built from the linear shock-speed relation $U_s=c_0+s\,u_p$, and $\rho\Gamma$ constant; with a constant reference pressure and constant $\Gamma$ it is the stiffened gas. Its waves have no closed form: the shock follows from a root find on the Hugoniot relation and the rarefaction from an integral along the isentrope (merlin_materials::fluid::numerical). The copper impact runs it.

The flux of direction $d$ is the flux of direction 0 with the velocity components permuted. This rotational invariance is why the code is written once for one direction and applied along each axis (see Numerical methods).

The exact Riemann solution

A Riemann problem is the initial value problem for two constant states with one jump. Its solution is self-similar, a function of $x/t$ alone: three waves leave the jump and separate four constant regions (Toro 2009, Chapter 4). The outer two waves are nonlinear, each a shock or a rarefaction fan; the middle one is a contact, across which pressure and velocity are continuous and only the density (and any transverse velocity) jumps. In the figure the left wave is a rarefaction, bounded by its head and tail, and the right wave a shock, as in Sod’s tube.

The wave fan of the Euler equations in the $(x,t)$ plane. The regions $\mathbf{u}_*^L$ and $\mathbf{u}_*^R$ share the pressure $p_*$ and velocity $u_*$.

For the Euler equations the solution is closed-form up to one scalar root. The pressure $p_*$ in the region between the nonlinear waves solves

f_L(p_*)+f_R(p_*)+u_R-u_L=0,\qquad
f_K(p)=\begin{cases}
(p-p_K)\sqrt{\dfrac{A_K}{p+B_K}} & p>p_K\ \text{(shock)}\\[2ex]
\dfrac{2c_K}{\gamma-1}\left[\left(\dfrac{p}{p_K}\right)^{\frac{\gamma-1}{2\gamma}}-1\right] & p\le p_K\ \text{(rarefaction)}
\end{cases}

with $A_K=\dfrac{2}{(\gamma+1)\rho_K}$ and $B_K=\dfrac{\gamma-1}{\gamma+1}p_K$, and the velocity there is $u_*=\tfrac12(u_L+u_R)+\tfrac12\bigl(f_R(p_*)-f_L(p_*)\bigr)$. The solution sampled at $x/t$ is the reference for every shock tube in this book. It is implemented twice, independently: in Rust ([merlin_finite_volume::riemann::ExactRiemann]) and in Python (src/pymerlin/exact_solutions/euler.py), and the two are checked against Toro’s published star values (Tables 4.1 and 4.2).

Exact solutions used for smooth and radial flows

• The isentropic vortex (Shu 1998): a vortex of strength $\beta$ carried by a uniform flow is an exact solution of the 2D equations for all time. It is smooth, so it measures the order of accuracy in two dimensions.
• Planar shocks, whose speed and post-shock state follow from Rankine-Hugoniot.
• Linear sound waves, Noh’s implosion, the Sedov-Taylor blast wave and solid-body rotation, in the curvilinear section.

Curvilinear coordinates

With the radius $r$ as the first coordinate, the Euler equations in cylindrical ($\alpha=1$) or spherical ($\alpha=2$) symmetry are

\begin{aligned}
\partial_t\rho + r^{-\alpha}\partial_r\!\left(r^{\alpha}\rho u\right) &= 0,\\
\partial_t(\rho u) + r^{-\alpha}\partial_r\!\left(r^{\alpha}(\rho u^2+p)\right) &= \frac{\alpha\,p}{r},\\
\partial_t E + r^{-\alpha}\partial_r\!\left(r^{\alpha}(E+p)u\right) &= 0 .
\end{aligned}

In the axisymmetric $(r,z)$ form the $z$ direction is added with the Cartesian flux. The only term that is not a divergence is the pressure source $\alpha p/r$ in the radial momentum, and the finite-volume form handles it by removing the pressure part of the flux before weighting with face areas (see Numerical methods).

Swirl. A cylindrical flow can also turn about the axis with the azimuthal velocity $w$. The equations gain a variable and two terms,

\begin{aligned}
\partial_t(\rho u)+r^{-1}\partial_r\!\left(r(\rho u^2+p)\right)+\partial_z(\rho uv) &= \frac{p+\rho w^2}{r},\\
\partial_t(\rho w)+r^{-1}\partial_r\!\left(r\rho uw\right)+\partial_z(\rho vw) &= -\frac{\rho uw}{r},
\end{aligned}

the centrifugal force $\rho w^2/r$ joining the pressure in the radial momentum and the Coriolis term $-\rho uw/r$ acting on the swirl. Mass, axial momentum and energy are unchanged. The second equation is the conservation of the angular momentum,

\partial_t(r\rho w)+r^{-1}\partial_r\!\left(r\cdot r\rho uw\right)+\partial_z(\rho vw\,r)=0,

in which the Coriolis term has disappeared, and this is the form the scheme discretises (see the methods page). Across the axis both $u$ and $w$ change sign, because the directions $\hat r$ and $\hat\theta$ both turn round. Solid-body rotation $w=\Omega r$ of a gas of constant density is an exact steady state when the pressure balances the centrifugal force, $\partial_rp=\rho w^2/r$, that is $p=p_0+\tfrac12\rho\Omega^2r^2$.

Sound waves. For a small perturbation of a gas at rest the equations reduce to $p'_{tt}=c^2\bigl(p'_{rr}+\tfrac{n}{r}p'_r\bigr)$, with $n=1$ for the cylinder and $n=2$ for the sphere. The sphere has the outgoing solution $p' = \rho_0 c\, f'(r-ct)/r$ from the potential $\varphi=f(r-ct)/r$ (Landau and Lifshitz, Fluid Mechanics, Section 69). The cylinder has no such single-profile solution, so the reference is a Hankel transform of the initial pulse, $p'(r,t)=\int_0^\infty k\hat P(k)J_0(kr)\cos(kct)\,dk$, evaluated numerically in Python (src/pymerlin/exact_solutions/acoustic.py).

Noh’s implosion (Noh 1987). A cold gas flows towards the origin at unit speed. For $\nu=\alpha+1$ (2 for the cylinder, 3 for the sphere), a shock leaves the origin at speed $s=(\gamma-1)/2$. Ahead of it $\rho=\rho_0(1+t/r)^{\nu-1}$, $u=-1$, $p=0$; behind it the gas is at rest,

\rho=\rho_0\left(\frac{\gamma+1}{\gamma-1}\right)^{\nu},\qquad u=0,\qquad p=\tfrac12(\gamma-1)\rho .

For $\gamma=5/3$ these are the densities 16 (cylinder) and 64 (sphere) and the shock speed $1/3$ seen in the cylindrical and spherical pages.

The Sedov-Taylor blast wave

Sedov (1959; Landau and Lifshitz, Fluid Mechanics, Section 106): an energy $E$ released at a point (a line, a plane) into a cold gas of density $\rho_0$ drives a strong shock at $R(t)=\bigl(Et^2/(\alpha\rho_0)\bigr)^{1/(\nu+2)}$, where $\nu=1,2,3$ for planar, cylindrical and spherical symmetry. Behind it the flow depends on $\xi=r/R(t)$ only. With $\lambda=2/(\nu+2)$, write

u=\lambda\frac rt\,V(\xi),\qquad c^2=\lambda^2\Bigl(\frac rt\Bigr)^2Z(\xi),\qquad \rho=\rho_0\,G(\xi),

and with $D=\xi\,d/d\xi$ the equations of continuity, entropy and momentum become the linear system

\begin{aligned}
(V-1)\,D\ln G+DV+\nu V &= 0,\\
D\ln Z+(1-\gamma)\,D\ln G &= \frac{2(1-\lambda V)}{\lambda(V-1)},\\
-V+\lambda V^2+\lambda(V-1)\,DV+\frac{\lambda}{\gamma}Z\,(2+D\ln G+D\ln Z) &= 0,
\end{aligned}

for $DV$, $D\ln Z$ and $D\ln G$. They are integrated from the shock $\xi=1$, where the strong-shock jump conditions give $V=2/(\gamma+1)$, $Z=2\gamma(\gamma-1)/(\gamma+1)^2$ and $G=(\gamma+1)/(\gamma-1)$, to the centre. The constant $\alpha$ follows from the total energy,

\alpha=\Omega_\nu\lambda^2\int_0^1G\Bigl(\tfrac12V^2+\frac{Z}{\gamma(\gamma-1)}\Bigr)\xi^{\nu+1}\,d\xi,

with $\Omega_\nu=2,2\pi,4\pi$ the angle of the geometry. For $\gamma=1.4$ it is 1.0774 (planar), 0.98408 (cylinder) and 0.85107 (sphere, Sedov’s 0.851). It is a limit: it describes the flow once the shock is far outside the region that held the energy and the gas ahead of it is cold, which the simulations approach from a ball of finite radius (the cylindrical and spherical pages).

The level set

An interface is the zero level set of a scalar function $\phi(\mathbf{x},t)$, negative on one side and positive on the other (Osher and Sethian 1988; Osher and Fedkiw 2003). Moving it with a velocity field $\mathbf{u}$ means advecting $\phi$:

\phi_t+\mathbf{u}\cdot\nabla\phi=0 .

This is not a conservation law in $\phi$: it is a Hamilton-Jacobi equation with the Hamiltonian $H=\mathbf{u}\cdot\mathbf{p}$, linear in the gradient $\mathbf{p}$, and nothing is conserved (see the note on area). Information travels along $\mathbf{u}$, so the derivative is taken on the upwind side. The velocity is prescribed in this book (uniform, rotation, a reversing vortex, or an expression in $\mathbf{x}$ and $t$); coupling it to the gas is planned. Solutions stay smooth as long as $\phi$ does, so the schemes are judged by an order of convergence and by returning a shape to where it started.

Distance to the interface of an ellipse. Characteristics of $\lvert\nabla\phi\rvert=1$ run along the normals; where those from different parts of the interface meet, the nearest point jumps and $\phi$ has a kink. That kink is why the pseudo-time reinitialisation uses Weno5 whatever rule advects $\phi$.

Signed distance. Advection stretches $\phi$, so one may want to rebuild a function with the same zero level set and $\lvert\nabla\phi\rvert=1$, the signed distance $d(\mathbf{x})$ to the interface. It solves the Eikonal equation $\lvert\nabla\phi\rvert=1$, reached either as the steady state of the pseudo-time equation $\phi_\tau+S(\phi_0)(\lvert\nabla\phi\rvert-1)=0$ with $S$ a smoothed sign (Sussman, Smereka and Osher 1994), or directly by sweeping (Zhao 2005). The exact signed distance of a circle, $\lvert\mathbf{x}-\mathbf{c}\rvert-R$, is the reference (high-order page).

What an answer must satisfy

These are the properties the tests check beyond a comparison with an exact solution:

Property	Statement
Conservation	On a periodic or closed domain, $\int\mathbf{u}\,dV$ is constant to round-off; with open ends it changes by exactly the boundary flux
Positivity	$\rho>0$ and $p>0$ throughout, including the near-vacuum and $10^5$ pressure-ratio problems
Rankine-Hugoniot	Shocks sit where conservation places them
Entropy	Expansion shocks do not appear
No new extrema	Limited schemes create no new extrema in scalar problems
Symmetry	A problem with a symmetry keeps it, in some cases to the last bit
Equilibrium	A gas at rest stays at rest, exactly, in any geometry

References

• LeVeque, R. J. (2002), Finite Volume Methods for Hyperbolic Problems, Cambridge.
• Toro, E. F. (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer.
• Shu, C.-W. (1998), Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, Lecture Notes in Mathematics 1697, Springer (the isentropic vortex).
• Noh, W. F. (1987), Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux, J. Comput. Phys. 72, 78-120.
• Landau, L. D. and Lifshitz, E. M., Fluid Mechanics, Sections 69 and 106.
• Osher, S. and Sethian, J. A. (1988), Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79.
• Osher, S. and Fedkiw, R. (2003), Level Set Methods and Dynamic Implicit Surfaces, Springer.
• Sussman, M., Smereka, P. and Osher, S. (1994), A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114.
• Zhao, H. (2005), A fast sweeping method for Eikonal equations, Math. Comp. 74.
• Sedov, L. I. (1959), Similarity and Dimensional Methods in Mechanics, Academic Press, Chapter IV.