Restore pT equilibrium when pTg relaxation uses up a phase

[aliencaocao]

Fixes #1820.

Change

In s_infinite_relaxation_k (src/common/m_phase_change.fpp), right after the pTg solve, the code checks whether either reacting mass is exactly zero. If it is, the existing pT solver is called again for the final masses. The call happens before the volume fractions and phasic energies are rebuilt:

call s_infinite_ptg_relaxation_k(j, k, l, pS, rhoe, q_cons_vf, TS)

! A depleted phase leaves Gibbs equality unconstrained; solve pT equilibrium again for the final masses
if (q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
    & l) == 0.0_wp .or. q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) == 0.0_wp) then
    call s_infinite_pt_relaxation_k(j, k, l, 2, pS, p_infpT, q_cons_vf, rhoe, TS)
end if

The test is exact because the pTg solver projects the liquid mass with min(max(ml, 0), mT). A depleted phase is therefore exactly 0. The vapor mass mT − ml is exactly 0 when the liquid mass reaches mT. States where both phases are present take the same path as in the original code.

Why this is correct

The pTg system solves Gibbs equality g_l(p, T) = g_v(p, T) together with the energy constraint for the liquid mass and pressure. The formulation follows Wilfong et al. (2025), Section 4.1.3 and MFC equations.md, Section 8.2. Gibbs equality is the condition for coexisting phases. When the projection removes one phase entirely, it cannot hold. The (p, T) where the Newton iteration stops is then not an equilibrium of the fluids that remain. Their equilibrium is the pT problem (Pelanti and Shyue 2014, Section 4, Step 3a; Wilfong et al. 2025, Section 4.1.3):

Σₖ mₖ / ρₖ(p, T) = 1,      Σₖ mₖ eₖ(p, T) = ρe

That is exactly what s_infinite_pt_relaxation_k solves (Section 8.1), with the masses that the pTg step leaves and the unchanged ρe. The rebuild of αₖ = mₖ/ρₖ and αₖρₖeₖ from that (p, T) therefore satisfies Σαₖ = 1 and Σαₖρₖeₖ = ρe by construction.

Before and after

The figure shows the air and vapor bubble from the issue, from the initial state that pre_process writes to 20 steps (1.8 ns):

• Before: 1,969 depleted cells sit at 103.45% filling with −0.43% energy closure. The 552 cells depleted at step 20 still sit at 103.45%.
• After: the filling error is exactly 0, and the energy closure is 8.5 × 10⁻¹⁶ in the initial state. After 20 steps, they are 4.6 × 10⁻¹¹ and 5.7 × 10⁻¹², within the tolerance of the pT solver.

Tests and goldens

The Phase Change model 6 cases and the phasechange_bubble example pass through depleted cells. Their goldens therefore change (for the example already at step 0, because pre_process relaxes the initial condition with the same routine). The goldens were regenerated on CPU with GCC 16.2 in double precision without case optimization:

• 02BD9B5E 2D -> Phase Change model 6 -> 3 Fluid(s) -> model equation -> 3
• 76C663B7 1D -> Phase Change model 6 -> 2 Fluid(s) -> model equation -> 2
• 8D6E051E 1D -> Phase Change model 6 -> 3 Fluid(s) -> model equation -> 3
• AAE5C27B 2D -> Phase Change model 6 -> 3 Fluid(s) -> model equation -> 2
• B97D4966 2D -> Example -> phasechange_bubble
• D21F4F38 1D -> Phase Change model 6 -> 2 Fluid(s) -> model equation -> 3
• D8136392 1D -> Phase Change model 6 -> 3 Fluid(s) -> model equation -> 2

Related: #1927 concerns when pTg relaxation is activated. This PR concerns what happens after it has run.

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

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[github-actions]

Claude Code Review

Head SHA: bc23e06

Files changed:

• 1
• src/common/m_phase_change.fpp

Findings:
No issues found. Checked src/common/m_phase_change.fpp (the new re-solve in s_infinite_relaxation_k against s_infinite_pt_relaxation_k and s_infinite_ptg_relaxation_k).

[codecov]

Codecov Report

❌ Patch coverage is 60.00000% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.89%. Comparing base (47568a0) to head (917c161).

Files with missing lines	Patch %	Lines
src/common/m_phase_change.fpp	60.00%	0 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2052      +/-   ##
==========================================
- Coverage   61.90%   61.89%   -0.01%     
==========================================
  Files          86       86              
  Lines       22885    22887       +2     
  Branches     3369     3370       +1     
==========================================
+ Hits        14166    14167       +1     
  Misses       6220     6220              
- Partials     2499     2500       +1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[github-actions]

Lines of Code

File	Lines	Diff
src/common/m_phase_change.fpp	289	+4

