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..b9b7e1986b2
--- /dev/null
+++ b/opm/material/constraintsolvers/SaturationPressure.hpp
@@ -0,0 +1,1502 @@
+// -*- 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
+#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.
+ * 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
+ * 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.
+ * 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.
+ * 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.
+ */
+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. 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,
+ 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});
+
+ // 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,
+ // 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::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);
+ }
+
+ /*!
+ * \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);
+ return positiveFinite_(vmL) && positiveFinite_(vmV)
+ && std::min(vmL, vmV) > 2.0 * clampedMolarVolume_
+ && std::abs(vmL - vmV)
+ > rootVolumeTolerance_ * std::max(vmL, vmV);
+ }
+
+ /*!
+ * \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, maxAccelerationRatio_);
+ const Scalar remaining = ratio / (1.0 - ratio);
+ 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.
+ 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) {
+ // 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;
+ }
+ 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 % accelerationInterval_ == accelerationInterval_ - 1) {
+ 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 = [&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) {
+ 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 = jacobianStep_;
+ 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] -= limitLogStep_(r[c]);
+ }
+ }
+ }
+ 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. 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,
+ * \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, 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
+ * 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 = [&fs, &z, eosType, &Kp](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 > distinctCompositionDistance_ && 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 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;
+ 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.0 * clampedMolarVolume_) {
+ 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;
+ };
+
+ // 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 = initialWalkStep;
+ bool inside = false;
+ // Halve the seed offset until a trial lands inside the envelope.
+ constexpr int maxSeedRefinements = 20;
+ // 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;
+ 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_) {
+ 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 * gapProbeSpacing * 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 = initialWalkStep;
+ break;
+ }
+ if (bracketed) {
+ break;
+ }
+ if (recovered) {
+ continue;
+ }
+ return false;
+ }
+ continue;
+ }
+ if (fHi <= 0.0) {
+ bracketed = true;
+ break;
+ }
+ lo = hi;
+ fLo = fHi;
+ yLo = yHi;
+ step = std::min(maxWalkStride, Scalar{2} * step);
+ }
+ if (!bracketed) {
+ 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.
+ Scalar fraction = std::clamp(fLo / (fLo - fHi),
+ minSecantFraction, maxSecantFraction);
+ Scalar mid{};
+ CompVec Y{};
+ Scalar f{};
+ 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) {
+ if (z[c] > 0.0) {
+ Y[c] = std::exp((1.0 - fraction) * std::log(yLo[c])
+ + fraction * std::log(yHi[c]));
+ }
+ }
+ 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(std::ldexp(static_cast(numerator), -level))) {
+ evaluated = true;
+ break;
+ }
+ }
+ }
+ 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;
+ 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$. 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
+ * \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, 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.
+ *
+ * \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.
+ constexpr Scalar scanStep = 0.9;
+ constexpr int maxPressureIterations = 200;
+ constexpr int maxSubstitutionIterations = 500;
+
+ // 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);
+ 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 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);
+
+ 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) < trivialKDistance_);
+ 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 % accelerationInterval_ == accelerationInterval_ - 1) {
+ accelerate_(K, d, dPrev);
+ }
+ 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);
+ 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) < trivialKDistance_);
+ K[c] = newK;
+ }
+ substitutionConverged = true;
+ }
+ }
+
+ const bool rootsDistinct = substitutionConverged
+ && rootsDistinctAtComposition_(fs, z, eosType);
+
+ const auto markOutside = [&pOut, &haveOut, &haveFOut, &fOut, &p = std::as_const(p)](
+ 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, &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);
+ };
+
+ // 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;
+ 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.
+ // 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 > distinctCompositionDistance_ || rootsDistinct)
+ && directionMatches) {
+ const auto stability =
+ knownPhaseStability_(fs, z, knownPhaseIdx,
+ candidateWilsonK, 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 + limitLogStep_(lnpFixed - lnp)));
+ 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 + limitLogStep_(lnpNext - lnp);
+ 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..3d7406188c1
--- /dev/null
+++ b/tests/material/test_saturation_pressure.cpp
@@ -0,0 +1,980 @@
+// -*- 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(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.
+ // 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(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
+ // 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(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.
+ // 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 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,
+ 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);
+}