14. Numerical Methods

Every astrophysical simulation solves a problem that is not the one written down, and this module says which one. It derives the Courant condition that Module 10 used, the viscosity that every grid code adds to the equations it solves, and the exact Riemann solver on which Godunov's method is built. It then runs three one-dimensional codes against closed forms that this book already owns. Four checks: a stability boundary reproduced, a convergence rate predicted with no fitted parameter, a blast wave in the right place with its peak too low, and a particle method whose error at a contact does not decrease with resolution.

Contents · ← Module 13

Module 10 §3 ended with a number: at a Reynolds number of $10^{6}$, one outer turnover of a directly simulated flow costs $10^{18}$ point-updates, and the $10$ pc molecular cloud of Module 10 §6 needs $6.24\times10^{17}$ grid points. Its conclusion was that "every astrophysical code therefore solves a different problem from the one written down here", and it named this module to say which problem that is. Module 10's count also used a result it did not derive: that an explicit code is limited to a time step of about one cell crossing, the Courant condition. Module 5 §8 named this module for how a gravitational collapse is followed once linear theory stops.

The originating plan of the book added one more debt: to say which results of Modules 1 to 13 came from a simulation and which from a closed-form solution. §11 pays it with a census of every module, and the census finds eight places where a shipped page calls the output of a simulation or a numerical model "measured" or "observed".

1. What a code solves

The equations of Module 2 are conservation laws. In one dimension, with $U = (\rho, \rho u, E)$ the densities of mass, momentum and total energy and $F(U) = (\rho u, \rho u^{2} + P, u(E + P))$ their fluxes, they read $\partial_t U + \partial_x F(U) = 0$, where $E = P/(\gamma-1) + \rho u^{2}/2$ for an ideal gas. A code replaces the continuous $U(x,t)$ by a finite list of numbers. The method used for every grid code in this module keeps the conservation exact.

Part (b) is why this module uses conservative schemes only. A scheme that is not conservative can converge to a discontinuity that moves at the wrong speed, because the jump conditions of Module 8 are the integral form of conservation and nothing else. The run checks part (a) on the spherical blast of §6: over $80{,}935$ steps the total energy decreases by $3.9\times10^{-16}$ of itself, which is rounding. That code updates $U_iV_i$, with $V_i$ the cell volume, through faces of area $r^{2}$, and the same telescoping conserves its total mass and energy; its momentum equation has a geometric source, and radial momentum is not conserved in spherical symmetry.

2. The Courant condition

The simplest conservation law is the advection of a scalar at a constant speed $a \gt 0$: $\partial_t u + a\,\partial_x u = 0$, with exact solution $u(x,t) = u_0(x - at)$. Two schemes for it carry the whole of this section. The upwind scheme takes the flux from the cell the flow comes from, $F_{i+1/2} = a u_i$, so (1.1) becomes $$ u_i^{n+1} = u_i^n - C\,(u_i^n - u_{i-1}^n), \qquad C = \frac{a\Delta t}{\Delta x}. \tag{2.1}$$ The forward-time centred-space scheme (FTCS) takes the average of the two neighbours, $F_{i+1/2} = a(u_i + u_{i+1})/2$, and gives $u_i^{n+1} = u_i^n - \tfrac{C}{2}(u_{i+1}^n - u_{i-1}^n)$. Here $C$ is written bare for the Courant number of the linear problem; it is the $C_{\rm CFL}$ of Definition 1 with $S_{\max} = a$.

Courant, Friedrichs and Lewy gave this argument in 1928. In the English translation of 1967 they write the time step as $h$ and the space step as $kh$, and conclude that the difference solution "in general cannot converge" when $k \lt 1$ (p. 228) and "If now $k \ge 1$ as we must require" (p. 230). For a unit wave speed $k \ge 1$ is $\Delta t/\Delta x \le 1$; the conversion to (2.2) is this book's. The argument proves necessity only, and the next proposition shows why that matters.

The trivial limit. At $C = 1$ the upwind scheme is $u_i^{n+1} = u_{i-1}^n$: it moves the data exactly one cell per step, which is the exact solution, and (2.3) gives $|G| = 1$ for every mode. Both schemes use only the three points of Proposition 2, so both satisfy it for $C \le 1$, and FTCS is still unstable. The domain-of-dependence condition is necessary; stability is a separate property of each scheme.

