Add the five-equation two-phase model by gouarin · Pull Request #13 · hpc-maths/samurai-euler

gouarin · 2026-09-13T18:39:20Z

Lot 4 of the plan: two compressible fluids sharing the mesh, in pressure and
velocity equilibrium, with an interface between them. It opens the cases of the
article a monofluid solver cannot reach, and two of them are here.

d_t (alpha_i rho_i) + div(alpha_i rho_i u) = 0      i = 0, 1
d_t (rho u)         + div(rho u @ u + p I) = 0
d_t E               + div((E + p) u)       = 0
d_t alpha_0         + u . grad(alpha_0)    = 0

dim + 4 components per cell, two new binaries (two_phase_1d,
two_phase_2d), and the monofluid solver untouched.

The model

Each phase is a stiffened gas; the mixture is one whose coefficients depend on
the volume fraction through two sums that are linear in alpha. That
linearity is the whole model: it is what lets a uniform pressure survive the
averaging of two fluids in one cell.

The four conservation laws go through the Riemann solvers the monofluid ones do.
The fifth does not: written as div(alpha u) - alpha div(u), it is discretized
with the contact velocity u* the Riemann solver already computes, and its
alpha_i term is the cell's own volume fraction — so the two cells an interface
separates receive different contributions. samurai calls that a non-conservative
flux and takes a pair of values per face, which is what the two-phase scheme
returns; the conservation laws take the conservative pair.

What it is held to

Check	Result
`alpha = 1` everywhere is the Euler system	Sod's tube run by both solvers, 3 Riemann solvers x 2 orders, agrees to 1e-12
Interface carried by a uniform flow	pressure uniform to 1e-11 relative, velocity to 1e-14
Periodic box	two masses, momentum and energy conserved to 1e-12
Water-air tube vs the exact solution	L1 on rho 1.5e-3 at level 10, contact within one cell of x = 1.1428

The exact solution is python/exact_two_phase_riemann.py, a two-material
stiffened-gas Riemann solver: p* = 4.796906e5 Pa, u* = 491.97 m/s,
rho*_water = 800.33, rho*_air = 2.758. The tube is measured against that
rather than against a figure read by eye.

The absent phase in the pure-fluid test is water, so anything leaking it into
the answer leaks a gigapascal.

Test cases

water_air_shock_tube (§6.1.1), triple_point (§6.1.2, with the two gases the
article gives rather than the single gamma the monofluid version had to settle
for), advected_interface (the interface condition), sod_x_pure (one fluid, to
compare against the monofluid solver).

The article gives neither the domain of the tube nor its diaphragm; its figure
gives both — [-2, 2] and 0.7 — and the exact solution confirms it, the contact
standing at 1.14 and the rarefaction spanning -1.69 to -0.49 at t_f. In 2D the
domain is [-2,2] x [-0.4,0.4], where level 6 is exactly the 320 x 64
equivalent resolution the article quotes; run there, the star density comes
out at 800.32 against 800.33 exact.

What the shared code had to give

Three changes, all in the direction of saying what was already true:

the slope limiter is templated on the number of components rather than on the
dimension — a limiter has no opinion on which component is which;
the solid wall no longer names the Euler layout: both models put the momentum
at 2 + d, each layout holds itself to it with a static_assert, and the wall
mirrors detail::momentum(normal);
a two-phase run writes the two partial densities alongside the mixture
density, alpha * rho being the mass of neither phase in a mixed cell.

THINC interface sharpening (second commit)

--thinc reconstructs the volume fraction inside a mixed cell as a hyperbolic
tangent instead of a limited straight line, with the interface position fixed in
closed form by conservation of the cell average. It replaces the reconstruction
of the volume fraction and of nothing else; the energy is rebuilt from the
sharpened alpha, and the mixture law being linear in alpha is what keeps the
interface condition exact.

Water-air tube, level 10, second order:

	interface	contact	L1 alpha	L1 rho
