Bug report: unbounded Newton loops in m_phase_change.fpp hang the solver (max_iter declared but never used)

[tao13146500073]

Summary

The Newton iteration loops in s_infinite_pt_relaxation_k, s_infinite_ptg_relaxation_k, and the TSat solver have no iteration limit. The module-level parameter

fortran

integer, parameter :: max_iter = 1e8_wp   ! src/common/m_phase_change.fpp:27

is declared but never referenced anywhere in the file. When a cell fails to converge, the do while loop never exits and the run hangs indefinitely (on GPU: the kernel never returns and the host blocks forever in cuStreamSynchronize).

Reproduced with the official example, unmodified: examples/2D_phasechange_bubble.

Reproduction

bash

./mfc.sh run examples/2D_phasechange_bubble/case.py -t pre_process simulation --gpu -n 1

Hangs at t_step = 8 every time.

Variation | Result -- | -- GPU build (nvfortran 25.9, OpenACC) | hangs at t_step = 8 CPU build (gfortran 13.3, --no-gpu --no-mpi) | hangs at t_step = 8 Nx = Ny = 1600 (example default) | hangs at t_step = 8 Nx = Ny = 400 | hangs at t_step = 8

Independent of hardware, compiler, MPI, and grid resolution — i.e. a deterministic non-convergence at the same step, not a race or a memory issue.

Evidence

gdb backtrace of the hung process (GPU build):

#8  cuStreamSynchronize () from /lib/x86_64-linux-gnu/libcuda.so.1
#9  __pgi_uacc_cuda_wait (lineno=-99, async=-1, dindex=1) at ../../src/cuda_wait.c:77
#10 __pgi_uacc_computedone2 (async=-1, ...) at ../../src/computeexitdone.c:59
#11 m_phase_change::s_infinite_relaxation_k (q_cons_vf=...)
        at /home/ubuntu/MFC/src/common/m_phase_change.fpp:250
#12 m_start_up::s_perform_time_step (t_step=8, time_avg=0.0166...)
        at /home/ubuntu/MFC/src/simulation/m_start_up.fpp:631
#13 p_main () at /home/ubuntu/MFC/src/simulation/p_main.fpp:78

The process shows GPU-Util 100% at low power draw (~106 W / 450 W on an RTX 4090), consistent with a small number of threads spinning in a scalar loop rather than doing useful parallel work.

A --debug build additionally reports, before the hang:

FATAL ERROR: DEVICE FORTRAN AUTO ALLOCATION FAILED (USE NVCOMPILER_ACC_CUDA_HEAPSIZE)
Accelerator Fatal Error: call to cuStreamSynchronize returned error 719
  (CUDA_ERROR_LAUNCH_FAILED): Launch failed (often invalid pointer dereference)
  File: src/simulation/m_data_output.fpp

Raising NVCOMPILER_ACC_CUDA_HEAPSIZE to 4 GiB did not change the behaviour (still hangs at t_step = 8), so the device-heap message appears to be a secondary symptom rather than the cause.

Affected code

Three loops, none of which check ns against any bound:

fortran

! s_infinite_pt_relaxation_k, ~line 320
do while ((abs(pS - pO) > palpha_eps) .and. (abs((pS - pO)/pO) > palpha_eps/1.e4_wp) .or. (ns == 0))

! s_infinite_ptg_relaxation_k, ~line 393
do while (((sqrt(R2D(1)**2 + R2D(2)**2) > ptgalpha_eps) .and. ((sqrt(R2D(1)**2 + R2D(2)**2)/rhoe) > (ptgalpha_eps/1.e6_wp) &
          & )) .or. (ns == 0))

! TSat solver, ~line 591
do while ((abs(FT) > ptgalpha_eps) .or. (ns == 0))

Also note max_iter is declared integer but initialised from a real literal (1e8_wp), which is at best misleading.

Suggested minimal fix

Bound each loop with the existing max_iter parameter (patch attached / below). This does not fix the underlying non-convergence, but it turns a hard hang into a bounded, observable failure — and lets the run proceed.

A proper fix would additionally need to decide what to do with cells that exhaust max_iter (skip relaxation and retain the pre-relaxation state, flag them, and optionally report a count per step).

Environment

• MFC master @ 2b1087cc
• nvfortran 25.9 (NVIDIA HPC SDK), OpenACC, CUDA on RTX 4090
• gfortran 13.3.0 (CPU build, --no-mpi)
• Ubuntu 24.04

[sbryngelson]

Will look into this

[sbryngelson]

Thanks for the unusually thorough report — the diagnosis is exactly right, including that the --debug heapsize/error-719 message is a secondary symptom. Fix up in #1653. Here is what we found.

Provenance of the hang. max_iter has been dead since the very first phase-change PR (#179): it was declared and added to the !$acc declare list, but the three do while loops were written convergence-only and never referenced it. A later cleanup (#999) removed it from the (by then macro-migrated) GPU declaration while tidying, leaving the fully-orphaned declaration you found. So it's a latent birth defect rather than a regression — the loops shipped without the cap the author had scaffolded.

Which loop actually spins. Instrumenting the -O3 build (small cap + state dump) on your unmodified 2D_phasechange_bubble shows the non-convergers are s_TSat (~4250 cell-visits) and s_infinite_ptg_relaxation_k (~730); the pT loop converges fine (0). The stuck cells are the low-pressure vapor bubble cells (pSat ≈ 2.3–5.5 kPa). Both loops compare a dimensional Gibbs-free-energy residual (natural magnitude ~10²–10³) against the fixed ptgalpha_eps = 1e-2 under Om = 1e-3 underrelaxation, so the residual stalls at a floor far above tolerance and the exit condition is effectively unreachable. s_TSat is worst — it has no relative/normalized fallback at all.

A -O3 vs -O0 Heisenbug worth flagging for anyone debugging this: under -O3 it's a true infinite loop (finite oscillation, abs(pS-pO) > eps never goes false). Under -O0/--debug the same cell instead goes NaN, which exits the loop (NaN comparisons are false), poisons the field, and WENO then crawls on pow(NaN) — looks hung but is inching forward. So --debug masks the real mechanism.

Fix (#1653): exactly your suggested minimal fix — bound all three loops by max_iter and accept the last iterate on exhaustion; also cleaned up the misleading integer :: max_iter = 1e8_wp to a plain integer. All 18 phase-change regression tests pass byte-identical (the cap sits above the legitimate convergence ceiling, so it truncates no converging cell), and your bubble case now advances past t_step = 8.

As you noted, this bounds the hang but doesn't fix the underlying non-convergence for those low-pressure cells (normalizing the Gibbs residual / adding a relative criterion to s_TSat and the pTg solver). That's a numerics change requiring physics validation, so it'll be tracked as a separate follow-up.

[sbryngelson]

@JRChreim you may want to observe the changed I've been making