What this does to Module 10's count. Module 10 Proposition 4 took the step of a simulation as $\delta t \sim \eta/U$: one smallest cell, crossed at the flow speed. For an explicit code that solves the compressible equations, the signal speed in Definition 1 is $|u| + c_{\rm s}$, and the step is $\Delta x/(U + c_{\rm s})$. That takes $1 + 1/M$ times as many steps, with $M = U/c_{\rm s}$ the Mach number. For Module 10's laboratory air, the thermal Mach number $U/\bar v = 0.0218$ (module10.html:148) and $c_{\rm s}/\bar v = (\pi\gamma/8)^{1/2} = 0.7415$ at $\gamma = 1.4$ give $M = 0.0294$ and a factor of $35.01$. For the supersonic cloud of Module 10 §6, $U = 2.64\ \mathrm{km\,s^{-1}}$ against $c_{\rm s} = 0.2439\ \mathrm{km\,s^{-1}}$ gives $M = 10.825$ and a factor of $1.0924$. Module 10's $\delta t \sim \eta/U$ is the step of a code that has removed sound waves, which is what an incompressible code does. It is right to within $9.24$ per cent for the cloud, and it undercounts a compressible code's work on the wind tunnel by a factor of $35.01$. Module 10 is not changed. Its $\mathrm{Re}^{3}$ counts the steps of a code without sound waves, so for a compressible explicit code it is a lower bound on the work, short by the factor $1 + 1/M$. Module 10 calls only its snapshot column a floor (module10.html:283); this result shows that its work column is one too.

3. The viscosity every grid carries

A stable scheme is not an exact one. The question that module10.html:285 left open is what equation the upwind scheme actually solves, and the answer is found by expanding (2.1) in a Taylor series about the grid point.

Two limits. At $C = 1$, $D_{\rm num} = 0$, which agrees with the exact shift found in §2. At $C = 0$ the time step is zero and nothing moves. $D_{\rm num}$ has its largest value, $a\Delta x/2$, but it acts for no time: the variance growth per step, $2D_{\rm num}\Delta t$, is zero. (3.1) is an expansion; the next proposition shows that, for the variance of a profile, it is exact.

The run advects a Gaussian 200 steps on 400 cells at $C = 0.25$, $0.5$ and $0.9$ and finds the ratio of measured to predicted variance growth to be $1.000000000$ each time. A code that runs at a flow speed $U$ on $N$ cells across a length $L$ therefore has a numerical Reynolds number $$ \mathrm{Re}_{\rm num} \;=\; \frac{UL}{D_{\rm num}} \;=\; \frac{2N}{1 - C}. \tag{3.3}$$ At $N = 1024$ this is $4096$ at $C = 0.5$ and $2.048\times10^{4}$ at $C = 0.9$. That is the answer to module10.html:285 for the upwind scheme: the code solves Euler's equations plus a viscosity that the grid sets, and the Reynolds number it runs at is a property of $N$ and $C$ and not of the gas. Better schemes than upwind reduce $D_{\rm num}$ in smooth flow, and §10 gives the bound that holds for any scheme.

4. The Riemann problem and Godunov's method

Upwinding a scalar is easy because the scalar has one speed and its sign says where the information comes from. The Euler equations have three signal speeds, $u - c_{\rm s}$, $u$ and $u + c_{\rm s}$, and they can point in different directions. Godunov's answer, in 1959, was to solve the problem exactly at every cell face.

The trivial limit. If the two states are equal, $f_{\rm L}(P) + f_{\rm R}(P) = 0$ at $P = P_{\rm L}$, so $P^{*} = P_{\rm L}$ and $u^{*} = u_{\rm L}$: nothing happens, as it must. If (4.2) fails, the two rarefactions pull the gas apart faster than sound can fill the gap, and a vacuum opens between them; the solver of this module stops with an error rather than return a pressure. The run solves (4.1) by Newton's method from a linearised first guess and samples the solution at any $x/t$ by the standard construction of each wave; the whole solver is riemann_star and riemann_sample in afd/figs/m14_solvers.py.

5. The anchor: Sod's shock tube

In 1977 Gary Sod compared the difference schemes of his day on one problem: a tube of gas at rest with $\rho = 1$, $P = 1$ on the left of a diaphragm and $\rho = 0.125$, $P = 0.1$ on the right, $\gamma = 1.4$ (his report COO-3077-145, p. 17, on a grid of $\Delta x = 0.01$ at Courant number $0.9$). The problem has been the first test of every code since, because its exact solution contains all three waves of Definition 2: a rarefaction moving left, a contact, and a shock moving right. This book evaluates it at $t = 0.2$ with the diaphragm at $x = 0.5$ on $[0, 1]$.