diffuse	9 cells	1.1426	1.8e-3	1.53e-3
`--thinc`	3 cells	1.1426	6.9e-4	1.54e-3

Three cells is the number the article quotes. The contact does not move and the
density and pressure errors do not change. On the triple point the mixed cells
fall from 1912 to 1147.

The steepness is weighted by the interface normal (beta_d = beta |n_d| + 0.001)
so an oblique interface is not carved into steps; the normal needs the gradient
in every direction, which a flux stencil cannot see, so it is computed once per
step over the mesh and read from a field.

Held to: the interface condition with the sharpening (1D and 2D, both
orders), conservation to 1e-12, alpha in [0, 1] through the shock tube, and a
run with no interface untouched to the bit — Sod's tube gives the same answer
with and without --thinc, which is what says the sharpening keys on the volume
fraction and not on steepness.

166 fast tests, 27 slow.

Shock-bubble interaction (fourth commit)

Section 6.1.3, the one case of the article validated against a laboratory rather
than against another code. It was listed here as out of reach, and that was about
the wrong thing: what a workstation cannot do is the article's 5120 x 1024 for a
millisecond, which its schlieren pictures need. The five wavefront velocities of
its table 3 are wave speeds, and they are inside the experimental margin at a
sixteenth of that resolution.

m/s	level 8	level 9	article	Haas & Sturtevant
shock	417.2	422.4 ± 1.7	423.2 ± 0.6	410 ± 41
refracted	946.6	951.6 ± 2.7	953 ± 7	900 ± 90
transmitted	382.3	381.6 ± 0.2	381.2 ± 0.7	393 ± 39
downstream	141.5	135.4 ± 1.8	141.5 ± 1.9	145 ± 15
jet	221.1	225.5 ± 1.4	222.9 ± 2.2	230 ± 23

Level 9 is a quarter of the article's resolution and fifteen minutes; three of
the five agree with it to two parts in a thousand. The shock reaches the bubble
at 58.8 µs against the "about 58" the article quotes.

python/shock_bubble_waves.py measures them as the article does, with one
simplification that takes an eye out of the loop: the incident shock, the
refracted wave and the transmitted wave are one measurement and not three, the
leading pressure front on the centreline being each of them in turn, so the
segments are separated at the edges of the bubble.

Two things about the initial state, both in the header. Equation (8) of the
article, taken literally, is not a shock — the pressure ratio is exactly Mach
1.22 but the density behind it is 1.6571 where Rankine-Hugoniot asks for 1.6295,
which would put the front at 401 m/s. The solver settles it in the first
microseconds and 422 comes out, which is why the article's own 58 µs and its
423 m/s agree with each other and not with 401. And the boundaries are named side
by side: samurai attaches boundary conditions rather than replacing them, so a
wall everywhere plus an inflow on one side would leave that side carrying both.

Not here yet

A positivity-preserving prediction for the adapted two-phase mesh, which leans on
the admissibility floor for now, and the multiple-bubble case of section 6.1.4 —
the one above at an equivalent 32768² over 3.3e5 time steps, where the model is
there and the machine is not.

