Module cylindrical 

Cylindrical coordinates

The Euler equations in cylindrical symmetry: a one-dimensional radial problem (geometry = "cylindrical", one axis, the radius $r$) and the axisymmetric $(r,z)$ problem (two axes, radius first). The radial sweep weights the fluxes with face areas and divides by the cell volume, with the pressure source in the form that cancels for a gas at rest (Numerical methods); the $z$ sweep is the Cartesian one. Swirl (azimuthal velocity) has its own section. Problem files are in problems/cylindrical/. Test paths are relative to crates/merlin-finite-volume/tests/ unless they start with crates/. The spherical counterpart is spherical, and the source terms are what the two pages test.

Case	Reference	Headline result
Gas at rest	exact equilibrium	$L(u)$ exactly zero
Sound wave	exact linear solution (Hankel transform, in Python)	order 1.97, 1.97, 2.06 (fit 2.00)
Noh’s implosion	exact	shock at 0.19875 (exact 0.2); plateau 15.93 (exact 16)
Sedov-Taylor blast	exact self-similar solution (in Python)	shock radius 0.71375 (exact 0.71310); error 2.5e-3 at 800 cells
Blast	a fine run in the same geometry	order 1.08 on the finest pair
Axisymmetric can	conservation	drift 0 / 0
Swirl	exact steady state and invariants	angular momentum drift 0; rotation error order 3.3, 3.1, 2.9
Against 3D	the 3D Cartesian ball	1.46 % relative $L_1$ difference

A gas at rest

A uniform gas at rest must stay at rest. With a geometric source term written as a separate approximation, the radial flux difference and the source cancel only to truncation error, and a gas at rest slowly acquires motion. Here the update is exactly zero: every value of $L(u)$ equals 0.0, for $r_{\min}=0$ and 0.5, first order and MUSCL, and with an axial flow in 2D.

Tests that enforce this

Enforced by geometry/operator.rs: a_gas_at_rest_gets_exactly_zero_rhs_with_the_first_order_scheme, a_gas_at_rest_gets_exactly_zero_rhs_with_muscl_too, a_gas_at_rest_with_axial_flow_stays_as_it_is_in_two_dimensions.

Outgoing sound wave

A pulse of relative amplitude $10^{-7}$ started at rest, small enough that the nonlinear terms are negligible, so the exact linear solution is the reference. A cylinder has no single-profile solution, so the reference is a Hankel transform computed in Python (src/pymerlin/exact_solutions/acoustic.py, with the Bessel functions from their integral form; see Equations). MUSCL, HLLC, SSP-RK3, 100 to 800 cells, $t=0.3$:

Cells	100	200	400	800
$L_1$ density error (units of the pulse)	9.00e-3	2.29e-3	5.86e-4	1.41e-4
Observed order		1.97	1.97	2.06

Fitted order 2.00. The values at cell centres are used for the initial data and the reference.

This study runs from Python, not cargo test, and runs in continuous integration: uv run pymerlin convergence cylindrical/euler/acoustic --expect-order=2 --tolerance=0.4. Alongside it, uv run pymerlin check --self-test acoustic checks the reference itself: the wave-equation residual of the solution and the initial data of the shipped problem files, so the files and the reference cannot drift apart. The Rust test of the same source terms is the spherical one, next page.

Noh’s implosion

A cold gas flows towards the axis at unit speed and a shock leaves it at $s=(\gamma-1)/2$, with the gas at rest behind. The solution is exact in this geometry (Noh 1987; see Equations): for $\gamma=5/3$ the shock is at $r=0.2$ at $t=0.6$ with density 16 behind it. 400 cells (cylindrical/euler/noh.toml).

The exact density at $t=0.6$ (cylinder and sphere). Behind the shock at $r=0.2$ the gas is at rest at the plateau density; ahead of it the inflow compresses the gas as $(1+t/r)^{\nu-1}$.

• Shock at 0.19875 against 0.2: half a cell.
• Mean density behind the shock 15.93 against 16.
• $L_1$ density error over $0.05<r<0.35$ at 100, 200 and 400 cells: 0.145, 0.0402, 0.0261; observed orders 1.85, then 0.62.
• MUSCL beats first order by more than 20 % at 200 cells.

Known error. The first few cells at the axis are wrong (a density spike and a dip). This “wall heating” is the standard Noh artefact of Godunov-type schemes, and it shrinks with the cell size without vanishing; the shock and the wall-heated centre limit the $L_1$ error to about first order, and the coarse pair is not yet in that range. The test therefore compares away from the centre.

Tests that enforce this

Enforced by geometry/noh.rs: the_density_behind_the_shock_is_the_exact_plateau, the_shock_is_at_the_exact_radius, the_error_falls_as_the_grid_is_refined, muscl_is_more_accurate_than_first_order_on_the_implosion; the shipped file by crates/merlin-cli/tests/run/curvilinear_problems.rs: noh_files_show_the_exact_gas_behind_the_shock.