What the exact solution makes possible is a measurement of the error of a code, region by region. The run applies Godunov's method at $C_{\rm CFL} = 0.9$ on $N = 100$ to $6400$ cells, measures the $L_1$ error of the density, $\sum_i|\rho_i - \rho_{\rm exact}(x_i)|\Delta x$, and splits it at the midpoints between the waves ($x = 0.585718$ and $0.767961$) into the fan, the contact and the shock.

The contact is predicted in advance. Inside the star region $P$ and $u$ are uniform. The Riemann problem between two cells with the same $P$ and $u$ and different $\rho$ has as its solution the contact alone, moving at $u^{*}$, so Godunov's flux of mass there is exactly $\rho_i u^{*}$: the upwind flux. The density across the contact is therefore advected by (2.1) at the Courant number $C_{\rm c} = u^{*}\Delta t/\Delta x$, which the run records step by step (its mean is $0.3810$ at $N = 6400$). Because $u^{*} \gt 0$, the upwind cell is the left one. A step has no variance, so Proposition 5 is applied to its differences. Let $d_i = \rho_{i-1} - \rho_i$. Because (2.1) is linear and its coefficients do not depend on $i$, the differences obey (2.1) too, with the Courant number $C_{\rm c}$ of each step. Initially $d_i$ is a single spike of height $\Delta\rho_{\rm c}$ at the contact. Each step replaces $d_i$ by the weighted average $(1 - C_{\rm c})d_i + C_{\rm c}d_{i-1}$, so for $0 \le C_{\rm c} \le 1$ the spike stays non-negative and Proposition 5 applies at every step: the total stays $\Delta\rho_{\rm c}$, and the variance grows to $\sigma^{2} = \sum_n C_{\rm c}(1 - C_{\rm c})\Delta x^{2}$. After many steps the sum of these independent two-point moves has a Gaussian profile, by the central limit theorem, whose conditions hold because each move is at most one cell and the variances sum without limit. The density is $\rho^{*}_{\rm L}$ minus the running sum of the $d_i$, so it is the step smoothed by that Gaussian. A step of height $\Delta\rho_{\rm c}$ smoothed by a Gaussian of width $\sigma$ has an $L_1$ error $2\Delta\rho_{\rm c}\int_0^\infty[1 - \Phi(y/\sigma)]\,\mathrm{d}y$, with $\Phi$ the standard normal distribution function, which is $$ E_{\rm c} \;=\; \Delta\rho_{\rm c}\,\sigma\left(\frac{2}{\pi}\right)^{1/2}, \qquad \sigma \propto \Delta x^{1/2}. \tag{5.1}$$ With $\Delta\rho_{\rm c} = 0.160746$ this has no free parameter. Because $\Delta t \propto \Delta x$, the number of steps is proportional to $1/\Delta x$, so $\sigma^{2} \propto \Delta x$ and the contact error falls as $\Delta x^{1/2}$. A shock does not smear this way: its characteristics converge on it, the width stays a fixed number of cells, and its error falls as $\Delta x$.

The fan, with order $0.752$, has no prediction here. Its solution is continuous, but its head and tail are corners where the derivative jumps, and a first-order scheme smears a corner too; this module does not derive the rate at which it does so. The ratio of $1.0659$ rather than $1$ is error that (5.1) does not model. The contact window also contains both star plateaus, and this module does not identify the source of the excess. The excess falls with $N$ faster than $E_{\rm c}$ does: it is $1.1293 - 1$ of $E_{\rm c}$ at $N = 100$ and $1.0659 - 1$ at $N = 6400$, and the ratio falls at every doubling.

6. A closed form tests a code: the Sedov blast

Module 8 §4 solved the blast wave of a point explosion exactly: the shock is at $R = \xi_0(Et^{2}/\rho_0)^{1/5}$ with $\xi_0 = 1.032777$ for $\gamma = 7/5$. Kamm and Timmes (2007) turned that into a standard problem for codes: $\rho_0 = 1$, $\gamma = 1.4$ and $E = 0.851072$ in spherical geometry put the shock at $r = 1.0$ at $t = 1.0$. The run imports Module 8's integrator unchanged (the routines are copied, not retyped), recomputes $\xi_0 = 1.032777$ and so an exact radius of $1.000000$ at $t = 1$.