Two compressible fluids sharing the mesh, in pressure and velocity equilibrium, with an interface between them: four conservation laws and one advection equation, dim + 4 components per cell against dim + 2 for the monofluid solver. d_t (alpha_i rho_i) + div(alpha_i rho_i u) = 0 i = 0, 1 d_t (rho u) + div(rho u @ u + p I) = 0 d_t E + div((E + p) u) = 0 d_t alpha_0 + u . grad(alpha_0) = 0 It opens the cases of the article that a monofluid solver cannot reach. Two of them are here: the water-air shock tube of section 6.1.1 and the triple point of section 6.1.2 as the article poses it, with two gases rather than the single gamma `triple_point_single_gamma` had to settle for. Each phase is a stiffened gas and the mixture is one whose coefficients depend on the volume fraction, through two sums that are linear in alpha: 1 / (gamma_m - 1) = sum_i alpha_i / (gamma_i - 1) gamma_m pi_m / (gamma_m - 1) = sum_i alpha_i gamma_i pi_i / (gamma_i - 1) That linearity is the whole model. It is what lets a uniform pressure survive the averaging of two fluids in one cell, and the interface test below measures exactly it. The four conservation laws go through the Riemann solvers the monofluid ones do, with the mixture sound speed and the partial densities carried through the contact. The fifth does not: written as div(alpha u) - alpha div(u), it is discretized with the contact velocity the Riemann solver already computes, and the term in alpha_i is the cell's own volume fraction, so the two cells an interface separates receive different contributions. samurai calls that a non-conservative flux and takes a pair of values per face, which is what the two-phase scheme returns. The conservation laws take the conservative pair and are conserved to the last bit. The volume fraction appears in the reconstruction and in the Hancock predictor as two more equations of the primitive system, and everything else -- the registry, the boundary conditions, the time integrators, the performance metrics, the limiter -- is the machinery already in the repository, used as it is. Three things had to give a little to make that possible, and all three are in the direction of saying what was already true: - the slope limiter is templated on the number of components rather than on the dimension, a limiter having no opinion on which component is which; - the solid wall boundary condition no longer names the Euler layout. Both models put the momentum at 2 + d, each layout holds itself to it with a static_assert, and the wall mirrors `detail::momentum(normal)`; - the output of a two-phase run carries the two partial densities alongside the mixture density, because alpha * rho is not the mass of either phase in a mixed cell and anything checking conservation needs the real ones. What the model is held to, from the sharpest statement down: - With alpha = 1 everywhere it IS the Euler system, and the two solvers must agree to round-off. Sod's tube run by both, at three Riemann solvers and two orders, agrees to twelve digits. The absent phase is water, so anything leaking it into the answer leaks a gigapascal. - An interface carried by a uniform flow leaves the pressure uniform to 1e-11 relative and the velocity to 1e-14. This is the property the model exists for; a conservative scheme on two materials fails it at any resolution. - A periodic box conserves the two masses, the momentum and the energy to 1e-12, which is what says the non-conservative term reaches the right cells with the right sign. - The water-air tube against its exact solution, which python/exact_two_phase_riemann.py computes for two stiffened gases: p* = 4.796906e5 Pa, u* = 491.97 m/s, rho*_water = 800.33, rho*_air = 2.758. At level 10 and second order the L1 error on the density is 1.5e-3, the contact lands within one cell of x = 1.1428 and the interface sits on about nine cells. The article gives neither the domain of the tube nor the position of its diaphragm; its figure gives both. The tube is [-2, 2] with the diaphragm at 0.7, and the exact solution confirms it: at t_f the contact stands at 1.14, the head of the rarefaction at -1.69 and its tail at -0.49, which is what the figure shows. In two dimensions the domain is [-2, 2] x [-0.4, 0.4], where level 6 is exactly the 320 x 64 equivalent resolution the article quotes; run there, the star density comes out at 800.32 against 800.33 exact. Not here yet, and named as such in the headers: THINC interface sharpening, which is what takes the interface from nine cells to two or three and which the article publishes both with and without; and a positivity-preserving prediction for the adapted two-phase mesh, which for now leans on the admissibility floor instead. One time loop serves both dimensions, instantiated by two three-line files, rather than one main per dimension as the monofluid solver has.

