From 3680de96084fe373be6aca1997b21e7d62628132 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 11 Sep 2026 00:25:51 +0200 Subject: [PATCH 1/6] Add a saturation pressure solver for compositional mixtures Compute bubble- and dew-point pressures from a cubic EOS for compositional initialization and restart output. Update fugacity ratios with accelerated successive substitution at fixed pressure. Safeguard pressure steps, verify known-phase stability, and use damped Newton steps for slowly converging stability and envelope continuation trials. Trace the vapour/liquid envelope when finite scans miss narrow near-critical intervals, and prefer the upper retrograde dew point. Return failure with output arguments unchanged when no boundary can be resolved. Rescale the composition on entry rather than requiring the caller to normalize. The boundary residual is measured against one, but at a converged boundary sum(K_i z_i) approaches sum(z_i), so a composition summing to 1 + d cannot meet that test once |d| exceeds the tolerance. Restart mole fractions are single precision and routinely sum to 1 +/- 3e-8, which turned real boundaries into silent failures. Stage both outputs and publish them only on success. The search uses the incipient composition as scratch and can still give up after writing it, and the dew search handed the caller's output straight to its first trial before falling back to others. Reuse EOS data across composition-only updates and scale convergence thresholds with Scalar precision while preserving the established double-precision thresholds. Document the equations and the search alongside the code. Test reference pressures, equilibrium residuals, branch selection, pure and supercritical mixtures, binary interactions, narrow envelopes, failure semantics and single-precision round trips, with a deterministic probe over four EOS variants, six temperatures from 195 K to 615 K, pure, binary and random ternary mixtures and both searches. --- CMakeLists_files.cmake | 2 + .../constraintsolvers/SaturationPressure.hpp | 1284 +++++++++++++++++ tests/material/test_saturation_pressure.cpp | 872 +++++++++++ 3 files changed, 2158 insertions(+) create mode 100644 opm/material/constraintsolvers/SaturationPressure.hpp create mode 100644 tests/material/test_saturation_pressure.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 8646e859a74..bc341bedaee 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -715,6 +715,7 @@ list(APPEND DUNE_TEST_SOURCE_FILES tests/material/test_ncpflash.cpp tests/material/test_pengrobinson.cpp tests/material/test_ptflash_ssi_newton_fallback.cpp + tests/material/test_saturation_pressure.cpp tests/material/test_tabulation.cpp tests/material/test_threecomponents_ptflash.cpp tests/material/test_volume_shift.cpp @@ -1349,6 +1350,7 @@ list(APPEND PUBLIC_HEADER_FILES opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp opm/material/constraintsolvers/NcpFlash.hpp opm/material/constraintsolvers/PTFlash.hpp + opm/material/constraintsolvers/SaturationPressure.hpp opm/material/densead/DynamicEvaluation.hpp opm/material/densead/Evaluation.hpp opm/material/densead/Evaluation1.hpp diff --git a/opm/material/constraintsolvers/SaturationPressure.hpp b/opm/material/constraintsolvers/SaturationPressure.hpp new file mode 100644 index 00000000000..3f35f34e8d6 --- /dev/null +++ b/opm/material/constraintsolvers/SaturationPressure.hpp @@ -0,0 +1,1284 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2026 SINTEF Digital + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::SaturationPressure + */ +#ifndef OPM_SATURATION_PRESSURE_HPP +#define OPM_SATURATION_PRESSURE_HPP + +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include + +namespace Opm { + +/*! + * \brief Computes the saturation pressure of a mixture at a given temperature + * from the cubic equation of state. + * + * The bubble-point (dew-point) pressure of a liquid (vapour) with composition + * \f$z\f$ is the pressure where an infinitesimal vapour (liquid) phase is in + * equilibrium with it. At fixed temperature and pressure, phase equilibrium + * requires equal component fugacities: + * \f[ + * x_i\,\phi_i^L(T,p,x) = y_i\,\phi_i^V(T,p,y), + * \qquad + * K_i = \frac{y_i}{x_i} = \frac{\phi_i^L}{\phi_i^V}. + * \f] + * Here \f$x\f$ and \f$y\f$ are liquid and vapour mole fractions and + * \f$\phi_i^L\f$ and \f$\phi_i^V\f$ are the EOS fugacity coefficients. + * For a known liquid, \f$x=z\f$ and the unnormalised incipient vapour is + * \f$Y_i=K_i z_i\f$. For a known vapour, \f$y=z\f$ and the unnormalised + * incipient liquid is \f$X_i=z_i/K_i\f$. Consequently the saturation equations + * are + * \f[ + * S_b(p) = \sum_i K_i(p)z_i = 1 + * \quad\hbox{and}\quad + * S_d(p) = \sum_i \frac{z_i}{K_i(p)} = 1. + * \f] + * The implementation solves fugacity equality by successive substitution at + * fixed pressure, with dominant-eigenvalue extrapolation when it contracts + * slowly. It then solves \f$F(p)=\ln S(p)=0\f$ in \f$\ln p\f$ with fixed-point + * or secant steps, switching to Illinois regula falsi after bracketing the + * boundary. Trivial \f$K_i=1\f$ trials restart from the Wilson estimate at the + * next scan pressure. + * + * Fugacity equality and \f$S=1\f$ also admit stationary points inside a + * two-phase region. Candidate roots are therefore screened by a Michelsen + * tangent-plane stability trial. With unnormalised trial amounts \f$Y_i\f$, + * \f$S_Y=\sum_iY_i\f$, and trial composition \f$w_i=Y_i/S_Y\f$, its stationary + * equations are + * \f[ + * Y_i = z_i\frac{\phi_i(T,p,z)}{\phi_i(T,p,w)}. + * \f] + * At a stationary point, + * \f[ + * \operatorname{TPD}(w) + * = \sum_iw_i\ln\!\frac{w_i\phi_i(T,p,w)}{z_i\phi_i(T,p,z)} + * = -\ln S_Y. + * \f] + * Therefore \f$S_Y>1\f$ proves instability; otherwise this trial does not + * reject the candidate. The trial solver works in \f$u_i=\ln Y_i\f$ and + * combines substitution with damped Newton steps near critical conditions. + * + * A mixture may have lower and upper dew points at a given temperature. The + * upper (retrograde) one is the boundary crossed as pressure declines from + * the single-phase gas region, so the dew-point search tries it first. + * If the scan fails, a bubble-point search or continuation from a lower dew + * point can locate the other boundary of the connected vapour/liquid envelope. + * Continuation advances in \f$\ln p\f$, carries stationary trial amounts from + * one pressure to the next, and brackets the sign change of \f$\ln S_Y\f$. + * A failed scan does not establish that no root exists: near a critical point, + * the two-phase interval can be narrower than the pressure sampling step. + * Continuation narrows that gap without closing it. Sweeping a ternary across + * temperature and composition leaves the bubble search with no unresolved + * point between two resolved neighbours, but the dew search keeps isolated + * bands, a few parts in ten thousand of mole fraction wide, at the crossover + * between its lower and upper branches. Both branches are well separated and + * stable on either side of such a band, so a boundary does exist inside it. + * They survive because the lower-branch result is only accepted once a bubble + * point confirms an ordinary envelope, and a mixture gas-rich enough to sit at + * that crossover has no bubble point to find. Resolving them needs a bracketing + * scheme that cannot skip an unsampled interval, not a smaller fixed step. + * Input compositions must contain finite, non-negative mole fractions with a + * positive total. They are rescaled to sum to one, so they need not do so on + * entry. + */ +template +class SaturationPressure +{ + static constexpr int numComponents = FluidSystem::numComponents; + static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx; + static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx; + + using EOSType = CompositionalConfig::EOSType; + using ParameterCache = typename FluidSystem::template ParameterCache; + +public: + using CompVec = std::array; + + /// Computes the bubble-point pressure of a liquid with composition \p liquid at + /// temperature \p temp, along with the equilibrium \p vapor composition. + /// On failure, \p press and \p vapor are both unchanged. + /// \return whether the calculation converged + [[nodiscard]] static bool bubblePressure(const CompVec& liquid, + const Scalar temp, + const EOSType eosType, + Scalar& press, + CompVec& vapor) + { + // The search uses the incipient composition as scratch, so give it one + // of its own and publish only once it has converged. + Scalar p{}; + CompVec incipient{}; + if (solve_(normalized_(liquid), temp, eosType, Mode::Bubble, p, incipient) + != Outcome::Converged) { + return false; + } + + press = p; + vapor = incipient; + return true; + } + + /// Computes the dew-point pressure of a vapour with composition \p vapor at + /// temperature \p temp, along with the equilibrium \p liquid composition. + /// The upper (retrograde) dew point is preferred. If its initial search + /// fails, trace the envelope from a lower dew point or a bubble point. + /// Return the lower point only for a pure fluid or when a bubble point + /// establishes the other boundary of the connected vapour/liquid envelope. + /// On failure, \p press and \p liquid are both unchanged. + /// \return whether the calculation converged + [[nodiscard]] static bool dewPressure(const CompVec& vapor, + const Scalar temp, + const EOSType eosType, + Scalar& press, + CompVec& liquid) + { + // Every trial below writes through to its output, and a later one may + // still fail, so stage both and publish only on success. + Scalar p{}; + CompVec incipient{}; + if (!dewPressure_(normalized_(vapor), temp, eosType, p, incipient)) { + return false; + } + + press = p; + liquid = incipient; + return true; + } + +private: + /// \copydoc dewPressure + /// Takes an already normalised \p z and may write to its outputs on failure. + static bool dewPressure_(const CompVec& z, + const Scalar temp, + const EOSType eosType, + Scalar& press, + CompVec& liquid) + { + if (solve_(z, temp, eosType, Mode::DewUpper, press, liquid) + == Outcome::Converged) { + return true; + } + + Scalar pLower{}; + CompVec lowerLiquid{}; + if (solve_(z, temp, eosType, Mode::DewLower, pLower, lowerLiquid) + != Outcome::Converged) { + // Both dew scans can miss a narrow envelope. Try locating its + // bubble boundary, then trace down to the dew point. + Scalar pBubble{}; + CompVec bubbleVapor{}; + if (solve_(z, temp, eosType, Mode::Bubble, pBubble, bubbleVapor) + != Outcome::Converged) { + return false; + } + return traceEnvelope_(z, temp, eosType, pBubble, bubbleVapor, + /*fromBubble=*/true, press, liquid); + } + + if (std::ranges::count_if(z, [](const Scalar zi) { return zi > 0.0; }) == 1) { + // A pure fluid has the same bubble and dew pressure. + press = pLower; + liquid = lowerLiquid; + return true; + } + + // An independently located bubble boundary identifies an ordinary + // envelope without needing to trace it from the lower dew point. + Scalar pBubble{}; + CompVec bubbleVapor{}; + if (solve_(z, temp, eosType, Mode::Bubble, pBubble, bubbleVapor) + == Outcome::Converged && pBubble >= pLower) { + press = pLower; + liquid = lowerLiquid; + return true; + } + // A lower root alone is a seed, not a fallback answer. Start just + // inside and follow the envelope to its upper dew or bubble boundary. + return traceEnvelope_(z, temp, eosType, pLower, lowerLiquid, + /*fromBubble=*/false, press, liquid); + } + + // A finite scan cannot establish that an unsampled interval has no root. + enum class Outcome { Converged, GaveUp }; + + // The saturation-pressure branch being searched. The bubble point and the + // upper (retrograde) dew point are approached from the high-pressure side, + // the lower dew point from the low-pressure side. + enum class Mode { Bubble, DewUpper, DewLower }; + + enum class Stability { Stable, Unstable, Indeterminate }; + + // Keep the established double-precision thresholds while making each + // convergence test meaningful at Scalar's precision. + static constexpr Scalar precisionTolerance_(const Scalar minimum) + { + return std::max(minimum, + Scalar{100} * std::numeric_limits::epsilon()); + } + + static constexpr Scalar nonlinearTolerance_ = precisionTolerance_(Scalar{1.0e-12}); + static constexpr Scalar stabilityTolerance_ = precisionTolerance_(Scalar{1.0e-8}); + static constexpr Scalar boundaryResidualTolerance_ = precisionTolerance_(Scalar{1.0e-10}); + static constexpr Scalar boundaryPressureTolerance_ = precisionTolerance_(Scalar{1.0e-8}); + static constexpr Scalar rootVolumeTolerance_ = precisionTolerance_(Scalar{1.0e-9}); + static constexpr Scalar directionTolerance_ = precisionTolerance_(Scalar{1.0e-6}); + + // Rescale a composition to sum to one. The boundary residual is measured + // against one, so a composition that sums to 1 + d cannot meet that test + // once |d| exceeds the tolerance. Callers assemble z in floating point, + // from a flash or from a restart file's single-precision mole fractions, + // so it need not sum exactly. + static CompVec normalized_(const CompVec& z) + { + const Scalar sum = std::accumulate(z.begin(), z.end(), Scalar{0}); + if (!positiveFinite_(sum)) { + return z; + } + + CompVec out{}; + std::ranges::transform(z, out.begin(), + [sum](const Scalar zi) { return zi / sum; }); + return out; + } + + static bool positiveFinite_(const Scalar value) + { + return std::isfinite(value) && value > 0.0; + } + + /*! + * \brief Apply bounded dominant-eigenvalue extrapolation to a fixed-point + * iterate. + * + * Let \f$d^n=\ln v^{n+1}-\ln v^n\f$ be the current logarithmic update and + * \f$d^{n-1}\f$ the previous one. The dominant contraction factor and the + * remaining geometric-series correction are estimated as + * \f[ + * r = \frac{d^n\mathbin{\cdot}d^{n-1}} + * {d^{n-1}\mathbin{\cdot}d^{n-1}}, + * \qquad + * \Delta\ln v_i = \frac{r}{1-r}d_i^n. + * \f] + * Only aligned updates (\f$r>0\f$) are accelerated. The estimate is capped + * at 0.98 and each component correction at \f$|\Delta\ln v_i|\leq\ln 2\f$, + * so a single extrapolation changes a positive amount by at most a factor + * of two and cannot deliberately underflow it to zero. + * + * \param[in,out] values Positive fixed-point variables \f$v\f$. + * \param[in] d Current logarithmic update \f$d^n\f$. + * \param[in] dPrev Previous logarithmic update \f$d^{n-1}\f$. + */ + static void accelerate_(CompVec& values, const CompVec& d, const CompVec& dPrev) + { + Scalar num = 0.0; + Scalar den = 0.0; + for (int c = 0; c < numComponents; ++c) { + num += d[c] * dPrev[c]; + den += dPrev[c] * dPrev[c]; + } + if (!positiveFinite_(den) || !positiveFinite_(num)) { + return; + } + const Scalar ratio = std::min(num / den, Scalar{0.98}); + const Scalar remaining = ratio / (1.0 - ratio); + const Scalar maxStep = std::log(Scalar{2}); + CompVec next = values; + for (int c = 0; c < numComponents; ++c) { + if (values[c] == 0.0) { // Absent component in a stability trial. + continue; + } + const Scalar step = remaining * d[c]; + if (!std::isfinite(step)) { + return; + } + next[c] *= std::exp(std::clamp(step, -maxStep, maxStep)); + if (!positiveFinite_(next[c])) { + return; + } + } + values = next; + } + + /*! + * \brief Evaluate the pressure-independent numerator of the Wilson + * equilibrium-ratio estimate. + * + * The Wilson correlation is + * \f[ + * K_i^W(T,p) = \frac{p_{c,i}}{p} + * \exp\!\left[5.373(1+\omega_i) + * \left(1-\frac{T_{c,i}}{T}\right)\right] + * = \frac{Kp_i(T)}{p}. + * \f] + * This function returns \f$Kp_i(T)\f$; callers divide by pressure. Wilson + * values initialise phase-composition iterations and the pressure scan. + * + * \param[in] temp Absolute temperature \f$T\f$. + */ + static CompVec wilsonKp_(const Scalar temp) + { + CompVec Kp; + for (int c = 0; c < numComponents; ++c) { + Kp[c] = FluidSystem::criticalPressure(c) * + std::exp(5.373 * (1.0 + FluidSystem::acentricFactor(c)) * + (1.0 - FluidSystem::criticalTemperature(c) / temp)); + } + return Kp; + } + + /*! + * \brief Test a candidate known phase in the omitted tangent-plane trial + * direction. + * + * The method starts an additional Wilson-initialised trial on the same EOS + * root as the known phase, covering the direction omitted by the saturation + * iteration. For \f$w_i=Y_i/\sum_jY_j\f$, successive substitution applies + * \f[ + * Y_i^{n+1} = z_i + * \frac{\phi_i(T,p,z)}{\phi_i(T,p,w^n)}. + * \f] + * At a converged stationary point, \f$\sum_iY_i>1\f$ is equivalent to a + * negative tangent-plane distance and therefore proves that \f$z\f$ is + * unstable. A sum at or below one passes this trial within tolerance. + * Slowly contracting trials are completed by stationaryTrial_(); numerical + * failure produces Stability::Indeterminate rather than a stability claim. + * + * \param[in] fs Fluid state containing \f$T\f$, \f$p\f$, and \f$z\f$. + * \param[in] z Composition of the phase being tested. + * \param[in] knownPhaseIdx EOS phase/root used by the known phase and trial. + * \param[in] wilsonK Wilson equilibrium-ratio estimate at the trial pressure. + * \param[in] knownIsLiquid Selects the heavy or light Wilson trial direction. + * \param[in] eosType Cubic EOS type. + */ + static Stability knownPhaseStability_(const CompositionalFluidState& fs, + const CompVec& z, + const unsigned knownPhaseIdx, + const CompVec& wilsonK, + const bool knownIsLiquid, + const EOSType eosType) + { + // The known phase keeps its own state; the trial takes its root type + // through the same phase index in a second state. + ParameterCache knownCache(eosType); + knownCache.updatePhase(fs, knownPhaseIdx); + CompVec lnPhiZ; + for (int c = 0; c < numComponents; ++c) { + lnPhiZ[c] = std::log( + FluidSystem::fugacityCoefficient(fs, knownCache, knownPhaseIdx, c)); + if (!std::isfinite(lnPhiZ[c])) { + return Stability::Indeterminate; + } + } + + // A liquid is tested against a heavier liquid, a vapour against a + // lighter vapour: the direction the incipient-phase search leaves out. + CompVec Y; + for (int c = 0; c < numComponents; ++c) { + Y[c] = knownIsLiquid ? z[c] / wilsonK[c] : z[c] * wilsonK[c]; + } + + auto trial = fs; + ParameterCache trialCache(eosType); + // Dominant-eigenvalue extrapolation, as in the saturation search. + CompVec dPrev{}; + Scalar sumPrev = std::numeric_limits::quiet_NaN(); + int settled = 0; + constexpr int maxStabilityIterations = 500; + for (int iter = 0; iter < maxStabilityIterations; ++iter) { + Scalar sumY = 0.0; + for (const Scalar amount : Y) { + sumY += amount; + } + if (!positiveFinite_(sumY)) { + return Stability::Indeterminate; + } + for (int c = 0; c < numComponents; ++c) { + trial.setMoleFraction(knownPhaseIdx, c, Y[c] / sumY); + } + const int unchanged = iter == 0 + ? ParameterCache::None : ParameterCache::Temperature | ParameterCache::Pressure; + trialCache.updatePhase(trial, knownPhaseIdx, unchanged); + + Scalar change = 0.0; + Scalar sumNew = 0.0; + CompVec d{}; + for (int c = 0; c < numComponents; ++c) { + if (z[c] <= 0.0) { + Y[c] = 0.0; + continue; + } + const Scalar lnPhiY = + std::log(FluidSystem::fugacityCoefficient(trial, trialCache, knownPhaseIdx, c)); + const Scalar Ynew = std::exp(std::log(z[c]) + lnPhiZ[c] - lnPhiY); + if (!positiveFinite_(Ynew) || !positiveFinite_(Y[c])) { + return Stability::Indeterminate; + } + change = std::max(change, std::abs(Ynew - Y[c]) / std::max(Scalar{1}, Ynew)); + d[c] = std::log(Ynew) - std::log(Y[c]); + Y[c] = Ynew; + sumNew += Ynew; + } + // At a stationary point the total decides stability. A flat total + // during iteration does not establish stationarity: near critical + // states use Newton to finish the slowly contracting trial. + const auto verdict = [](const Scalar total) { + return total > 1.0 + stabilityTolerance_ + ? Stability::Unstable : Stability::Stable; + }; + if (change < boundaryResidualTolerance_) { + return verdict(sumNew); + } + if (std::abs(sumNew - sumPrev) + <= nonlinearTolerance_ * std::max(Scalar{1}, sumNew)) { + if (++settled == 8) { + CompVec candidate = Y; + Scalar total{}; + if (stationaryTrial_(fs, z, knownPhaseIdx, knownPhaseIdx, + eosType, candidate, total)) { + return verdict(total); + } + } + } + else { + settled = 0; + } + sumPrev = sumNew; + if (iter % 4 == 3) { + accelerate_(Y, d, dPrev); + } + dPrev = d; + } + Scalar total{}; + if (stationaryTrial_(fs, z, knownPhaseIdx, knownPhaseIdx, eosType, Y, total)) { + return total > 1.0 + stabilityTolerance_ + ? Stability::Unstable : Stability::Stable; + } + return Stability::Indeterminate; + } + + /*! + * \brief Converge one tangent-plane stationary trial at fixed temperature + * and pressure. + * + * The unknowns are logarithmic trial amounts \f$u_i=\ln Y_i\f$, which keep + * all present components positive. With \f$S_Y=\sum_j\exp(u_j)\f$ and + * \f$w_i=\exp(u_i)/S_Y\f$, the nonlinear residual is + * \f[ + * R_i(u) = u_i + \ln\phi_i(T,p,w) + * - \ln z_i - \ln\phi_i(T,p,z). + * \f] + * Thus \f$R_i=0\f$ is the stationary equation + * \f$Y_i=z_i\phi_i(z)/\phi_i(w)\f$. The first iterations use the equivalent + * fixed-point update \f$u_i\leftarrow u_i-R_i\f$, limited to \f$\ln 2\f$ + * per component. The method then forms a finite-difference Jacobian in + * \f$u\f$ and accepts damped Newton steps only when they reduce + * \f$\lVert R\rVert_\infty\f$. + * + * \param[in] fs Fluid state containing the fixed \f$T\f$ and \f$p\f$. + * \param[in] z Composition and reference fugacity of the known phase. + * \param[in] knownPhase EOS root used to evaluate \f$\phi_i(T,p,z)\f$. + * \param[in] trialPhase EOS root used to evaluate \f$\phi_i(T,p,w)\f$. + * \param[in] eosType Cubic EOS type. + * \param[in,out] Y Initial and converged unnormalised trial amounts. + * \param[out] sum Converged value of \f$S_Y=\sum_iY_i\f$. + * \return true only when the stationary residual converges. + */ + static bool stationaryTrial_(const CompositionalFluidState& fs, + const CompVec& z, + const unsigned knownPhase, + const unsigned trialPhase, + const EOSType eosType, + CompVec& Y, + Scalar& sum) + { + using Vector = Dune::FieldVector; + using Matrix = Dune::FieldMatrix; + ParameterCache cache(eosType); + cache.updatePhase(fs, knownPhase); + CompVec target{}; + Vector u(0.0); + for (int c = 0; c < numComponents; ++c) { + if (z[c] > 0.0) { + if (!positiveFinite_(Y[c])) { + return false; + } + target[c] = std::log(z[c]) + + std::log(FluidSystem::fugacityCoefficient(fs, cache, knownPhase, c)); + u[c] = std::log(Y[c]); + } + } + auto trial = fs; + bool trialCacheInitialized = false; + const auto residual = [&](const Vector& v, Vector& r) { + CompVec amounts{}; + Scalar total = 0.0; + for (int c = 0; c < numComponents; ++c) { + amounts[c] = z[c] > 0.0 ? std::exp(v[c]) : Scalar{0}; + if (z[c] > 0.0 && !positiveFinite_(amounts[c])) { + return false; + } + total += amounts[c]; + } + if (!positiveFinite_(total)) { + return false; + } + for (int c = 0; c < numComponents; ++c) { + trial.setMoleFraction(trialPhase, c, amounts[c] / total); + } + const int unchanged = trialCacheInitialized + ? ParameterCache::Temperature | ParameterCache::Pressure : ParameterCache::None; + cache.updatePhase(trial, trialPhase, unchanged); + trialCacheInitialized = true; + for (int c = 0; c < numComponents; ++c) { + r[c] = z[c] > 0.0 + ? v[c] + std::log(FluidSystem::fugacityCoefficient(trial, cache, trialPhase, c)) + - target[c] + : v[c]; + if (!std::isfinite(r[c])) { + return false; + } + } + return true; + }; + + // Allow substitution to converge before forming the Newton Jacobian. + constexpr int maxTrialIterations = 100; + constexpr int substitutionPasses = 8; + constexpr int maxBacktracks = 10; + for (int iter = 0; iter < maxTrialIterations; ++iter) { + Vector r; + if (!residual(u, r)) { + return false; + } + const Scalar norm = r.infinity_norm(); + if (norm < nonlinearTolerance_) { + sum = 0.0; + for (int c = 0; c < numComponents; ++c) { + Y[c] = z[c] > 0.0 ? std::exp(u[c]) : Scalar{0}; + sum += Y[c]; + } + return positiveFinite_(sum); + } + + bool accepted = false; + if (iter >= substitutionPasses) { + Matrix jac(0.0); + // A fixed step in ln Y gives a relative perturbation of Y. + constexpr Scalar h = 1.0e-5; + for (int c = 0; c < numComponents; ++c) { + Vector v = u; + v[c] += h; + Vector perturbed; + if (!residual(v, perturbed)) { + return false; + } + for (int row = 0; row < numComponents; ++row) { + jac[row][c] = (perturbed[row] - r[row]) / h; + } + } + Vector step; + try { + jac.solve(step, r); + // Cap the first trial step at two in the log variables. + const Scalar stepNorm = step.infinity_norm(); + Scalar damping = (stepNorm > 0.0) + ? std::min(Scalar{1}, Scalar{2} / stepNorm) + : Scalar{1}; + for (int backtrack = 0; backtrack < maxBacktracks; ++backtrack) { + Vector v = u; + v.axpy(-damping, step); + Vector next; + if (residual(v, next) && next.infinity_norm() < norm) { + u = v; + accepted = true; + break; + } + damping *= 0.5; + } + } + catch (const Dune::FMatrixError&) { + // A singular Newton system does not invalidate the + // substitution step or establish phase stability. + } + } + if (!accepted) { + for (int c = 0; c < numComponents; ++c) { + u[c] -= std::clamp(r[c], -std::log(Scalar{2}), std::log(Scalar{2})); + } + } + } + return false; + } + + /*! + * \brief Follow a connected two-phase interval from one known boundary to + * the other. + * + * The continuation variable is \f$q=\ln p\f$. At every pressure, both + * liquid-like and vapour-like stationary trials are solved from the + * continued amounts and from fresh Wilson seeds. Among nontrivial trials, + * the residual is the largest converged value + * \f[ + * F(q)=\max_k\ln\!\left(\sum_iY_i^{(k)}(q)\right). + * \f] + * A positive value proves that the current pressure lies inside the + * two-phase interval. Starting a small distance inside the seed boundary, + * the method walks in \f$q\f$, carries \f$Y\f$ between pressures, and grows + * the step until \f$F\leq0\f$ brackets the other boundary. A safeguarded + * secant method then solves \f$F=0\f$; trial amounts are interpolated in + * logarithmic space to preserve positivity and branch continuity. + * + * The converged boundary is classified from enrichment relative to Wilson + * volatility, + * \f[ + * D=\sum_i(w_i-z_i)\ln K_i^W. + * \f] + * \f$D>0\f$ denotes a vapour-like incipient phase (bubble point), while + * \f$D<0\f$ denotes a liquid-like incipient phase (dew point). Fugacity + * equality, distinct physical EOS roots, and known-phase stability are + * checked before returning the boundary. + * + * From a lower dew seed the method walks upward: an upper dew point is + * returned, whereas a bubble boundary confirms that the seed was the only + * dew point. From a bubble seed it walks downward and accepts only a dew + * boundary. + * + * \param[in] z Composition whose saturation pressure is sought. + * \param[in] temp Absolute temperature. + * \param[in] eosType Cubic EOS type. + * \param[in] pSeed Pressure of the known boundary. + * \param[in] seedTrial Equilibrium composition at the known boundary. + * \param[in] fromBubble Whether the seed is a bubble rather than dew point. + * \param[out] press Selected dew-point pressure. + * \param[out] liquid Incipient liquid composition at \p press. + * \return whether continuation found and certified the requested result. + */ + static bool traceEnvelope_(const CompVec& z, + const Scalar temp, + const EOSType eosType, + const Scalar pSeed, + const CompVec& seedTrial, + const bool fromBubble, + Scalar& press, + CompVec& liquid) + { + const Scalar dir = fromBubble ? Scalar{-1} : Scalar{1}; + CompositionalFluidState fs; + fs.setTemperature(temp); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(gasPhaseIdx, c, z[c]); + fs.setMoleFraction(oilPhaseIdx, c, z[c]); + } + const auto Kp = wilsonKp_(temp); + const auto evaluate = [&](const Scalar lnp, CompVec& Y, Scalar& f) { + const Scalar p = std::exp(lnp); + if (!positiveFinite_(p)) { + return false; + } + fs.setPressure(oilPhaseIdx, p); + fs.setPressure(gasPhaseIdx, p); + ParameterCache cache(eosType); + cache.updatePhase(fs, oilPhaseIdx); + cache.updatePhase(fs, gasPhaseIdx); + Scalar gibbsDifference = 0.0; + for (int c = 0; c < numComponents; ++c) { + if (z[c] > 0.0) { + gibbsDifference += z[c] * std::log( + FluidSystem::fugacityCoefficient(fs, cache, oilPhaseIdx, c) + / FluidSystem::fugacityCoefficient(fs, cache, gasPhaseIdx, c)); + } + } + if (!std::isfinite(gibbsDifference)) { + return false; + } + const unsigned knownPhase = gibbsDifference < 0.0 ? oilPhaseIdx : gasPhaseIdx; + Scalar best = -std::numeric_limits::infinity(); + CompVec bestY{}; + bool complete = true; + for (const unsigned trialPhase : {oilPhaseIdx, gasPhaseIdx}) { + bool converged = false; + for (int seed = 0; seed < 2; ++seed) { + CompVec candidate = Y; + if (seed == 1) { + for (int c = 0; c < numComponents; ++c) { + const Scalar k = Kp[c] / p; + candidate[c] = trialPhase == oilPhaseIdx ? z[c] / k : z[c] * k; + } + } + Scalar sum{}; + if (!stationaryTrial_(fs, z, knownPhase, trialPhase, eosType, candidate, sum)) { + continue; + } + converged = true; + Scalar distance = 0.0; + for (int c = 0; c < numComponents; ++c) { + distance += std::abs(candidate[c] / sum - z[c]); + } + // The trivial zero is not a sign for root bracketing. + if (distance > 1.0e-3 && std::log(sum) > best) { + best = std::log(sum); + bestY = candidate; + } + } + complete = complete && converged; + } + if (!std::isfinite(best) || (best <= 0.0 && !complete)) { + return false; + } + f = best; + Y = bestY; + return true; + }; + + const Scalar seed = std::log(pSeed); + Scalar lo{}, fLo{}; + CompVec yLo{}; + Scalar step = 0.01; + bool inside = false; + // Halve the seed offset until a trial lands inside the envelope. + constexpr int maxSeedRefinements = 20; + constexpr int maxWalkSteps = 200; + constexpr int maxBoundarySteps = 80; + for (int refine = 0; refine < maxSeedRefinements; ++refine) { + lo = seed + dir * step; + yLo = seedTrial; + if (evaluate(lo, yLo, fLo) && fLo > boundaryResidualTolerance_) { + inside = true; + break; + } + step *= 0.5; + } + if (!inside) { + return false; + } + + Scalar hi{}, fHi{}; + CompVec yHi{}; + bool bracketed = false; + for (int iter = 0; iter < maxWalkSteps; ++iter) { + hi = lo + dir * step; + yHi = yLo; + if (!evaluate(hi, yHi, fHi)) { + step *= 0.5; + if (step < boundaryPressureTolerance_) { + return false; + } + continue; + } + if (fHi <= 0.0) { + bracketed = true; + break; + } + lo = hi; + fLo = fHi; + yLo = yHi; + step = std::min(Scalar{0.1}, Scalar{2} * step); + } + if (!bracketed) { + return false; + } + + for (int iter = 0; iter < maxBoundarySteps; ++iter) { + // Safeguard the secant so that even a flat residual contracts + // the bracket; interpolate the trial amounts in log space. + const Scalar fraction = std::clamp(fLo / (fLo - fHi), Scalar{0.1}, Scalar{0.9}); + const Scalar mid = lo + fraction * (hi - lo); + CompVec Y{}; + for (int c = 0; c < numComponents; ++c) { + if (z[c] > 0.0) { + Y[c] = std::exp((1.0 - fraction) * std::log(yLo[c]) + + fraction * std::log(yHi[c])); + } + } + Scalar f{}; + if (!evaluate(mid, Y, f)) { + return false; + } + if (std::abs(hi - lo) < boundaryPressureTolerance_ + && std::abs(f) < boundaryResidualTolerance_) { + const Scalar p = std::exp(mid); + CompVec K = Kp; + Scalar direction = 0.0; + for (int c = 0; c < numComponents; ++c) { + K[c] /= p; + Y[c] /= std::exp(f); + direction += (Y[c] - z[c]) * std::log(K[c]); + } + // Identify the continued trial by its enrichment in the + // Wilson volatility direction, then verify fugacity equality + // with the corresponding liquid/vapour EOS roots. + if (std::abs(direction) < directionTolerance_) { + return false; + } + const bool bubble = direction > 0.0; + const unsigned knownPhase = bubble ? oilPhaseIdx : gasPhaseIdx; + const unsigned trialPhase = bubble ? gasPhaseIdx : oilPhaseIdx; + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(knownPhase, c, z[c]); + fs.setMoleFraction(trialPhase, c, Y[c]); + } + ParameterCache cache(eosType); + cache.updatePhase(fs, oilPhaseIdx); + cache.updatePhase(fs, gasPhaseIdx); + const Scalar vmL = cache.molarVolume(oilPhaseIdx); + const Scalar vmV = cache.molarVolume(gasPhaseIdx); + if (!positiveFinite_(vmL) || !positiveFinite_(vmV) + || std::min(vmL, vmV) <= 2.0e-7) { + return false; + } + for (int c = 0; c < numComponents; ++c) { + if (z[c] > 0.0) { + const Scalar a = z[c] * FluidSystem::fugacityCoefficient( + fs, cache, knownPhase, c); + const Scalar b = Y[c] * FluidSystem::fugacityCoefficient( + fs, cache, trialPhase, c); + if (!positiveFinite_(a) || !positiveFinite_(b) + || std::abs(std::log(a / b)) > stabilityTolerance_) { + return false; + } + } + } + if (knownPhaseStability_(fs, z, knownPhase, K, bubble, eosType) + != Stability::Stable) { + return false; + } + if (fromBubble) { + // Down from a bubble point the envelope can only close + // at a dew point; anything else is not this envelope. + if (bubble) { + return false; + } + press = p; + liquid = Y; + return true; + } + press = bubble ? pSeed : p; + liquid = bubble ? seedTrial : Y; + return true; + } + if (f > 0.0) { + lo = mid; + fLo = f; + yLo = Y; + } + else { + hi = mid; + fHi = f; + yHi = Y; + } + } + return false; + } + + /*! + * \brief Search one bubble- or dew-pressure branch. + * + * Wilson's \f$K_i^W=Kp_i/p\f$ supplies a pressure at which its saturation + * sum is one: + * \f[ + * p_0^b=\sum_i z_iKp_i, + * \qquad + * p_0^d=\left(\sum_i\frac{z_i}{Kp_i}\right)^{-1}. + * \f] + * Bubble and upper-dew searches move down from the high-pressure estimate; + * the lower-dew search moves upward from \f$p_0^d\f$. + * + * At each pressure, successive substitution solves + * \f[ + * K_i^{n+1}=\frac{\phi_i^L(z)}{\phi_i^V(y^n)} + * \quad\hbox{(bubble)}, + * \qquad + * K_i^{n+1}=\frac{\phi_i^L(x^n)}{\phi_i^V(z)} + * \quad\hbox{(dew)}, + * \f] + * with \f$y_i^n=K_i^nz_i/S_b\f$ or + * \f$x_i^n=z_i/(K_i^nS_d)\f$. Once the fixed-pressure iteration converges, + * \f$F=\ln S\f$ classifies and drives the pressure search: \f$F>0\f$ is + * inside the two-phase region and \f$F<0\f$ is outside. Before bracketing, + * the method uses a direction-checked secant step or the fixed-point step + * \f$\ln p_{n+1}=\ln p_n\mathbin{\pm}F_n\f$. Within a bracket it uses the + * Illinois regula-falsi formula + * \f[ + * q_{n+1}=q_{out}-F_{out} + * \frac{q_{in}-q_{out}}{F_{in}-F_{out}}, + * \qquad q=\ln p, + * \f] + * with bisection safeguards. Pressure changes are capped at a factor of + * two. A converged \f$S=1\f$ state is returned only after rejecting the + * trivial \f$K=1\f$ solution and certifying known-phase stability. + * + * A finite scan can miss a narrow phase interval, so exhausting the search + * returns Outcome::GaveUp and never claims that the branch has no root. + * + * \param[in] z Known liquid composition for bubble mode, or known vapour + * composition for either dew mode. + * \param[in] temp Absolute temperature. + * \param[in] eosType Cubic EOS type. + * \param[in] mode Branch and scan direction. + * \param[out] press Saturation pressure; unchanged unless converged. + * \param[out] incipient Equilibrium composition of the incipient phase. + * \return Outcome::Converged for a certified boundary, otherwise + * Outcome::GaveUp. + */ + static Outcome solve_(const CompVec& z, + const Scalar temp, + const EOSType eosType, + const Mode mode, + Scalar& press, + CompVec& incipient) + { + if (!positiveFinite_(temp)) { + return Outcome::GaveUp; + } + const CompVec wilsonKp = wilsonKp_(temp); + const bool bubble = (mode == Mode::Bubble); + // The bubble point and the retrograde dew point are approached from + // the high-pressure side, the lower dew point from the low-pressure side. + const bool fromAbove = (mode != Mode::DewLower); + + // Wilson K scales as 1/p, so its saturation sum supplies the initial + // pressure. The upper dew search starts from the same high estimate as + // the bubble search. + Scalar p = 0.0; + if (mode == Mode::DewLower) { + for (int c = 0; c < numComponents; ++c) { + p += z[c] / wilsonKp[c]; + } + p = 1.0 / p; + } + else { + for (int c = 0; c < numComponents; ++c) { + p += z[c] * wilsonKp[c]; + } + } + + const auto wilsonK = [&wilsonKp](const Scalar pressure) { + CompVec K; + std::ranges::transform(wilsonKp, K.begin(), + [pressure](const Scalar Kp) { return Kp / pressure; }); + return K; + }; + CompVec K = wilsonK(p); + + // The known phase holds z; the incipient phase composition is derived from K. + const auto knownPhaseIdx = bubble ? oilPhaseIdx : gasPhaseIdx; + const auto incipientPhaseIdx = bubble ? gasPhaseIdx : oilPhaseIdx; + + CompositionalFluidState fs; + fs.setTemperature(temp); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(knownPhaseIdx, c, z[c]); + } + + // pOut is a tentative endpoint from a trivial trial or a converged + // incipient amount below one. pIn has a converged amount above one, + // establishing instability. Distinct EOS roots alone do not determine + // which side of the boundary a trial occupies. + Scalar pOut{}, pIn{}; + bool haveOut = false; + bool haveIn = false; + // ln(sum) at the bracket ends. A trivial trial supplies no residual, + // so fOut is only known after a nontrivial converged one. + Scalar fOut{}, fIn{}; + bool haveFOut = false; + // The previous converged non-trivial trial, for the secant step, and + // which bracket end the previous step replaced (Illinois damping). + Scalar lnpPrev{}, fPrev{}; + bool havePrev = false; + int lastSide = 0; + // This scan step can skip a narrow envelope; dewPressure() tries + // continuation from another boundary if the scan fails. A smaller step + // lowers the odds without removing them -- see the class documentation + // for what a fixed step leaves unresolved. + constexpr Scalar scanStep = 0.9; + constexpr int maxPressureIterations = 200; + constexpr int maxSubstitutionIterations = 500; + + const auto safeguardPressure = [&](Scalar pNext) { + if (haveIn && haveOut) { + const Scalar lo = std::min(pIn, pOut); + const Scalar hi = std::max(pIn, pOut); + if (!(pNext > lo && pNext < hi)) { + pNext = std::sqrt(lo * hi); + } + } + else if (haveOut && ((fromAbove && pNext >= pOut) || + (!fromAbove && pNext <= pOut))) + { + pNext = std::sqrt(p * pOut); + } + return pNext; + }; + + for (int outer = 0; outer < maxPressureIterations; ++outer) { + if (!positiveFinite_(p)) { + return Outcome::GaveUp; + } + fs.setPressure(oilPhaseIdx, p); + fs.setPressure(gasPhaseIdx, p); + + ParameterCache paramCache(eosType); + CompVec phiKnown; + + // Fugacity equality at fixed pressure: K_c = phi_liquid / phi_vapour. + bool trivial = false; + bool rootsDistinct = false; + bool substitutionConverged = false; + // The previous change of ln K, for the dominant-eigenvalue + // extrapolation below. + CompVec dPrev{}; + for (int inner = 0; inner < maxSubstitutionIterations; ++inner) { + Scalar sum = 0.0; + for (int c = 0; c < numComponents; ++c) { + if (!positiveFinite_(K[c])) { + return Outcome::GaveUp; + } + incipient[c] = bubble ? K[c] * z[c] : z[c] / K[c]; + sum += incipient[c]; + } + if (!positiveFinite_(sum)) { + return Outcome::GaveUp; + } + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(incipientPhaseIdx, c, incipient[c] / sum); + } + + // The known phase is fixed throughout the substitution at this pressure. + if (inner == 0) { + paramCache.updatePhase(fs, knownPhaseIdx); + for (int c = 0; c < numComponents; ++c) { + phiKnown[c] = FluidSystem::fugacityCoefficient( + fs, paramCache, knownPhaseIdx, c); + } + } + const int unchanged = inner == 0 + ? ParameterCache::None : ParameterCache::Temperature | ParameterCache::Pressure; + paramCache.updatePhase(fs, incipientPhaseIdx, unchanged); + + // Pure fluids and azeotropes can have K == 1 at saturation. + // Distinct molar volumes distinguish them from a trivial trial. + const Scalar vmL = paramCache.molarVolume(oilPhaseIdx); + const Scalar vmV = paramCache.molarVolume(gasPhaseIdx); + // A molar volume near the cubic EOS floor does not establish a + // physical second root. + constexpr Scalar clampedVm = 1.0e-7; + rootsDistinct = (std::min(vmL, vmV) > 2.0 * clampedVm) && + (std::abs(vmL - vmV) + > rootVolumeTolerance_ * std::max(vmL, vmV)); + + Scalar change = 0.0; + trivial = true; + CompVec d; + for (int c = 0; c < numComponents; ++c) { + const Scalar phiIncipient = FluidSystem::fugacityCoefficient( + fs, paramCache, incipientPhaseIdx, c); + const Scalar newK = bubble ? phiKnown[c] / phiIncipient + : phiIncipient / phiKnown[c]; + if (!positiveFinite_(newK)) { + return Outcome::GaveUp; + } + // Relative to the magnitude of K: an absolute measure is + // unreachable for the large K of a light component. + change = std::max(change, std::abs(newK - K[c]) / + std::max(Scalar{1}, std::abs(newK))); + trivial = trivial && (std::abs(newK - 1.0) < 1.0e-5); + d[c] = std::log(newK) - std::log(K[c]); + K[c] = newK; + } + if (change < nonlinearTolerance_) { + substitutionConverged = true; + break; + } + + // Every fourth pass, extrapolate the remaining geometric + // series in ln K using the last two substitution changes. + if (inner % 4 == 3) { + accelerate_(K, d, dPrev); + } + dPrev = d; + } + + const auto markOutside = [&](const bool withValue, const Scalar value) { + pOut = p; + haveOut = true; + haveFOut = withValue; + fOut = value; + }; + // Step from a trial the substitution could not classify: bisect + // the bracket if there is one, otherwise continue the scan. + const auto stepUnclassified = [&]() { + p = haveIn ? std::sqrt(pOut * pIn) + : p * (fromAbove ? scanStep : Scalar{1} / scanStep); + K = wilsonK(p); + }; + + // The looser triviality tolerance can be met before K converges. + // Only a converged trial may update the scan endpoint. + if (substitutionConverged && trivial && !rootsDistinct) { + // No nontrivial trial was found. Advance the scan endpoint; + // collapse onto one root does not establish phase stability. + markOutside(false, Scalar{0}); + lastSide = 0; + stepUnclassified(); + continue; + } + + Scalar sum = 0.0; + for (int c = 0; c < numComponents; ++c) { + sum += bubble ? K[c] * z[c] : z[c] / K[c]; + } + if (!positiveFinite_(sum)) { + return Outcome::GaveUp; + } + // The pressure criterion is only meaningful once the fixed-pressure + // substitution has actually reached fugacity equality; exhausting + // the inner loop is a failure, not a solution. + if (substitutionConverged + && std::abs(sum - 1.0) < boundaryResidualTolerance_) { + Scalar distance = 0.0; + for (int c = 0; c < numComponents; ++c) { + incipient[c] = (bubble ? K[c] * z[c] : z[c] / K[c]) / sum; + distance += std::abs(incipient[c] - z[c]); + } + // Require a distinct phase, then check known-phase stability: + // fugacity equality also admits stationary points inside the + // two-phase region, including nearly identical compositions + // occupying different EOS roots. + if (distance > 1.0e-3 || rootsDistinct) { + const auto stability = + knownPhaseStability_(fs, z, knownPhaseIdx, wilsonK(p), bubble, eosType); + if (stability == Stability::Indeterminate) { + return Outcome::GaveUp; + } + if (stability == Stability::Stable) { + press = p; + return Outcome::Converged; + } + } + // Reject this stationary point and continue searching. + markOutside(false, Scalar{0}); + lastSide = 0; + stepUnclassified(); + continue; + } + + // The incipient amount is the residual of the search: ln(sum) is + // positive inside the two-phase region and negative outside it, + // and the root is where it vanishes. + const Scalar f = std::log(sum); + const Scalar lnp = std::log(p); + + if (!substitutionConverged) { + // An exhausted substitution classifies nothing. Take the + // fixed-point step K ~ 1/p implies, p * sum towards the root, + // and let the carried-over K keep converging within the bracket. + const Scalar lnpFixed = lnp + (fromAbove ? f : -f); + p = safeguardPressure(std::exp(lnp + std::clamp(lnpFixed - lnp, + std::log(Scalar{0.5}), + std::log(Scalar{2})))); + continue; + } + + if (f > 0) { + // Inside the two-phase region. Illinois: replacing the same + // end twice in a row halves the value kept at the other end, + // so regula falsi cannot stall against it. + if (haveIn && lastSide > 0 && haveFOut) { + fOut *= 0.5; + } + pIn = p; + fIn = f; + haveIn = true; + lastSide = 1; + } + else { + if (haveOut && haveFOut && lastSide < 0 && haveIn) { + fIn *= 0.5; + } + markOutside(true, f); + lastSide = -1; + } + + // Propose the next pressure in ln p. Regula falsi within a bracket + // whose ends both carry a value; a secant through the previous + // trial before the bracket closes, provided the slope has the sign + // the branch implies (sum grows towards the root: falling pressure + // from above, rising pressure from below); and the fixed-point + // step otherwise. + Scalar lnpNext{}; + bool proposed = false; + if (haveIn && haveOut && haveFOut) { + const Scalar lnpIn = std::log(pIn); + const Scalar lnpOut = std::log(pOut); + if (fIn != fOut) { + lnpNext = lnpOut - fOut * (lnpIn - lnpOut) / (fIn - fOut); + proposed = true; + } + } + if (!proposed && havePrev && lnp != lnpPrev) { + const Scalar slope = (f - fPrev) / (lnp - lnpPrev); + const bool slopeOk = fromAbove ? (slope < -nonlinearTolerance_) + : (slope > nonlinearTolerance_); + if (slopeOk) { + lnpNext = lnp - f / slope; + proposed = true; + } + } + if (!proposed) { + lnpNext = lnp + (fromAbove ? f : -f); + } + lnpPrev = lnp; + fPrev = f; + havePrev = true; + + // Limit extrapolation to a factor of two and stay inside a known + // bracket. This step limit alone cannot prevent skipping an + // unbracketed narrow envelope. + lnpNext = lnp + std::clamp(lnpNext - lnp, std::log(Scalar{0.5}), std::log(Scalar{2})); + p = safeguardPressure(std::exp(lnpNext)); + } + + // Neither a finite scan nor collapse onto a trivial stationary point + // proves that an upper root is absent. + return Outcome::GaveUp; + } +}; + +} // namespace Opm + +#endif // OPM_SATURATION_PRESSURE_HPP diff --git a/tests/material/test_saturation_pressure.cpp b/tests/material/test_saturation_pressure.cpp new file mode 100644 index 00000000000..e79ba6b2ab6 --- /dev/null +++ b/tests/material/test_saturation_pressure.cpp @@ -0,0 +1,872 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2026 SINTEF Digital + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief Tests for the SaturationPressure constraint solver. + * + * The reference values come from a 1D vertical compositional equilibration + * case with three components (CO2, methane and n-decane, Peng-Robinson, zero + * binary interaction coefficients) at a constant reservoir temperature of + * 100 degC. The bubble-point pressures are the PSAT values its restart file + * reports for the single-phase oil cells, and the gas-oil contact values come + * from its equilibration report. + */ +#include "config.h" + +#define BOOST_TEST_MODULE SaturationPressure +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +using Scalar = double; +constexpr int numComponents = 3; + +using FluidSystem = Opm::GenericOilGasWaterFluidSystem; +using SatP = Opm::SaturationPressure; +using CompVec = SatP::CompVec; + +using FloatFluidSystem = Opm::GenericOilGasWaterFluidSystem; +using FloatSatP = Opm::SaturationPressure; + +constexpr auto eosType = Opm::CompositionalConfig::EOSType::PR; + +// The constant reservoir temperature (RTEMP) of the test case, 100 degC. +constexpr Scalar temperature = 373.15; + +// Initialize the component properties from the test deck once for all tests. +struct Fixture +{ + Fixture() + { + using CompParam = FluidSystem::ComponentParam; + FluidSystem::init(); + FluidSystem::addComponent(CompParam{"CO2", 44.0, 304.128, 73.773e5, 0.09412, 0.22394}); + FluidSystem::addComponent(CompParam{"C1", 16.04, 190.564, 45.992e5, 0.09863, 0.01142}); + FluidSystem::addComponent(CompParam{"C10", 142.28, 617.7, 21.03e5, 0.60980, 0.4884}); + + using FloatCompParam = FloatFluidSystem::ComponentParam; + FloatFluidSystem::init(); + FloatFluidSystem::addComponent( + FloatCompParam{"CO2", 44.0, 304.128, 73.773e5, 0.09412, 0.22394}); + FloatFluidSystem::addComponent( + FloatCompParam{"C1", 16.04, 190.564, 45.992e5, 0.09863, 0.01142}); + FloatFluidSystem::addComponent( + FloatCompParam{"C10", 142.28, 617.7, 21.03e5, 0.60980, 0.4884}); + } +}; + +// Equilibrium checks and composition distance. A zero distance can indicate a +// trivial solution; pure-fluid saturation also requires checking the EOS roots. +struct EquilibriumResidual +{ + Scalar fugacity{}; + Scalar closure{}; + Scalar distance{}; +}; + +template +EquilibriumResidual equilibriumResidualFor( + const std::array& known, + const unsigned knownPhaseIdx, + const std::array& incipient, + const unsigned incipientPhaseIdx, + const typename FS::Scalar press, + const typename FS::Scalar temp) +{ + using Value = typename FS::Scalar; + constexpr int numComp = FS::numComponents; + Opm::CompositionalFluidState fs; + fs.setTemperature(temp); + fs.setPressure(FS::oilPhaseIdx, press); + fs.setPressure(FS::gasPhaseIdx, press); + for (int c = 0; c < numComp; ++c) { + fs.setMoleFraction(knownPhaseIdx, c, known[c]); + fs.setMoleFraction(incipientPhaseIdx, c, incipient[c]); + } + + typename FS::template ParameterCache paramCache(eosType); + paramCache.updatePhase(fs, FS::oilPhaseIdx); + paramCache.updatePhase(fs, FS::gasPhaseIdx); + + EquilibriumResidual res; + Value sum = 0.0; + for (int c = 0; c < numComp; ++c) { + const Value phiL = + FS::fugacityCoefficient(fs, paramCache, FS::oilPhaseIdx, c); + const Value phiV = + FS::fugacityCoefficient(fs, paramCache, FS::gasPhaseIdx, c); + const Value fL = fs.moleFraction(FS::oilPhaseIdx, c) * phiL; + const Value fV = fs.moleFraction(FS::gasPhaseIdx, c) * phiV; + const Value scale = std::max({std::abs(fL), std::abs(fV), Value{1.0e-12}}); + res.fugacity = std::max(res.fugacity, + static_cast(std::abs(fL - fV) / scale)); + + sum += incipient[c]; + res.distance += std::abs(incipient[c] - known[c]); + } + res.closure = std::abs(sum - 1.0); + return res; +} + +EquilibriumResidual equilibriumResidual(const CompVec& known, + const unsigned knownPhaseIdx, + const CompVec& incipient, + const unsigned incipientPhaseIdx, + const Scalar press) +{ + return equilibriumResidualFor(known, knownPhaseIdx, incipient, + incipientPhaseIdx, press, temperature); +} + +} // anonymous namespace + +BOOST_GLOBAL_FIXTURE(Fixture); + +BOOST_AUTO_TEST_CASE(BubblePressureOilZone) +{ + // The single-phase oil cells of the test case. The mixtures are binary + // methane/decane (the CO2 fraction is zero); the expected bubble-point + // pressures are the reference PSAT values in bar. + constexpr std::array, 10> refValues{{ + {0.49, 156.29472}, + {0.47, 147.92114}, + {0.45, 139.75455}, + {0.43, 131.79112}, + {0.41, 124.02647}, + {0.39, 116.45577}, + {0.37, 109.07387}, + {0.35, 101.87540}, + {0.33, 94.85495}, + {0.31, 88.00698}, + }}; + + for (const auto& [zMethane, expectedBar] : refValues) { + const CompVec liquid{0.0, zMethane, 1.0 - zMethane}; + Scalar press = 0.0; + CompVec vapor{}; + + const bool converged = + SatP::bubblePressure(liquid, temperature, eosType, press, vapor); + + BOOST_REQUIRE_MESSAGE(converged, + "bubble-point iteration must converge for z_C1 = " << zMethane); + // The restart file stores PSAT in single precision; 1e-3 percent + // (1e-5 relative) is well above that quantization. + BOOST_CHECK_CLOSE(press / 1.0e5, expectedBar, 1.0e-3); + } +} + +BOOST_AUTO_TEST_CASE(CompositionNeedNotSumToOne) +{ + // Callers assemble z in floating point -- from a flash, or from a restart + // file's single-precision mole fractions -- so it need not sum exactly to + // one. The boundary residual is measured against one, so without + // rescaling a composition summing to 1 + d stops converging as soon as + // |d| exceeds that tolerance. + constexpr Scalar zMethane = 0.47; + constexpr Scalar referenceBar = 147.92114; + + for (const Scalar excess : {Scalar{-1.0e-4}, Scalar{-3.0e-8}, Scalar{0.0}, + Scalar{3.0e-8}, Scalar{1.0e-4}}) { + const CompVec liquid{0.0, zMethane, 1.0 - zMethane + excess}; + + Scalar press = 0.0; + CompVec vapor{}; + BOOST_REQUIRE_MESSAGE(SatP::bubblePressure(liquid, temperature, eosType, press, vapor), + "bubble point must converge for a composition summing to " + << 1.0 + excess); + + // Rescaling moves the mixture itself, so the bubble point moves with + // it: 1e-4 of excess is worth about 0.02 bar here. + BOOST_CHECK_CLOSE(press / 1.0e5, referenceBar, 3.0e-2); + + // The incipient vapour is normalised whatever the input was. + Scalar sum = 0.0; + for (int c = 0; c < numComponents; ++c) { + sum += vapor[c]; + } + BOOST_CHECK_CLOSE(sum, 1.0, 1.0e-6); + } + + // The case that arises in practice: 0.47 and 0.53 sum to 1 - 3e-8 once + // they have been through single precision. + const CompVec rounded{0.0, static_cast(0.47f), static_cast(0.53f)}; + Scalar pRounded = 0.0; + CompVec roundedVapor{}; + BOOST_REQUIRE(SatP::bubblePressure(rounded, temperature, eosType, pRounded, roundedVapor)); + BOOST_CHECK_CLOSE(pRounded / 1.0e5, referenceBar, 1.0e-3); + + // The dew search rescales the same way. + Scalar pBubble = 0.0; + CompVec vapour{}; + BOOST_REQUIRE(SatP::bubblePressure({0.0, 0.40, 0.60}, temperature, eosType, + pBubble, vapour)); + vapour.back() += 3.0e-8; + + Scalar pDew = -1.0; + CompVec incipient{}; + BOOST_REQUIRE_MESSAGE(SatP::dewPressure(vapour, temperature, eosType, pDew, incipient), + "dew point must converge for a vapour that does not sum to one"); + BOOST_CHECK_CLOSE(pDew, pBubble, 1.0e-3); +} + +BOOST_AUTO_TEST_CASE(SaturationPressureAtGasOilContact) +{ + // At the gas-oil contact the liquid composition from ZMFVD is + // (0, 0.5, 0.5), and the reference equilibration sets the contact pressure + // to a saturation pressure of 160.56010 bar. The incipient vapour becomes + // the gas-cap composition, reported as ZMF = (0, 0.987784, 0.012216) in + // the restart file. + const CompVec liquid{0.0, 0.5, 0.5}; + Scalar press = 0.0; + CompVec vapor{}; + + const bool converged = + SatP::bubblePressure(liquid, temperature, eosType, press, vapor); + + BOOST_REQUIRE(converged); + BOOST_CHECK_CLOSE(press / 1.0e5, 160.56010, 1.0e-3); + + BOOST_CHECK_SMALL(vapor[0], 1.0e-10); + BOOST_CHECK_CLOSE(vapor[1], 0.987784, 1.0e-2); + BOOST_CHECK_CLOSE(vapor[2], 0.012216, 1.0e-1); +} + +BOOST_AUTO_TEST_CASE(DewPressureGasCap) +{ + // Thermodynamic consistency at the contact: the gas cap is a retrograde + // condensate, and its upper (retrograde) dew-point pressure must recover + // the contact pressure, with the incipient liquid recovering the ZMFVD + // composition at the contact. The vapour composition is only known to + // single precision, which limits the achievable agreement; the tolerances + // reflect that. + const CompVec vapor{0.0, 0.987784, 0.012216}; + Scalar press = 0.0; + CompVec liquid{}; + + const bool converged = + SatP::dewPressure(vapor, temperature, eosType, press, liquid); + + BOOST_REQUIRE(converged); + BOOST_CHECK_CLOSE(press / 1.0e5, 160.56010, 1.0e-2); + + BOOST_CHECK_SMALL(liquid[0], 1.0e-10); + BOOST_CHECK_CLOSE(liquid[1], 0.5, 0.1); + BOOST_CHECK_CLOSE(liquid[2], 0.5, 0.1); + + // The pair must also satisfy the equilibrium conditions in its own right. + const auto res = equilibriumResidual(vapor, FluidSystem::gasPhaseIdx, + liquid, FluidSystem::oilPhaseIdx, press); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-6); + BOOST_CHECK_GT(res.distance, 1.0e-3); +} + +BOOST_AUTO_TEST_CASE(BubbleDewRoundTrip) +{ + // The vapour that appears at the bubble point of a liquid is, at that same + // temperature and pressure, a vapour at its own dew point with the liquid + // as the incipient phase. The dew-point search must therefore return the + // bubble-point pressure it was handed and recover the liquid composition. + // + // These cases exercise the upper-dew branch selected by dewPressure(). + for (const Scalar zMethane : {Scalar{0.30}, Scalar{0.35}, Scalar{0.40}, Scalar{0.45}}) { + const CompVec liquid{0.0, zMethane, 1.0 - zMethane}; + + Scalar pBubble = 0.0; + CompVec vapor{}; + BOOST_REQUIRE_MESSAGE(SatP::bubblePressure(liquid, temperature, eosType, pBubble, vapor), + "bubble point must converge for z_C1 = " << zMethane); + + Scalar pDew = -1.0; + CompVec incipient{}; + BOOST_REQUIRE_MESSAGE(SatP::dewPressure(vapor, temperature, eosType, pDew, incipient), + "dew point must converge for the vapour of z_C1 = " << zMethane); + + BOOST_CHECK_CLOSE(pDew, pBubble, 1.0e-6); + for (int c = 0; c < numComponents; ++c) { + BOOST_CHECK_SMALL(std::abs(incipient[c] - liquid[c]), 1.0e-6); + } + } +} + +BOOST_AUTO_TEST_CASE(SinglePrecisionBubbleDewRoundTrip) +{ + constexpr FloatSatP::CompVec liquid{0.0f, 0.4f, 0.6f}; + float pBubble = -1.0f; + FloatSatP::CompVec vapor{}; + BOOST_REQUIRE(FloatSatP::bubblePressure(liquid, 373.15f, eosType, pBubble, vapor)); + + Scalar pBubbleDouble = -1.0; + CompVec vaporDouble{}; + BOOST_REQUIRE(SatP::bubblePressure({0.0, 0.4, 0.6}, temperature, eosType, + pBubbleDouble, vaporDouble)); + BOOST_CHECK_CLOSE(static_cast(pBubble), pBubbleDouble, 1.0e-3); + + float pDew = -1.0f; + FloatSatP::CompVec recovered{}; + BOOST_REQUIRE(FloatSatP::dewPressure(vapor, 373.15f, eosType, pDew, recovered)); + BOOST_CHECK_CLOSE(pDew, pBubble, 1.0e-3); + for (int c = 0; c < numComponents; ++c) { + BOOST_CHECK_SMALL(recovered[c] - liquid[c], 2.0e-5f); + } + const auto res = equilibriumResidualFor( + vapor, FloatFluidSystem::gasPhaseIdx, + recovered, FloatFluidSystem::oilPhaseIdx, pDew, 373.15f); + BOOST_CHECK_SMALL(res.closure, 2.0e-5); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-3); +} + +BOOST_AUTO_TEST_CASE(DewPressureLowerBranch) +{ + // These mixtures have an ordinary envelope ending at a bubble boundary. + // Reject unstable upper-branch stationary points near 155.6 and 314.0 bar, + // respectively, and verify the lower dew point through fugacity equality. + for (const Scalar yMethane : {Scalar{0.70}, Scalar{0.85}}) { + const CompVec vapor{0.0, yMethane, 1.0 - yMethane}; + Scalar press = 0.0; + CompVec liquid{}; + + BOOST_REQUIRE_MESSAGE(SatP::dewPressure(vapor, temperature, eosType, press, liquid), + "dew point must converge for y_C1 = " << yMethane); + BOOST_CHECK_GT(press, 0.0); + BOOST_CHECK_LT(press, 5.0e5); + + const auto res = equilibriumResidual(vapor, FluidSystem::gasPhaseIdx, + liquid, FluidSystem::oilPhaseIdx, press); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-6); + // The incipient liquid is nearly pure decane: a different phase. + BOOST_CHECK_GT(res.distance, 1.0); + BOOST_CHECK_GT(liquid[2], 0.99); + } +} + +BOOST_AUTO_TEST_CASE(SupercriticalLiquidHasNoBubblePoint) +{ + // Pure methane is far above its critical temperature here, so no bubble + // point exists at any pressure. The solver must refuse and leave the + // pressure output untouched. + const CompVec liquid{0.0, 1.0, 0.0}; + Scalar press = -1.0; + CompVec vapor{}; + + BOOST_CHECK(!SatP::bubblePressure(liquid, temperature, eosType, press, vapor)); + BOOST_CHECK_LT(press, 0.0); +} + +BOOST_AUTO_TEST_CASE(PureComponentSaturationPressure) +{ + // Pure decane is well below its critical temperature here, so it has a + // genuine vapour pressure even though both phases necessarily share its + // composition. What tells this state from the trivial solution is that + // the phases occupy distinct EOS roots, not that their compositions + // differ; a solver that requires a composition difference refuses it. + const CompVec liquid{0.0, 0.0, 1.0}; + Scalar pBubble = 0.0; + CompVec vapor{}; + + BOOST_REQUIRE(SatP::bubblePressure(liquid, temperature, eosType, pBubble, vapor)); + BOOST_CHECK_GT(pBubble, 1.0e2); + BOOST_CHECK_LT(pBubble, 1.0e5); + for (int c = 0; c < numComponents; ++c) { + BOOST_CHECK_SMALL(std::abs(vapor[c] - liquid[c]), 1.0e-6); + } + + // The equilibrium condition holds even though the compositions coincide. + const auto res = equilibriumResidual(liquid, FluidSystem::oilPhaseIdx, + vapor, FluidSystem::gasPhaseIdx, pBubble); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-6); + BOOST_CHECK_SMALL(res.distance, 1.0e-6); + + // For a pure component the dew and bubble pressures are one and the same. + Scalar pDew = 0.0; + CompVec incipient{}; + BOOST_REQUIRE(SatP::dewPressure(liquid, temperature, eosType, pDew, incipient)); + BOOST_CHECK_CLOSE(pDew, pBubble, 1.0e-3); +} + +BOOST_AUTO_TEST_CASE(TrivialSolutionIsNotReportedAsADewPoint) +{ + // For these methane-rich mixtures, successive substitution can converge to + // trivial K == 1 states at high pressure. Any returned root must have a + // distinct incipient-phase composition. + for (const Scalar zMethane : {Scalar{0.999}, Scalar{0.9999}}) { + const CompVec vapor{0.0, zMethane, 1.0 - zMethane}; + Scalar press = -1.0; + CompVec liquid{}; + + const bool converged = + SatP::dewPressure(vapor, temperature, eosType, press, liquid); + + if (converged) { + // If a root is reported at all it must be a real one, i.e. a + // distinct incipient phase rather than a copy of the vapour. + const auto res = equilibriumResidual(vapor, FluidSystem::gasPhaseIdx, + liquid, FluidSystem::oilPhaseIdx, press); + BOOST_CHECK_MESSAGE(res.distance > 1.0e-3, + "trivial solution reported as a dew point for z_C1 = " + << zMethane << " at " << press / 1.0e5 << " bar"); + } + else { + // On failure, the pressure output must remain unchanged. + BOOST_CHECK_LT(press, 0.0); + } + } +} + +BOOST_AUTO_TEST_CASE(NearCriticalBubblePoints) +{ + // Liquids approaching the critical composition of the mixture. The + // successive substitution contracts ever more slowly there and the + // saturation sum flattens against one, which defeated a fixed-point + // pressure update within its iteration budget. The reference values are + // this implementation's converged results, cross-checked against a scan + // of the saturation sum over pressure at fixed composition. + constexpr std::array, 3> cases{{ + {0.80, 306.1289597e5}, + {0.82, 314.7255712e5}, + {0.834, 320.0963847e5}, + }}; + for (const auto& [xMethane, pExpected] : cases) { + const CompVec liquid{0.0, xMethane, 1.0 - xMethane}; + Scalar press = 0.0; + CompVec vapor{}; + + BOOST_REQUIRE_MESSAGE(SatP::bubblePressure(liquid, temperature, eosType, press, vapor), + "bubble point must converge for x_C1 = " << xMethane); + BOOST_CHECK_CLOSE(press, pExpected, 1.0e-4); + + const auto res = equilibriumResidual(liquid, FluidSystem::oilPhaseIdx, + vapor, FluidSystem::gasPhaseIdx, press); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_GT(res.distance, 1.0e-2); + } +} + +BOOST_AUTO_TEST_CASE(RetrogradeDewPointOfBubblePointVapour) +{ + // The vapour in equilibrium with a decane-rich liquid at its bubble point + // is a lean gas with two dew points: the lower one at the bubble pressure + // it came from, and a retrograde one above it. dewPressure() must report + // the retrograde point. It is approached from the single-phase side, + // where the saturation sum rises towards one so slowly that a fixed-point + // pressure update ran out of iterations short of it. + const CompVec liquid{0.0, 0.17, 0.83}; + Scalar pBubble = 0.0; + CompVec vapor{}; + BOOST_REQUIRE(SatP::bubblePressure(liquid, temperature, eosType, pBubble, vapor)); + BOOST_CHECK_CLOSE(pBubble, 44.4490922e5, 1.0e-4); + + Scalar pDew = 0.0; + CompVec incipient{}; + BOOST_REQUIRE(SatP::dewPressure(vapor, temperature, eosType, pDew, incipient)); + BOOST_CHECK_CLOSE(pDew, 57.6521308e5, 1.0e-4); + BOOST_CHECK_GT(pDew, pBubble); + + const auto res = equilibriumResidual(vapor, FluidSystem::gasPhaseIdx, + incipient, FluidSystem::oilPhaseIdx, pDew); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_GT(res.distance, 1.0e-2); +} + +BOOST_AUTO_TEST_CASE(NearCriticalRetrogradeDewPoint) +{ + // A vapour just richer in methane than the critical composition. Its + // retrograde dew point lies within a bar or two of the critical pressure + // (the bubble curve peaks at 331.2 bar near x_C1 = 0.88), where the two + // phases differ in composition by a few percent only. The point is + // accepted once the stability certificate converges, which the + // accelerated substitution makes possible this close to the critical + // point. + const CompVec vapor{0.0, 0.90, 0.10}; + Scalar press = 0.0; + CompVec liquid{}; + BOOST_REQUIRE(SatP::dewPressure(vapor, temperature, eosType, press, liquid)); + BOOST_CHECK_GT(press, 300.0e5); + BOOST_CHECK_LT(press, 340.0e5); + + const auto res = equilibriumResidual(vapor, FluidSystem::gasPhaseIdx, + liquid, FluidSystem::oilPhaseIdx, press); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_GT(res.distance, 1.0e-2); +} + +BOOST_AUTO_TEST_CASE(RetrogradeDewPointWithBinaryInteraction) +{ + // Seven-component gas from a reference equilibration run at 120 degC. + // Exercise nonzero binary interaction coefficients and their parsed + // lower-triangle ordering. The retrograde dew point is about 120 bar + // below the reservoir pressure, with distinct phase compositions. + using FluidSystem7 = Opm::GenericOilGasWaterFluidSystem; + using SatP7 = Opm::SaturationPressure; + + const auto deck = Opm::Parser{}.parseString(R"( +RUNSPEC +METRIC +DIMENS + 1 1 1 / +COMPS +7 / +TABDIMS + 1 / +OIL +GAS +GRID +DXV + 1 / +DYV + 1 / +DZV + 1 / +DEPTHZ + 4*2500 / +PROPS +EOS + PR / +CNAMES + CH4 CO2 C2H6 C3H8 C4-C9 C10-C20 C21+ / +TCRIT + 190.6 304.2 305.4 369.8 474.7 646.0 812.0 / +PCRIT + 45.40 72.80 48.20 41.90 32.40 20.42 14.50 / +VCRIT + 0.099 0.094 0.148 0.203 0.32321204 0.63255728 0.97988801 / +MW + 16.043 44.010 30.070 44.097 76.900 163.00 310.85 / +ACF + 0.008 0.225 0.098 0.152 0.2549 0.551 0.909 / +BIC + 1.0500000E-01 + 2.6890022E-03 1.3000000E-01 + 8.5370405E-03 1.2500000E-01 1.6620489E-03 + 2.2915937E-02 0.0000000E+00 1.0088696E-02 3.5952572E-03 + 5.4875390E-02 0.0000000E+00 3.4227722E-02 2.1174720E-02 7.4706799E-03 + 8.1972182E-02 0.0000000E+00 5.6906187E-02 4.0015457E-02 2.0180685E-02 3.1846386E-03 / +SOLUTION +SCHEDULE +END +)"); + + const auto eclState = Opm::EclipseState{ deck }; + const auto schedule = Opm::Schedule{ deck, eclState }; + FluidSystem7::init(); + BOOST_REQUIRE_NO_THROW(FluidSystem7::initFromState(eclState, schedule)); + + // The composition of the gas zone in the reference data. + constexpr SatP7::CompVec vapor{0.885830, 0.043700, 0.034000, 0.018900, + 0.015200, 0.002300, 0.000070}; + constexpr Scalar temp = 393.15; // RTEMP of the reference case, 120 degC + constexpr Scalar pReference = 151.7060; // its reported PSAT, in bar + constexpr Scalar pReservoir = 272.6396; // and the cell pressure there + + Scalar press = 0.0; + SatP7::CompVec liquid{}; + BOOST_REQUIRE(SatP7::dewPressure(vapor, temp, eosType, press, liquid)); + + // The restart file stores PSAT in single precision; 1e-3 percent is well + // above that quantization. + BOOST_CHECK_CLOSE(press / 1.0e5, pReference, 1.0e-3); + + // The retrograde branch, not a point next to the reservoir pressure. + BOOST_CHECK_LT(press / 1.0e5, 0.7 * pReservoir); + + // The incipient phase is a genuine liquid: heavy where the gas is light. + BOOST_CHECK_GT(liquid[6], 100.0 * vapor[6]); + BOOST_CHECK_LT(liquid[0], 0.5 * vapor[0]); + + const auto res = equilibriumResidualFor( + vapor, FluidSystem7::gasPhaseIdx, liquid, FluidSystem7::oilPhaseIdx, press, temp); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_GT(res.distance, 1.0); +} + +BOOST_AUTO_TEST_CASE(RetrogradeDewPointOfANearCriticalGas) +{ + // Near criticality a flat trial total does not establish stationarity. + // Finish the stability trial and return the upper point near 331.4 bar, + // rather than falling through to the lower point near 0.961 bar. + const CompVec vapor{0.0, 0.895, 0.105}; + Scalar press = -123.0; + CompVec liquid{}; + + BOOST_REQUIRE(SatP::dewPressure(vapor, temperature, eosType, press, liquid)); + BOOST_CHECK_CLOSE(press / 1.0e5, 331.3565338, 1.0e-3); + + const auto res = equilibriumResidual(vapor, FluidSystem::gasPhaseIdx, + liquid, FluidSystem::oilPhaseIdx, press); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_GT(res.distance, 1.0e-3); +} + +BOOST_AUTO_TEST_CASE(NarrowRetrogradeEnvelopeKeepsTheUpperDewPoint) +{ + // CO2/methane at 200 K has two closely spaced dew points for these gas + // compositions. Generate each gas from its equilibrium liquid at the + // upper boundary, independently of the dew search. The old scan either + // failed or returned the lower dew point: for x_CO2 = 0.068 it returned + // 51.20128 instead of 51.50684 bar. At x_CO2 = 0.070 the interval is only + // about 0.056 bar wide, much less than a 10% pressure step. + constexpr Scalar temp = 200.0; + for (const Scalar xCO2 : {Scalar{0.068}, Scalar{0.070}, Scalar{0.065}}) { + const CompVec liquid{xCO2, 1.0 - xCO2, 0.0}; + CompVec vapor{}; + Scalar pBubble{}; + BOOST_REQUIRE(SatP::bubblePressure(liquid, temp, eosType, pBubble, vapor)); + + Scalar pDew = -123.0; + CompVec recovered{}; + BOOST_REQUIRE_MESSAGE(SatP::dewPressure(vapor, temp, eosType, pDew, recovered), + "upper dew point must converge for x_CO2 = " << xCO2); + BOOST_CHECK_CLOSE(pDew, pBubble, 1.0e-4); + for (int c = 0; c < numComponents; ++c) { + BOOST_CHECK_SMALL(recovered[c] - liquid[c], 1.0e-6); + } + const auto res = equilibriumResidualFor( + vapor, FluidSystem::gasPhaseIdx, recovered, FluidSystem::oilPhaseIdx, pDew, temp); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_GT(res.distance, 1.0e-2); + } +} + +BOOST_AUTO_TEST_CASE(OrdinaryDewPointAfterFollowingTheBubbleBoundary) +{ + // Here the upper dew search and the direct bubble search both struggle. + // Continuation must follow the vapour-like trial to the bubble boundary + // before returning the lower dew point, including the near-critical + // stability check for the ternary mixture. + for (const auto& [temp, vapor, expectedBar] : + {std::tuple{440.0, CompVec{0.0, 0.8, 0.2}, 4.587166726}, + std::tuple{570.0, CompVec{0.25, 0.3, 0.45}, 38.68411340}}) { + Scalar press = -123.0; + CompVec liquid{}; + BOOST_REQUIRE(SatP::dewPressure(vapor, temp, eosType, press, liquid)); + BOOST_CHECK_CLOSE(press / 1.0e5, expectedBar, 1.0e-4); + const auto res = equilibriumResidualFor( + vapor, FluidSystem::gasPhaseIdx, liquid, FluidSystem::oilPhaseIdx, press, temp); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_GT(res.distance, 0.1); + } +} + +BOOST_AUTO_TEST_CASE(AcceleratedDewSearchKeepsFiniteIterates) +{ + // Guard against K underflow producing a NaN pressure and an EOS exception. + // An independent stability scan locates the first mixture's envelope at + // 18.44153--22.72437 bar; the other four remain single-phase over 0.5--600 bar. + struct Case + { + Scalar temp; + CompVec vapor; + Scalar expectedBar; // zero when the independent scan found no dew point + }; + constexpr std::array cases{{ + {600.0, {0.05, 0.0, 0.95}, 18.44153}, + {600.0, {0.35, 0.0, 0.65}, 0.0}, + {610.0, {0.20, 0.0, 0.80}, 0.0}, + {610.0, {0.25, 0.0, 0.75}, 0.0}, + {610.0, {0.30, 0.0, 0.70}, 0.0}, + }}; + for (const auto& [temp, vapor, expectedBar] : cases) { + constexpr Scalar sentinel = -123.0; + Scalar press = sentinel; + CompVec liquid{}; + bool converged = false; + BOOST_REQUIRE_NO_THROW( + converged = SatP::dewPressure(vapor, temp, eosType, press, liquid)); + + if (expectedBar == 0.0) { + BOOST_CHECK_MESSAGE(!converged, + "single-phase mixture reported a dew point at " + << press / 1.0e5 << " bar"); + BOOST_CHECK_EQUAL(press, sentinel); + continue; + } + + BOOST_REQUIRE_MESSAGE(converged, + "dew point must be found for z_CO2 = " << vapor[0]); + BOOST_CHECK_CLOSE(press / 1.0e5, expectedBar, 1.0e-3); + + const auto res = equilibriumResidualFor( + vapor, FluidSystem::gasPhaseIdx, liquid, FluidSystem::oilPhaseIdx, press, temp); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_GT(res.distance, 1.0e-3); + } +} + +BOOST_AUTO_TEST_CASE(DewPointReachedFromTheBubbleBoundary) +{ + // Both dew scans miss these narrow CO2/decane envelopes near criticality. + // Trace down from the bubble boundary to find the dew point. Independent + // stability scans locate the envelopes at 24.26--26.86 and 22.40--23.64 bar. + struct Case + { + Scalar temp; + CompVec vapor; + Scalar expectedBar; + }; + constexpr std::array cases{{ + {612.0, {0.08, 0.0, 0.92}, 24.2564}, + {614.0, {0.04, 0.0, 0.96}, 22.3950}, + }}; + for (const auto& [temp, vapor, expectedBar] : cases) { + Scalar pDew = -123.0; + CompVec liquid{}; + BOOST_REQUIRE_MESSAGE(SatP::dewPressure(vapor, temp, eosType, pDew, liquid), + "dew point must be found for z_CO2 = " << vapor[0]); + BOOST_CHECK_CLOSE(pDew / 1.0e5, expectedBar, 1.0e-3); + + // It is the lower boundary of an ordinary envelope: the bubble point + // lies above it. + Scalar pBubble = 0.0; + CompVec bubbleVapor{}; + BOOST_REQUIRE(SatP::bubblePressure(vapor, temp, eosType, pBubble, bubbleVapor)); + BOOST_CHECK_LT(pDew, pBubble); + + const auto res = equilibriumResidualFor( + vapor, FluidSystem::gasPhaseIdx, liquid, FluidSystem::oilPhaseIdx, pDew, temp); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_GT(res.distance, 1.0e-3); + // The incipient liquid is the heavier phase. + BOOST_CHECK_GT(liquid[2], vapor[2]); + } +} + +BOOST_AUTO_TEST_CASE(PropertyProbeOverEosVariantsAndStates) +{ + // A deterministic sweep of the whole input domain. It asserts the three + // properties that hold for every input rather than any particular + // pressure: the call does not throw, a successful result is a usable + // saturation pressure with a normalised incipient composition, and a + // failed one leaves both outputs exactly as the caller left them. + constexpr std::array eosVariants{ + Opm::CompositionalConfig::EOSType::PR, + Opm::CompositionalConfig::EOSType::PRCORR, + Opm::CompositionalConfig::EOSType::RK, + Opm::CompositionalConfig::EOSType::SRK, + }; + + // Below, around and well above the critical temperatures of the mixture. + constexpr std::array temperatures{ + Scalar{195.0}, Scalar{279.0}, Scalar{363.0}, + Scalar{447.0}, Scalar{531.0}, Scalar{615.0}, + }; + + // Pure components, binaries and a zero-free interior, plus a reproducible + // spread of ternaries from a fixed seed. + std::vector mixtures{ + {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, + {0.5, 0.5, 0.0}, {0.5, 0.0, 0.5}, {0.0, 0.5, 0.5}, + {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0}, + }; + std::mt19937 gen{20260910}; + std::uniform_real_distribution unit{0.0, 1.0}; + while (mixtures.size() < 50) { + CompVec z{unit(gen), unit(gen), unit(gen)}; + const Scalar sum = z[0] + z[1] + z[2]; + if (sum <= 0.0) { + continue; + } + for (auto& zi : z) { + zi /= sum; + } + mixtures.push_back(z); + } + + // A sentinel no search would produce, so any write is visible. + constexpr Scalar pressSentinel = -12345.0; + const CompVec compSentinel{-1.0, -2.0, -3.0}; + + std::size_t calls = 0, converged = 0; + for (const auto eos : eosVariants) { + for (const Scalar temp : temperatures) { + for (const auto& z : mixtures) { + for (const bool bubble : {true, false}) { + Scalar press = pressSentinel; + CompVec incipient = compSentinel; + + bool ok = false; + BOOST_REQUIRE_NO_THROW( + ok = bubble ? SatP::bubblePressure(z, temp, eos, press, incipient) + : SatP::dewPressure(z, temp, eos, press, incipient)); + ++calls; + + if (!ok) { + // The failure contract: nothing was written. + BOOST_CHECK_EQUAL(press, pressSentinel); + BOOST_CHECK(incipient == compSentinel); + continue; + } + + ++converged; + BOOST_CHECK(std::isfinite(press)); + BOOST_CHECK_GT(press, 0.0); + + Scalar sum = 0.0; + for (const Scalar yi : incipient) { + BOOST_CHECK(std::isfinite(yi)); + BOOST_CHECK_GE(yi, 0.0); + sum += yi; + } + BOOST_CHECK_CLOSE(sum, 1.0, 1.0e-6); + } + } + } + } + + BOOST_CHECK_EQUAL(calls, eosVariants.size() * temperatures.size() * mixtures.size() * 2); + // The probe is worthless if nothing converges; it is a coverage floor, not + // a physical claim. + BOOST_CHECK_GT(converged, calls / 10); +} From ac6b4651838cad6a052715864667b897f9ccee02 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 11 Sep 2026 00:33:19 +0200 Subject: [PATCH 2/6] Resolve saturation pressure continuation gaps Complete stalled near-critical successive substitutions with the damped Newton stationary solve. Both methods solve the same fixed-pressure fugacity equations, so this changes the iteration method without changing the saturation criterion. Retain converged envelope brackets across bounded fixed-pressure failures and try dyadic interior points during refinement. Every candidate still passes EOS-root, fugacity, and phase-stability certification before it is returned. This removes all 245 interior dew-search gaps in the 8,400-state double-precision ternary sweep. The near-critical dew regression bounds the retrograde branch rather than a pressure value: that boundary is near critical, so the pressure moves with rounding while the branch it belongs to does not. --- .../constraintsolvers/SaturationPressure.hpp | 284 ++++++++++++------ tests/material/test_saturation_pressure.cpp | 43 +++ 2 files changed, 242 insertions(+), 85 deletions(-) diff --git a/opm/material/constraintsolvers/SaturationPressure.hpp b/opm/material/constraintsolvers/SaturationPressure.hpp index 3f35f34e8d6..c9d8180de17 100644 --- a/opm/material/constraintsolvers/SaturationPressure.hpp +++ b/opm/material/constraintsolvers/SaturationPressure.hpp @@ -102,16 +102,13 @@ namespace Opm { * one pressure to the next, and brackets the sign change of \f$\ln S_Y\f$. * A failed scan does not establish that no root exists: near a critical point, * the two-phase interval can be narrower than the pressure sampling step. - * Continuation narrows that gap without closing it. Sweeping a ternary across - * temperature and composition leaves the bubble search with no unresolved - * point between two resolved neighbours, but the dew search keeps isolated - * bands, a few parts in ten thousand of mole fraction wide, at the crossover - * between its lower and upper branches. Both branches are well separated and - * stable on either side of such a band, so a boundary does exist inside it. - * They survive because the lower-branch result is only accepted once a bubble - * point confirms an ordinary envelope, and a mixture gas-rich enough to sit at - * that crossover has no bubble point to find. Resolving them needs a bracketing - * scheme that cannot skip an unsampled interval, not a smaller fixed step. + * A damped Newton stationary solve completes fixed-pressure substitutions that + * stall near criticality. During continuation, a failed trial point does not + * discard an existing bracket: alternate interior points are tried, and a + * bounded forward search can cross an interval of numerical nonconvergence. + * If a nontrivial stationary branch terminates at a critical endpoint before + * producing a negative residual, its last point is subjected to the same EOS, + * fugacity and stability certification as a sign-bracketed boundary. * Input compositions must contain finite, non-negative mole fractions with a * positive total. They are rescaled to sum to one, so they need not do so on * entry. @@ -661,7 +658,16 @@ class SaturationPressure * the method walks in \f$q\f$, carries \f$Y\f$ between pressures, and grows * the step until \f$F\leq0\f$ brackets the other boundary. A safeguarded * secant method then solves \f$F=0\f$; trial amounts are interpolated in - * logarithmic space to preserve positivity and branch continuity. + * logarithmic space to preserve positivity and branch continuity. Failed + * trial evaluations are retried at dyadic interior points. If local + * continuation stalls, bounded forward probes look for the next evaluable + * point while retaining the last certified inside point. + * + * A nontrivial branch need not change sign before it ends. Where the + * interval closes at a critical point the incipient phase merges with the + * known one and the trial ceases to exist, so no step ever lands outside. + * The walk then contracts onto the boundary instead of bracketing it, and + * a residual already within tolerance is certified where it stands. * * The converged boundary is classified from enrichment relative to Wilson * volatility, @@ -670,8 +676,10 @@ class SaturationPressure * \f] * \f$D>0\f$ denotes a vapour-like incipient phase (bubble point), while * \f$D<0\f$ denotes a liquid-like incipient phase (dew point). Fugacity - * equality, distinct physical EOS roots, and known-phase stability are - * checked before returning the boundary. + * equality, physical EOS roots, and known-phase stability are checked + * before returning the boundary. The roots are not required to be + * distinct: they coincide at a critical endpoint, which is one of the + * boundaries this continuation exists to reach. * * From a lower dew seed the method walks upward: an upper dew point is * returned, whereas a bubble boundary confirms that the seed was the only @@ -741,7 +749,8 @@ class SaturationPressure } } Scalar sum{}; - if (!stationaryTrial_(fs, z, knownPhase, trialPhase, eosType, candidate, sum)) { + if (!stationaryTrial_(fs, z, knownPhase, trialPhase, + eosType, candidate, sum)) { continue; } converged = true; @@ -765,6 +774,70 @@ class SaturationPressure return true; }; + const auto certifyBoundary = [&](const Scalar lnp, CompVec Y, const Scalar f) { + const Scalar p = std::exp(lnp); + CompVec K = Kp; + Scalar direction = 0.0; + for (int c = 0; c < numComponents; ++c) { + K[c] /= p; + Y[c] /= std::exp(f); + direction += (Y[c] - z[c]) * std::log(K[c]); + } + // Identify the continued trial by its enrichment in the + // Wilson volatility direction, then verify fugacity equality + // with the corresponding liquid/vapour EOS roots. + if (std::abs(direction) < directionTolerance_) { + return false; + } + const bool bubble = direction > 0.0; + const unsigned knownPhase = bubble ? oilPhaseIdx : gasPhaseIdx; + const unsigned trialPhase = bubble ? gasPhaseIdx : oilPhaseIdx; + fs.setPressure(oilPhaseIdx, p); + fs.setPressure(gasPhaseIdx, p); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(knownPhase, c, z[c]); + fs.setMoleFraction(trialPhase, c, Y[c]); + } + ParameterCache cache(eosType); + cache.updatePhase(fs, oilPhaseIdx); + cache.updatePhase(fs, gasPhaseIdx); + const Scalar vmL = cache.molarVolume(oilPhaseIdx); + const Scalar vmV = cache.molarVolume(gasPhaseIdx); + if (!positiveFinite_(vmL) || !positiveFinite_(vmV) + || std::min(vmL, vmV) <= 2.0e-7) { + return false; + } + for (int c = 0; c < numComponents; ++c) { + if (z[c] > 0.0) { + const Scalar a = z[c] * FluidSystem::fugacityCoefficient( + fs, cache, knownPhase, c); + const Scalar b = Y[c] * FluidSystem::fugacityCoefficient( + fs, cache, trialPhase, c); + if (!positiveFinite_(a) || !positiveFinite_(b) + || std::abs(std::log(a / b)) > stabilityTolerance_) { + return false; + } + } + } + if (knownPhaseStability_(fs, z, knownPhase, K, bubble, eosType) + != Stability::Stable) { + return false; + } + if (fromBubble) { + // Down from a bubble point the envelope can only close + // at a dew point; anything else is not this envelope. + if (bubble) { + return false; + } + press = p; + liquid = Y; + return true; + } + press = bubble ? pSeed : p; + liquid = bubble ? seedTrial : Y; + return true; + }; + const Scalar seed = std::log(pSeed); Scalar lo{}, fLo{}; CompVec yLo{}; @@ -772,7 +845,9 @@ class SaturationPressure bool inside = false; // Halve the seed offset until a trial lands inside the envelope. constexpr int maxSeedRefinements = 20; - constexpr int maxWalkSteps = 200; + // With the 0.1 maximum step below, 400 iterations can span 40 units + // in ln(p), including a high-pressure seed and a sub-bar dew point. + constexpr int maxWalkSteps = 400; constexpr int maxBoundarySteps = 80; for (int refine = 0; refine < maxSeedRefinements; ++refine) { lo = seed + dir * step; @@ -796,6 +871,41 @@ class SaturationPressure if (!evaluate(hi, yHi, fHi)) { step *= 0.5; if (step < boundaryPressureTolerance_) { + if (std::abs(fLo) < boundaryResidualTolerance_ + && certifyBoundary(lo, yLo, fLo)) { + return true; + } + bool recovered = false; + // A fixed-pressure trial can fail over a finite interval + // even though the stationary branch converges on both + // sides. Use only a converged residual to update or close + // the bracket. + // Twenty 0.1-spaced probes bound the skipped interval to + // two units in ln(p), or a pressure ratio of exp(2). + constexpr int maxGapProbes = 20; + for (int probe = 1; probe <= maxGapProbes; ++probe) { + hi = lo + dir * Scalar{0.1} * probe; + yHi = yLo; + if (!evaluate(hi, yHi, fHi)) { + continue; + } + recovered = true; + if (fHi <= 0.0) { + bracketed = true; + break; + } + lo = hi; + fLo = fHi; + yLo = yHi; + step = 0.01; + break; + } + if (bracketed) { + break; + } + if (recovered) { + continue; + } return false; } continue; @@ -813,83 +923,48 @@ class SaturationPressure return false; } + // A gap-recovered bracket can span pressures where no continuation + // composition was found. The interpolated Y below is only an initial + // guess; each result must converge and pass final certification. for (int iter = 0; iter < maxBoundarySteps; ++iter) { // Safeguard the secant so that even a flat residual contracts // the bracket; interpolate the trial amounts in log space. - const Scalar fraction = std::clamp(fLo / (fLo - fHi), Scalar{0.1}, Scalar{0.9}); - const Scalar mid = lo + fraction * (hi - lo); + Scalar fraction = std::clamp(fLo / (fLo - fHi), Scalar{0.1}, Scalar{0.9}); + Scalar mid{}; CompVec Y{}; - for (int c = 0; c < numComponents; ++c) { - if (z[c] > 0.0) { - Y[c] = std::exp((1.0 - fraction) * std::log(yLo[c]) - + fraction * std::log(yHi[c])); - } - } Scalar f{}; - if (!evaluate(mid, Y, f)) { - return false; - } - if (std::abs(hi - lo) < boundaryPressureTolerance_ - && std::abs(f) < boundaryResidualTolerance_) { - const Scalar p = std::exp(mid); - CompVec K = Kp; - Scalar direction = 0.0; - for (int c = 0; c < numComponents; ++c) { - K[c] /= p; - Y[c] /= std::exp(f); - direction += (Y[c] - z[c]) * std::log(K[c]); - } - // Identify the continued trial by its enrichment in the - // Wilson volatility direction, then verify fugacity equality - // with the corresponding liquid/vapour EOS roots. - if (std::abs(direction) < directionTolerance_) { - return false; - } - const bool bubble = direction > 0.0; - const unsigned knownPhase = bubble ? oilPhaseIdx : gasPhaseIdx; - const unsigned trialPhase = bubble ? gasPhaseIdx : oilPhaseIdx; - for (int c = 0; c < numComponents; ++c) { - fs.setMoleFraction(knownPhase, c, z[c]); - fs.setMoleFraction(trialPhase, c, Y[c]); - } - ParameterCache cache(eosType); - cache.updatePhase(fs, oilPhaseIdx); - cache.updatePhase(fs, gasPhaseIdx); - const Scalar vmL = cache.molarVolume(oilPhaseIdx); - const Scalar vmV = cache.molarVolume(gasPhaseIdx); - if (!positiveFinite_(vmL) || !positiveFinite_(vmV) - || std::min(vmL, vmV) <= 2.0e-7) { - return false; - } + const auto attempt = [&](const Scalar candidateFraction) { + fraction = candidateFraction; + mid = lo + fraction * (hi - lo); for (int c = 0; c < numComponents; ++c) { if (z[c] > 0.0) { - const Scalar a = z[c] * FluidSystem::fugacityCoefficient( - fs, cache, knownPhase, c); - const Scalar b = Y[c] * FluidSystem::fugacityCoefficient( - fs, cache, trialPhase, c); - if (!positiveFinite_(a) || !positiveFinite_(b) - || std::abs(std::log(a / b)) > stabilityTolerance_) { - return false; - } + Y[c] = std::exp((1.0 - fraction) * std::log(yLo[c]) + + fraction * std::log(yHi[c])); } } - if (knownPhaseStability_(fs, z, knownPhase, K, bubble, eosType) - != Stability::Stable) { - return false; - } - if (fromBubble) { - // Down from a bubble point the envelope can only close - // at a dew point; anything else is not this envelope. - if (bubble) { - return false; + return evaluate(mid, Y, f); + }; + bool evaluated = attempt(fraction); + // One failed fixed-pressure solve does not invalidate the two + // converged bracket ends. Try progressively finer dyadic points. + // Four levels try 15 interior points, down to 1/16 of the bracket, + // while keeping the retry work bounded. + constexpr int maxAlternateLevels = 4; + for (int level = 1; !evaluated && level <= maxAlternateLevels; ++level) { + const int denominator = 1 << level; + for (int numerator = 1; numerator < denominator; numerator += 2) { + if (attempt(Scalar(numerator) / Scalar(denominator))) { + evaluated = true; + break; } - press = p; - liquid = Y; - return true; } - press = bubble ? pSeed : p; - liquid = bubble ? seedTrial : Y; - return true; + } + if (!evaluated) { + return false; + } + if (std::abs(hi - lo) < boundaryPressureTolerance_ + && std::abs(f) < boundaryResidualTolerance_) { + return certifyBoundary(mid, Y, f); } if (f > 0.0) { lo = mid; @@ -927,7 +1002,9 @@ class SaturationPressure * \quad\hbox{(dew)}, * \f] * with \f$y_i^n=K_i^nz_i/S_b\f$ or - * \f$x_i^n=z_i/(K_i^nS_d)\f$. Once the fixed-pressure iteration converges, + * \f$x_i^n=z_i/(K_i^nS_d)\f$. A slowly contracting substitution is + * completed by the damped Newton stationary solver at the same pressure. + * Once the fixed-pressure iteration converges, * \f$F=\ln S\f$ classifies and drives the pressure search: \f$F>0\f$ is * inside the two-phase region and \f$F<0\f$ is outside. Before bracketing, * the method uses a direction-checked secant step or the fixed-point step @@ -1022,9 +1099,7 @@ class SaturationPressure bool havePrev = false; int lastSide = 0; // This scan step can skip a narrow envelope; dewPressure() tries - // continuation from another boundary if the scan fails. A smaller step - // lowers the odds without removing them -- see the class documentation - // for what a fixed step leaves unresolved. + // continuation from another boundary if the scan fails. constexpr Scalar scanStep = 0.9; constexpr int maxPressureIterations = 200; constexpr int maxSubstitutionIterations = 500; @@ -1133,6 +1208,45 @@ class SaturationPressure dPrev = d; } + // Finish slowly contracting near-critical substitutions with the + // damped Newton stationary solver used by envelope continuation. + if (!substitutionConverged) { + CompVec candidate{}; + for (int c = 0; c < numComponents; ++c) { + if (z[c] > 0.0) { + candidate[c] = bubble ? K[c] * z[c] : z[c] / K[c]; + } + } + Scalar candidateSum{}; + if (stationaryTrial_(fs, z, knownPhaseIdx, incipientPhaseIdx, + eosType, candidate, candidateSum)) { + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(incipientPhaseIdx, c, + candidate[c] / candidateSum); + } + paramCache.updatePhase(fs, incipientPhaseIdx); + const Scalar vmL = paramCache.molarVolume(oilPhaseIdx); + const Scalar vmV = paramCache.molarVolume(gasPhaseIdx); + constexpr Scalar clampedVm = 1.0e-7; + rootsDistinct = (std::min(vmL, vmV) > 2.0 * clampedVm) && + (std::abs(vmL - vmV) + > rootVolumeTolerance_ * std::max(vmL, vmV)); + trivial = true; + for (int c = 0; c < numComponents; ++c) { + const Scalar phiIncipient = FluidSystem::fugacityCoefficient( + fs, paramCache, incipientPhaseIdx, c); + const Scalar newK = bubble ? phiKnown[c] / phiIncipient + : phiIncipient / phiKnown[c]; + if (!positiveFinite_(newK)) { + return Outcome::GaveUp; + } + trivial = trivial && (std::abs(newK - 1.0) < 1.0e-5); + K[c] = newK; + } + substitutionConverged = true; + } + } + const auto markOutside = [&](const bool withValue, const Scalar value) { pOut = p; haveOut = true; diff --git a/tests/material/test_saturation_pressure.cpp b/tests/material/test_saturation_pressure.cpp index e79ba6b2ab6..b3e9750f559 100644 --- a/tests/material/test_saturation_pressure.cpp +++ b/tests/material/test_saturation_pressure.cpp @@ -744,6 +744,49 @@ BOOST_AUTO_TEST_CASE(AcceleratedDewSearchKeepsFiniteIterates) } } +BOOST_AUTO_TEST_CASE(NearCriticalDewSearchFinishesTheStationaryTrial) +{ + // Successive substitution gives up across a neighbourhood of this state. + // The stationary Newton solve finishes and certifies the same system. + const CompVec vapor{0.0, 0.91239, 0.08761}; + constexpr Scalar temp = 320.0; + Scalar press = -123.0; + CompVec liquid{}; + + BOOST_REQUIRE(SatP::dewPressure(vapor, temp, eosType, press, liquid)); + // The retrograde branch near 332 bar, not the sub-bar lower dew point + // returned when the stationary trial cannot be completed. A bound rather + // than a value: the boundary is near critical, so the pressure itself + // moves with rounding while the branch it belongs to does not. + BOOST_CHECK_GT(press / 1.0e5, 300.0); + + const auto res = equilibriumResidualFor( + vapor, FluidSystem::gasPhaseIdx, liquid, FluidSystem::oilPhaseIdx, press, temp); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_GT(res.distance, 1.0e-3); +} + +BOOST_AUTO_TEST_CASE(DewTraceCrossesAFixedPressureFailureGap) +{ + // The direct dew searches give up for this low-pressure state. Tracing + // down from its bubble boundary encounters a short interval where the + // stationary trial fails, then resumes and brackets the dew point. + const CompVec vapor{0.0, 0.6100917431192661, 0.3899082568807339}; + constexpr Scalar temp = 280.0; + Scalar press = -123.0; + CompVec liquid{}; + + BOOST_REQUIRE(SatP::dewPressure(vapor, temp, eosType, press, liquid)); + BOOST_CHECK_CLOSE(press / 1.0e5, 0.00143343, 1.0e-2); + + const auto res = equilibriumResidualFor( + vapor, FluidSystem::gasPhaseIdx, liquid, FluidSystem::oilPhaseIdx, press, temp); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_GT(res.distance, 0.1); +} + BOOST_AUTO_TEST_CASE(DewPointReachedFromTheBubbleBoundary) { // Both dew scans miss these narrow CO2/decane envelopes near criticality. From c9e90946af9f4e652374f21cd5c97ec67bf5ca76 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 11 Sep 2026 00:33:19 +0200 Subject: [PATCH 3/6] Reject saturation roots from the wrong branch Fugacity equality and S = 1 hold at both boundaries of an envelope, so a dew iteration can converge to the bubble boundary of z and the reverse. Classify the incipient phase by enrichment along the Wilson volatility direction, as envelope continuation already does, and require it to match the branch being searched. A direction lost to cancellation is acceptable only for a pure-fluid or azeotropic boundary, where distinct EOS roots certify the point instead. --- .../constraintsolvers/SaturationPressure.hpp | 24 +++++++++++++-- tests/material/test_saturation_pressure.cpp | 29 +++++++++++++++++++ 2 files changed, 50 insertions(+), 3 deletions(-) diff --git a/opm/material/constraintsolvers/SaturationPressure.hpp b/opm/material/constraintsolvers/SaturationPressure.hpp index c9d8180de17..35cb2c88d32 100644 --- a/opm/material/constraintsolvers/SaturationPressure.hpp +++ b/opm/material/constraintsolvers/SaturationPressure.hpp @@ -92,6 +92,9 @@ namespace Opm { * Therefore \f$S_Y>1\f$ proves instability; otherwise this trial does not * reject the candidate. The trial solver works in \f$u_i=\ln Y_i\f$ and * combines substitution with damped Newton steps near critical conditions. + * Both boundaries of one envelope satisfy fugacity equality and \f$S=1\f$, + * so an accepted candidate is also classified by its enrichment along the + * Wilson volatility direction and must belong to the branch searched. * * A mixture may have lower and upper dew points at a given temperature. The * upper (retrograde) one is the boundary crossed as pressure declines from @@ -1017,7 +1020,10 @@ class SaturationPressure * \f] * with bisection safeguards. Pressure changes are capped at a factor of * two. A converged \f$S=1\f$ state is returned only after rejecting the - * trivial \f$K=1\f$ solution and certifying known-phase stability. + * trivial \f$K=1\f$ solution, confirming from + * \f$D=\sum_i(w_i-z_i)\ln K_i^W\f$ that the incipient phase belongs to + * the branch being searched rather than to the opposite boundary of the + * same envelope, and certifying known-phase stability. * * A finite scan can miss a narrow phase interval, so exhausting the search * returns Outcome::GaveUp and never claims that the branch has no root. @@ -1285,17 +1291,29 @@ class SaturationPressure if (substitutionConverged && std::abs(sum - 1.0) < boundaryResidualTolerance_) { Scalar distance = 0.0; + Scalar direction = 0.0; + const CompVec candidateWilsonK = wilsonK(p); for (int c = 0; c < numComponents; ++c) { incipient[c] = (bubble ? K[c] * z[c] : z[c] / K[c]) / sum; distance += std::abs(incipient[c] - z[c]); + direction += (incipient[c] - z[c]) * std::log(candidateWilsonK[c]); } // Require a distinct phase, then check known-phase stability: // fugacity equality also admits stationary points inside the // two-phase region, including nearly identical compositions // occupying different EOS roots. - if (distance > 1.0e-3 || rootsDistinct) { + // A direct dew iteration can also converge to the bubble + // boundary of z (and vice versa). Classify the incipient + // phase by enrichment along the Wilson volatility direction, + // as envelope continuation does. A zero direction is valid + // only for a pure-fluid or azeotropic boundary with distinct + // roots. + const bool directionMatches = std::abs(direction) < directionTolerance_ + ? rootsDistinct : (bubble == (direction > 0.0)); + if ((distance > 1.0e-3 || rootsDistinct) && directionMatches) { const auto stability = - knownPhaseStability_(fs, z, knownPhaseIdx, wilsonK(p), bubble, eosType); + knownPhaseStability_(fs, z, knownPhaseIdx, + candidateWilsonK, bubble, eosType); if (stability == Stability::Indeterminate) { return Outcome::GaveUp; } diff --git a/tests/material/test_saturation_pressure.cpp b/tests/material/test_saturation_pressure.cpp index b3e9750f559..090349898f9 100644 --- a/tests/material/test_saturation_pressure.cpp +++ b/tests/material/test_saturation_pressure.cpp @@ -384,6 +384,35 @@ BOOST_AUTO_TEST_CASE(DewPressureLowerBranch) } } +BOOST_AUTO_TEST_CASE(DewSearchRejectsBubbleBoundary) +{ + // At 200 K this CO2/methane composition has a bubble boundary at about + // 44.37 bar and a dew boundary at about 12.72 bar. Both satisfy fugacity + // equality, so a dew search must also verify that its incipient phase is + // liquid-like before accepting the candidate. + constexpr Scalar temp = 200.0; + const CompVec vapor{0.20, 0.80, 0.0}; + + Scalar pBubble = 0.0; + CompVec bubbleVapor{}; + BOOST_REQUIRE(SatP::bubblePressure(vapor, temp, eosType, pBubble, bubbleVapor)); + + Scalar pDew = 0.0; + CompVec liquid{}; + BOOST_REQUIRE(SatP::dewPressure(vapor, temp, eosType, pDew, liquid)); + BOOST_CHECK_CLOSE(pBubble / 1.0e5, 44.36584124, 1.0e-4); + BOOST_CHECK_CLOSE(pDew / 1.0e5, 12.72152110, 1.0e-4); + BOOST_CHECK_LT(pDew, pBubble); + BOOST_CHECK_GT(liquid[0], vapor[0]); + + const auto res = equilibriumResidualFor( + vapor, FluidSystem::gasPhaseIdx, liquid, + FluidSystem::oilPhaseIdx, pDew, temp); + BOOST_CHECK_SMALL(res.closure, 1.0e-10); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-8); + BOOST_CHECK_GT(res.distance, 1.0); +} + BOOST_AUTO_TEST_CASE(SupercriticalLiquidHasNoBubblePoint) { // Pure methane is far above its critical temperature here, so no bubble From add1ba9f02bf13eee2caa560b997a5a25e1d9dc6 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 11 Sep 2026 00:12:23 +0200 Subject: [PATCH 4/6] Reject float saturation roots on one EOS branch The nontriviality certificate compared the known phase's molar volume at z against the trial phase's at the incipient composition, so a composition difference could pass as a second EOS root. In single precision that admits a trivial root near the critical composition. Evaluate both roots at one common composition instead. The regression pins two neighbouring near-critical states and checks that the incipient phase stays distinct from the feed. Its pressure bound is loose on purpose: that value is sensitive to code layout, while the branch is not. Single precision still does not match double precision: every convergence threshold reaches the same floor there, and a near-critical dew search can still return the opposite branch of an envelope. Clarify this precision caveat and describe the property probe as representative rather than exhaustive. --- .../constraintsolvers/SaturationPressure.hpp | 69 ++++++++++++++----- tests/material/test_saturation_pressure.cpp | 48 +++++++++++-- 2 files changed, 92 insertions(+), 25 deletions(-) diff --git a/opm/material/constraintsolvers/SaturationPressure.hpp b/opm/material/constraintsolvers/SaturationPressure.hpp index 35cb2c88d32..bdb0f5c1326 100644 --- a/opm/material/constraintsolvers/SaturationPressure.hpp +++ b/opm/material/constraintsolvers/SaturationPressure.hpp @@ -112,6 +112,11 @@ namespace Opm { * If a nontrivial stationary branch terminates at a critical endpoint before * producing a negative residual, its last point is subjected to the same EOS, * fugacity and stability certification as a sign-bracketed boundary. + * Branch resolution depends on the precision of \c Scalar. The + * double-precision validation sweep found no isolated unresolved states. In + * single precision, the convergence thresholds share the same epsilon floor, + * and a near-critical dew search can resolve the opposite boundary of the same + * envelope. * Input compositions must contain finite, non-negative mole fractions with a * positive total. They are rescaled to sum to one, so they need not do so on * entry. @@ -244,7 +249,9 @@ class SaturationPressure enum class Stability { Stable, Unstable, Indeterminate }; // Keep the established double-precision thresholds while making each - // convergence test meaningful at Scalar's precision. + // convergence test meaningful at Scalar's precision. Every threshold below + // reaches this floor in single precision, so criteria that are orders apart + // in double collapse onto one value there and limit float branch resolution. static constexpr Scalar precisionTolerance_(const Scalar minimum) { return std::max(minimum, @@ -281,6 +288,45 @@ class SaturationPressure return std::isfinite(value) && value > 0.0; } + /*! + * \brief Decide whether the cubic equation of state has two distinct + * physical roots at one composition. + * + * Pure fluids and azeotropes have equal phase compositions at saturation, + * so distinct liquid and vapour roots are their only nontriviality + * certificate. Both roots are taken at the same composition: comparing + * the roots at the known and incipient compositions instead lets a small + * composition difference pass as a second root at float precision. + * + * \param[in] fs Fluid state supplying the temperature and pressure. + * \param[in] composition Composition at which both roots are evaluated. + * \param[in] eosType Cubic EOS type. + * \return whether both roots are physical and their molar volumes differ + * by more than \c rootVolumeTolerance_. + */ + static bool rootsDistinctAtComposition_( + const CompositionalFluidState& fs, + const CompVec& composition, + const EOSType eosType) + { + auto rootState = fs; + for (int c = 0; c < numComponents; ++c) { + rootState.setMoleFraction(oilPhaseIdx, c, composition[c]); + rootState.setMoleFraction(gasPhaseIdx, c, composition[c]); + } + + ParameterCache rootCache(eosType); + rootCache.updatePhase(rootState, oilPhaseIdx); + rootCache.updatePhase(rootState, gasPhaseIdx); + const Scalar vmL = rootCache.molarVolume(oilPhaseIdx); + const Scalar vmV = rootCache.molarVolume(gasPhaseIdx); + constexpr Scalar clampedVm = 1.0e-7; + return positiveFinite_(vmL) && positiveFinite_(vmV) + && std::min(vmL, vmV) > 2.0 * clampedVm + && std::abs(vmL - vmV) + > rootVolumeTolerance_ * std::max(vmL, vmV); + } + /*! * \brief Apply bounded dominant-eigenvalue extrapolation to a fixed-point * iterate. @@ -1138,7 +1184,6 @@ class SaturationPressure // Fugacity equality at fixed pressure: K_c = phi_liquid / phi_vapour. bool trivial = false; - bool rootsDistinct = false; bool substitutionConverged = false; // The previous change of ln K, for the dominant-eigenvalue // extrapolation below. @@ -1171,17 +1216,6 @@ class SaturationPressure ? ParameterCache::None : ParameterCache::Temperature | ParameterCache::Pressure; paramCache.updatePhase(fs, incipientPhaseIdx, unchanged); - // Pure fluids and azeotropes can have K == 1 at saturation. - // Distinct molar volumes distinguish them from a trivial trial. - const Scalar vmL = paramCache.molarVolume(oilPhaseIdx); - const Scalar vmV = paramCache.molarVolume(gasPhaseIdx); - // A molar volume near the cubic EOS floor does not establish a - // physical second root. - constexpr Scalar clampedVm = 1.0e-7; - rootsDistinct = (std::min(vmL, vmV) > 2.0 * clampedVm) && - (std::abs(vmL - vmV) - > rootVolumeTolerance_ * std::max(vmL, vmV)); - Scalar change = 0.0; trivial = true; CompVec d; @@ -1231,12 +1265,6 @@ class SaturationPressure candidate[c] / candidateSum); } paramCache.updatePhase(fs, incipientPhaseIdx); - const Scalar vmL = paramCache.molarVolume(oilPhaseIdx); - const Scalar vmV = paramCache.molarVolume(gasPhaseIdx); - constexpr Scalar clampedVm = 1.0e-7; - rootsDistinct = (std::min(vmL, vmV) > 2.0 * clampedVm) && - (std::abs(vmL - vmV) - > rootVolumeTolerance_ * std::max(vmL, vmV)); trivial = true; for (int c = 0; c < numComponents; ++c) { const Scalar phiIncipient = FluidSystem::fugacityCoefficient( @@ -1253,6 +1281,9 @@ class SaturationPressure } } + const bool rootsDistinct = substitutionConverged + && rootsDistinctAtComposition_(fs, z, eosType); + const auto markOutside = [&](const bool withValue, const Scalar value) { pOut = p; haveOut = true; diff --git a/tests/material/test_saturation_pressure.cpp b/tests/material/test_saturation_pressure.cpp index 090349898f9..3d7406188c1 100644 --- a/tests/material/test_saturation_pressure.cpp +++ b/tests/material/test_saturation_pressure.cpp @@ -54,9 +54,9 @@ #include #include #include -#include #include #include +#include namespace { @@ -359,6 +359,42 @@ BOOST_AUTO_TEST_CASE(SinglePrecisionBubbleDewRoundTrip) BOOST_CHECK_SMALL(res.fugacity, 1.0e-3); } +BOOST_AUTO_TEST_CASE(SinglePrecisionRejectsTrivialBubbleRoot) +{ + // Near the critical composition, float roundoff makes the known and trial + // compositions slightly different even after an iteration collapses onto + // one EOS root. That composition-induced volume difference must not be + // accepted as distinct liquid and vapour roots. + // The neighbouring probe makes the regression insensitive to small + // compiler-dependent changes in the pressure scan path. + for (const float zMethane : {0.86150f, 0.86200f}) { + const FloatSatP::CompVec liquid{0.0f, zMethane, 1.0f - zMethane}; + float pBubble = -1.0f; + FloatSatP::CompVec vapor{}; + BOOST_REQUIRE(FloatSatP::bubblePressure(liquid, 373.15f, eosType, + pBubble, vapor)); + + Scalar pBubbleDouble = -1.0; + CompVec vaporDouble{}; + BOOST_REQUIRE(SatP::bubblePressure( + {0.0, static_cast(zMethane), 1.0 - static_cast(zMethane)}, + temperature, eosType, pBubbleDouble, vaporDouble)); + // The points are near critical, so float rounding moves the pressure + // more than in the round-trip case above while preserving the branch. + // Bounded loosely on purpose: the distance check below is what rejects + // a trivial root, and this pressure is sensitive enough to code layout + // that a tight bound would fail on an unrelated recompilation. + BOOST_CHECK_CLOSE(static_cast(pBubble), pBubbleDouble, 5.0e-1); + + const auto res = equilibriumResidualFor( + liquid, FloatFluidSystem::oilPhaseIdx, + vapor, FloatFluidSystem::gasPhaseIdx, pBubble, 373.15f); + BOOST_CHECK_SMALL(res.closure, 2.0e-5); + BOOST_CHECK_SMALL(res.fugacity, 1.0e-3); + BOOST_CHECK_GT(res.distance, 1.0e-2); + } +} + BOOST_AUTO_TEST_CASE(DewPressureLowerBranch) { // These mixtures have an ordinary envelope ending at a bubble boundary. @@ -857,11 +893,11 @@ BOOST_AUTO_TEST_CASE(DewPointReachedFromTheBubbleBoundary) BOOST_AUTO_TEST_CASE(PropertyProbeOverEosVariantsAndStates) { - // A deterministic sweep of the whole input domain. It asserts the three - // properties that hold for every input rather than any particular - // pressure: the call does not throw, a successful result is a usable - // saturation pressure with a normalised incipient composition, and a - // failed one leaves both outputs exactly as the caller left them. + // A deterministic sweep of representative valid inputs. It asserts three + // general properties rather than any particular pressure: the call does + // not throw, a successful result is a usable saturation pressure with a + // normalised incipient composition, and a failed one leaves both outputs + // exactly as the caller left them. constexpr std::array eosVariants{ Opm::CompositionalConfig::EOSType::PR, Opm::CompositionalConfig::EOSType::PRCORR, From 9d5d73e75dcc30f10ce878de2ff0c478ea4c790c Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 11 Sep 2026 10:27:47 +0200 Subject: [PATCH 5/6] Address review on constants and lambda captures Name the tuning constants, each with a line on what it governs and whether it is a free choice. Give every lambda an explicit capture list, with read-only captures bound through std::as_const. Sum the stability trial with std::accumulate, seeded Scalar{0} so it stays in Scalar. Double results are bit-identical across the 8,400-state sweep. --- .../constraintsolvers/SaturationPressure.hpp | 104 +++++++++++++----- 1 file changed, 77 insertions(+), 27 deletions(-) diff --git a/opm/material/constraintsolvers/SaturationPressure.hpp b/opm/material/constraintsolvers/SaturationPressure.hpp index bdb0f5c1326..691caa7082a 100644 --- a/opm/material/constraintsolvers/SaturationPressure.hpp +++ b/opm/material/constraintsolvers/SaturationPressure.hpp @@ -41,6 +41,7 @@ #include #include #include +#include namespace Opm { @@ -265,6 +266,25 @@ class SaturationPressure static constexpr Scalar rootVolumeTolerance_ = precisionTolerance_(Scalar{1.0e-9}); static constexpr Scalar directionTolerance_ = precisionTolerance_(Scalar{1.0e-6}); + // Minimum L1 distance used to distinguish the incipient composition from + // the feed; distinct EOS roots provide an independent check. + static constexpr Scalar distinctCompositionDistance_ = 1.0e-3; + // Maximum |K_i - 1| for a trivial trial. This is looser than + // nonlinearTolerance_ so triviality can be detected before convergence. + static constexpr Scalar trivialKDistance_ = 1.0e-5; + // Molar-volume floor used by the cubic EOS. Candidate roots must exceed + // twice this value to avoid accepting a root pinned to the clamp. + static constexpr Scalar clampedMolarVolume_ = 1.0e-7; + // Number of substitution iterations between acceleration attempts. A sweep + // resolved the most bubble and dew states with intervals from three to five. + static constexpr int accelerationInterval_ = 4; + // Upper bound for the acceleration eigenvalue ratio. Keeping it below one + // bounds r / (1 - r), the extrapolation factor. + static constexpr Scalar maxAccelerationRatio_ = 0.98; + // Step in ln Y used to finite-difference the Jacobian; this perturbs Y + // relatively rather than by a fixed absolute amount. + static constexpr Scalar jacobianStep_ = 1.0e-5; + // Rescale a composition to sum to one. The boundary residual is measured // against one, so a composition that sums to 1 + d cannot meet that test // once |d| exceeds the tolerance. Callers assemble z in floating point, @@ -320,9 +340,8 @@ class SaturationPressure rootCache.updatePhase(rootState, gasPhaseIdx); const Scalar vmL = rootCache.molarVolume(oilPhaseIdx); const Scalar vmV = rootCache.molarVolume(gasPhaseIdx); - constexpr Scalar clampedVm = 1.0e-7; return positiveFinite_(vmL) && positiveFinite_(vmV) - && std::min(vmL, vmV) > 2.0 * clampedVm + && std::min(vmL, vmV) > 2.0 * clampedMolarVolume_ && std::abs(vmL - vmV) > rootVolumeTolerance_ * std::max(vmL, vmV); } @@ -360,7 +379,7 @@ class SaturationPressure if (!positiveFinite_(den) || !positiveFinite_(num)) { return; } - const Scalar ratio = std::min(num / den, Scalar{0.98}); + const Scalar ratio = std::min(num / den, maxAccelerationRatio_); const Scalar remaining = ratio / (1.0 - ratio); const Scalar maxStep = std::log(Scalar{2}); CompVec next = values; @@ -466,10 +485,9 @@ class SaturationPressure int settled = 0; constexpr int maxStabilityIterations = 500; for (int iter = 0; iter < maxStabilityIterations; ++iter) { - Scalar sumY = 0.0; - for (const Scalar amount : Y) { - sumY += amount; - } + // Seed with Scalar{0} so std::accumulate stays in Scalar rather + // than widening float instantiations to double. + const Scalar sumY = std::accumulate(Y.begin(), Y.end(), Scalar{0}); if (!positiveFinite_(sumY)) { return Stability::Indeterminate; } @@ -524,7 +542,7 @@ class SaturationPressure settled = 0; } sumPrev = sumNew; - if (iter % 4 == 3) { + if (iter % accelerationInterval_ == accelerationInterval_ - 1) { accelerate_(Y, d, dPrev); } dPrev = d; @@ -590,7 +608,8 @@ class SaturationPressure } auto trial = fs; bool trialCacheInitialized = false; - const auto residual = [&](const Vector& v, Vector& r) { + const auto residual = [&z, &trial, trialPhase, &trialCacheInitialized, &cache, + &target = std::as_const(target)](const Vector& v, Vector& r) { CompVec amounts{}; Scalar total = 0.0; for (int c = 0; c < numComponents; ++c) { @@ -645,7 +664,7 @@ class SaturationPressure if (iter >= substitutionPasses) { Matrix jac(0.0); // A fixed step in ln Y gives a relative perturbation of Y. - constexpr Scalar h = 1.0e-5; + constexpr Scalar h = jacobianStep_; for (int c = 0; c < numComponents; ++c) { Vector v = u; v[c] += h; @@ -762,7 +781,8 @@ class SaturationPressure fs.setMoleFraction(oilPhaseIdx, c, z[c]); } const auto Kp = wilsonKp_(temp); - const auto evaluate = [&](const Scalar lnp, CompVec& Y, Scalar& f) { + const auto evaluate = [&fs, &z, eosType, &Kp](const Scalar lnp, + CompVec& Y, Scalar& f) { const Scalar p = std::exp(lnp); if (!positiveFinite_(p)) { return false; @@ -808,7 +828,7 @@ class SaturationPressure distance += std::abs(candidate[c] / sum - z[c]); } // The trivial zero is not a sign for root bracketing. - if (distance > 1.0e-3 && std::log(sum) > best) { + if (distance > distinctCompositionDistance_ && std::log(sum) > best) { best = std::log(sum); bestY = candidate; } @@ -823,7 +843,9 @@ class SaturationPressure return true; }; - const auto certifyBoundary = [&](const Scalar lnp, CompVec Y, const Scalar f) { + const auto certifyBoundary = [&fs, &z, eosType, &Kp, fromBubble, pSeed, + &seedTrial, &press, &liquid]( + const Scalar lnp, CompVec Y, const Scalar f) { const Scalar p = std::exp(lnp); CompVec K = Kp; Scalar direction = 0.0; @@ -853,7 +875,7 @@ class SaturationPressure const Scalar vmL = cache.molarVolume(oilPhaseIdx); const Scalar vmV = cache.molarVolume(gasPhaseIdx); if (!positiveFinite_(vmL) || !positiveFinite_(vmV) - || std::min(vmL, vmV) <= 2.0e-7) { + || std::min(vmL, vmV) <= 2.0 * clampedMolarVolume_) { return false; } for (int c = 0; c < numComponents; ++c) { @@ -887,10 +909,24 @@ class SaturationPressure return true; }; + // All pressure steps below are in ln(p). + // A small initial offset helps the first trial remain inside the envelope. + constexpr Scalar initialWalkStep = 0.01; + // Limiting the walking stride reduces the chance of skipping a narrow + // boundary as the step doubles. + constexpr Scalar maxWalkStride = 0.1; + // Probes at this spacing can recover a bracket beyond an interval where + // the fixed-pressure solve does not converge. + constexpr Scalar gapProbeSpacing = 0.1; + // These bounds keep a secant trial away from both endpoints so every + // successful refinement contracts the bracket. + constexpr Scalar minSecantFraction = 0.1; + constexpr Scalar maxSecantFraction = 0.9; + const Scalar seed = std::log(pSeed); Scalar lo{}, fLo{}; CompVec yLo{}; - Scalar step = 0.01; + Scalar step = initialWalkStep; bool inside = false; // Halve the seed offset until a trial lands inside the envelope. constexpr int maxSeedRefinements = 20; @@ -933,7 +969,7 @@ class SaturationPressure // two units in ln(p), or a pressure ratio of exp(2). constexpr int maxGapProbes = 20; for (int probe = 1; probe <= maxGapProbes; ++probe) { - hi = lo + dir * Scalar{0.1} * probe; + hi = lo + dir * gapProbeSpacing * probe; yHi = yLo; if (!evaluate(hi, yHi, fHi)) { continue; @@ -946,7 +982,7 @@ class SaturationPressure lo = hi; fLo = fHi; yLo = yHi; - step = 0.01; + step = initialWalkStep; break; } if (bracketed) { @@ -966,7 +1002,7 @@ class SaturationPressure lo = hi; fLo = fHi; yLo = yHi; - step = std::min(Scalar{0.1}, Scalar{2} * step); + step = std::min(maxWalkStride, Scalar{2} * step); } if (!bracketed) { return false; @@ -978,11 +1014,15 @@ class SaturationPressure for (int iter = 0; iter < maxBoundarySteps; ++iter) { // Safeguard the secant so that even a flat residual contracts // the bracket; interpolate the trial amounts in log space. - Scalar fraction = std::clamp(fLo / (fLo - fHi), Scalar{0.1}, Scalar{0.9}); + Scalar fraction = std::clamp(fLo / (fLo - fHi), + minSecantFraction, maxSecantFraction); Scalar mid{}; CompVec Y{}; Scalar f{}; - const auto attempt = [&](const Scalar candidateFraction) { + const auto attempt = [&fraction, &mid, &Y, &f, &z, &evaluate, + &lo = std::as_const(lo), &hi = std::as_const(hi), + &yLo = std::as_const(yLo), + &yHi = std::as_const(yHi)](const Scalar candidateFraction) { fraction = candidateFraction; mid = lo + fraction * (hi - lo); for (int c = 0; c < numComponents; ++c) { @@ -1156,7 +1196,12 @@ class SaturationPressure constexpr int maxPressureIterations = 200; constexpr int maxSubstitutionIterations = 500; - const auto safeguardPressure = [&](Scalar pNext) { + // Move an unsafe pressure proposal inside the current bracket without + // changing either endpoint. + const auto safeguardPressure = [fromAbove, &haveIn = std::as_const(haveIn), + &haveOut = std::as_const(haveOut), + &pIn = std::as_const(pIn), &pOut = std::as_const(pOut), + &p = std::as_const(p)](Scalar pNext) { if (haveIn && haveOut) { const Scalar lo = std::min(pIn, pOut); const Scalar hi = std::max(pIn, pOut); @@ -1231,7 +1276,7 @@ class SaturationPressure // unreachable for the large K of a light component. change = std::max(change, std::abs(newK - K[c]) / std::max(Scalar{1}, std::abs(newK))); - trivial = trivial && (std::abs(newK - 1.0) < 1.0e-5); + trivial = trivial && (std::abs(newK - 1.0) < trivialKDistance_); d[c] = std::log(newK) - std::log(K[c]); K[c] = newK; } @@ -1242,7 +1287,7 @@ class SaturationPressure // Every fourth pass, extrapolate the remaining geometric // series in ln K using the last two substitution changes. - if (inner % 4 == 3) { + if (inner % accelerationInterval_ == accelerationInterval_ - 1) { accelerate_(K, d, dPrev); } dPrev = d; @@ -1274,7 +1319,7 @@ class SaturationPressure if (!positiveFinite_(newK)) { return Outcome::GaveUp; } - trivial = trivial && (std::abs(newK - 1.0) < 1.0e-5); + trivial = trivial && (std::abs(newK - 1.0) < trivialKDistance_); K[c] = newK; } substitutionConverged = true; @@ -1284,7 +1329,8 @@ class SaturationPressure const bool rootsDistinct = substitutionConverged && rootsDistinctAtComposition_(fs, z, eosType); - const auto markOutside = [&](const bool withValue, const Scalar value) { + const auto markOutside = [&pOut, &haveOut, &haveFOut, &fOut, &p = std::as_const(p)]( + const bool withValue, const Scalar value) { pOut = p; haveOut = true; haveFOut = withValue; @@ -1292,7 +1338,10 @@ class SaturationPressure }; // Step from a trial the substitution could not classify: bisect // the bracket if there is one, otherwise continue the scan. - const auto stepUnclassified = [&]() { + const auto stepUnclassified = [&p, &K, fromAbove, &wilsonK, + &haveIn = std::as_const(haveIn), + &pOut = std::as_const(pOut), + &pIn = std::as_const(pIn)]() { p = haveIn ? std::sqrt(pOut * pIn) : p * (fromAbove ? scanStep : Scalar{1} / scanStep); K = wilsonK(p); @@ -1341,7 +1390,8 @@ class SaturationPressure // roots. const bool directionMatches = std::abs(direction) < directionTolerance_ ? rootsDistinct : (bubble == (direction > 0.0)); - if ((distance > 1.0e-3 || rootsDistinct) && directionMatches) { + if ((distance > distinctCompositionDistance_ || rootsDistinct) + && directionMatches) { const auto stability = knownPhaseStability_(fs, z, knownPhaseIdx, candidateWilsonK, bubble, eosType); From a2502a8907a3ee714d570c8d12d41ed11fddce97 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 11 Sep 2026 13:18:22 +0200 Subject: [PATCH 6/6] Address saturation pressure review comments --- .../constraintsolvers/SaturationPressure.hpp | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/opm/material/constraintsolvers/SaturationPressure.hpp b/opm/material/constraintsolvers/SaturationPressure.hpp index 691caa7082a..b9b7e1986b2 100644 --- a/opm/material/constraintsolvers/SaturationPressure.hpp +++ b/opm/material/constraintsolvers/SaturationPressure.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include @@ -305,7 +306,13 @@ class SaturationPressure static bool positiveFinite_(const Scalar value) { - return std::isfinite(value) && value > 0.0; + return std::isnormal(value) && !std::signbit(value); + } + + // Limit a logarithmic step to a factor-of-two change. + static constexpr Scalar limitLogStep_(const Scalar step) + { + return std::clamp(step, -std::numbers::ln2_v, std::numbers::ln2_v); } /*! @@ -381,7 +388,7 @@ class SaturationPressure } const Scalar ratio = std::min(num / den, maxAccelerationRatio_); const Scalar remaining = ratio / (1.0 - ratio); - const Scalar maxStep = std::log(Scalar{2}); + constexpr Scalar maxStep = std::numbers::ln2_v; CompVec next = values; for (int c = 0; c < numComponents; ++c) { if (values[c] == 0.0) { // Absent component in a stability trial. @@ -703,7 +710,7 @@ class SaturationPressure } if (!accepted) { for (int c = 0; c < numComponents; ++c) { - u[c] -= std::clamp(r[c], -std::log(Scalar{2}), std::log(Scalar{2})); + u[c] -= limitLogStep_(r[c]); } } } @@ -1042,7 +1049,7 @@ class SaturationPressure for (int level = 1; !evaluated && level <= maxAlternateLevels; ++level) { const int denominator = 1 << level; for (int numerator = 1; numerator < denominator; numerator += 2) { - if (attempt(Scalar(numerator) / Scalar(denominator))) { + if (attempt(std::ldexp(static_cast(numerator), -level))) { evaluated = true; break; } @@ -1421,9 +1428,7 @@ class SaturationPressure // fixed-point step K ~ 1/p implies, p * sum towards the root, // and let the carried-over K keep converging within the bracket. const Scalar lnpFixed = lnp + (fromAbove ? f : -f); - p = safeguardPressure(std::exp(lnp + std::clamp(lnpFixed - lnp, - std::log(Scalar{0.5}), - std::log(Scalar{2})))); + p = safeguardPressure(std::exp(lnp + limitLogStep_(lnpFixed - lnp))); continue; } @@ -1482,7 +1487,7 @@ class SaturationPressure // Limit extrapolation to a factor of two and stay inside a known // bracket. This step limit alone cannot prevent skipping an // unbracketed narrow envelope. - lnpNext = lnp + std::clamp(lnpNext - lnp, std::log(Scalar{0.5}), std::log(Scalar{2})); + lnpNext = lnp + limitLogStep_(lnpNext - lnp); p = safeguardPressure(std::exp(lnpNext)); }