Fix the 6-equation phase energies at characteristic boundaries and the wall contact angle with surface tension

[wilfonba]

Description

A liquid drop touching a characteristic wall (bc = -5) blew up, and a drop touching a ghost-cell no-slip wall (-16) came out with a Laplace jump 39 % too low. The cause is two independent defects. The first also affects characteristic boundaries that are not walls.

1. The 6-equation model's phase internal energies did not follow the characteristic boundary update. In s_cbc, the boundary flux is rewritten only for equations up to adv%end, so the phase internal energies alpha_k rho_k e_k kept their Riemann flux built from extrapolated ghost cells. Their non-conservative source -alpha_k p div(u) uses the face interface velocity, which s_cbc sets at the boundary face equal to the interior face's, so it sees no boundary-normal compression in the boundary cell. The mixture energy, meanwhile, gets the characteristic rate, which does. Every characteristic boundary cell therefore drifted out of sum_k alpha_k rho_k e_k = E - KE each step. Pressure relaxation absorbs that where a cell holds one phase or a uniform mixture; where an interface crosses the boundary, a trace phase is driven to a negative density ("Pressure relaxation produced a non-physical phasic density"). This affects any 6-equation characteristic boundary that an interface crosses, not only walls: an interface leaving through a -8 outflow made an 18 % pressure error (gas/gas) or stopped the run (gas/liquid).

2. Walls did not give the color function a 90 degree contact angle. The color function c and its gradient c_divs were copied into the ghost cells of characteristic and ghost-cell walls (zeroth-order extrapolation). The central-difference gradient in the boundary cell then has a wall-normal component, so the capillary stress sigma n1 n2 acts along the wall at every contact line. Symmetry planes (-2) already reflected both, which is why they were right. Where the mirror goes depends on where the boundary condition acts: a ghost-cell wall lies on the domain face, while a characteristic boundary imposes its conditions on the boundary cell's state (the characteristic amplitudes are one-sided derivatives at that cell's center, and the face flux is chosen so that cell evolves by them).

Changes