The diffuse interface of the five-equation model spreads over about ten cells and keeps spreading. The sharpening the article uses reconstructs the volume fraction inside a mixed cell as a hyperbolic tangent rather than as a limited straight line: alpha_i(X) = 1/2 [ 1 + tanh( beta (sigma X + x_c) ) ], X in [0, 1] sigma says which way the interface faces, beta how steep it is, and x_c where it sits inside the cell. Only x_c is unknown, and it is fixed by the one thing the reconstruction must not break: the profile has to average to the cell average it came from. Integrating the tanh and solving for tanh(beta x_c) gives it in closed form, so conservation is exact rather than iterated towards. It replaces the reconstruction of the volume fraction and of nothing else. Pressure, velocity and the two partial densities keep their MUSCL slopes, and the total energy is rebuilt from the sharpened alpha by prim2cons. That is what makes it safe: the mixture law is linear in alpha, so the pressure of a uniform flow survives whatever the volume fraction does, and the interface test still measures 1e-11 with the sharpening on. The predictor is handed the jump across the cell as the slope of alpha, which is the only place that slope appears. The steepness is weighted by the interface normal, beta_d = beta |n_d| + 0.001, so a face the interface runs parallel to is not sharpened at all; without it the scheme carves steps into an interface that crosses the mesh at an angle. The normal needs the gradient in every direction and a flux stencil is a line -- the same wall the transverse terms of the Hancock predictor run into -- so it is computed once per time step over the whole mesh, into a field the flux function reads through a pointer. Once per step and not once per stage: the interface moves by less than a cell in a step, and a normal one stage old is a smaller error than the one the sharpening is correcting. Measured on the water-air tube at level 10, second order: interface contact L1 alpha L1 rho diffuse 9 cells 1.1426 1.8e-3 1.53e-3 --thinc 3 cells 1.1426 6.9e-4 1.54e-3 Three cells is the number the article quotes. The contact does not move, the density and the pressure keep the errors they had against the exact solution, and the volume fraction is two and a half times more accurate. On the triple point the mixed cells fall from 1912 to 1147. Four things are asserted on, beyond the width itself. The interface condition survives the sharpening, in one and two dimensions and at both orders. So do the two masses, the momentum and the energy, to 1e-12. The volume fraction stays in [0, 1] through the shock tube. And a run with no interface at all is untouched to the bit: Sod's tube through the two-phase solver gives the same answer with and without --thinc, which is what says the sharpening keys on the volume fraction and not on steepness -- a shock is steep, and nothing here reaches for it. --thinc needs --order 2; there is no reconstruction to replace at first order.

The three water-air reference comparisons failed on the CI runner and passed locally, on one cell out of a hundred and twenty-eight, by 3e-12 relative against a tolerance of 1e-12. Nothing is wrong with either run. The tube spans a gigapascal against an atmosphere, and four hundred time steps of a nonlinear scheme over a jump of four decades turn the last bit of a sum -- contracted into a fused multiply-add here and not there -- into a difference a second architecture does not reproduce. The monofluid references hold 1e-12 across the same two machines because their fields are of order one and their dynamic range is small; this one cannot, and pretending otherwise would only mean regenerating it on whichever machine ran it last. So the tolerance says what two machines actually agree on: 1e-9 relative, three hundred times the observed difference and still nine digits of every cell. A perturbation of 1e-8 on a single cell -- smaller than any change of the scheme could be -- is still caught, which is what the number has to be worth.

…lation The three water-air reference comparisons failed on the CI runner and passed here. The previous commit read that as round-off and loosened the tolerance to 1e-9; the runner then failed by 3e-4, which is not round-off, and the loosening was the wrong answer to the wrong diagnosis. What it actually is, measured rather than guessed: the same source compiled with -ffp-contract=off -- the one thing that differs between an ARM machine that fuses a multiply and an add and an x86 runner that does not -- reproduces the runner's numbers to the digit, 3.1e-4 relative on the density and 2.3e-3 on the pressure. The amplifier is the problem, not the model. The rarefaction of the tube runs into a liquid at a gigapascal, and ahead of its analytic head the scheme leaves a foot where the density is 1000 minus something tiny, a number built entirely by cancellation; two hundred time steps turn one bit there into three parts in ten thousand. No tolerance both accepts that and asserts anything. So the field comparison moves to `sod_x_pure`, whose fields are of order one and which the same -ffp-contract experiment leaves inside atol 1e-10 and rtol 1e-12 -- the tolerance every other reference in the repository is held to, and now a tolerance that has been checked against the thing that broke rather than against this machine alone. It exercises the same kernel over a shock, a contact and a rarefaction, and the three solvers still disagree enough for the test that tells them apart. The water-air tube keeps a fast test, on the quantities that are not built by cancellation: the density of the plateau between the rarefaction and the contact, which the scheme reaches rather than computes its way through, and the position of the contact. Its exact-solution comparison is unchanged and stays in the slow suite, where the resolution is high enough for the errors to mean something. A uniform mesh is necessary for a field comparison and not sufficient. The readme says so now, next to the rule it qualifies.

