fix(phase_change): globalize pTg/TSat Newton solvers so they converge (esp. GPU)#1654
Draft
sbryngelson wants to merge 4 commits into
Draft
fix(phase_change): globalize pTg/TSat Newton solvers so they converge (esp. GPU)#1654sbryngelson wants to merge 4 commits into
sbryngelson wants to merge 4 commits into
Conversation
…de#1646) The three Newton solvers in m_phase_change (s_infinite_pt_relaxation_k, s_infinite_ptg_relaxation_k, s_TSat) looped on a convergence-only condition with no iteration limit. The module parameter max_iter was declared but never referenced (dead since MFlowCode#179), so a cell that fails to converge hung the solver forever (on GPU: kernel never returns, host blocks in cuStreamSynchronize). Reproduced with the stock examples/2D_phasechange_bubble. Bound each loop by max_iter and accept the last iterate on exhaustion, turning an unbounded hang into a bounded, finite step. Reduce the misleading real-literal init (max_iter = 1e8_wp) to a clean integer (100000). Golden-safe: all 18 phase-change regression tests pass byte-identical, confirming the cap sits above the legitimate convergence ceiling and truncates no converging cell.
…ally converge (MFlowCode#1646) Beyond bounding the loops, the pTg-equilibrium and saturation-temperature Newton solvers never reached their tolerance for low-pressure cells: a fixed 1e-3 underrelaxation stalled the iterates far from the root (Gibbs residual frozen ~1e3). This reworks both to converge robustly with bounded, uniform iteration counts (no GPU warp divergence). s_TSat: reformulated onto the dimensionless, strictly-monotone saturation relation G(T)=A+B/T+C*ln(T)+D*ln(p+p_inf,l)-ln(p+p_inf,v)=0 (the A/B/C/D coefficients computed in s_initialize_phasechange_module were previously dead), solved with a safeguarded Newton (rtsafe: Newton step when it stays bracketed and makes progress, bisection otherwise). Cannot diverge, needs no accurate guess. s_infinite_ptg_relaxation_k: damped Newton with backtracking line search on the same analytic 2x2 Jacobian, every step projected onto the physical bounds 0<=ml<=mT and pS>pmin, with a singular-Jacobian guard. Residual computation extracted into a device helper s_compute_ptg_residual so trial states can be evaluated without mutating q_cons. Also removed dead locals (p_infpTg, TJac) and the now-unused p_infpT argument, and rewrote the pT relative-convergence test in multiply form to avoid dividing by pO (pO=0 on the first pass). Validation: bubble cells now converge in <=25 (pTg) / <=60 (TSat) iterations instead of stalling. All 18 phase-change golden tests pass on CPU (gfortran) and GPU (nvfortran/acc, H200) - CPU/GPU agree to tolerance. The 6 model-6 goldens shift by ~1e-6 relative (same equilibrium root, consistent non-lagged auxiliaries); model-5/pT goldens unchanged. The max_iter cap from the parent commit is retained as a backstop.
2 tasks
Contributor
There was a problem hiding this comment.
Pull request overview
Note
Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.
This PR addresses phase-change Newton non-convergence (especially problematic on GPU) by replacing fixed underrelaxation with safeguarded Newton methods and by updating phase-change golden artifacts to reflect the new equilibrium solves.
Changes:
- Reworks
s_TSatinto a bracketed/safeguarded Newton root find (rtsafe-style) using precomputed saturation-curve coefficients. - Replaces fixed-underrelaxed pTg Newton with a damped Newton + backtracking line search and introduces a residual helper for trial evaluations.
- Regenerates phase-change golden metadata files for the affected model-6 tests.
Reviewed changes
Copilot reviewed 7 out of 13 changed files in this pull request and generated 7 comments.
Show a summary per file
| File | Description |
|---|---|
| src/common/m_phase_change.fpp | Implements new safeguarded TSat solver and damped/line-searched pTg Newton; adds iteration/backstop parameters. |
| tests/D8136392/golden-metadata.txt | Regenerated golden metadata reflecting updated equilibrium behavior/environment. |
| tests/D21F4F38/golden-metadata.txt | Regenerated golden metadata reflecting updated equilibrium behavior/environment. |
| tests/AAE5C27B/golden-metadata.txt | Regenerated golden metadata reflecting updated equilibrium behavior/environment. |
| tests/8D6E051E/golden-metadata.txt | Regenerated golden metadata reflecting updated equilibrium behavior/environment. |
| tests/76C663B7/golden-metadata.txt | Regenerated golden metadata reflecting updated equilibrium behavior/environment. |
| tests/02BD9B5E/golden-metadata.txt | Regenerated golden metadata reflecting updated equilibrium behavior/environment. |
Comment on lines
+565
to
+573
| ! log(p_sat + p_inf,v) is undefined for p_sat + p_inf,v <= 0; no valid saturation temperature | ||
| ! exists there, so leave the incoming temperature untouched. | ||
| if (pSat + ps_inf(vp) <= 0.0_wp) then | ||
| TSat = TSIn | ||
| return | ||
| end if | ||
|
|
||
| ! calculating the jacobian | ||
| dFdT = -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*log(TSat) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) - 1) & | ||
| & *log(pSat + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(pSat + ps_inf(vp)) | ||
| ! Pressure-dependent part of G(T), constant across the iteration. | ||
| pterm = D*log(pSat + ps_inf(lp)) - log(pSat + ps_inf(vp)) |
Comment on lines
+432
to
+435
| ! pressure floor (stiffened gas requires pS + ps_inf > 0 for both phases) | ||
| pmin = -min(ps_inf(lp), ps_inf(vp)) + 1.0_wp | ||
|
|
||
| ! calculating the (2D) Jacobian Matrix used in the solution of the pTg-quilibrium model | ||
| call s_compute_ptg_residual(ml, mT, pS, j, k, l, q_cons_vf, rhoe, R2D, TS, mCP, mQ, mCVGP, mCVGP2, mCPD) |
Comment on lines
+464
to
+471
| detJ = Jac(1, 1)*Jac(2, 2) - Jac(1, 2)*Jac(2, 1) | ||
| ! singular Jacobian: no usable Newton direction, accept the current (best) state | ||
| if (detJ == 0.0_wp) exit | ||
|
|
||
| InvJac(1, 1) = Jac(2, 2)/detJ | ||
| InvJac(1, 2) = -Jac(1, 2)/detJ | ||
| InvJac(2, 1) = -Jac(2, 1)/detJ | ||
| InvJac(2, 2) = Jac(1, 1)/detJ |
Comment on lines
+27
to
+32
| integer, parameter :: max_iter = 100000 !< max Newton iterations before accepting the last iterate | ||
| real(wp), parameter :: pCr = 4.94e7_wp !< Critical pressure of water [Pa] | ||
| real(wp), parameter :: TCr = 385.05_wp + 273.15_wp !< Critical temperature of water [K] | ||
| real(wp), parameter :: Tsat_min = 1.0_wp !< lower bracket for the saturation-temperature root [K] | ||
| real(wp), parameter :: tsat_tol = 1.0e-12_wp !< relative convergence tolerance for the s_TSat root find | ||
| integer, parameter :: ptg_ls_max = 30 !< max backtracking-line-search halvings in the pTg solver |
Comment on lines
+355
to
+358
| !> Evaluate the pTg-equilibrium residual R2D and temperature TS at a trial state (ml, pS) WITHOUT mutating q_cons_vf, so the | ||
| !! Newton driver can line-search. The total reacting mass mT is conserved, so the reacting masses are (ml, mT - ml) and only the | ||
| !! inert fluids are read from q_cons_vf. Also returns the mixture sums the Jacobian and the final temperature need. | ||
| subroutine s_compute_ptg_residual(ml, mT, pS, j, k, l, q_cons_vf, rhoe, R2D, TS, mCP, mQ, mCVGP, mCVGP2, mCPD) |
Comment on lines
+13
to
+15
| CMake v3.26.5 on atl1-1-02-008-32-0.pace.gatech.edu | ||
|
|
||
| C : GNU v13.3.0 (/usr/bin/cc) | ||
| Fortran : GNU v13.3.0 (/usr/bin/gfortran) | ||
| C : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc) |
| OpenMP : OFF | ||
|
|
||
| Fypp : /home/user/Documents/GitHub/MFC-JRChreim/build/venv/bin/fypp | ||
| Fypp : /storage/project/r-sbryngelson3-0/sbryngelson3/MFC-phasechange/build/venv/bin/fypp |
tao13146500073
approved these changes
Jul 17, 2026
tao13146500073
left a comment
tao13146500073
left a comment
There was a problem hiding this comment.
Confirmed on my end — the stock 2D_phasechange_bubble that hung at t_step 8 now runs past 200 steps on NVHPC/OpenACC. Also verified convergence independently: dropping max_iter from 100000 to 100 leaves Time/step unchanged (1.422 s both), so no cell is hitting the backstop. My case is 3-fluid (liquid water / vapor / air) with cyl_coord = T, which I think exercises the inert-fluid path in s_compute_ptg_residual that Copilot flagged on pmin — no issue observed.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1654 +/- ##
==========================================
- Coverage 60.95% 59.53% -1.42%
==========================================
Files 83 83
Lines 20003 21141 +1138
Branches 2983 3137 +154
==========================================
+ Hits 12193 12587 +394
- Misses 5782 6444 +662
- Partials 2028 2110 +82 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
sbryngelson
marked this pull request as draft
July 17, 2026 13:47
…PU divergence) The overheated-vapor/subcooled-liquid pT shortcut produced states O(1) different from the pTg equilibrium, selected by a sub-ULP near-tie comparison (TSSL < TSatSL*(1-eps)). Backend rounding differences flipped the branch, inverting the phase and diverging CPU vs GPU on all model-6 cases. pTg's Newton projects the liquid mass onto [0,mT], so it recovers the single-phase limits itself; the shortcuts are removed and pTg is solved unconditionally. Deletes now-dead s_TSat and the saturation-curve coefficients. All 18 phase-change goldens pass CPU and GPU (machine-zero agreement); 6 model-6 goldens regenerated, model-5 unchanged.
…ewton-convergence # Conflicts: # src/common/m_phase_change.fpp
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
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.
Description
Follow-up to #1653 (fyi to @JRChreim) That PR bounded the phase-change Newton loops so a non-convergent cell can no longer hang the solver. This PR fixes the underlying non-convergence so those cells actually reach equilibrium — and does so with bounded, uniform iteration counts, which matters most on GPU.
Root cause of the non-convergence. Instrumenting the stock
examples/2D_phasechange_bubbleshowed the non-convergers weres_TSatands_infinite_ptg_relaxation_k(the pT solver was fine). Both applied a fixedOm = 1e-3underrelaxation that stalled the iterates far from the root — e.g. the pTg Gibbs residual froze at ~1e3 for the low-pressure vapor cells and never moved. The convergence criterion was never the problem; the solver was.s_TSat— reformulated + safeguardedThe saturation temperature is the root of Gibbs equality
g_l = g_v, which for stiffened gases is the dimensionless, strictly-monotone relationThe
A/B/C/Dcoefficients were already computed ins_initialize_phasechange_module— and previously dead (referenced nowhere). This PR solvesG(T)=0with a safeguarded Newton (rtsafe: Newton step when it stays bracketed and reduces the interval, bisection otherwise).Gis strictly increasing over[Tsat_min, TCr]for a liquid–vapor pair, so the root is unique; the solver cannot diverge and needs no accurate initial guess. Out-of-window pressures clamp gracefully.s_infinite_ptg_relaxation_k— globalizedSame analytic 2×2 Jacobian, but the fixed underrelaxation is replaced by a damped Newton with backtracking line search, every step projected onto the physical bounds
0 ≤ ml ≤ mTandpS > pmin, plus a singular-Jacobian guard. The residual evaluation is extracted into a device helpers_compute_ptg_residualso trial states can be evaluated without mutatingq_cons_vf(the total reacting massmTis conserved by construction, somlis the only mass unknown).Incidental (in the routines being rewritten)
p_infpTgandTJac, and the now-unusedp_infpTargument.pO(pO = 0on the first pass) — addresses @copilot's note on fix(phase_change): bound Newton relaxation loops by max_iter #1653.Why this matters on GPU
These solvers run inside
$:GPU_PARALLEL_LOOPcalling[seq]device routines, so a warp cannot retire until its slowest cell's loop exits — a single stalled cell stalls 31 idle lanes. The reformulated solvers converge in a bounded, uniform handful of iterations, which is the real GPU win beyond "no hang." New device code uses theGPU_*macros only.Testing
max_iterbackstop.model 6goldens (which exercise pTg/TSat) are regenerated — they shift by ~1e-6 relative (same equilibrium root, now reached with consistent non-lagged auxiliaries); the 12model 5/ pT-only goldens are unchanged.2D_phasechange_bubblecase that hung att_step = 8on GPU now runs cleanly (reachedt_step > 400).max_itercap from fix(phase_change): bound Newton relaxation loops by max_iter #1653 is retained as a permanent backstop.Checklist
model 6golden files)GPU changes
Note on stacking
This branch builds on #1653, so until that merges the diff here shows both commits (the
max_itercap + this convergence rewrite). Happy to rebase once #1653 lands, or to fold #1653 into this PR if reviewers prefer the complete fix in one place.