Second order: MUSCL reconstruction, Hancock and SSP-RK2 - #10
Merged
Merged
Conversation
The floor has been applied at every step of the 3D run since it was written, and nobody knows whether it ever fires. It now returns the number of cells it changed and the run prints the total, which answers the question the code has been carrying a comment about: a floor that never fires can go, one that fires away from level interfaces is a statement about the scheme above it, and one that fires at all has moved mass and energy, the projection not being conservative.
The three solvers were each written inside their own scheme builder, and the three builders differed only by the twenty lines in the middle. A second-order scheme needs those same twenty lines fed with reconstructed face values instead of cell averages, which under the old shape meant three more copies. They are now free functions taking the two states an interface separates, nothing more: riemann::rusanov, riemann::hll, riemann::hllc. The Godunov scheme that hands them the two cell averages is written once, in schemes/godunov.hpp, and templated on which solver to call, so the three still produce one scheme type and the dispatch still returns it. Each solver takes its states twice, as primitives and as the matching conservative vector: the wave speeds come from one, the jumps from the other, and the caller has both in hand. No result moves: the expressions are the same, in the same order, and the reference fields of all 17 cases are unchanged.
The solver was first order in space and in time: cell averages handed straight to the Riemann solver, explicit Euler on top. Two options now open it up. --order 2 reconstructs the values the solution takes on each face, from a limited slope per cell, and hands the solver those instead. Four cells per interface rather than two. The limiter is --slope-limiter, minmod, vanleer or moncen, moncen by default, with none available for measuring an order on a smooth solution, where clipping an extremum costs a fraction of an order and says nothing. --time-integrator picks how the step is taken. `euler` with a reconstructed flux is MUSCL-Hancock: the face values are advanced half a step inside the flux function, which makes one explicit step second order for one flux evaluation. `ssprk2` is Heun, two stages. `auto`, the default, means Euler at order 1 and SSP-RK2 at order 2, which is the pairing that reaches second order in any dimension. Measured: on the advected pulse, in one dimension, MUSCL-Hancock converges at 2.01 and SSP-RK2 at 2.00, Hancock being twice as accurate for half the work; on the vortex, in two dimensions, SSP-RK2 holds 2.13 on a uniform mesh and 2.04 on an adapted one, and Hancock drops to 1.05. The drop is the transverse terms the predictor cannot see: it advances a face value with the equations normal to that face, and a flux stencil is a line. At CFL 0.05 it climbs back to 1.77, which is how we know the reconstruction itself is sound. Boundary conditions follow: Imposed, Reflective and a new Outflow are written for any even stencil width, and bc::ghost_layers, set once from the order, says whether they fill one layer of ghost cells or the two a reconstruction reads. The free stream stays exact to the last bit at order 2, including across level jumps, which is what says the wide conditions and the prediction operator agree with the wider stencil. Order 1 is untouched: same default, same integrator, same reference fields.
A bump of density carried by a uniform flow at uniform pressure. With p and u constant the system collapses onto rho_t + u rho_x = 0, so the exact solution is the initial profile translated, undeformed, on a contact wave. It fills the hole the vortex leaves. Measuring an order needs a solution that is smooth and known, the vortex is the only one the repository had, and it is two-dimensional by nature. One dimension is exactly where the Hancock predictor is expected to be second order, and where it can be told apart from a reconstruction that is merely plausible. The pulse is narrow enough and starts far enough from the boundaries that the imposed ambient state at either end is the exact solution there to better than 1e-12 for the whole run. Registered in one, two and three dimensions, and pinned in the non-regression suite in one.
The invariants now run at order 1 and order 2. The second is where they bite: a reconstruction reads a second layer of ghost cells and takes slopes across level jumps, which is two more ways to break a uniform state. The free stream on a mesh pinned by --refine-boundary is the sharp one, and it comes back at exactly zero. Four validation runs keep the order once it is gained: the vortex at order 2 on a uniform mesh and on an adapted one with epsilon scaled as h^2, both above 1.9; the pulse in one dimension for each integrator, which is the test that says the Hancock predictor is right; and the same pulse under the default limiter, held above 1.8, so that a limiter which has stopped limiting and one which clips everything both fail. 80 fast tests in 36 s, 14 slow ones in 31 s.
The bottom boundary was a value function: samurai calls one of those once per stencil, with the cell just inside the boundary, and the value it returns is repeated into every ghost layer. For the half of that boundary which is an inflow the value is a constant and repeating it is right. For the half which is a wall it is not: a wall is a mirror, the k-th ghost reflects the k-th cell inside, and at order 2 the second layer was carrying the reflection of the first cell instead of the second. The MUSCL slope in the first ghost is then taken over (w, m(w), m(w)) and comes out at zero, so the reconstruction at the wall drops back to first order and loses the symmetry the wall is made of, right where the jet runs along it. The condition is now a condition rather than a value function, which gives it every cell of the stencil and lets it fill each layer with what belongs there. At order 1 nothing moves, to the bit: the three double Mach references are unchanged. At order 2, on a 1024 x 256 mesh at t = 0.2, the density differs by up to 7.2 from what the old boundary produced, all of it within five rows of the wall, and the peak in the jet goes from 22.0 to 23.1.
gouarin
force-pushed
the
second-order-muscl
branch
from
September 13, 2026 06:05
0cae2d3 to
a457cdb
Compare
The Hancock predictor is complete in one dimension and incomplete in more,
where the transverse terms are out of reach of a line stencil. SSP-RK2 works
around that by adding a stage. Strang removes the difficulty instead: sweep one
direction at a time, and every sweep is a one-dimensional problem for which the
normal-only predictor is exactly right.
`--time-integrator strang` runs X(dt/2) Y(dt) X(dt/2) in two dimensions, the
symmetric sequence that cancels the splitting error to second order, and
X(dt/2) Y(dt/2) Z(dt) Y(dt/2) X(dt/2) in three. A sweep is the same scheme with
a flux function on one direction only; samurai skips a direction whose function
is not set, so a sweep costs one pass and not dim of them.
Measured on the isentropic vortex at t = 0.2, uniform mesh, unlimited slopes:
integrator L1 order L1 at level 7
ssprk2 2.13 1.41e-5
strang 2.23 8.00e-6
so Strang is second order and 1.8 times more accurate at every resolution
measured. On an adapted mesh with epsilon scaled as h^2 it holds 2.13. Counted
in flux passes it is the cheaper of the two, 2*dim - 1 against 2*dim, but each
sweep pays for a full ghost update: measured wall time comes out even in two
dimensions and about 15% against it in three.
In one dimension it reduces to the single Hancock step, and the test asserts
that the two agree bit for bit rather than merely closely.
The free stream stays exact under it, in two and three dimensions, and the
closed box conserves to round-off, so the sweeps leave the invariants where
they were.
Of the integrators that reach second order in more than one dimension, it is the accurate one: 2.23 against 2.13 on the vortex, and an L1 error 1.8 times smaller at every resolution measured, for the same wall time in two dimensions. Order 1 keeps explicit Euler, so nothing that was measured before moves. SSP-RK2 stays a named option, and it is the one to reach for as a counter-check, owing nothing to a predictor or to the order of the sweeps.
…a shock Three things the review asked for before merge. The readme listed --scheme and --cfl and knew nothing of --order, --slope-limiter or --time-integrator. It now carries the three, and a section on what the integrators are not interchangeable about, since at order 2 in more than one dimension the choice decides whether the scheme is second order at all. Two comments in the test suite still told a reader to tighten the thresholds "when the MUSCL-Hancock update lands", which is this branch. Rewritten to say what holds now. The gap they were pointing at was real: the order was measured on smooth solutions and with no limiter, so a limiter gone wrong at a discontinuity would have passed the whole suite. Both shock tubes now run at both orders, each against its own thresholds, and the second order earns them: on the 123 problem the L1 error on the density falls from 3.5e-3 to 3.9e-4 and the near-vacuum pressure from 1.27 times the exact star value to 1.004; on Sod from 6.7e-3 to 1.6e-3. That measurement immediately shows something the gap was hiding: the transverse velocity of the rotated Sod, which is what the rotation was introduced to measure, gets four times WORSE at order 2, 2.6e-3 against 9.6e-3. A reconstruction taken direction by direction has no reason to treat a discontinuity at 45 degrees as well as one along an axis. Every other error on that case falls by four or more. The thresholds record both. Three of the latent risks are answered here too, none of which moves a result: an assert on the time step the Hancock predictor reads, a note on get_max_lambda that the bound it returns is the one-dimensional one and that an unsplit integrator in 3D at --cfl 0.4 is already above its stability limit, and a note that a time-dependent boundary is read at t^n by every stage and every sweep.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Second order in space and in time, behind three options. Order 1 is untouched:
same default, same results, same reference fields.
At order 2 each cell gets a limited slope and the Riemann solver is handed the
two face values, on a stencil of four.
euleralso advances those values half astep inside the flux, which is MUSCL-Hancock;
ssprk2is Heun;strangsweepsone direction at a time.
autois Euler at order 1, Strang at order 2.The orders, measured
The three are not interchangeable. The Hancock predictor advances a face value
with the equations taken normal to that face, the transverse terms it drops are
of the same order as the ones it keeps, and a line stencil cannot reach them.
Three measurements say so rather than assume it: at CFL 0.05 the order climbs
back to 1.77, in one dimension it reaches 2.01, and the other two integrators
reach 2.1 on the same mesh. Strang removes the difficulty instead of working
around it, a sweep having no transverse direction, and is 1.8 times more
accurate than SSP-RK2 at equal resolution for the same wall time in 2D.
What came with it
advected_pulse, a smooth one-dimensional solution: a density bump on acontact wave, exact. The vortex is two-dimensional by nature, so there was no
way to measure an order in one dimension, which is where a correct Hancock
predictor can be told from a plausible one.
Boundary conditions of the right width.
Imposed,Reflectiveand a newOutflowfor any even stencil,bc::ghost_layerspicking one layer or two fromthe order. Two are needed whichever integrator drives the scheme, since a sweep
along x reads the same four cells along x. Dropping the second costs nothing
where the solution is flat against the boundary and a great deal where it is
not: on the double Mach it moves the density by 13 within five rows of the wall.
A wall that survives the wider stencil. The double Mach bottom was a value
function, and samurai repeats what one of those returns into every ghost layer.
Right for the inflow half, wrong for the wall half: the second layer carried the
reflection of the first cell, the MUSCL slope in the ghost came out at zero, and
the wall fell back to first order under the jet. Order 1 unchanged to the bit,
order 2 moved by up to 7.2 in density near the wall.
The Riemann solvers, extracted into free functions of the two states an
interface separates, called by the Godunov scheme and the MUSCL one alike, and
a count on the positivity floor.
Tests
81 fast tests in 37 s, invariants at both orders. Both shock tubes run at both
orders now, each against its own thresholds, so that a limiter gone wrong at a
discontinuity cannot pass on the strength of the smooth measurements, which use
no limiter at all. The second order earns them, L1 on the density of the 123
problem falling from 3.5e-3 to 3.9e-4, and makes one number worse: the
transverse velocity of the rotated Sod, the isotropy measure the rotation exists
for, goes from 2.6e-3 to 9.6e-3. A reconstruction taken direction by direction
has no reason to treat a discontinuity at 45 degrees as well as one along an
axis.
Five slow runs keep the order: the vortex under each multi-dimensional
integrator and on an adapted mesh, the pulse under each one-dimensional
integrator, and the pulse under the default limiter held above 1.8.
Noticed, not fixed
The 3D Sedov at level 4 simulates nothing: the cell is 0.125 wide and the
nearest cell centre sits at r = 0.108 against
r_blast = 0.1, so threereferences and a positivity test are pinning a uniform state. Level 5 fixes it.
Three risks are documented where they live, each of which would move results.
limit_positivityfloors the state after a whole step and not after each stageor sweep, so a near-vacuum can reach the next one unfloored.
get_max_lambdareturns the one-dimensional wave speed bound, which an unsplit integrator in
three dimensions at
--cfl 0.4is already above. A time-dependent boundary isread at
t^nby every stage and every sweep.The
NOTEinmain_2d.cppabout the two different predictions atinitialisation stands; reconciling it moves 2D results.