• m_cbc.fpp: for the 6-equation model the characteristic update also covers the phase internal energies (flux_cbc_index = int_en%end). Each gets d/dt [alpha_k (Gamma_k p + Pi_k) + alpha_k rho_k qv_k] from the same dp/dt, d alpha_k/dt and d(alpha_k rho_k)/dt as the mixture energy. Gamma_k, Pi_k and their density derivatives come from s_phase_coefficients, so every equation of state is covered (constant for ideal and stiffened gases, state dependent for Mie-Gruneisen, JWL and Vinet), exactly as the existing mixture-energy update already does through s_compute_mixture_coefficients_dt. The boundary face's interface velocity stays the interior one, so -alpha_k p div(u) adds nothing there.
• m_variables_conversion.fpp: s_convert_primitive_to_flux_variables (used only by s_cbc's boundary extrapolation) now fills the phase internal energy fluxes for the 6-equation model, through s_phase_internal_energy.
• m_boundary_primitives.fpp: s_color_function_reflective now loops over the gradient components with one routine, s_mirror_into_ghost_cells, negating the normal one, and mirrors either about the domain face (off = 1, unchanged for -2) or about the boundary cell's center (off = 0).
• m_boundary_common.fpp: with surface tension, the color function and its gradient are mirrored at -5 walls (about the boundary cell's center) and at -15 / -16 walls (about the face). Other characteristic boundaries keep extrapolation: an interface passes through them.
• cases.py: seven new regression tests (below).

The 5-equation model's characteristic boundaries are unchanged, bit for bit, as are symmetry planes.

Verification

Every case was run on CPU (macOS, gfortran 15, Open MPI) with this branch ("after") and with its base 3dc5b2f2 ("before"), from the same inputs. Unless stated otherwise: WENO5 (mapped), HLLC, third-order Runge-Kutta, double precision.

An interface leaving through a characteristic outflow (fix 1, away from walls)

Setup. A uniform flow, p = 1 and u = U = 0.5, on x in [0, 1] with 200 cells: fluid 1 (ideal gas, gamma = 1.4, rho = 1) for x < 0.6 and fluid 2 for x > 0.6. The left boundary extrapolates (-3); the right one is the non-reflecting characteristic outflow (-8), which the interface reaches at t = 0.8. mpp_lim on, dt = 0.2 dx / (c_max + U), run to t = 1.2. The exact solution is the uniform state, so any pressure or velocity disturbance is made by the outflow. The 5-equation model, which this PR does not change, is the reference. Fluid 2 is

• gas: ideal gas, gamma = 5/3, rho = 0.5;
• stiffened liquid: gamma = 4.4, pi_inf = 6, rho = 1000;
• Mie-Gruneisen: rho0 = 0.9, c0 = 1, s = 1.5, Gamma0 = 0.4, at rho = 1.1, away from its reference density so the density-dependent Gamma(rho) and Pi(rho) matter;
• JWL: A = 6, B = 0.15, R1 = 4, R2 = 1, omega = 0.3, rho0 = 0.9, at rho = 1.1.

The Mie-Gruneisen and JWL cases run on a strip 25 cells high, periodic in y (the same 1D problem), because MFC writes the 6-equation pressure only in 2D and 3D. For the other two the pressure is computed from the conservative variables.

Fluid 2		6-eq before	6-eq after	5-eq
gas (rho 2:1)	max |p − 1|	1.8e-1	8.6e-5	2.4e-5
	max |u − U|	1.4e-1	9.2e-5	2.3e-5
stiffened liquid (rho 1:1000)	max |p − 1|	stops at t = 0.8	4.4e-3	1.4e-4
Mie-Gruneisen	max |p − 1|	2.7e-1	6.4e-4	1.8e-4
	max |u − U|	2.6e-1	4.3e-4	9.7e-5
JWL	max |p − 1|	1.7e-1	3.1e-4	5.4e-5
	max |u − U|	1.6e-1	2.5e-4	3.9e-5

On the base code the gas/gas outflow sends a 17 % pressure jump back into the domain as soon as the interface reaches it; after the fix the pressure stays flat (the remaining 1e-4 is not visible at this scale):

With periodic boundaries instead of the outflow (no characteristic boundary), the 6-equation stiffened-liquid case keeps the pressure uniform to 5e-10, so what remains after the fix is the boundary's (see Known open issues).

A drop on one wall: two contact lines (fixes 1 and 2)

Setup. examples/2D_laplace_pressure_jump: a liquid drop (stiffened gas, gamma = 2.1, pi_inf = 1e6, rho = 1000) of radius R = 0.15 in an ideal gas (gamma = 1.4, rho = 1), sigma = 8, so sigma/R = 53.3 Pa; p = 1e5 initially in both, no velocity, inviscid. 50 x 50 cells on [0, 0.375]^2, dt = 5e-6, 100 000 steps (t = 0.5), 6- and 5-equation models. The drop is centered on the x-beg boundary at mid height, so it meets that wall in two contact lines; the other boundaries extrapolate. With a symmetry plane (-2) on x-beg this is half of a free drop, so every wall with a 90 degree contact angle must reproduce it: the Laplace jump sigma/R and spurious currents that decay. The drop is centered on the wall, which for the characteristic wall means on the boundary cells' centers. The jump is the mean pressure within R/2 of the drop's center minus the mean beyond 3R/2.

Pressure at the end of each 6-equation run (or at its last save, if it stopped):

Δp / (σ/R) at t = 0.5, spurious max |u|	6-eq before	6-eq after	5-eq before	5-eq after
symmetry plane (-2), reference		0.9990, 0.031		0.9994, 0.035
characteristic wall (-5)	stops before t = 0.05	1.0077, 0.029	0.89–0.93, spurious |u| growing from 0.04 to 0.21; stops after t = 0.45	1.0112, 0.028
no-slip wall (-16)	0.868, 0.083	0.9991, 0.031	0.869, 0.089	0.9994, 0.034

The -16 wall now reproduces the symmetry plane. The characteristic wall gives the right jump and decaying spurious currents and no longer blows up, but it is not yet steady: the wall cell at each contact line carries a pressure excess that grows through the run (6-equation: 53 Pa at t = 0.1, 77 Pa at 0.3, 120 Pa at 0.5, against 53.5 Pa in the drop; 5-equation similar), visible as the two bright cells on the wall in the field above. On the symmetry plane and the -16 wall the same cells stay at 53.8 Pa. See

Known open issues.

A drop in a wall corner (fix 2 for -16; an open case for -5)

Setup. The same drop, centered on the corner of the x-beg and y-beg boundaries, which are both walls of the type under test; the corner lies inside the drop. With symmetry planes it is a quarter of a free drop.

Δp / (σ/R) at t = 0.5	6-eq before	6-eq after	5-eq before	5-eq after
symmetry planes (-2), reference	0.9992	(same)	0.9995	(same)
no-slip walls (-16)	0.61	0.9992	0.61	0.9995
characteristic walls (-5)	stops before t = 0.05	stops before t = 0.45	stops before t = 0.35	stops at t = 0.48

The -16 drop now matches the symmetry planes to 3e-5. With two characteristic walls meeting inside the drop the runs last much longer but still fail; see Known open issues.

Tests

• Seven new tests: model_eqns=3 -> bc=-8 in 1D, 2D and 3D (A70E837F, B77B013C, ECD339A3) pin the 6-equation characteristic update; capillary=T -> model_eqns=3 -> bc=-5 and bc=-16 in 2D and 3D (048B5618, A414808B, 9EC4E2FA, D8E467A9) pin contact lines on both wall types. No existing golden changed.
• Full suite on macOS: 726 passed, 0 failed, 39 skipped.

Known open issues

Neither is introduced by this PR; both cases failed sooner, or outright, before it.

1. Contact lines on characteristic walls are not yet steady. On a single characteristic wall the contact-line wall cell's pressure excess grows through the run (one-wall drop above) without stopping it by t = 0.5. Where two characteristic walls meet inside the drop, the spurious capillary currents also feed a tangential velocity mode along the walls next to the corner that grows by roughly tenfold every 0.05 in time and stops the run (corner drop above). Both need the capillary forcing: the corner drop with sigma = 0 stays at rest to 1e-11. Contact lines on the ghost-cell wall (-16) and on symmetry planes show neither.
2. A 6-equation characteristic outflow emits a small acoustic wave when a liquid interface leaves through it, about 1e-3 of the pressure for a 1000:1 interface (5-equation: 4e-6), traveling back into the domain (top-right panel of the outflow figure). Before this PR the same run stopped as the interface reached the boundary.

AI disclosure

This PR was prepared with Claude Code (Anthropic), which carried out the diagnosis, the code changes, the tests and the verification runs above.

---

Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

---

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

PR template credit: junegunn

🤖 Generated with Claude Code

[sbryngelson]

a github (non-self-hosted) runner update caused an issue that i fixed (for your debug jobs). syncing that pr from master should fix your failing tests

[github-actions]

Lines of Code

File	Lines	Diff
src/common/m_boundary_primitives.fpp	1439	-33
src/simulation/m_cbc.fpp	1132	+18
src/common/m_variables_conversion.fpp	970	+10
src/common/m_boundary_common.fpp	455	+7

Directory	Lines	Diff
common	10510	-16
simulation	28544	+18
total	47587	+2