Expand description
§One dimension
Cartesian, one space dimension, on $[0,1]$ unless stated, split into 4 blocks, one thread.
Problem files are in problems/1d/. Test paths are relative to crates/merlin-finite-volume/tests/
unless they start with crates/. “Drift” is the relative change of an integrated quantity
between the first and last snapshot.
| Case | Reference | Headline result |
|---|---|---|
| Linear advection | exact translation | upwind order 0.99, Lax-Wendroff 2.00, MUSCL 2.03 in $L_1$ |
| Burgers shock | Rankine-Hugoniot | shock at 0.50001 (exact 0.5) |
| Burgers rarefaction | exact fan | within bounds, $L_1$ 5.9e-4 |
| Transonic fan | exact fan | opens correctly; the unfixed Roe flux does not |
| Smooth data forming a shock | characteristics | shock values 0.847570 vs 0.847576; order 1.99 to 2.07 in $L_1$ |
| Toro’s star values | published table | five of five |
| Sod’s tube | exact Riemann solution | plateau $p_*$ 0.30313 vs 0.30313 |
| Stationary contact | exact (the initial data) | density error exactly 0 |
| Toro’s tests 2, 3, 4 | exact Riemann solution | positive, right plateaus |
| Smooth density wave | exact translation | MUSCL order 2.02 to 2.04 |
§Linear advection
Problem. A square wave and a smooth bump are carried once around a periodic unit domain at speed 1
(1000 cells, $t=1$). After one period the exact solution is the initial data.
Result.
| Problem file | Scheme | Steps | $L_1$ error | Min / max of $u$ | Drift |
|---|---|---|---|---|---|
advection/upwind.toml | first-order upwind, forward Euler, CFL 0.8 | 1256 | 2.67e-2 | 1.2e-28 / 1.0000 | 2.70e-15 |
advection/lax_wendroff.toml | Lax-Wendroff, CFL 0.8 | 1256 | 1.36e-2 | -0.222 / 1.222 | 1.11e-15 |
advection/muscl.toml | MUSCL, van Leer, SSP-RK3, CFL 0.5 | 2000 | 8.82e-3 | -1.8e-78 / 1.0000 | 7.71e-14 |
Observed order on smooth data $\sin 2\pi x$ (128 to 2048 cells, one period, CFL 0.5; fitted slope):
| Scheme | $L_1$ | $L_2$ | $L_\infty$ |
|---|---|---|---|
| upwind | 0.988 | 0.988 | 0.987 |
| Lax-Wendroff | 2.000 | 2.000 | 2.000 |
| MUSCL (van Leer) + SSP-RK3 | 2.032 | 1.756 | 1.392 |
Reading it. Upwind is first order and Lax-Wendroff second order in every norm. MUSCL is
second order in $L_1$ and lower in $L_2$ and $L_\infty$ because the limiter flattens smooth
extrema, as designed (Numerical methods). All three integrate the
solution to round-off (drift 1e-15 to 8e-14). The overshoots of Lax-Wendroff (1.222 and -0.222)
and the smearing of upwind are the expected costs of a linear scheme of each order; see
Known limitations.
Tests that enforce this
Enforced by schemes/upwind_convergence.rs: upwind_is_first_order_in_every_norm,
upwind_is_first_order_for_leftward_waves; schemes/lax_wendroff.rs: lax_wendroff_is_second_order_in_every_norm,
courant_number_one_translates_data_exactly, lax_wendroff_overshoots_a_square_wave_while_upwind_stays_bounded;
schemes/muscl.rs: van_leer_is_second_order_in_l1, minmod_is_nearly_second_order_in_l1,
monotonized_central_is_nearly_second_order_in_l1, limited_schemes_create_no_new_extrema_and_do_not_increase_total_variation.
The results do not depend on the block count (errors_do_not_depend_on_the_number_of_blocks), and the HDF5 run
holds exactly the numbers of the CSV run (pymerlin check: 9 snapshots of 1000 cells and the integral history, bit-identical).
§Burgers shock and rarefaction
$u_t+(u^2/2)_x=0$, 400 cells on $[0,1]$, open ends. The errors are against the exact Riemann solution.
Shock (burgers/shock.toml): $u=1$ left of $x=0.3$ and 0 right of it; first-order Godunov flux, CFL 0.8,
$t=0.4$. The Rankine-Hugoniot speed is $(1+0)/2$, so the shock should be at $x=0.5$.
- Shock position, from the crossing of
$u=0.5$: 0.50001 (exact 0.5). - The jump occupies two cells;
$L_1$error 8.81e-4. - The integral grows from 0.3 to 0.5, exactly the inflow through the open left end: the flux
$1/2$for time 0.4 puts in$0.2$. The change is what the boundary flux accounts for.
Rarefaction (burgers/rarefaction.toml): $u=0$ left, 1 right; MUSCL (van Leer), Godunov flux, SSP-RK3, CFL 0.5,
$t=0.4$. The jump opens into the fan $u=(x-0.3)/t$.
$L_1$error 5.92e-4,$L_\infty$8.08e-3; the solution stays inside the initial bounds$[0,1]$(to 1e-6).- The integral falls from 0.7 to 0.5 because
$0.5\times0.4=0.2$leaves through the right end.
Tests that enforce this
Enforced by equations/burgers.rs: every_scheme_puts_the_shock_where_rankine_hugoniot_says,
a_shock_converges_at_first_order_for_every_scheme, a_rarefaction_fan_converges,
the_shock_is_sharp_and_stays_inside_the_initial_bounds, monotone_schemes_obey_the_maximum_principle_and_are_tvd_through_a_shock,
local_lax_friedrichs_smears_a_shock_more_than_godunov; the exact solutions themselves by equations/burgers_reference.rs:
a_shock_moves_at_the_rankine_hugoniot_speed, a_rarefaction_fan_is_linear_between_its_edges,
the_averages_integrate_to_the_boundary_flux; end to end from the files by crates/merlin-cli/tests/run/burgers_problems.rs:
the_shock_sits_at_the_rankine_hugoniot_position, the_rarefaction_is_the_linear_fan.
§Transonic rarefaction
burgers/transonic.toml: $u=-1$ left and +1 right of $x=0$ on $[-1,1]$, first-order Godunov, $t=0.5$. The
fan passes through the sonic point $u=0$, which is where a flux without an entropy fix fails.
$L_1$error 1.32e-2,$L_\infty$3.62e-2 (first-order smoothing of the two fan corners).- The two cells at the origin read -0.019 and +0.019, symmetric about zero.
- With
riemann_solver = "roe_no_fix"the same problem keeps a frozen jump from -1 to +1 ($L_1$error 0.500), an expansion shock that is not the solution. This is shown on purpose in Known limitations, and the same test suite proves the entropy-satisfying fluxes (Godunov, Engquist-Osher, local Lax-Friedrichs) do not fail.
Tests that enforce this
Enforced by equations/burgers.rs: the_entropy_satisfying_fluxes_open_a_transonic_jump_into_a_fan,
the_roe_flux_without_a_fix_keeps_an_entropy_violating_stationary_jump, the_roe_flux_without_a_fix_is_right_at_a_shock;
equations/burgers_reference.rs: the_transonic_fan_is_odd_about_the_origin;
crates/merlin-cli/tests/run/burgers_problems.rs: the_transonic_problem_opens_into_a_fan_with_a_correct_flux,
the_unfixed_flux_leaves_the_transonic_jump_frozen.
§Smooth data forming a shock
burgers/sine.toml: $u=\sin 2\pi x$, periodic; MUSCL (van Leer), Godunov, SSP-RK3, CFL 0.5. The characteristics cross at
$t_b=1/(2\pi)\approx0.159$, so the exact solution is known before that time, and the shock position and values are
known after it from the characteristic equation.
$t$ | $L_1$ error | $L_\infty$ error |
|---|---|---|
| 0.05 | 1.16e-5 | 3.39e-4 |
| 0.10 | 2.35e-5 | 5.19e-4 |
| 0.15 (just before the shock) | 8.04e-5 | 2.50e-3 |
- At
$t=0.4$the shock sits at$x=0.5$(by symmetry) with values$\pm0.847570$. Solving$\xi+t\sin2\pi\xi=\tfrac12$gives 0.847576. $\max|u|\le0.99997$at every snapshot: no new extrema.- The mean of
$u$stays at zero (7e-17 at$t=0.4$).
Observed order, 100 to 800 cells, $t=0.1$, CFL 0.4 (between successive resolutions):
| Scheme | $L_1$ orders | $L_2$ orders |
|---|---|---|
| MUSCL + SSP-RK3, Godunov | 1.99, 2.02, 2.07 | 1.75, 1.76, 1.76 |
| first order + forward Euler, Godunov | 0.94, 0.96, 0.98 | 0.87, 0.93, 0.96 |
Tests that enforce this
Enforced by equations/burgers.rs: muscl_is_second_order_on_a_smooth_solution,
first_order_schemes_converge_at_first_order_on_a_smooth_solution, limited_muscl_creates_no_new_extrema_through_a_shock;
equations/burgers_reference.rs: the_smooth_solution_satisfies_its_characteristic_relation,
the_smooth_solution_is_periodic_and_conserves_its_mean; crates/merlin-cli/tests/run/burgers_problems.rs:
smooth_data_forms_a_shock_and_the_integral_stays_zero.
§Toro’s star values
The exact Euler solver is the reference for everything below, so it is checked first. Toro (2009, Tables 4.1 and 4.2) gives
five tests with their star-region pressure and velocity. pymerlin check --self-test euler runs the independent Python solver
on them with a relative tolerance of 2e-5:
| Test | $p_*$ (script) | $u_*$ (script) | Published: $p_*$, $u_*$ |
|---|---|---|---|
| 1, Sod | 0.30313 | 0.927453 | 0.30313, 0.92745 |
| 2, “123” | 0.00189387 | 0 | 0.00189, 0 |
| 3, strong left blast | 460.894 | 19.5975 | 460.894, 19.5975 |
| 4, colliding shocks | 1691.65 | 8.68977 | 1691.64, 8.68975 |
| 5, stationary contact (moving frame) | 460.894 | 1.4e-6 | 460.894, 0 |
(1691.65 against 1691.64 is the last printed digit; the script’s value is 1691.647 to more digits.) The Rust solver is tested against the same table, and mass, momentum and energy are conserved by the sampled solution in each.
Tests that enforce this
Enforced by riemann/euler_exact.rs: the_star_region_matches_the_published_values_for_all_five_tests,
the_solution_conserves_mass_momentum_and_energy, pressure_and_velocity_are_continuous_across_the_contact,
mirroring_the_problem_mirrors_the_solution.
§Sod’s shock tube
euler/sod.toml: $\gamma=1.4$, 400 cells, MUSCL (van Leer), SSP-RK2, CFL 0.4, HLLC, outflow ends, $t=0.25$
(550 steps). Left $(\rho,u,p)=(1,0,1)$, right $(0.125,0,0.1)$, diaphragm at 0.5.
- The plateau between the rarefaction and the contact (
$x$from 0.55 to 0.70) has$p=0.30313$,$u=0.92746$; exact$p_*=0.30313$,$u_*=0.927453$. $L_1$errors: density 1.65e-3, velocity 2.87e-3, pressure 1.01e-3.- Mass drift 1.97e-16, energy drift 1.61e-16. Momentum rises by exactly 0.225, which is
$(1-0.1)\times0.25$: the pressure difference at the two open ends times the time.
Tests that enforce this
Enforced by crates/merlin-cli/tests/run/euler_problems.rs: sods_problem_has_the_published_plateau_values,
a_gas_run_reports_mass_momentum_and_energy; equations/euler.rs: every_riemann_flux_converges_on_sods_problem,
a_periodic_wave_conserves_mass_momentum_and_energy_to_round_off, a_closed_tube_conserves_mass_and_energy_through_many_reflections,
a_reflecting_wall_acts_like_a_mirror_image.
§Stationary contact
euler/contact.toml, $t=0.5$ (1676 steps): only the density jumps (1 to 0.125) at velocity 0 and pressure 1, so the exact
solution is the initial data for all time.
The density error is exactly 0, $u$ and $p$ stay exactly uniform, and the drift of mass, momentum and energy is 0.
With HLL in place of HLLC the contact is smeared over 43 cells (Known limitations); that is the
contrast that makes this a meaningful test of the flux.
Tests that enforce this
Enforced by riemann/euler_fluxes.rs: a_stationary_contact_is_kept_exactly_by_hllc_and_the_exact_flux_only,
hllc_carries_a_moving_contact_without_diffusion; equations/euler.rs:
a_stationary_contact_is_kept_by_hllc_and_the_exact_flux_but_smeared_by_the_others;
crates/merlin-cli/tests/run/euler_problems.rs: hllc_keeps_the_stationary_contact_and_hll_smears_it.
§Water tube (a stiffened gas)
euler/water_tube.toml, $t=180\,\mu\text{s}$: water at $10^9$ Pa against water at $10^5$ Pa, a stiffened gas
$p=(\gamma-1)\rho e-\gamma p_\infty$ with $\gamma=4.4$, $p_\infty=6\times10^8$ Pa (Ivings, Causon and Toro 1998), on
400 cells with MUSCL (van Leer), SSP-RK2 and HLLC; HLLC is the only flux the program offers for a fluid other than the
ideal gas (Known limitations).
The reference is the Riemann solution of the ideal gas in the shifted pressure $\tilde p=p+p_\infty$: with the energy
$E-p_\infty$ the stiffened-gas equations are the ideal-gas equations, because a constant drops out of $p_x$ and
$u(E+p)=u\big((E-p_\infty)+\tilde p\big)$. The solver of that solution is itself checked against Toro’s published values
(above). The $L_1$ error as a fraction of the range of each variable, with the observed order of the
fit over 100 to 800 cells:
| Variable | 100 | 200 | 400 | 800 | Order |
|---|---|---|---|---|---|
| density | 1.55e-2 | 8.31e-3 | 4.35e-3 | 2.76e-3 | 0.84 |
| pressure | 7.94e-3 | 4.19e-3 | 1.90e-3 | 1.04e-3 | 0.99 |
The order is that of Toro’s tests: the discontinuities limit the $L_1$ convergence to about one, and the density
(which carries the contact) falls a little slower than the pressure. Every wave stays inside the tube at the output time, so the open ends play no part
(mass and energy are conserved to round-off, and the momentum changes by the pressure difference times the time, as it must). The same scheme run on the stiffened gas and on the
ideal gas in the shifted pressure gives the same flow to $10^{-10}$ of the largest value, and the result does not depend on
the number of blocks. The run through the program (merlin water_tube.toml) is within 1e-2 of the density range of the exact solution, and its
plateau pressure is within 2 % of the exact star pressure.
Tests that enforce this
Enforced by equations/euler_stiffened.rs: a_water_shock_tube_converges_to_its_exact_solution,
a_stiffened_gas_is_an_ideal_gas_in_the_shifted_pressure, the_stiffened_runs_do_not_depend_on_the_block_layout;
crates/merlin-cli/tests/run/euler_problems.rs: the_water_tube_matches_the_exact_stiffened_gas_solution.
§Copper impact (a Mie-Gruneisen solid)
euler/copper_impact.toml, $t=6\,\mu\text{s}$: two copper plates approach at 1000 m/s each (density 8930 kg/m^3, $c_0=3940$ m/s,
$s=1.489$, $\rho\Gamma=\Gamma_0\rho_0$ with $\Gamma_0=2$; Marsh 1980 as given by Menikoff and Plohr 1989, Section VI), on 400 cells
with MUSCL (van Leer), SSP-RK2 and HLLC. Where they meet each is brought to rest by a shock. The shock-speed relation of the material,
$U_s=c_0+s\,u_p$, is the speed of the shock relative to the plate it runs into, so with $u_p=1000$ m/s it gives $U_s=5429$ m/s, the
shocked copper at rest at $p^*=\rho_0U_su_p=48.48$ GPa and density $\rho_0U_s/(U_s-u_p)=10\,946$ kg/m^3, and, because the plate moves towards
the shock, a shock that runs at $U_s-u_p=4429$ m/s in the tube, $0.0266$ m from the middle at the output time. That is the whole
solution, exact and from the relation alone: the model is built from it, but the scheme, the conversions between conserved and primitive
variables and the wave speeds are not, and have to reproduce it.
| Variable | 100 | 200 | 400 | 800 | Order |
|---|---|---|---|---|---|
| density | 1.11e-2 | 7.85e-3 | 3.43e-3 | 1.41e-3 | 1.01 |
| pressure | 7.09e-3 | 6.32e-3 | 2.54e-3 | 8.96e-4 | 1.03 |
($L_1$ error as a fraction of the jump, observed order of the fit.) At 400 cells the pressure between the shocks is $p^*$ to 1e-4
(measured spread 2e-5), the velocity there is zero to 0.005 m/s, and each shock is within 0.1 of a cell of $U_s-u_p$ from the middle. The
density on that plateau is about 5e-4 low, the “wall heating” that Godunov-type schemes show at a symmetry plane that starts as a jump. The
plates keep entering through the open ends, so mass and energy grow by exactly what flows in ($2\rho_0u_pt$ and $\rho_0u_p^3t$), to round-off, and the
momentum stays zero; the solution is mirror-symmetric and the same for any block layout, bit for bit. pymerlin plot draws no reference for
this problem (a Mie-Gruneisen Riemann problem has no solution in the Python tools).
Tests that enforce this
Enforced by equations/euler_mie_gruneisen.rs: a_copper_impact_converges_to_the_shock_speed_relation,
the_shocked_copper_is_at_the_pressure_the_relation_gives_and_the_shocks_are_in_place,
a_copper_impact_conserves_mass_momentum_and_energy_through_its_open_ends,
the_impact_is_mirror_symmetric_and_independent_of_the_block_layout;
crates/merlin-cli/tests/run/euler_problems.rs: the_copper_impact_matches_the_shock_speed_relation.
§Toro’s other tests
Test 2, the “123” problem (euler/123.toml, $t=0.15$, 416 steps): two gases rush apart at speed 2, leaving a
near-vacuum (exact $p_*=0.00189$).
- Density and pressure stay positive (minimum 0.0200 and 0.0033 at
$t=0.15$); momentum stays at 7e-17 (a symmetric problem). $L_1$errors: density 2.34e-3, velocity 1.21e-2, pressure 1.05e-3.- Visible weakness. In the central cells the pressure is 0.00335 against the exact 0.00189 and the density 0.0200 against
0.0219. With 1600 cells the pressure there is 0.00216 and the
$L_1$errors fall (density 9.50e-4, pressure 2.59e-4). The integrated errors converge; the central values of this problem converge slowly.
Test 3, strong blast (euler/blast.toml, $t=0.012$, 643 steps): the pressure jumps by $10^5$.
- Plateau between the fan and the contact:
$p=460.93$(range 460.78 to 460.96) and$u=19.5955$; exact$p_*=460.894$,$u_*=19.5975$. - Nothing went negative: minimum density 0.568, minimum pressure 0.01 (the undisturbed right state).
- Mass drift 0, energy drift 3.6e-16, and the momentum gain 11.99988 equals
$(1000-0.01)\times0.012$. - Visible weakness. A 2 % overshoot in
$u$at the foot of the fan (20.008 against 19.5975), and the contact is smeared over 7 cells (10 to 90 %).
Test 4, colliding shocks (euler/colliding_shocks.toml, $t=0.035$, 1052 steps):
- Central plateau:
$p$mean 1691.06 (range 1682.8 to 1700.2) and$u$mean 8.687; exact$p_*=1691.647$,$u_*=8.68977$. The ripples are about 0.5 % of the pressure; the contact is 7 cells wide (10 to 90 %). - The ends are open, so the integrals change; no conservation claim is made for this problem.
Tests that enforce this
Enforced by equations/euler.rs: toro_tests_converge_to_the_exact_solution_with_hllc,
the_extreme_problems_stay_positive_and_finite_with_every_flux; crates/merlin-cli/tests/run/euler_problems.rs:
the_colliding_shocks_make_the_published_central_state, the_extreme_problems_stay_physical.
§Smooth density wave
$\rho=1+0.2\sin2\pi x$, $u=1$, $p=1$, periodic, one period (the exact solution is the initial data), HLLC, 100 to 800
cells, CFL 0.4:
| Scheme | $L_1$ order of density, between successive resolutions |
|---|---|
| MUSCL + SSP-RK3 | 2.02, 2.03, 2.04 |
| first order + forward Euler | 0.94, 0.97, 0.99 |
Tests that enforce this
Enforced by equations/euler.rs: a_smooth_density_wave_converges_at_second_order_with_muscl,
a_smooth_density_wave_converges_at_first_order_without_reconstruction; the one-dimensional results do not depend on the block count
(results_do_not_depend_on_the_block_layout, results_do_not_depend_on_the_number_of_blocks).
§What holds in every 1D problem
- Conservation: periodic and closed problems conserve mass, momentum and energy to round-off, and open problems change by exactly the boundary flux, as the Sod momentum shows.
- Block independence: any number of blocks and threads gives the same bits.
- Positivity: the near-vacuum and the
$10^5$pressure-ratio problems stay positive with every flux. - File against library: a run from the TOML file equals the library-driven run bit for bit
(
crates/merlin-cli/tests/run/run_from_file.rs:a_file_run_matches_the_library_driven_run_bit_for_bit,a_muscl_run_matches_the_library_driven_run_bit_for_bit).
§References
- Toro, E. F. (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer, Chapter 4 and Tables 4.1 and 4.2 (the five tests).
- LeVeque, R. J. (2002), Finite Volume Methods for Hyperbolic Problems, Cambridge, Chapters 11 and 12.
- Menikoff, R. and Plohr, B. J. (1989), The Riemann problem for fluid flow of real materials, Rev. Mod. Phys. 61, the Mie-Gruneisen law (Section VI) and the wave curves of a general equation of state.
- Marsh, S. P. (ed.) (1980), LASL Shock Hugoniot Data, University of California Press, the shock-speed relation of copper.
- Ivings, M. J., Causon, D. M. and Toro, E. F. (1998), On Riemann solvers for compressible liquids, Int. J. Numer. Meth. Fluids 28, the stiffened gas and its exact Riemann solution.