Sedov-Taylor blast

A hot disc of radius 0.02 ($p=10^4$ against $10^{-3}$, energy per unit length $E=31.4$) in a cold gas at rest drives a strong cylindrical shock. Once the shock is several radii of the disc away, the flow forgets how the energy was deposited and becomes the self-similar solution of Sedov (1959), with $R(t)=(Et^2/\alpha\rho_0)^{1/4}$ and $\alpha=0.98408$ for $\gamma=1.4$, computed in Python (src/pymerlin/exact_solutions/sedov.py; the solution is derived on the Equations page). 400 cells, $t=0.09$ (cylindrical/euler/sedov.toml):

• Shock radius 0.71375 against 0.71310: 0.1 %, a quarter of a cell.
• The density peak is 5.29 against the exact 6: a shock smeared over about three cells cannot reach it.
• Mass and energy drift 1.2e-13 and 5.0e-13 (the shock has not reached the open outer end).
• $L_1$ density error against the exact solution (units of the density range) at 100, 200, 400, 800 cells: 1.38e-2, 6.01e-3, 4.11e-3, 2.54e-3, observed orders 1.19, 0.55, 0.70 (fit 0.79).

Why this is a bound and not an order. The solution is the point-source limit, and the simulation starts from a disc of finite radius, so the error does not converge to zero however fine the grid: it falls to a few parts in a thousand and stays there. The check is therefore that the finest error is below 0.005. The reference is only drawn once the shock is 6 radii out and the pressure ratio is at least $10^3$, so the blast problem above, whose ball is a quarter of the domain, gets none.

Both checks run from Python, in continuous integration, not in cargo test: uv run pymerlin check --self-test sedov and uv run pymerlin convergence cylindrical/euler/sedov --max-error=0.005. The self-test is what makes the reference trustworthy: $\alpha$ for the sphere with $\gamma=1.4$ is Sedov’s 0.851, the mass inside the shock equals the mass swept up for all three geometries and $\gamma=1.2,1.4,5/3$ (relative error below 1e-4), and the sampled solution satisfies the continuity, momentum and entropy equations by finite differences in $r$ and $t$ (residuals below 1e-3), independently of the derivation.

Blast against a fine run

A hot core $r<0.25$ ($p=10$ against 0.1) in a closed tube, $t=0.2$ (cylindrical/euler/blast.toml). There is no exact solution, so the reference is a 1600-cell run in the same geometry, and the coarse runs (50 to 400 cells) converge to it:

Cells	50	100	200	400
$L_1$ density error	8.4e-2	7.4e-2	4.4e-2	2.1e-2

The order on the finest pair is 1.08, as expected for a flow containing a shock.

Tests that enforce this

Enforced by geometry/blast.rs: a_blast_converges_to_the_fine_run_in_the_same_geometry.

The axisymmetric can

The $(r,z)$ problem: a ball of high pressure in a closed can (cylindrical/euler/axisymmetric_blast.toml, 128 by 192 cells, $t=0.4$), with total masses given by the volumes: $3\pi/2$ for the can. Relative drift of mass and energy over 1475 steps: 0 and 0 (to the last bit).

Tests that enforce this

Enforced by geometry/conservation.rs: a_closed_axisymmetric_can_conserves_mass_and_energy_with_both_methods, a_closed_cylinder_or_ball_conserves_mass_and_energy (cylinder: drift 9.5e-14 and 1.5e-13 over 2879 steps).

Swirl

With swirl = true the gas has an azimuthal velocity $w$ about the axis: the variables are rho, u (radial), v (axial), w and p, on a one- or two-axis cylindrical domain (spherical domains refuse it). The conserved quantity is the angular momentum $r\rho w$, the centrifugal force $\rho w^2/r$ acts on the radial momentum, and the symmetry axis reverses both $u$ and $w$; how the scheme does this is on the methods page. The angular momentum is reported as its own integral, $\int r\rho w\,dV$, because the plain integral of $\rho w$ is not constant.

Check	Result
$w=0$ against the gas without swirl	density, both momenta and energy equal to the last bit at $t=0.05$ on 24 by 24 cells (same number of steps); the swirl stays exactly 0
Closed can with a swirling ring, 32 by 32 cells, $t=0.15$	relative drift of the angular momentum 0, mass 0, energy 1.8e-16 (unsplit); all three 0 (Strang)
Solid-body rotation $w=\Omega r$, $p=p_0+\rho\Omega^2r^2/2$: an exact steady state	radial momentum, $L_1$ over $r<0.6$ at $t=0.2$, 32 to 256 cells: 2.71e-5, 2.70e-6, 3.17e-7, 4.35e-8; orders 3.33, 3.09, 2.86 (MUSCL)
The same, first order	2.40e-4, 8.57e-5, 4.19e-5, 2.14e-5; orders 1.48, 1.03, 0.97
The shipped swirling_blast, 96 by 144 cells, 847 steps	mass and energy drift 0; the angular momentum is 3.53416 against the exact initial $2\pi\cdot1.5\cdot\tfrac14\cdot1.5=3.53429$ (4e-5, from the cell-centre values) and drifts by 1.3e-16
Threads	1, 2 and 5 threads give bit-identical fields