The code is Godunov's method in spherical geometry: cells of volume $(r_+^{3} - r_-^{3})/3$ per steradian, face areas $r^{2}$, and in the momentum equation the geometric source $P(r_+^{2} - r_-^{2})/V$, which is the cell integral of $2Pr$ and makes a uniform pressure an exact steady state. The energy is deposited as heat in the innermost cell, the ambient pressure is $10^{-6}$, and the grid covers $0 \le r \le 1.2$. The run uses $n = 100$, $200$ and $400$ cells at $C_{\rm CFL} = 0.5$. The largest run takes $80{,}935$ steps, because the gas at the centre becomes hot and rarefied and its sound speed sets the step.

The two halves of CHECK 3 point opposite ways for one reason. In the exact solution almost all the swept-up mass sits in a thin shell behind the shock (Module 8 §4). A scheme that smears by Proposition 4 spreads that shell over several cells and lowers its peak, while the position of the shock is fixed by conservation of energy and mass, which the scheme keeps exactly by Proposition 1. A code judged only on its peak density at $n = 400$ would appear to violate Module 8's jump condition by $26.65$ per cent. It does not: the peak is the quantity a first-order code converges on most slowly.

7. Smoothed particle hydrodynamics

A grid code fixes the cells and moves the fluid through them. Smoothed particle hydrodynamics (SPH) does the opposite: the fluid is a set of particles of fixed mass that move with the flow, and every field is an average over neighbours. It needs no grid, it follows collapsing gas to high density without refinement (its spatial resolution follows the density; its mass resolution does not, §9), and it conserves momentum exactly.

Here bare $h$ is the smoothing length, a module-local meaning: Module 7's $h_B$, $h_S$ and Module 8's mixing width $h$ do not occur in this module. With $h = 1.2\,\Delta x$ the kernel reaches $2.4$ particle spacings on each side.

What is stated. Shocks need dissipation that (7.2) does not have. Monaghan's artificial viscosity (1992, p. 550, printed without an equation number) is $\Pi_{ab} = (-\alpha\bar c_{ab}\mu_{ab} + \beta\mu_{ab}^{2})/\bar\rho_{ab}$ for approaching particles and zero otherwise, with $\mu_{ab} = \bar h_{ab}v_{ab}x_{ab}/(x_{ab}^{2} + 0.01\bar h_{ab}^{2})$ (his eq. 4.2), bars denoting pair averages, and "$\alpha = 1$ and $\beta = 2$" (p. 551), the values used here. The thermal energy per unit mass evolves by $\mathrm{d}u_a/\mathrm{d}t = \sum_b m_b(P_a/\rho_a^{2} + \Pi_{ab}/2)\,v_{ab}\,\partial_{x_a}\overline{W}_{ab}$, which with (7.2) conserves the total energy; Problem D3 shows it. The integrator is kick–drift–kick with a step of $0.2\,h/(c_{\rm s} + |v|)$.

8. Grid against particles at a contact

The same Sod problem, run with SPH, separates the two methods at one place. The particles have equal mass, so the right-hand gas, eight times less dense, has particles eight times further apart. The run places 200, 400 or 800 particles per unit length on the left, freezes the six outermost particles at each end of $[-0.5, 1.5]$ as walls, and measures the pressure near the contact at $t = 0.2$, where the exact pressure is $P^{*} = 0.303130$ everywhere between the fan and the shock.

Price describes this as "a large 'blip' in the pressure profile around the contact discontinuity" (arXiv preprint, p. 15). That run takes the density from the continuity equation, and he attributes the blip to "the unresolved gradient in density and a slight overshoot in the thermal energy at the same location". With the density computed by the sum of Definition 4, as here, he finds that "a spurious blip in the pressure is still apparent due to the overshoot in thermal energy" (p. 16). This run shows that overshoot: at 800 particles per unit length the thermal energy within $0.05$ of the contact reaches $3.0943$ against the exact $2.8535$, and $P = (\gamma-1)\rho u$ carries it into the pressure. His own tube uses $\gamma = 5/3$ with 450 particles; this module uses Sod's $\gamma = 1.4$, so the numbers are not his, only the finding. A grid code has the opposite behaviour: Godunov's method smears the contact in density by Proposition 5, but it never produces a pressure spike there, because both $\rho$ and $E$ are smeared by the same flux.

