Skip to content

Tutorial 1 — Your first run: Sod shock tube

Time: ~10 minutes. Prerequisites: aeolus1d installed, a working terminal. Anchored on: examples/01_single_pipe_sod.toml.

The Sod problem is the textbook benchmark for compressible-flow solvers. A long pipe with a notional diaphragm at its centre separates high-pressure gas on the left from low-pressure gas on the right. At t = 0 the diaphragm bursts, sending a rarefaction wave left and a shock + contact discontinuity right. The exact solution is known analytically (Sod 1978), so we use it to validate every numerical scheme aeolus1d ships.

After this tutorial you'll have run aeolus1d, opened the HDF5 output, plotted it, and compared against the analytical answer — all the basic verbs.

Step 1 — Look at the case

Open examples/01_single_pipe_sod.toml in any text editor. The whole case is ~30 lines. Read top to bottom:

[case]
name = "sod_shock_tube"
t_end = 0.2          # simulate to t = 0.2 (Sod's reference)
scheme = "muscl"     # 2nd-order spatial reconstruction

[gas]
gamma = 1.4          # ratio of specific heats
R_gas = 1.0          # Sod is non-dimensional

The [case] block sets globals: end time and which numerical scheme. It is written in the compact form: every CaseHeader knob omitted from [case] fills from its schema default — here that means cfl = 0.5 (the Courant number that controls the time step) and limiter = "minmod" (the TVD limiter for the MUSCL slopes), among ~60 others. aeolus1d expand case.toml --full rewrites the block with every field explicit if you want to see them all; run-settings.md documents each. muscl + minmod is the standard 2nd-order TVD combination; the other scheme options are first_order (simpler, more diffusive) and mp5 (5th-order, much sharper).

[[pipe]]
id = "tube"
length = 1.0
n_cells = 400
area = 1.0

  [pipe.init]
  kind = "riemann"
  x0 = 0.5
    [pipe.init.left]
    rho = 1.0
    u = 0.0
    p = 1.0
    [pipe.init.right]
    rho = 0.125
    u = 0.0
    p = 0.1

One pipe, length 1.0, divided into 400 cells. The riemann init kind places a discontinuity at x0 = 0.5 with the left state on one side and the right state on the other — that's the diaphragm.

[[bc]]
pipe = "tube"
end = "left"
kind = "transmissive"

[[bc]]
pipe = "tube"
end = "right"
kind = "transmissive"

Both ends are transmissive (zero-gradient): waves leave without reflecting. For a Sod test on this time scale the waves don't reach the boundaries anyway, but it's good hygiene.

Step 2 — Run it

There are two equivalent paths. The terminal is faster for batch / scripted use; the GUI is the one to pick when you want to see the simulation tube paint up live or click around the result. Pick whichever:

Option A — Terminal

aeolus1d run examples/01_single_pipe_sod.toml --out /tmp/sod.h5

Expected output (abridged — the solver-settings banner lists every active knob, flagging non-defaults with *):

Running case 'sod_shock_tube' to t_end=0.2 s ...
Aeolus solver settings
  driver                   : explicit single-pipe (ssp_rk2)
  topology                 : 1 pipe, 0 junctions
  ...
  ( * = non-default vs schema )
  finished at t = 0.2 s
  wall time   = 0.36 s  (1.8x realtime)
  pipe 'tube'      : rho=[0.125, 1] u=[0, 0.9294] p=[0.1, 1]
  -> wrote /tmp/sod.h5
  -> wrote /tmp/sod.manifest.json

The solver chose its own time step from the CFL constraint (default cfl = 0.5); no setup beyond what you saw in the TOML. (On a fresh environment the very first run pays a one-time Numba kernel compile — run aeolus1d warmup once to get it out of the way; on an interactive terminal you'll also see periodic t = <now> / <t_end> progress lines.)

Option B — GUI

aeolus1d gui

(or python -m aeolus1d_gui from a repo checkout, or the aeolus1d_gui/run_gui.sh helper.)

Then:

  1. File → Open → examples/01_single_pipe_sod.toml. The tube appears as a single rectangle on the canvas with a transmissive BC marker on each end.
  2. Click Run in the bottom Runner dock (or press F5). The GUI saves the case to a temp dir, launches aeolus1d run as a subprocess, and streams the output into the dark console strip below the canvas.
  3. On a clean exit (code 0), the post-processing tab automatically loads the output HDF5 and renders the pipe as a coloured tube. Use the field-selector combo to flip between ρ, u, p, T.
  4. The 2D Plot dock on the right shows the cell-resolved profile of whatever field you've selected; drag the time slider to scrub through the run if you saved a time series (this case doesn't — go to Step 4 for the analytical comparison).

Step 3 — Look at the output

import h5py, numpy as np
with h5py.File("/tmp/sod.h5") as f:
    pipe = f["pipes/tube"]
    x = pipe["x"][:]      # cell-centre positions
    rho = pipe["rho"][:]  # density at t_final
    u = pipe["u"][:]      # velocity
    p = pipe["p"][:]      # pressure

print(f"rho range: {rho.min():.3f} → {rho.max():.3f}")
print(f"u   range: {u.min():.3f} → {u.max():.3f}")
print(f"p   range: {p.min():.3f} → {p.max():.3f}")

Expected:

rho range: 0.125 → 1.000
u   range: 0.000 → 0.929
p   range: 0.100 → 1.000

The ranges match the initial high / low states because the rarefaction's leading edge hasn't reached the left boundary, the shock's leading edge hasn't reached the right, and the high-density gas hasn't compressed any further than its initial value (it's expanding, not compressing).

Step 4 — Plot it

The repo ships a plot script that does the comparison against the analytical solution:

python examples/plots/plot_01_sod.py

That writes /tmp/aeolus_example_01.png. Open it. You'll see three panels (ρ, u, p), each comparing the numerical run to the analytical Riemann solution. The numerical curves should overlay the analytical to within a cell or two on every wave (rarefaction fan, contact discontinuity, shock).

The L1 error on density is reported on the figure — for the default 400-cell mesh with MUSCL+minmod it should be around 0.4%. aeolus1d's CI gate for this case is L1(ρ) < 2%.

Step 5 — Modify and rerun

Knobs to try (each is a one-line edit + rerun):

  • Mesh refinement. Bump n_cells = 400 to 800 and re-run. The L1 error should roughly halve (1st-order convergence near discontinuities is the textbook result for TVD MUSCL).

  • Drop to first order. Change scheme = "muscl" to scheme = "first_order". The shock and contact will smear noticeably; L1 will roughly triple. (The limiter setting only applies to MUSCL — it's harmless to leave at its default.)

  • Bump to 5th order. Change scheme = "mp5" and re-plot. The contact discontinuity becomes much sharper; L1 drops by roughly 2–3× on the same mesh. MP5 needs Ng = 3 ghost cells, which aeolus1d configures automatically based on the scheme.

  • Different test. Change the left/right rho, u, p values to a different shock-tube problem (Lax, Shu-Osher) and see the same machinery handle it.

What you've done

  • Loaded a TOML case via the CLI.
  • Read the HDF5 output programmatically with h5py.
  • Compared a numerical run to an analytical reference.
  • Probed the schemes by toggling one knob.

That's the basic loop: edit TOML → run CLI → read HDF5 → plot. Every more-complex case in aeolus1d is the same loop with more blocks in the TOML.

Next: Networks and junctions →

Now that one pipe works, learn how to wire several pipes into a network through 0-D plena.