Section 6.1.3 of the article, and the one case in it validated against a laboratory rather than against another code. A Mach 1.22 shock runs down a 445 x 89 mm tube into a 25 mm helium bubble; the bubble carries sound at three times the speed of the air around it, so the refracted wave outruns the incident shock, the bubble caves in on its upstream side and drives a jet through itself. Five fronts come out of it and Haas and Sturtevant measured their speeds in 1987, which is what table 3 of the article compares itself with. It was listed here as out of reach, and that was about the wrong thing. What is out of reach on a workstation is the article's 5120 x 1024 for a millisecond, which is what its schlieren pictures need. The five velocities do not: they are wave speeds, and they are inside the experimental margin at a sixteenth of that resolution. python/shock_bubble_waves.py measures them the way the article does -- a snapshot every ten microseconds, each front read along the axis of the tube, a straight line through its positions -- with one simplification that removes an eye from the loop. The incident shock, the refracted wave and the transmitted wave are not three measurements but one: the leading pressure front on the centreline IS the incident shock while it is right of the bubble, the refracted wave while it is inside it, and the transmitted wave once it is out. Its speed changes at each crossing, and the segments are separated at the edges of the bubble. level 8 level 9 article experiment shock 417.2 422.4±1.7 423.2±0.6 410±41 refracted 946.6 951.6±2.7 953±7 900±90 transmitted 382.3 381.6±0.2 381.2±0.7 393±39 downstream 141.5 135.4±1.8 141.5±1.9 145±15 jet 221.1 225.5±1.4 222.9±2.2 230±23 Level 9 is a quarter of the article's resolution and fifteen minutes; three of the five agree with it to two parts in a thousand, and all five are inside the 10% the experiment quotes. The shock reaches the bubble at 58.8 microseconds against the "about 58" the article gives, which is the geometry and the initial state in a single number. Two things are worth recording about the initial state, both in the header. Its states are equation (8) of the article and taken literally they are not a shock: the pressure ratio is exactly Mach 1.22 but the density behind it is 1.6571 where Rankine-Hugoniot asks for 1.6295, which would put the front at 401 m/s. The solver settles that in the first microseconds and 422 comes out, which is also why the article's own 58 microseconds and its 423 m/s agree with each other and not with 401. The numbers are left exactly as the article gives them. And the boundaries are named side by side rather than as a wall everywhere plus exceptions: samurai attaches boundary conditions instead of replacing them, so a side carrying both a wall and an inflow would be decided by the order they were attached in. `bc::imposed_state` is the generic half of `bc::imposed`, which took a monofluid primitive state and could not serve a two-phase one.

gouarin added 5 commits September 13, 2026 20:38

Navigation Menu

Uh oh!

Add the five-equation two-phase model - #13

Add the five-equation two-phase model#13
gouarin wants to merge 5 commits into
mainfrom
five-equation-model

gouarin commented Sep 13, 2026 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Uh oh!

Conversation

gouarin commented Sep 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

The model

What it is held to

Test cases

What the shared code had to give

THINC interface sharpening (second commit)

Shock-bubble interaction (fourth commit)

Not here yet

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

gouarin commented Sep 13, 2026 •

edited

Loading