Solid-body rotation is the check of the centrifugal term: the pressure gradient balances $\rho w^2/r$ only if the term is right, and without it the gas would not stay put. The profiles are polynomials of degree two, which the linear reconstruction nearly captures, so the error falls faster than second order (about 3); the test requires a fitted order above 1.8 and an error below 1e-6 on the finest grid. The exact steady state of the first-order scheme falls at first order, as the scheme does.

cylindrical/euler/swirling_blast.toml is a hot ball on the axis of a gas in solid-body rotation, $w=1.5r$, with the pressure that balances it. The blast expands, and the swirl inside it spins down because the same angular momentum is carried to larger radii:

Tests that enforce this

Enforced by geometry/swirl.rs: without_swirl_the_gas_with_three_velocities_is_the_gas_with_two, angular_momentum_mass_and_energy_are_conserved_to_round_off, a_swirling_run_is_identical_on_any_number_of_threads, solid_body_rotation_is_a_steady_state_to_at_least_second_order, solid_body_rotation_is_first_order_without_reconstruction; crates/merlin-cli/tests/run/curvilinear_problems.rs: the_swirling_blast_reports_and_conserves_its_angular_momentum (and the shipped file in closed_problems_conserve_the_volume_integrals_of_mass_and_energy); the integral itself by crates/merlin-core/src/blocks.rs: the_angular_momentum_integral_weights_by_the_radius; the file format by crates/merlin-cli/src/config/tests.rs: a_swirling_gas_has_an_azimuthal_velocity_and_an_angular_momentum, the_symmetry_axis_reverses_the_swirl_too_but_a_wall_does_not, swirl_needs_a_cylindrical_domain.

Axisymmetric ball against the 3D ball

A ball of radius 0.25 with $p=10$ in the corner of the octant, 24 cells per side, $t=0.1$, density along the axis: the axisymmetric run (an $(r,z)$ discretisation) against the 3D Cartesian run. They are two discretisations of one flow, so they cannot agree exactly; the relative $L_1$ difference is 1.46 %, and the test requires it to stay below 3 %.

Tests that enforce this

Enforced by geometry/blast.rs: an_axisymmetric_ball_matches_the_cartesian_ball_on_its_axis.

Orders at a glance

Error against refinement for the cylindrical cases, each series divided by its own coarsest error. The shock problems (Noh, Sedov) flatten towards first order or below; the smooth ones follow their design order.

Equivalence and independence

• A 2D flow without $z$ dependence equals the 1D cylinder bit for bit (density, radial momentum, energy; axial momentum exactly 0), on 24 by 10 cells (geometry/operator.rs: a_two_dimensional_flow_without_height_dependence_is_the_one_dimensional_cylinder).
• The right-hand side is independent of the block layout and the thread count, and the mass and energy sums equal the boundary flux (the_rhs_does_not_depend_on_the_block_layout_or_the_threads, the_rhs_of_the_conserved_mass_and_energy_sums_to_the_boundary_flux); whole runs on 1 to 8 threads are bit-identical (geometry/conservation.rs: a_blast_on_a_curvilinear_grid_is_identical_on_any_number_of_threads, the_answer_does_not_depend_on_the_block_layout; crates/merlin-cli/tests/run/curvilinear_problems.rs: a_curvilinear_run_is_identical_on_any_number_of_threads).
• A file run equals the library run bit for bit (a_file_run_matches_the_library_driven_run_bit_for_bit).
• The time step accounts for the small cells next to the axis, where the Courant number is up to 2 times the nominal one, and a strong outward flow at the axis keeps a positive density (geometry/conservation.rs: the_time_step_accounts_for_the_small_cells_next_to_the_axis, a_strong_outward_flow_at_the_axis_keeps_a_positive_density).
• A scheme that has its own time step (Lax-Wendroff) is refused on a curvilinear grid (a_scheme_with_its_own_time_step_is_refused_on_a_curvilinear_grid).

References

• 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.
• Morse, P. M. and Ingard, K. U., Theoretical Acoustics, Chapter 7 (the Hankel-transform solution).
• Toro, E. F. (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics, Chapter 1 (the Euler equations with geometric source terms).
• Sedov, L. I. (1959), Similarity and Dimensional Methods in Mechanics, Academic Press, Chapter IV; Landau, L. D. and Lifshitz, E. M., Fluid Mechanics, Section 106 (the point explosion).