Why this matters beyond one dimension. A contact between gases of different entropy is exactly the interface of Module 7 §6. Price argues that the same pressure blip acts on a shearing interface as an artificial surface tension, and surface tension is what Module 7's dispersion relation shows stabilising short wavelengths. Agertz et al. (2007) ran the Kelvin–Helmholtz instability with density contrasts of 8 and 10, with a grid code (Enzo) and an SPH code (Gasoline), and found that in SPH "Perturbations are damped out very quickly" (arXiv preprint, p. 13), while the grid code grew the rolls. They attribute the failure to density gradients; Price (arXiv preprint, p. 28) attributes it to entropy gradients, "not density gradients as has often been assumed", and the interface named above follows Price. This module does not reproduce the two-dimensional test; it reproduces its one-dimensional cause.

9. Following a collapse

Module 5 §3 derived the Jeans length $\lambda_{\rm J} = (\pi c_{\rm s}^{2}/G\rho_0)^{1/2}$ and showed that perturbations longer than it grow. A code has its own shortest length, the cell. If the cell is comparable with $\lambda_{\rm J}$, rounding errors on the grid scale can grow under gravity where the gas should be stable, and the code fragments a cloud that should not fragment.

For Module 5's $10$ K core at $n(\mathrm{H_2}) = 10^{4}\ \mathrm{cm^{-3}}$, $\lambda_{\rm J} = 0.19479$ pc $= 40{,}178$ AU (the run converts Module 5's rounded $0.19479$ pc; module05.html:170 prints $40177$ AU from the unrounded value) and $M_{\rm J} = 2.658\,M_\odot$, so the cell must be at most $1.004\times10^{4}$ AU. At $10^{6}$, $10^{8}$ and $10^{10}\ \mathrm{cm^{-3}}$ the limit falls to $1004$, $100.4$ and $10.04$ AU. A uniform grid that covers four initial Jeans lengths, $0.77916$ pc, at the resolution needed at $10^{10}\ \mathrm{cm^{-3}}$ needs $1.6\times10^{4}$ cells per side and $4.096\times10^{12}$ in three dimensions. That is why collapse is followed with adaptive mesh refinement, which halves the cell only where $J$ would exceed the limit: a factor $10^{3}$ in $\lambda_{\rm J}$ is $9.966$ levels of factor-2 refinement. The whole of this paragraph assumes the gas stays at $10$ K; it stops being true when the gas becomes opaque to its own cooling radiation, which this book does not compute.

The SPH counterpart is stated, not derived. Bate and Burkert (1997) find that SPH resolves gravitational fragmentation only when the Jeans mass exceeds the mass of about twice the number of neighbours, $M_{\rm res} \approx 2N_{\rm neigh}m$ (their eq. 6, p. 1068), with $N_{\rm neigh} = 50$. At $10^{10}\ \mathrm{cm^{-3}}$ the core's $M_{\rm J}$ is $0.002658\,M_\odot$, so each particle may carry at most $2.66\times10^{-5}\,M_\odot$. The condition does not take care of itself as the gas collapses. The particles crowd together, so the smoothing length shrinks and the spatial resolution follows the density without refinement. The particle mass $m$ does not change, and $M_{\rm J} \propto \rho^{-1/2}$ falls. A calculation therefore resolves fragmentation only up to the density at which $M_{\rm J}$ has fallen to $2N_{\rm neigh}m$, which Bate and Burkert give as a maximum resolvable density (their eq. 7, p. 1068). The particle mass must be chosen at the start for the highest density the run is to follow; the limit above is the one for this core followed to $10^{10}\ \mathrm{cm^{-3}}$. A grid meets the Truelove condition by adding cells where they are needed; an SPH run of fixed particle mass meets the Bate and Burkert condition only if its particle mass was chosen small enough in advance.

10. The price of resolution

Module 10 §3 priced resolving the dissipation scale. This module adds two results. The first is a bound on what any scheme can represent, and it returns Module 10's own count. The second is the cost of a contact, and it makes the count worse.

