Diverging Newton (water equation) with saturation hysteresis on a large black-oil model — Flow 2026.04

[mayur-nawal]

Non-convergence (diverging Newton, water equation) with saturation hysteresis on a large black-oil model — Flow 2026.04

Describe the bug

On a large ECLIPSE (E100) black-oil model, OPM Flow fails to converge on the first timestep whenever saturation hysteresis (SATOPTS HYSTER) is enabled. The Newton iteration diverges: the water residual CNV(W) starts O(1) and explodes to ~1e3–1e4 within 1–2 iterations, with saturation updates DS ≈ 9 (physically impossible for a [0,1] saturation), while CNV(O)/CNV(G) stay O(1–100) and mass-balance residuals are ~1e-6. Timesteps are cut to ~1e-6 days and it still diverges.

With hysteresis disabled, the identical model runs cleanly end-to-end (full field history, 0 convergence failures, physically sensible results).

Environment

• OPM Flow 2026.04 (release build), Linux, MPI (30 ranks).
• Model: ~1.5M active cells (124×200×200, 1,547,477 active), FIELD units, OIL WATER GAS DISGAS BRINE.
• Hysteresis: SATOPTS / HYSTER / (RUNSPEC); EHYSTR 0.1 4 4 1 'PC' 'RETR' 6* 1 / (PROPS, HYSTMOD=4 Killough); ENDSCALE.
• Sat functions: SWOF/SGOF (48 tables) with SATNUM (drainage) + IMBNUM (imbibition) regions (~1.75M cells have SATNUM≠IMBNUM).

Convergence evidence (first timestep, from .DBG)

Iter   MB(O)     MB(W)     MB(G)     CNV(O)    CNV(W)     CNV(G)     DS        CNVFLAG
  0  9.69e-07  5.90e-10  3.25e-07  5.41e-02  1.170e+00  2.50e-02      -       Relaxed
  1  2.10e-09  1.93e-06  7.75e-10  2.427e+0  8.409e+03  9.63e-01  1.07e-02   TuningDP
  2  5.22e-06  8.98e-06  1.24e-07  1.36e+02  1.091e+04  5.39e+01  8.94e+00   Strict

CNV(W) → 8409 → 10910; the wetting-phase update blows up. Full table attached (convergence_evidence.txt).

What we have ruled out (all reproduce the identical failure)

• Initialization — three independent methods: (a) enumerated/explicit PRESSURE/SWAT/SGAS/RS arrays; (b) a genuine RESTART (re-serialized the E100 restart via opm.io.ecl ERst/EclOutput; the hysteresis state arrays SGMAX/SOMAX/SWMAX/SWHY*/SGHY* are present and loaded); (c) EQUIL equilibration (EQUIL+RSVD+SWATINIT). All initialize cleanly (NOSIM passes); all diverge with hysteresis. → not an initialization / hysteresis-state problem.
• Hysteresis model — fails with HYSTMOD=4 and HYSTMOD=0.
• Brine — fails with BRINE on and off (PVTWSALT→PVTW).
• Solver options — no combination helps: --tolerance-cnv-relaxed=0.1..0.3, --tolerance-mb-relaxed=1e-5, --relaxed-max-pv-fraction=1.0, --ds-max=0.05, --newton-max-iterations=50, --time-step-after-event-in-days, --enable-tuning true/false, --project-saturations=true, --tolerance-max-ds=0.2. (Expected — relaxation cannot rescue a diverging Newton.)
• Saturation tables — all 48 SWOF/SGOF tables verified self-consistent (Sgmax = 1-Swco); per-cell SATNUM/IMBNUM endpoints consistent (0 violations). The printed Sgmax should not exceed 1-Swco warnings are floating-point borderline on the exact equality, not real violations.

Expected behavior

Hysteresis run converges (as it does with hysteresis off).

Additional context

• Runs with --check-satfunc-consistency=false (only the ~1e-17 borderline endpoint warnings, otherwise clean).
• The failure looks consistent with discontinuous wetting-phase relperm derivatives at drainage↔imbibition scanning-curve transitions destabilizing the water saturation equation.

convergence_evidence.txt
deck_hysteresis_config.txt

[atgeirr]

@totto82 have you seen similar trouble in hysteresis cases like this before?

[totto82]

Thank you for the report. It is not easy to assess the problem without looking more carefully at the model. Do you have a smaller model you could share that reproduces the problem? Either via GitHub or privately via email, if that is better for you.