Directory	Lines	Diff
common	10530	+4
total	47589	+4

[wilfonba]

@aliencaocao I changed the triggering condition from floating-point comparisons to 0.0_wp to "if PTG failed to converge". The 0.0_wp comparisons were safer here than usual, given the min/max logic in s_infinite_ptg_relaxation, but exact comparisons to floating-point values are still never good. The new condition is less specific than the old one, but solving the pt relaxation whenever the ptg relaxation fails to converge catches the depleted-phase case and other cases where the ptg relaxation can't be trusted. All the golden files passed without regeneration. Can you rerun your reproducer and post-process the results to see if it gives the right answer there? Once that's done, this should be good to go. The correction matches the Pelanti reference as far as I can tell.

[aliencaocao]

@wilfonba I have ran the reproducer using your fix and can confirm it fixes the issue in the same way.

[sbryngelson]

is this not an optional?

[sbryngelson]

prove this gives the right results, what are the costs if it does? tradeoffs?

[sbryngelson]

your PR shows that the solution changes and via an energy argument is better, but i don't see validation that this is the correct solution it's converging to? what about transient processes as well? is it robust under those? surely there's some cost to this as well?

[aliencaocao]

I think it would help to set the context if I go through what the bug is first so you can see how I am arriving at the fix, before we get to why it is correct and its costs.

Every time step, each cell runs two solves in a row. First, a pT solve, given these masses and this much internal energy, find the one pressure and temperature where the fluids exactly fill the cell and the energies add up. Then, if both reacting phases are present, a pTg solve, move mass between liquid and vapor until both would rather coexist. (that's the saturation condition).

Sometimes the right answer to the second solve is to convert all of it. The liquid fully evaporates. But the equation the pTg solver is trying to satisfy (liquid and vapor in equilibrium) is meaningless once one of them is gone. So the Newton just runs out of iterations

MFC/src/common/m_phase_change.fpp

Lines 369 to 372 in 47568a0

	do ns = 1, max_iter
	! converged on the absolute residual, or on the rhoe-relative residual (multiply form, rhoe > 0)
	if ((resnorm <= ptgalpha_eps) .or. (resnorm <= (ptgalpha_eps/1.e6_wp)*rhoe)) exit

and stops at some arbitrary pressure and temperature. The old code then built the cell from that leftover (p, T), giving wrong numeric results.

My fix is to ignore the leftover (p, T) and redo the first pT solve with the final masses. That problem always has an answer, and it just reuses an existing solver (relax_model = 5).

Why that is correct:
Since the pT solver only looks at (masses, energy) and ignores the incoming pressure, it gives the same (and correct) answer no matter where the failed solve happened to stop. This is also what the model's references prescribe for exactly this case: Pelanti & Shyue (2014, §4) run the (p,T) relaxation at the fixed masses after the Gibbs step (Step 3a), and the MFC paper Wilfong et al. (2025, §4.1.3) follows the same structure: once the Gibbs step empties a phase, the pT problem is the equilibrium of the fluids that remain. Thus if we trust that the existing pT solver is correct, then this must arrive at a correct solution.

To validate the solution independently of MFC's solvers, I took the output files and checked whether they indeed satisfy the equations of the pT problem. On D21F4F38 step 50, the two depleted cells are pure single phase, so the equilibrium has a closed form, p = (n−1)(ρe − m·q) − n·B. It matches the pressure in MFC's output to floating-point accuracy, with α = 1 exactly.

Cost:
This indeed introduces an extra, cheap pT solve on those cells, but that is the cost to be paid for correctness, otherwise the results are wrong to begin with and there is no meaning in discussing cost. The failing cells already run the full pTg Newton (2×2 solve with line search), which dominates the relaxation cost, and the added scalar solve is negligible compared to it. I actually have an optimisation for that dominant part and will PR it separately after this lands. But for this PR, I believe this is the lowest-cost correction we can do.

On transients: D21F4F38's depletion happens during the run, not just at the initial condition, and in the bubble reproducer 552 cells deplete between step 0 and step 20 and still close to 4.6e-11 / 5.7e-12. The re-solve only recomputes (p, T) so it adds no new feedback path into the time stepping. Thus it remains stable mid-run. If anything, the transients were the problem that the PR is trying to fix.

For the trade off we have 3 cases:

1. pTg converges: identical with current master branch
2. pTg fails, phase depleted (the reproducer I used): wrong results currently. The extra pT solve is necessary for correctness.
3. pTg fails but both phase still present (hypothetical case, I dont see any of that in the existing goldens): my original ==0.0 ignores this but @wilfonba 's trigger handles this. The latter is strictly better than the current code, which does nothing in this case and leaves the wrong result there.

[sbryngelson]

seems convincing