A bound for any scheme. No scheme can represent a structure narrower than one cell, so a code with $N$ cells across the outer scale $L$ has a dissipation scale of at least $\Delta x = L/N$. By Module 10's (3.4), $L/\eta = \mathrm{Re}^{3/4}$, so the largest Reynolds number such a code can represent is $\mathrm{Re} \le N^{4/3}$: $1.032\times10^{4}$ at $N = 1024$ and $6.554\times10^{4}$ at $N = 4096$. A wind-tunnel flow at $\mathrm{Re} = 10^{6}$ needs $3.162\times10^{4}$ cells per side, which is Module 10's number again. At $C = 0.5$ the upwind scheme falls short of this bound: by (3.3) it runs at $\mathrm{Re}_{\rm num} = 2N/(1-C)$, which at $N = 1024$ and $C = 0.5$ is $4096$, below $N^{4/3} = 1.032\times10^{4}$.

The cost of a contact. CHECK 2 says what resolution buys where the flow has a contact. Fitted over the last doubling, the contact error of §5 falls with order $0.5209$. To bring it from its value at $N = 6400$ to $10^{-4}$ needs $N = 1.864\times10^{5}$ cells. In three dimensions the work scales as $N^{3}$ cells times $N$ Courant-limited steps, so that is $7.19\times10^{5}$ times the work at $N = 6400$ for a factor $5.79$ in error. A code user who expects first-order convergence would budget $N^{4}$ against $\Delta x$, and would be short by the difference between an error falling as $\mathrm{work}^{-1/4}$ and one falling as $\mathrm{work}^{-1/8}$.

11. Which numbers in this book came from a simulation

The plan of this book asked this module to say which of its results came from a simulation and which from a closed form. Two readers went through Modules 1 to 13 and classified every number that a module compares against or takes as an input: (a) the output of a numerical simulation; (b) a numerical model built by solving equations on a computer, such as the BS05 solar model table or the ATLAS9 model atmosphere; (c) an ordinary differential equation with no closed form, integrated numerically, by this book's own script or, in 3 of the 30 rows, by a published integration (Chandrasekhar's in Module 3, Kandori's in Module 5), such as Lane–Emden, Bonnor–Ebert, Sedov, or the critical solutions of Bondi and Parker. The count is by table row, and one row can hold several numbers from one source.

The closed-form results that carry the book are the other column: the mean free paths of Module 1, the jump conditions and the $2/5$ exponent of Module 8, Kolmogorov's $5/3$ and $4/5$ in Module 10, the Eddington limit of Module 13. Class (c) is not simulation in the sense of this module: it is a quadrature whose error the book controls and reports, as Module 8 §4 does for $\xi_0$ to eight digits. Classes (a) and (b) are inputs whose error the book does not control.

The census also found places where a shipped page names a simulation or a numerical model as a measurement. They are listed here as findings. The pages are not edited.

The census did not re-open the verdicts of the modules that print these numbers. What changes is the kind of evidence: a check against a simulation tests whether two calculations agree, and a check against a measurement tests whether a calculation agrees with nature. The first seven rows are established from the shipped pages and their sources; the last is not, because the primary papers behind it were not read.

12. What this module captures, and what it misses

The module captures the three properties of a numerical scheme that decide whether its output can be trusted as a solution of the equations of this book: stability (the Courant condition), the equation it actually solves (the numerical viscosity), and its convergence on the discontinuities that astrophysical flows contain. It misses everything that makes a production code: more than one dimension, higher-order reconstruction, magnetic fields, self-gravity, radiation, and adaptive grids.

Of the debts the shipped pages and the plan printed against this module, six are paid: the Courant condition that Module 10 Proposition 4 used, by Propositions 2 and 3; which problem a code solves (module10.html:285, :918), by Propositions 4 and 5 and §10; how a collapse is followed (module05.html:767), by §9; Riemann solvers, by §4 and §5; grid against SPH, by §7 and §8; and the census, by §11. One is half paid, Pope's section, in the table above.

13. Problems

Notation

Sources

1. Courant, R., Friedrichs, K. & Lewy, H. (1928), "Über die partiellen Differenzengleichungen der mathematischen Physik", Math. Ann. 100, 32–74. Read in the English translation, "On the Partial Difference Equations of Mathematical Physics", IBM J. Res. Dev. 11 (1967) 215–234, translated by Phyllis Fox; the 1928 original was not seen. P. 228: the difference solution "in general cannot converge" for $k \lt 1$; p. 230: "If now $k \ge 1$ as we must require".
2. Sod, G. A. (1978), "A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws", J. Comput. Phys. 27, 1–31. Not obtained. Read instead: G. A. Sod, "A Survey of Numerical Methods for Compressible Fluids", Courant Institute report COO-3077-145 (June 1977), OSTI 7308221. P. 17 gives the initial states, $\gamma = 1.4$, $\Delta x = 0.01$ and Courant number $0.9$; Table II (p. 31) prints $0.426$, $0.927$, $0.303$, $0.266$, the constant states of Glimm's method on 9 interior points, not labelled as the exact solution. The pages read give no output time.
3. Marzouk, O. A. (2020), "The Sod gasdynamics problem as a tool for benchmarking face flux construction in the finite volume method", Scientific African 10, e00573, read as arXiv:2108.01709. Table 2, at $t = 0.2$, is CHECK 2 (a).
4. Lax, P. & Wendroff, B. (1960), "Systems of conservation laws", Comm. Pure Appl. Math. 13, 217–237. Not obtained. Read instead: the report of the same title, Los Alamos LA-2285 (1958), whose p. 9 states the theorem of Proposition 1 (b).
5. Godunov, S. K. (1959), Mat. Sb. 47(89), 271–306. Bibliographic record checked on Math-Net.Ru; the paper was not read, and nothing here depends on its text.
6. Banks, J. W., Aslam, T. & Rider, W. J. (2008), "On sub-linear convergence for linearly degenerate waves in capturing schemes", J. Comput. Phys. 227, 6985–7002, read as the authors' preprint (OSTI 945128). The smeared width $(4\nu t)^{1/2}$ with $\nu = O(\Delta x)$, hence $\Delta x^{1/2}$ at a contact.
7. Kamm, J. R. & Timmes, F. X. (2007), LA-UR-07-2849, their Sect. 6, as read for Module 8: $\rho_0 = 1$, $\gamma = 1.4$, $E = 0.851072$, shock at $r = 1.0$ at $t = 1.0$.
8. Monaghan, J. J. (1992), "Smoothed particle hydrodynamics", ARA&A 30, 543–574, read from the original scan. Density sum eq. 2.7 (p. 545); $\Pi_{ab}$ p. 550 (no equation number), $\mu_{ab}$ eq. 4.2; "$\alpha = 1$ and $\beta = 2$" and "$\eta^2 = 0.01h^2$", p. 551; the cubic spline and its normalisation $2/3$ in one dimension, pp. 554–555.
9. Price, D. J. (2008), "Modelling discontinuities and Kelvin–Helmholtz instabilities in SPH", J. Comput. Phys. 227, 10040–10057, read as arXiv:0709.2772v3; pages are the preprint's. The conductivity term eq. (27), its signal speed eq. (33); the "blip", p. 15 (density from the continuity equation) and p. 16 (density by summation); entropy against density gradients, p. 28; his shock tube at $\gamma = 5/3$ with 450 equal-mass particles.
10. Agertz, O. et al. (2007), "Fundamental differences between SPH and grid methods", MNRAS 380, 963–978, read as arXiv:astro-ph/0610051v1. The journal record lists 19 authors and the preprint 20; the preprint's text was read, and the quotation "Perturbations are damped out very quickly" is its p. 13.
11. Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, J. H., II, Howell, L. H. & Greenough, J. A. (1997), "The Jeans condition: a new constraint on spatial resolution in simulations of isothermal self-gravitational hydrodynamics", ApJ 489, L179–L183. $J \equiv \Delta x/\lambda_{\rm J}$, eq. (2), p. L180, and the value $0.25$; "necessary (but not, in general, sufficient)", p. L181.
12. Bate, M. R. & Burkert, A. (1997), "Resolution requirements for smoothed particle hydrodynamics calculations with self-gravity", MNRAS 288, 1060–1072, read from the ADS scan. $M_{\rm res} \approx M_{\rm tot}(2N_{\rm neigh}/N_{\rm tot})$, eq. (6), p. 1068; the maximum resolvable density, eq. (7), p. 1068; $N_{\rm neigh} = 50$.
13. Toro, E. F., Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer. Not read; cited for the name of the "123 problem" of Problem K1 only. Its values in K1 are this module's own run.
14. Pope, S. B. (2000), Turbulent Flows, Cambridge University Press. Not obtained. Only the caption of his Fig. 9.3 was read, with the asymptote $(R_\lambda/70)^6$.
15. The census of §11 was compiled by two readers of Modules 1–13 and is kept with the book's working records; every row carries the file:line it was read from.

Contents · ← Module 13