diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index b88aea638ea..db2f5a74067 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -104,6 +104,7 @@ list(APPEND MAIN_SOURCE_FILES
opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.cpp
opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.cpp
opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.cpp
+ opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.cpp
opm/input/eclipse/EclipseState/Geochemistry/SpeciesConfig.cpp
opm/input/eclipse/EclipseState/Geochemistry/MineralConfig.cpp
opm/input/eclipse/EclipseState/Geochemistry/IonExchangeConfig.cpp
@@ -714,8 +715,10 @@ 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
)
if(dune-common_FOUND)
@@ -970,6 +973,7 @@ list(APPEND PUBLIC_HEADER_FILES
opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp
opm/input/eclipse/EclipseState/Co2StoreConfig.hpp
opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp
+ opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.hpp
opm/input/eclipse/EclipseState/EclipseConfig.hpp
opm/input/eclipse/EclipseState/EclipseState.hpp
opm/input/eclipse/EclipseState/EndpointScaling.hpp
@@ -1346,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/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.cpp b/opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.cpp
new file mode 100644
index 00000000000..abe76916674
--- /dev/null
+++ b/opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.cpp
@@ -0,0 +1,78 @@
+/*
+ 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 3 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 .
+*/
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+namespace Opm {
+
+double moleFractionTolerance()
+{
+ return 1.0e-4;
+}
+
+double exactSumSlack(const std::size_t numValues)
+{
+ return 2.0 * numValues * std::numeric_limits::epsilon();
+}
+
+void normalizeMoleFractions(std::vector& fractions,
+ const std::string& what,
+ const KeywordLocation& location)
+{
+ const double sum = std::accumulate(fractions.begin(), fractions.end(), 0.0);
+
+ // A non-finite fraction makes the sum non-finite, and every comparison
+ // against a NaN is false: the checks below would all pass and the
+ // fractions would then be "normalized" by dividing through the NaN.
+ if (!std::isfinite(sum)) {
+ throw OpmInputError(fmt::format("The mole fractions of {} sum to {}, "
+ "which is not a finite number.", what, sum),
+ location);
+ }
+
+ const double deviation = std::abs(sum - 1.0);
+
+ if (deviation > moleFractionTolerance()) {
+ throw OpmInputError(fmt::format("The mole fractions of {} sum to {}, "
+ "which is not one.", what, sum),
+ location);
+ }
+
+ if (deviation > exactSumSlack(fractions.size())) {
+ // Printed round-trip: a deviation small enough to round away at a
+ // fixed precision is exactly the one worth naming.
+ OpmLog::warning(fmt::format("The mole fractions of {} sum to {}: they should "
+ "sum to unity and have been normalized.", what, sum));
+ }
+
+ for (auto& x : fractions) {
+ x /= sum;
+ }
+}
+
+} // namespace Opm
diff --git a/opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.hpp b/opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.hpp
new file mode 100644
index 00000000000..59183372a05
--- /dev/null
+++ b/opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.hpp
@@ -0,0 +1,53 @@
+/*
+ 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 3 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 .
+*/
+
+#ifndef OPM_NORMALIZE_MOLE_FRACTIONS_HPP
+#define OPM_NORMALIZE_MOLE_FRACTIONS_HPP
+
+#include
+#include
+#include
+
+namespace Opm {
+
+class KeywordLocation;
+
+/// How far a set of mole fractions may sum away from one before it is
+/// rejected. Writing a composition with a few digits costs this much.
+double moleFractionTolerance();
+
+/// Slack below which a sum counts as exactly one: summing n values of order
+/// one costs about n roundings, and representing them costs as many again.
+double exactSumSlack(std::size_t numValues);
+
+/// Scales \p fractions so that they sum to one, and says so when the scaling
+/// was more than the arithmetic of the sum.
+///
+/// \param what Names the input in the warning, e.g. "row 2 of COMPVD table 1"
+/// or "stream 'ISTR'".
+///
+/// \throw OpmInputError when the sum is too far from one to be the rounding of
+/// the values.
+void normalizeMoleFractions(std::vector& fractions,
+ const std::string& what,
+ const KeywordLocation& location);
+
+} // namespace Opm
+
+#endif // OPM_NORMALIZE_MOLE_FRACTIONS_HPP
diff --git a/opm/input/eclipse/EclipseState/InitConfig/Equil.cpp b/opm/input/eclipse/EclipseState/InitConfig/Equil.cpp
index 0ae470f1158..47401a6e12d 100644
--- a/opm/input/eclipse/EclipseState/InitConfig/Equil.cpp
+++ b/opm/input/eclipse/EclipseState/InitConfig/Equil.cpp
@@ -78,7 +78,10 @@ namespace Opm {
if (compositional) {
comp_init_type = record.getItem().get(0);
if (comp_init_type == 2 || comp_init_type == 3) {
- set_to_saturation_pressure = record.getItem().get(0) != 1;
+ // Item 11 has no default value; when left unset the pressure at
+ // the contact is set to the saturation pressure.
+ const auto& item = record.getItem();
+ set_to_saturation_pressure = !item.hasValue(0) || (item.get(0) != 1);
}
}
}
diff --git a/opm/input/eclipse/EclipseState/Tables/Tables.cpp b/opm/input/eclipse/EclipseState/Tables/Tables.cpp
index beeb1b3a3ee..70affa03846 100644
--- a/opm/input/eclipse/EclipseState/Tables/Tables.cpp
+++ b/opm/input/eclipse/EclipseState/Tables/Tables.cpp
@@ -95,6 +95,8 @@
#include
#include
+#include
+#include
#include
#include
@@ -114,6 +116,7 @@
#include
#include
#include
+#include
#include
#include
#include
@@ -2906,27 +2909,25 @@ ZmfvdTable::ZmfvdTable(const DeckItem& item, const int tableID, const int numCom
const auto nrows = item.data_size() / ncol;
const std::string tableName {"ZMFVD"};
+ std::vector moles(numComponents, 0.);
for (std::size_t row = 0; row < nrows; ++row) {
// Depth column
const std::size_t depthIdx = row * ncol;
const double siDepth = item.getSIDouble(depthIdx);
getColumn(0).addValue(siDepth, tableName);
- std::vector moles(numComponents, 0.);
// Component mole-fraction columns (dimensionless)
for (int c = 0; c < numComponents; ++c) {
const std::size_t compIdx = row * ncol + 1 + c;
- const auto mole_fraction = item.get(compIdx);
- moles[c] = mole_fraction;
- getColumn(1 + c).addValue(mole_fraction, tableName);
+ moles[c] = item.get(compIdx);
}
- // checking to make sure the sum of the mole fractions are 1.
- constexpr double epsilon = 1.e-5;
- const double sum_fractions = std::accumulate(moles.begin(), moles.end(), 0.);
- if (std::abs(sum_fractions - 1.) > epsilon) {
- const std::string reason = fmt::format("ZMFVD table {}: sum of mole fractions in row {} is not 1 (sum is {})",
- tableID + 1, row + 1, sum_fractions);
- throw OpmInputError(reason, location);
+
+ // Normalize, so the rounding never reaches the equilibration.
+ normalizeMoleFractions(moles,
+ fmt::format("row {} of ZMFVD table {}", row + 1, tableID + 1),
+ location);
+ for (int c = 0; c < numComponents; ++c) {
+ getColumn(1 + c).addValue(moles[c], tableName);
}
}
}
@@ -2986,6 +2987,7 @@ CompvdTable::CompvdTable(const DeckItem& item,
const auto& data = item.getData();
const std::string tableName{"COMPVD"};
+ std::vector moles(numComponents, 0.0);
for (std::size_t row = 0; row < nrows; ++row) {
const std::size_t rowStart = row * ncol;
@@ -2994,21 +2996,16 @@ CompvdTable::CompvdTable(const DeckItem& item,
getColumn(0).addValue(siDepth, tableName);
// Component mole-fraction columns (dimensionless).
- std::vector moles(numComponents, 0.0);
for (int c = 0; c < numComponents; ++c) {
- const auto z = data.at(rowStart + 1 + c);
- moles[c] = z;
- getColumn(1 + c).addValue(z, tableName);
+ moles[c] = data.at(rowStart + 1 + c);
}
- // Sum-to-one check, same epsilon is used in ZMFVD
- constexpr double epsilon = 1.e-5;
- const double sum_fractions = std::accumulate(moles.begin(), moles.end(), 0.);
- if (std::abs(sum_fractions - 1.) > epsilon) {
- const std::string reason = fmt::format(
- "COMPVD table {}: sum of mole fractions in row {} is not 1 (sum is {})",
- tableID + 1, row + 1, sum_fractions);
- throw OpmInputError(reason, location);
+ // Normalize, so the rounding never reaches the equilibration.
+ normalizeMoleFractions(moles,
+ fmt::format("row {} of COMPVD table {}", row + 1, tableID + 1),
+ location);
+ for (int c = 0; c < numComponents; ++c) {
+ getColumn(1 + c).addValue(moles[c], tableName);
}
// Phase flag: stored as a strong enum, validated to be exactly 0 or 1.
diff --git a/opm/input/eclipse/Schedule/Well/WellKeywordHandlers.cpp b/opm/input/eclipse/Schedule/Well/WellKeywordHandlers.cpp
index f08c95a1f27..406b36ce576 100644
--- a/opm/input/eclipse/Schedule/Well/WellKeywordHandlers.cpp
+++ b/opm/input/eclipse/Schedule/Well/WellKeywordHandlers.cpp
@@ -36,6 +36,7 @@
#include
#include
#include
+#include
#include
#include
#include
@@ -469,18 +470,19 @@ void handleWELLSTRE(HandlerContext& handlerContext)
auto& inj_streams = handlerContext.state().inj_streams;
for (const auto& record : handlerContext.keyword) {
const auto stream_name = record.getItem().getTrimmedString(0);
- const auto& composition = record.getItem().getSIDoubleData();
+ auto composition = record.getItem().getSIDoubleData();
const std::size_t num_comps = handlerContext.static_schedule().m_runspec.numComps();
if (composition.size() != num_comps) {
const std::string msg = fmt::format("The number of the composition values for stream '{}' is not the same as the number of components.", stream_name);
throw OpmInputError(msg, handlerContext.keyword.location());
}
- const double sum = std::accumulate(composition.begin(), composition.end(), 0.0);
- if (std::abs(sum - 1.0) > std::numeric_limits::epsilon()) {
- const std::string msg = fmt::format("The sum of the composition values for stream '{}' is not 1.0, but {}.", stream_name, sum);
- throw OpmInputError(msg, handlerContext.keyword.location());
- }
+ // A composition written with a few digits does not sum to one
+ // exactly; scale it rather than reject the deck over the rounding.
+ normalizeMoleFractions(composition,
+ fmt::format("stream '{}'", stream_name),
+ handlerContext.keyword.location());
+
auto composition_ptr = std::make_shared>(composition);
inj_streams.update(stream_name, std::move(composition_ptr));
}
diff --git a/opm/material/constraintsolvers/SaturationPressure.hpp b/opm/material/constraintsolvers/SaturationPressure.hpp
new file mode 100644
index 00000000000..09647ca3ef8
--- /dev/null
+++ b/opm/material/constraintsolvers/SaturationPressure.hpp
@@ -0,0 +1,368 @@
+// -*- 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
+
+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
+ * \c z is the pressure where the incipient vapour (liquid) phase appears. The
+ * equilibrium ratios K are obtained by successive substitution on the fugacity
+ * coefficient ratios at fixed pressure, and the pressure is updated to drive
+ * the total amount of the incipient phase, sum_c K_c z_c (or sum_c z_c / K_c),
+ * towards one. The pressure is approached from the single-phase side: when the
+ * substitution collapses onto the trivial solution K == 1, the pressure is
+ * moved into the two-phase region and the iteration restarted from the Wilson
+ * estimate. Pressures known to be single-phase and known to be two-phase are
+ * kept as a bracket around the saturation pressure, and the search bisects
+ * within it once both are known, so that a narrow phase envelope is not
+ * stepped over.
+ *
+ * A reservoir gas generally has two dew points at a given temperature. The
+ * upper (retrograde) one is the boundary crossed when the reservoir pressure
+ * declines, and it is the one the reference simulator reports as the
+ * saturation pressure of a gas. To match that reference behaviour the
+ * dew-point search targets the retrograde branch first and only falls back to
+ * the lower branch when no retrograde dew point is found. This remains to be
+ * revisited if the decision changes.
+ */
+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;
+
+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.
+ /// \return whether the calculation converged
+ static bool bubblePressure(const CompVec& liquid,
+ const Scalar temp,
+ const EOSType eosType,
+ Scalar& press,
+ CompVec& vapor)
+ { return solve_(liquid, temp, eosType, Mode::Bubble, press, vapor) == Outcome::Converged; }
+
+ /// 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; the lower one is only
+ /// searched when the upper branch provably has no root, so a plain
+ /// convergence failure reports false rather than the wrong branch.
+ /// On failure \p press is left untouched while \p liquid holds iteration
+ /// scratch and must not be read.
+ /// \return whether the calculation converged
+ static bool dewPressure(const CompVec& vapor,
+ const Scalar temp,
+ const EOSType eosType,
+ Scalar& press,
+ CompVec& liquid)
+ {
+ // The lower branch is only a fallback for a mixture that has no upper
+ // (retrograde) dew point. After a plain convergence failure it stays
+ // untried: reporting the lower branch then could return the wrong dew
+ // point of a mixture that does have both.
+ switch (solve_(vapor, temp, eosType, Mode::DewUpper, press, liquid)) {
+ case Outcome::Converged:
+ return true;
+ case Outcome::NoRoot:
+ return solve_(vapor, temp, eosType, Mode::DewLower, press, liquid)
+ == Outcome::Converged;
+ case Outcome::GaveUp:
+ return false;
+ }
+ return false;
+ }
+
+private:
+ // What a branch search established: a converged saturation point, positive
+ // evidence that the branch has no genuine root, or an exhausted iteration
+ // from which nothing can be concluded.
+ enum class Outcome { Converged, NoRoot, 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 };
+
+ // The Wilson correlation for K_c * press.
+ 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;
+ }
+
+ static Outcome solve_(const CompVec& z,
+ const Scalar temp,
+ const EOSType eosType,
+ const Mode mode,
+ Scalar& press,
+ CompVec& incipient)
+ {
+ 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. In every case the single-phase side is the one the scan
+ // starts on, and the pressure update moves away from it.
+ const bool fromAbove = (mode != Mode::DewLower);
+
+ // The Wilson estimate solves the saturation condition exactly since
+ // K ~ 1/p. The upper dew point is searched from the high-pressure
+ // side, so it starts from the same high estimate as the bubble point.
+ 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];
+ }
+ }
+
+ 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]);
+ }
+
+ // Bracket of the saturation pressure: "single" is a pressure known to
+ // lie on the single-phase side, "two" one known to lie inside the
+ // two-phase region. Once both are known the search bisects between
+ // them instead of stepping, so the boundary cannot be run past.
+ Scalar pSingle{}, pTwo{};
+ bool haveSingle = false;
+ bool haveTwo = false;
+ // Set whenever a scan point ends with the substitution exhausted. Such
+ // a point classifies nothing, so a scan that met one cannot conclude
+ // the branch has no root.
+ bool anyInconclusive = false;
+
+ // Step used to scan for the two-phase region before the bracket is
+ // closed. It is deliberately fine: a coarse step can cross a narrow
+ // phase envelope in one go and leave the mixture looking single-phase
+ // on both sides.
+ constexpr Scalar scanStep = 0.9;
+ constexpr int maxOuter = 200;
+ constexpr int maxInner = 100;
+
+ for (int outer = 0; outer < maxOuter; ++outer) {
+ fs.setPressure(oilPhaseIdx, p);
+ fs.setPressure(gasPhaseIdx, p);
+
+ // Fugacity equality at fixed pressure: K_c = phi_liquid / phi_vapour.
+ bool trivial = false;
+ bool rootsDistinct = false;
+ bool substitutionConverged = false;
+ for (int inner = 0; inner < maxInner; ++inner) {
+ Scalar sum = 0.0;
+ for (int c = 0; c < numComponents; ++c) {
+ incipient[c] = bubble ? K[c] * z[c] : z[c] / K[c];
+ sum += incipient[c];
+ }
+ for (int c = 0; c < numComponents; ++c) {
+ fs.setMoleFraction(incipientPhaseIdx, c, incipient[c] / sum);
+ }
+
+ typename FluidSystem::template ParameterCache paramCache(eosType);
+ paramCache.updatePhase(fs, oilPhaseIdx);
+ paramCache.updatePhase(fs, gasPhaseIdx);
+
+ // A pure component or an azeotrope has K = 1 at a genuine
+ // saturation point, where the two phases share a composition
+ // but occupy different EOS roots. The molar volumes tell that
+ // state apart from the trivial solution, whose phases are one
+ // and the same.
+ const Scalar vmL = paramCache.molarVolume(oilPhaseIdx);
+ const Scalar vmV = paramCache.molarVolume(gasPhaseIdx);
+ // The cubic EOS clamps an unphysical root to 1e-7 m^3/mol; a
+ // volume at the clamp is no real root, and treating it as a
+ // distinct phase would invent a saturation point for a
+ // supercritical mixture.
+ constexpr Scalar clampedVm = 1.0e-7;
+ rootsDistinct = (std::min(vmL, vmV) > 2.0 * clampedVm) &&
+ (std::abs(vmL - vmV) > 1.0e-9 * std::max(vmL, vmV));
+
+ Scalar change = 0.0;
+ trivial = true;
+ for (int c = 0; c < numComponents; ++c) {
+ const Scalar phiL = FluidSystem::fugacityCoefficient(
+ fs, paramCache, oilPhaseIdx, c);
+ const Scalar phiV = FluidSystem::fugacityCoefficient(
+ fs, paramCache, gasPhaseIdx, c);
+ const Scalar newK = phiL / phiV;
+ // 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);
+ K[c] = newK;
+ }
+ if (change < 1.0e-12) {
+ substitutionConverged = true;
+ break;
+ }
+ }
+
+ if (!substitutionConverged) {
+ anyInconclusive = true;
+ }
+
+ // The trivial test needs a converged K just as the pressure
+ // criterion below does: its threshold is 1e-5 while convergence
+ // is 1e-12, so an exhausted iteration can read as near-trivial
+ // while K is still moving, and would then close the bracket on
+ // a pressure it never classified.
+ if (substitutionConverged && trivial && !rootsDistinct) {
+ // The phases collapsed into one: the pressure lies on the
+ // single-phase side of the saturation pressure. Bisect
+ // towards a pressure already known to be two-phase, or scan
+ // on if the bracket is not closed yet.
+ pSingle = p;
+ haveSingle = true;
+ p = haveTwo ? std::sqrt(pSingle * pTwo)
+ : p * (fromAbove ? scanStep : Scalar{1} / scanStep);
+ K = wilsonK(p);
+ continue;
+ }
+
+ Scalar sum = 0.0;
+ for (int c = 0; c < numComponents; ++c) {
+ sum += bubble ? K[c] * z[c] : z[c] / K[c];
+ }
+ // 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) < 1.0e-10) {
+ 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]);
+ }
+ // A converged K that leaves the incipient phase with both the
+ // composition and the EOS root of the known one is the trivial
+ // solution wearing a disguise: it satisfies the saturation
+ // condition at an arbitrary pressure. Treat it as the
+ // single-phase point it is and keep searching. Mixtures where
+ // the roots have genuinely merged are refused for the same
+ // reason; returning nothing beats returning a false pressure.
+ if (distance > 1.0e-3 || rootsDistinct) {
+ press = p;
+ return Outcome::Converged;
+ }
+ pSingle = p;
+ haveSingle = true;
+ p = haveTwo ? std::sqrt(pSingle * pTwo)
+ : p * (fromAbove ? scanStep : Scalar{1} / scanStep);
+ K = wilsonK(p);
+ continue;
+ }
+
+ // Only a converged substitution certifies the pressure as lying
+ // inside the two-phase region; an exhausted one proves nothing
+ // and must not pollute the bracket.
+ if (substitutionConverged) {
+ pTwo = p;
+ haveTwo = true;
+ }
+ // Inside the two-phase region the incipient amount exceeds one and
+ // the pressure moves towards the saturation pressure: up on the
+ // bubble and retrograde branches, down on the lower dew branch.
+ const Scalar factor = (mode == Mode::DewLower) ? 1.0 / sum : sum;
+ Scalar pNext = p * std::clamp(factor, Scalar{0.5}, Scalar{2.0});
+ // Never step onto or past a pressure already known to be
+ // single-phase; bisect towards it instead.
+ if (haveSingle && ((fromAbove && pNext >= pSingle) ||
+ (!fromAbove && pNext <= pSingle)))
+ {
+ pNext = std::sqrt(p * pSingle);
+ }
+ p = pNext;
+ }
+
+ // The scan covered nine decades of pressure without meeting a
+ // two-phase state: the branch has no dew or bubble point to find.
+ // That conclusion only holds if every point along the way was
+ // classified; an exhausted substitution leaves the branch unproven.
+ if (!haveTwo) {
+ return anyInconclusive ? Outcome::GaveUp : Outcome::NoRoot;
+ }
+ // A bracket that collapsed without an accepted solution pinned the
+ // phase boundary down to a point where only the trivial solution
+ // lives; that too is positive evidence the branch has no genuine
+ // saturation point. A bracket still open is merely an unfinished
+ // search.
+ if (haveSingle &&
+ (std::abs(pSingle - pTwo) <= 1.0e-6 * std::max(pSingle, pTwo)))
+ {
+ return Outcome::NoRoot;
+ }
+ return Outcome::GaveUp;
+ }
+};
+
+} // namespace Opm
+
+#endif // OPM_SATURATION_PRESSURE_HPP
diff --git a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp
index 7d01b3f9a51..28512aa4c80 100644
--- a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp
+++ b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp
@@ -26,6 +26,7 @@
#ifndef OPM_GENERIC_OIL_GAS_WATER_FLUIDSYSTEM_HPP
#define OPM_GENERIC_OIL_GAS_WATER_FLUIDSYSTEM_HPP
+#include
#include
#include
@@ -93,15 +94,18 @@ namespace Opm {
Scalar critic_pres; // unit: parscal
Scalar critic_vol; // unit: m^3/kmol
Scalar acentric_factor; // unit: dimension less
+ Scalar volume_shift; // unit: dimension less (SSHIFT)
ComponentParam(const std::string_view name_, const Scalar molar_mass_, const Scalar critic_temp_,
- const Scalar critic_pres_, const Scalar critic_vol_, const Scalar acentric_factor_)
+ const Scalar critic_pres_, const Scalar critic_vol_, const Scalar acentric_factor_,
+ const Scalar volume_shift_ = 0.0)
: name(name_),
molar_mass(molar_mass_),
critic_temp(critic_temp_),
critic_pres(critic_pres_),
critic_vol(critic_vol_),
- acentric_factor(acentric_factor_)
+ acentric_factor(acentric_factor_),
+ volume_shift(volume_shift_)
{}
};
@@ -177,7 +181,10 @@ namespace Opm {
static_cast(eos_props.critical_temperature[c]),
static_cast(eos_props.critical_pressure[c]),
static_cast(eos_props.critical_volume[c] * 1.e3),
- static_cast(eos_props.acentric_factors[c])});
+ static_cast(eos_props.acentric_factors[c]),
+ c < eos_props.volume_shifts.size()
+ ? static_cast(eos_props.volume_shifts[c])
+ : Scalar{0}});
}
const auto& bic = eos_props.binary_interaction_coefficient;
@@ -220,6 +227,22 @@ namespace Opm {
return component_param_[compIdx].acentric_factor;
}
+
+ /*!
+ * \brief The volume shift of a component (SSHIFT) [].
+ *
+ * It corrects the molar volume of the equation of state, particularly
+ * for the liquid phase, and leaves the phase equilibrium untouched.
+ *
+ * \copydetails Doxygen::compIdxParam
+ */
+ static Scalar volumeShift(unsigned compIdx)
+ {
+ assert(isConsistent());
+ assert(compIdx < numComponents);
+
+ return component_param_[compIdx].volume_shift;
+ }
/*!
* \brief Critical temperature of a component [K].
*
@@ -325,7 +348,10 @@ namespace Opm {
assert(phaseIdx < numPhases);
if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) {
- return decay(fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx));
+ // The shift belongs here rather than in the cached volume: the
+ // fugacity coefficients are computed from the unshifted one.
+ const auto Vm = paramCache.correctedMolarVolume(fluidState, phaseIdx);
+ return decay(fluidState.averageMolarMass(phaseIdx) / Vm);
}
else {
const LhsEval& p = decay(fluidState.pressure(phaseIdx));
diff --git a/opm/material/fluidsystems/PTFlashParameterCache.hpp b/opm/material/fluidsystems/PTFlashParameterCache.hpp
index 85b2fb2a524..8bab0b1cc37 100644
--- a/opm/material/fluidsystems/PTFlashParameterCache.hpp
+++ b/opm/material/fluidsystems/PTFlashParameterCache.hpp
@@ -30,6 +30,8 @@
#ifndef OPM_PTFlash_PARAMETER_CACHE_HPP
#define OPM_PTFlash_PARAMETER_CACHE_HPP
+#include
+
#include
#include
#include
@@ -37,6 +39,8 @@
#include
+#include
+
#include
namespace Opm {
@@ -286,6 +290,60 @@ class PTFlashParameterCache
return Vm_[phaseIdx];
}
+ /*!
+ * \brief The volume shift of a phase, sum_c x_c s_c b_c [m^3/mol]
+ *
+ * It is not folded into molarVolume() because the fugacity coefficients
+ * are computed from the unshifted volume. Use correctedMolarVolume() to
+ * get the volume the fluid actually occupies.
+ *
+ * \param phaseIdx The fluid phase of interest
+ */
+ template
+ Scalar volumeShift(const FluidState& fluidState, unsigned phaseIdx) const
+ {
+ // b_c from the dimensionless B_c = b_c p / (R T).
+ const Scalar T = decay(fluidState.temperature(phaseIdx));
+ const Scalar p = decay(fluidState.pressure(phaseIdx));
+ const Scalar RT_p = Constants::R * T / p;
+
+ Scalar shift = 0;
+ for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
+ const Scalar b = decay(Bi(phaseIdx, compIdx)) * RT_p;
+ shift += decay(fluidState.moleFraction(phaseIdx, compIdx))
+ * FluidSystem::volumeShift(compIdx) * b;
+ }
+ return shift;
+ }
+
+ /*!
+ * \brief The molar volume the fluid occupies, shift included [m^3/mol]
+ *
+ * This is the physical volume of the phase: the density, the phase
+ * saturations and the transport properties all follow from it. Only the
+ * fugacity coefficients keep the unshifted molarVolume(), because the
+ * shift cancels from the equilibrium ratios and the two-parameter
+ * expression they use is derived for the unshifted root.
+ *
+ * \param phaseIdx The fluid phase of interest
+ */
+ template
+ Scalar correctedMolarVolume(const FluidState& fluidState, unsigned phaseIdx) const
+ {
+ const Scalar Vm = molarVolume(phaseIdx) - volumeShift(fluidState, phaseIdx);
+
+ // SSHIFT is unconstrained deck input. A shift larger than the molar
+ // volume leaves nothing behind, and every quantity derived from it
+ // would be meaningless rather than merely inaccurate.
+ if (!(scalarValue(Vm) > 0)) {
+ throw NumericalProblem(
+ fmt::format("The SSHIFT volume shift of phase {} leaves a corrected "
+ "molar volume of {}, which is not positive.",
+ phaseIdx, scalarValue(Vm)));
+ }
+ return Vm;
+ }
+
/*!
* \brief Returns the Peng-Robinson mixture parameters for the oil
diff --git a/tests/material/test_saturation_pressure.cpp b/tests/material/test_saturation_pressure.cpp
new file mode 100644
index 00000000000..c08eee318c5
--- /dev/null
+++ b/tests/material/test_saturation_pressure.cpp
@@ -0,0 +1,325 @@
+// -*- 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 are taken from a reference simulator run of 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 reported in the restart file for the
+ * single-phase oil cells, and the gas-oil contact values are taken from the
+ * equilibration report in the PRT file.
+ */
+#include "config.h"
+
+#define BOOST_TEST_MODULE SaturationPressure
+#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 = typename SatP::CompVec;
+
+constexpr auto eosType = Opm::CompositionalConfig::EOSType::PR;
+
+// The constant reservoir temperature (RTEMP) of the test case, 100 degC.
+constexpr Scalar temperature = 373.15;
+
+// The fluid system is initialized once with the component properties of the
+// test deck (TCRIT, PCRIT, ACF, MW, VCRIT). Only the critical properties and
+// the acentric factors enter the fugacity coefficients, but the full set is
+// provided for completeness.
+struct Fixture
+{
+ Fixture()
+ {
+ using CompParam = typename 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});
+ }
+};
+
+// How far a (known phase, incipient phase) pair at \p press is from being a
+// saturation point: the largest relative fugacity imbalance over the
+// components, how far the incipient composition is from summing to one, and
+// how far it is from the known composition (a vanishing distance is the
+// trivial solution, which satisfies the other two at any pressure).
+struct EquilibriumResidual
+{
+ Scalar fugacity{};
+ Scalar closure{};
+ Scalar distance{};
+};
+
+EquilibriumResidual equilibriumResidual(const CompVec& known,
+ const unsigned knownPhaseIdx,
+ const CompVec& incipient,
+ const unsigned incipientPhaseIdx,
+ const Scalar press)
+{
+ Opm::CompositionalFluidState fs;
+ fs.setTemperature(temperature);
+ fs.setPressure(FluidSystem::oilPhaseIdx, press);
+ fs.setPressure(FluidSystem::gasPhaseIdx, press);
+ for (int c = 0; c < numComponents; ++c) {
+ fs.setMoleFraction(knownPhaseIdx, c, known[c]);
+ fs.setMoleFraction(incipientPhaseIdx, c, incipient[c]);
+ }
+
+ typename FluidSystem::template ParameterCache paramCache(eosType);
+ paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
+ paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
+
+ EquilibriumResidual res;
+ Scalar sum = 0.0;
+ for (int c = 0; c < numComponents; ++c) {
+ const Scalar phiL =
+ FluidSystem::fugacityCoefficient(fs, paramCache, FluidSystem::oilPhaseIdx, c);
+ const Scalar phiV =
+ FluidSystem::fugacityCoefficient(fs, paramCache, FluidSystem::gasPhaseIdx, c);
+ const Scalar fL = fs.moleFraction(FluidSystem::oilPhaseIdx, c) * phiL;
+ const Scalar fV = fs.moleFraction(FluidSystem::gasPhaseIdx, c) * phiV;
+ const Scalar scale = std::max({std::abs(fL), std::abs(fV), Scalar{1.0e-12}});
+ res.fugacity = std::max(res.fugacity, 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;
+}
+
+} // 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 simulator PSAT values in bar.
+ const 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(SaturationPressureAtGasOilContact)
+{
+ // At the gas-oil contact the liquid composition from ZMFVD is
+ // (0, 0.5, 0.5) and the reference simulator reports "Pressure at gas-oil
+ // contact set to 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(DewPressureLowerBranch)
+{
+ // A methane-rich mixture that is lean enough to have no retrograde dew
+ // point at this temperature, so dewPressure() only succeeds through the
+ // lower-branch fallback. The assertion is the equilibrium condition
+ // itself rather than a previously recorded pressure, so the test states
+ // what a dew point is instead of what this solver happened to return.
+ const CompVec vapor{0.0, 0.90, 0.10};
+ Scalar press = 0.0;
+ CompVec liquid{};
+
+ BOOST_REQUIRE(SatP::dewPressure(vapor, temperature, eosType, press, liquid));
+ // A genuinely lower-branch pressure: the retrograde region of comparable
+ // mixtures sits above 100 bar, the lower dew point of this one near 1 bar.
+ BOOST_CHECK_GT(press, 0.0);
+ BOOST_CHECK_LT(press, 50.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 must be a genuinely different phase, and being the
+ // lower branch it is the heavy one: richer in decane than the vapour.
+ BOOST_CHECK_GT(res.distance, 1.0e-3);
+ BOOST_CHECK_GT(liquid[2], vapor[2]);
+}
+
+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)
+{
+ // Almost pure methane far above its critical temperature has no dew point
+ // here. The successive substitution can still settle onto K == 1 at an
+ // arbitrarily high pressure; that trivial solution satisfies the
+ // saturation condition numerically and must not be returned as an answer.
+ 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 {
+ // The expected outcome is an outright refusal, which must leave
+ // the pressure untouched rather than half-written.
+ BOOST_CHECK_LT(press, 0.0);
+ }
+ }
+}
diff --git a/tests/material/test_volume_shift.cpp b/tests/material/test_volume_shift.cpp
new file mode 100644
index 00000000000..73e00f0101c
--- /dev/null
+++ b/tests/material/test_volume_shift.cpp
@@ -0,0 +1,370 @@
+// -*- 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 that the SSHIFT volume shift reaches the density.
+ *
+ * The fluid is the seven-component one of a deck whose reference run reports
+ * a gas density of 165.87 kg/m3 at 272.599 bar and 393.15 K. Without the
+ * shift the equation of state gives 172.84.
+ */
+#include "config.h"
+
+#define BOOST_TEST_MODULE VolumeShift
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+#include
+#include
+
+namespace {
+
+using Scalar = double;
+constexpr int numComponents = 7;
+
+using FluidSystem = Opm::GenericOilGasWaterFluidSystem;
+// A distinct instantiation, so it carries its own static component parameters:
+// the same fluid with every shift set to zero. Comparing against it is what
+// makes the fugacity test below independent rather than self-referential.
+using UnshiftedSystem = Opm::GenericOilGasWaterFluidSystem;
+using CompVec = std::array;
+
+constexpr auto eosType = Opm::CompositionalConfig::EOSType::PR;
+
+// The state the reference reports the density at: the top of the column,
+// where RTEMP is 120 degC and the equilibrated pressure 272.599 bar.
+constexpr Scalar temperature = 393.15;
+constexpr Scalar pressure = 272.599e5;
+
+// A seven-component reservoir fluid: methane, CO2, ethane, propane and three
+// lumped fractions. The shifts are the SSHIFT values of the same deck.
+const CompVec z{0.88583, 0.0437, 0.034, 0.0189, 0.0152, 0.0023, 0.00007};
+const CompVec shift{-0.1595, -0.0817, -0.1134, -0.0863, -0.0243568,
+ 0.10784125400986, 0.206892761354429};
+
+struct Component
+{
+ const char* name;
+ Scalar molarMass; // kg/mol, as the parser hands it to the fluid system
+ Scalar criticalT; // K
+ Scalar criticalP; // Pa
+ Scalar criticalV; // m^3/kmol
+ Scalar acentric;
+};
+
+// MW, TCRIT, PCRIT, VCRIT and ACF as the deck gives them.
+const std::array components{{
+ {"CH4", 0.016043, 190.60, 45.40e5, 0.09900000, 0.00800},
+ {"CO2", 0.044010, 304.20, 72.80e5, 0.09400000, 0.22500},
+ {"C2H6", 0.030070, 305.40, 48.20e5, 0.14800000, 0.09800},
+ {"C3H8", 0.044097, 369.80, 41.90e5, 0.20300000, 0.15200},
+ {"C4-C9", 0.076900, 474.70, 32.40e5, 0.32321204, 0.25490},
+ {"C10-C20", 0.163000, 646.00, 20.42e5, 0.63255728, 0.55100},
+ {"C21+", 0.310850, 812.00, 14.50e5, 0.97988801, 0.90900},
+}};
+
+// The component parameters are static, so the fluid system carries one
+// configuration per process: it is initialized once, with the shifts.
+struct Fixture
+{
+ Fixture()
+ {
+ using CompParam = typename FluidSystem::ComponentParam;
+ FluidSystem::init();
+ for (int c = 0; c < numComponents; ++c) {
+ const auto& p = components[c];
+ FluidSystem::addComponent(CompParam{p.name, p.molarMass, p.criticalT,
+ p.criticalP, p.criticalV, p.acentric,
+ shift[c]});
+ }
+
+ using UnshiftedParam = typename UnshiftedSystem::ComponentParam;
+ UnshiftedSystem::init();
+ for (int c = 0; c < numComponents; ++c) {
+ const auto& p = components[c];
+ UnshiftedSystem::addComponent(UnshiftedParam{p.name, p.molarMass, p.criticalT,
+ p.criticalP, p.criticalV,
+ p.acentric, 0.0});
+ }
+ }
+};
+
+/// The gas density of the mixture, and the fugacity coefficients alongside it.
+struct PhaseState
+{
+ Scalar density{};
+ CompVec fugacityCoefficient{};
+};
+
+PhaseState gasState()
+{
+ Opm::CompositionalFluidState fs;
+ fs.setTemperature(temperature);
+ fs.setPressure(FluidSystem::oilPhaseIdx, pressure);
+ fs.setPressure(FluidSystem::gasPhaseIdx, pressure);
+ for (int c = 0; c < numComponents; ++c) {
+ fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, z[c]);
+ fs.setMoleFraction(FluidSystem::oilPhaseIdx, c, z[c]);
+ }
+
+ typename FluidSystem::template ParameterCache paramCache(eosType);
+ paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
+
+ PhaseState state;
+ state.density = FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx);
+ for (int c = 0; c < numComponents; ++c) {
+ state.fugacityCoefficient[c] =
+ FluidSystem::fugacityCoefficient(fs, paramCache, FluidSystem::gasPhaseIdx, c);
+ }
+ return state;
+}
+
+} // Anonymous namespace
+
+BOOST_GLOBAL_FIXTURE(Fixture);
+
+BOOST_AUTO_TEST_CASE(ShiftMovesTheDensityOntoTheReference)
+{
+ const auto gas = gasState();
+
+ // The reference reports 165.87 kg/m3 at the top of the column. The
+ // composition here is the one COMPVD states rather than the equilibrated
+ // one of that exact depth, which is worth a few tenths of a percent; the
+ // unshifted equation of state gives 172.84, so the tolerance is far
+ // tighter than the effect being tested.
+ BOOST_CHECK_CLOSE(gas.density, 165.87, 1.0);
+}
+
+BOOST_AUTO_TEST_CASE(ShiftLeavesTheEquilibriumRatiosAlone)
+{
+ // The invariant of a volume translation is the equilibrium ratio, not the
+ // fugacity coefficient itself. Peneloux scales the fugacity of component
+ // c by exp(-c_c p / (R T)), a factor set by the component and the state
+ // but not by the phase, so it cancels from K_c = phi_c^liquid /
+ // phi_c^vapour and leaves the phase split untouched. That ratio is what
+ // this checks, against a separate fluid system holding the same
+ // components with zero shifts.
+ Opm::CompositionalFluidState fs;
+ Opm::CompositionalFluidState fsRef;
+ fs.setTemperature(temperature);
+ fsRef.setTemperature(temperature);
+ for (int ph : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) {
+ fs.setPressure(ph, pressure);
+ fsRef.setPressure(ph, pressure);
+ }
+ for (int c = 0; c < numComponents; ++c) {
+ for (int ph : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) {
+ fs.setMoleFraction(ph, c, z[c]);
+ fsRef.setMoleFraction(ph, c, z[c]);
+ }
+ }
+
+ typename FluidSystem::template ParameterCache pc(eosType);
+ typename UnshiftedSystem::template ParameterCache pcRef(eosType);
+ for (int ph : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) {
+ pc.updatePhase(fs, ph);
+ pcRef.updatePhase(fsRef, ph);
+ }
+
+ // The shift is real, and the reference system carries none of it.
+ BOOST_CHECK_GT(std::abs(pc.volumeShift(fs, FluidSystem::gasPhaseIdx)), 0.0);
+ BOOST_CHECK_SMALL(pcRef.volumeShift(fsRef, UnshiftedSystem::gasPhaseIdx), 1.0e-30);
+
+ for (int c = 0; c < numComponents; ++c) {
+ const Scalar K =
+ FluidSystem::fugacityCoefficient(fs, pc, FluidSystem::oilPhaseIdx, c) /
+ FluidSystem::fugacityCoefficient(fs, pc, FluidSystem::gasPhaseIdx, c);
+ const Scalar KRef =
+ UnshiftedSystem::fugacityCoefficient(fsRef, pcRef, UnshiftedSystem::oilPhaseIdx, c) /
+ UnshiftedSystem::fugacityCoefficient(fsRef, pcRef, UnshiftedSystem::gasPhaseIdx, c);
+ BOOST_CHECK_CLOSE(K, KRef, 1.0e-10);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(ShiftDoesNotReachTheCachedVolume)
+{
+ // The equilibrium ratios above survive because the coefficients are built
+ // from the unshifted root. Pinning that here keeps a future change from
+ // folding the shift into the cache, where it would corrupt the
+ // two-parameter expression the coefficients use rather than translate it.
+ Opm::CompositionalFluidState fs;
+ fs.setTemperature(temperature);
+ fs.setPressure(FluidSystem::oilPhaseIdx, pressure);
+ fs.setPressure(FluidSystem::gasPhaseIdx, pressure);
+ for (int c = 0; c < numComponents; ++c) {
+ fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, z[c]);
+ fs.setMoleFraction(FluidSystem::oilPhaseIdx, c, z[c]);
+ }
+ typename FluidSystem::template ParameterCache pc(eosType);
+ pc.updatePhase(fs, FluidSystem::gasPhaseIdx);
+
+ const Scalar Vm = pc.molarVolume(FluidSystem::gasPhaseIdx);
+ const Scalar shift = pc.volumeShift(fs, FluidSystem::gasPhaseIdx);
+ BOOST_CHECK_CLOSE(pc.correctedMolarVolume(fs, FluidSystem::gasPhaseIdx),
+ Vm - shift, 1.0e-10);
+ // The density follows the corrected volume, not the cached one.
+ BOOST_CHECK_CLOSE(FluidSystem::density(fs, pc, FluidSystem::gasPhaseIdx),
+ fs.averageMolarMass(FluidSystem::gasPhaseIdx) / (Vm - shift),
+ 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(ShiftMovesTheLiquidDensityToo)
+{
+ // The liquid phase is where the two-parameter equation of state is worst
+ // and the shift matters most, so it gets its own check. The expected
+ // shift is rebuilt here from the deck's SSHIFT and b_c = Omega_b R Tc/pc
+ // rather than read back from the parameter cache, so the whole chain
+ // (SSHIFT -> b_c -> corrected volume -> density) is verified independently.
+ Opm::CompositionalFluidState fs;
+ fs.setTemperature(temperature);
+ fs.setPressure(FluidSystem::oilPhaseIdx, pressure);
+ fs.setPressure(FluidSystem::gasPhaseIdx, pressure);
+ for (int c = 0; c < numComponents; ++c) {
+ fs.setMoleFraction(FluidSystem::oilPhaseIdx, c, z[c]);
+ fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, z[c]);
+ }
+ typename FluidSystem::template ParameterCache pc(eosType);
+ pc.updatePhase(fs, FluidSystem::oilPhaseIdx);
+
+ constexpr Scalar OmegaB = 0.0777960739; // Peng-Robinson
+ constexpr Scalar R = 8.31446261815324;
+ Scalar expectedShift = 0.0;
+ for (int c = 0; c < numComponents; ++c) {
+ const auto& p = components[c];
+ const Scalar b = OmegaB * R * p.criticalT / p.criticalP;
+ expectedShift += z[c] * shift[c] * b;
+ }
+ // The constants here are written out independently of the library's, so
+ // the agreement is limited by their last digits rather than by round-off;
+ // the effect under test is a ten percent change in density.
+ BOOST_CHECK_CLOSE(pc.volumeShift(fs, FluidSystem::oilPhaseIdx), expectedShift, 1.0e-3);
+
+ const Scalar Vm = pc.molarVolume(FluidSystem::oilPhaseIdx);
+ const Scalar rho = FluidSystem::density(fs, pc, FluidSystem::oilPhaseIdx);
+ BOOST_CHECK_CLOSE(rho,
+ fs.averageMolarMass(FluidSystem::oilPhaseIdx) / (Vm - expectedShift),
+ 1.0e-4);
+ // The shifts are mostly negative, so the correction enlarges the volume
+ // and the shifted liquid is the lighter of the two.
+ BOOST_CHECK_LT(rho, fs.averageMolarMass(FluidSystem::oilPhaseIdx) / Vm);
+}
+
+// A three-component system of its own, so parsing a deck into it cannot
+// disturb the seven-component one the tests above share.
+using DeckSystem = Opm::GenericOilGasWaterFluidSystem;
+
+BOOST_AUTO_TEST_CASE(ParsedShiftReachesTheFluidSystem)
+{
+ // The tests above hand the shifts to addComponent() directly, which leaves
+ // the parsing path unproven. This one goes through the deck: SSHIFT is
+ // read into the compositional configuration and initFromState() must carry
+ // it onto the component parameters.
+ const auto deck = Opm::Parser{}.parseString(R"(
+RUNSPEC
+METRIC
+DIMENS
+ 1 1 1 /
+COMPS
+3 /
+TABDIMS
+ 1 /
+OIL
+GAS
+WATER
+GRID
+DXV
+ 1 /
+DYV
+ 1 /
+DZV
+ 1 /
+DEPTHZ
+ 4*2000 /
+PROPS
+CNAMES
+ C1 C10 CO2 /
+TCRIT
+ 190.6 617.7 304.2 /
+PCRIT
+ 46.0 21.1 73.8 /
+VCRIT
+ 0.0990 0.6240 0.0940 /
+MW
+ 16.043 142.285 44.010 /
+ACF
+ 0.008 0.4885 0.225 /
+SSHIFT
+ -0.1595 0.10784 -0.0817 /
+SOLUTION
+SCHEDULE
+END
+)");
+
+ const auto eclState = Opm::EclipseState{ deck };
+ const auto schedule = Opm::Schedule{ deck, eclState };
+
+ BOOST_REQUIRE_NO_THROW(DeckSystem::initFromState(eclState, schedule));
+
+ const std::array expected{-0.1595, 0.10784, -0.0817};
+ for (unsigned c = 0; c < 3; ++c) {
+ BOOST_CHECK_CLOSE(DeckSystem::volumeShift(c), expected[c], 1.0e-10);
+ }
+}
+
+BOOST_AUTO_TEST_CASE(AnOverlargeShiftIsRejectedRatherThanReturned)
+{
+ // SSHIFT is unconstrained input. A shift bigger than the molar volume
+ // makes the corrected volume non-positive, where the density would come
+ // back negative or infinite instead of failing.
+ using BigShift = Opm::GenericOilGasWaterFluidSystem;
+ using CompParam = typename BigShift::ComponentParam;
+ BigShift::init();
+ BigShift::addComponent(CompParam{"C1", 0.016043, 190.60, 45.40e5, 0.099, 0.008,
+ /*volume_shift=*/1.0e3});
+
+ Opm::CompositionalFluidState fs;
+ fs.setTemperature(Scalar{400});
+ fs.setPressure(BigShift::oilPhaseIdx, Scalar{100e5});
+ fs.setPressure(BigShift::gasPhaseIdx, Scalar{100e5});
+ fs.setMoleFraction(BigShift::gasPhaseIdx, 0, Scalar{1});
+ fs.setMoleFraction(BigShift::oilPhaseIdx, 0, Scalar{1});
+
+ typename BigShift::template ParameterCache pc(eosType);
+ pc.updatePhase(fs, BigShift::gasPhaseIdx);
+
+ BOOST_CHECK_THROW(BigShift::density(fs, pc, BigShift::gasPhaseIdx),
+ Opm::NumericalProblem);
+}
diff --git a/tests/parser/InitConfigTest.cpp b/tests/parser/InitConfigTest.cpp
index b4490703e61..99ca37a43a8 100644
--- a/tests/parser/InitConfigTest.cpp
+++ b/tests/parser/InitConfigTest.cpp
@@ -1,5 +1,6 @@
/*
Copyright 2015 Statoil ASA.
+ Copyright 2026 SINTEF Digital
This file is part of the Open Porous Media project (OPM).
@@ -251,6 +252,52 @@ SCHEDULE
)" };
}
+ std::string deckWithCompositionalEquilDefaultedItem11()
+ {
+ return { R"(RUNSPEC
+METRIC
+DIMENS
+ 10 10 10 /
+EQLDIMS
+1 100 20 1 1 /
+OIL
+WATER
+GAS
+COMPS
+3 /
+SOLUTION
+EQUIL
+ 2050 150 2300 0 2050 0 3* 3 /
+GRID
+START -- 0
+19 JUN 2007 /
+SCHEDULE
+)" };
+ }
+
+ std::string deckWithCompositionalEquilItem11One()
+ {
+ return { R"(RUNSPEC
+METRIC
+DIMENS
+ 10 10 10 /
+EQLDIMS
+1 100 20 1 1 /
+OIL
+WATER
+GAS
+COMPS
+3 /
+SOLUTION
+EQUIL
+ 2050 150 2300 0 2050 0 3* 3 1 /
+GRID
+START -- 0
+19 JUN 2007 /
+SCHEDULE
+)" };
+ }
+
std::string deckWithStrEquil()
{
return { R"(RUNSPEC
@@ -468,6 +515,37 @@ BOOST_AUTO_TEST_CASE(CompositionalEquilOperations)
BOOST_CHECK(record.setToSaturationPressure());
}
+BOOST_AUTO_TEST_CASE(CompositionalEquilDefaultedItem11)
+{
+ // EQUIL record:
+ // 2050 150 2300 0 2050 0 3* 3 /
+ // Item 10: COMP_INIT_TYPE = 3
+ // Item 11: defaulted => saturation pressure IS set
+ const auto deck = createDeck(deckWithCompositionalEquilDefaultedItem11());
+ const Runspec runspec(deck);
+ const InitConfig config(deck, runspec.phases(), runspec.compositionalMode());
+
+ const auto& record = config.getEquil().getRecord(0);
+
+ BOOST_CHECK_EQUAL(3, record.compositionalInitType());
+ BOOST_CHECK(record.setToSaturationPressure());
+}
+
+BOOST_AUTO_TEST_CASE(CompositionalEquilItem11One)
+{
+ // EQUIL record:
+ // 2050 150 2300 0 2050 0 3* 3 1 /
+ // Item 11: COMP_NOT_SET_SAT_PRESSURE = 1 => datum pressure is kept
+ const auto deck = createDeck(deckWithCompositionalEquilItem11One());
+ const Runspec runspec(deck);
+ const InitConfig config(deck, runspec.phases(), runspec.compositionalMode());
+
+ const auto& record = config.getEquil().getRecord(0);
+
+ BOOST_CHECK_EQUAL(3, record.compositionalInitType());
+ BOOST_CHECK(!record.setToSaturationPressure());
+}
+
BOOST_AUTO_TEST_CASE(StrEquilOperations)
{
const auto deck = createDeck(deckWithStrEquil());
diff --git a/tests/parser/ScheduleTests.cpp b/tests/parser/ScheduleTests.cpp
index 4ced860661c..fd8c81fb13a 100644
--- a/tests/parser/ScheduleTests.cpp
+++ b/tests/parser/ScheduleTests.cpp
@@ -5286,6 +5286,47 @@ SCHEDULE
}
}
+BOOST_AUTO_TEST_CASE(WELLSTRE_rounded_composition_is_normalized) {
+ // An injection stream written with a few digits sums to 0.999968 rather
+ // than one. It is accepted, scaled, and reaches the well through WINJGAS
+ // with its ratios intact.
+ const auto sched = make_schedule(gptable_deck(R"(
+WELSPECS
+ 'INJ' 'G1' 1 1 2000 'GAS' /
+/
+WELLSTRE
+ 'STR1' 0.749520 0.240448 0.010000 /
+/
+WINJGAS
+ 'INJ' 'STREAM' 'STR1' /
+/
+WCONINJE
+ 'INJ' 'GAS' 'OPEN' 'RATE' 100 /
+/
+TSTEP
+ 1 /
+)"));
+
+ const auto& composition = sched.getWell("INJ", 0).getInjectionProperties().gasInjComposition();
+ BOOST_REQUIRE_EQUAL(composition.size(), std::size_t{3});
+
+ const double sum = std::accumulate(composition.begin(), composition.end(), 0.0);
+ BOOST_CHECK_CLOSE(sum, 1.0, 1.0e-12);
+
+ // Scaling preserves the ratios of the stream.
+ BOOST_CHECK_CLOSE(composition[0] / composition[1], 0.749520 / 0.240448, 1.0e-10);
+ BOOST_CHECK_CLOSE(composition[0], 0.749520 / 0.999968, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(WELLSTRE_composition_far_from_one_is_rejected) {
+ // A sum too far from one to be the rounding of the values is still an error.
+ BOOST_CHECK_THROW(make_schedule(gptable_deck(R"(
+WELLSTRE
+ 'STR1' 0.70 0.20 0.05 /
+/
+)")), Opm::OpmInputError);
+}
+
BOOST_AUTO_TEST_CASE(GPTABLE_solution_seed_and_schedule_respec) {
const auto sched = make_schedule(R"(
RUNSPEC
diff --git a/tests/parser/TableManagerTests.cpp b/tests/parser/TableManagerTests.cpp
index 0ff83e82995..4e50b44b6a9 100644
--- a/tests/parser/TableManagerTests.cpp
+++ b/tests/parser/TableManagerTests.cpp
@@ -23,6 +23,9 @@
#include
+#include
+#include
+#include
#include
#include
@@ -3420,6 +3423,154 @@ END
}
+BOOST_AUTO_TEST_CASE(CompvdTable_NormalizationIsReported) {
+ // The row that was scaled must be named, with the sum it had before.
+ const auto deck = Opm::Parser{}.parseString(R"(
+RUNSPEC
+METRIC
+COMPS
+7 /
+EQLDIMS
+1 /
+PROPS
+COMPVD
+ 2573.5 0.73779563 0.04261 0.122439 0.04501 0.04351 0.007502 0.0011211 0 277.56 /
+END
+)");
+
+ std::ostringstream warnings;
+ Opm::OpmLog::addBackend("STREAM",
+ std::make_shared(warnings, Opm::Log::MessageType::Warning));
+ const auto tmgr = Opm::TableManager{ deck };
+ Opm::OpmLog::removeBackend("STREAM");
+
+ const std::string text = warnings.str();
+ BOOST_CHECK(text.find("row 1 of COMPVD table 1") != std::string::npos);
+ BOOST_CHECK(text.find("0.99998") != std::string::npos);
+ BOOST_CHECK(text.find("normalized") != std::string::npos);
+}
+
+BOOST_AUTO_TEST_CASE(ZmfvdTable_ExactRowIsNotReported) {
+ // A row that already sums to one is scaled by one, and says nothing.
+ const auto deck = Opm::Parser{}.parseString(R"(
+RUNSPEC
+METRIC
+COMPS
+3 /
+EQLDIMS
+1 /
+PROPS
+ZMFVD
+ 2000.0 0.3 0.3 0.4 /
+END
+)");
+
+ std::ostringstream warnings;
+ Opm::OpmLog::addBackend("STREAM",
+ std::make_shared(warnings, Opm::Log::MessageType::Warning));
+ const auto tmgr = Opm::TableManager{ deck };
+ Opm::OpmLog::removeBackend("STREAM");
+
+ BOOST_CHECK(warnings.str().find("normalized") == std::string::npos);
+}
+
+BOOST_AUTO_TEST_CASE(ZmfvdTable_TooFarFromOneIsRejected) {
+ // A row that misses one by more than rounding is still an error.
+ const auto deck = Opm::Parser{}.parseString(R"(
+RUNSPEC
+METRIC
+COMPS
+3 /
+EQLDIMS
+1 /
+PROPS
+ZMFVD
+ 2000.0 0.30 0.30 0.35 /
+END
+)");
+
+ BOOST_CHECK_THROW(Opm::TableManager{ deck }, Opm::OpmInputError);
+}
+
+BOOST_AUTO_TEST_CASE(CompvdTable_RoundedCompositionIsNormalized) {
+ // A lumped characterization quoted to a handful of digits sums to one only
+ // to within its rounding; the row is accepted and stored normalized. This
+ // row is taken from a customer deck and sums to 0.99998773.
+ const auto deck = Opm::Parser{}.parseString(R"(
+RUNSPEC
+METRIC
+COMPS
+7 /
+EQLDIMS
+1 /
+PROPS
+COMPVD
+ 2573.5 0.73779563 0.04261 0.122439 0.04501 0.04351 0.007502 0.0011211 0 277.56 /
+END
+)");
+
+ const auto tmgr = Opm::TableManager{ deck };
+ const auto& compvd = tmgr.getCompvdTables();
+ BOOST_REQUIRE_EQUAL(compvd.size(), std::size_t{1});
+
+ const auto& table = compvd.getTable(0);
+ double sum = 0.0;
+ for (int c = 0; c < 7; ++c) {
+ sum += table.getMoleFractionColumn(c)[0];
+ }
+ BOOST_CHECK_CLOSE(sum, 1.0, 1.0e-12);
+
+ // The normalization is a rescaling, so the component ratios are untouched.
+ BOOST_CHECK_CLOSE(table.getMoleFractionColumn(0)[0] / table.getMoleFractionColumn(1)[0],
+ 0.73779563 / 0.04261, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(ZmfvdTable_RoundedCompositionIsNormalized) {
+ const auto deck = Opm::Parser{}.parseString(R"(
+RUNSPEC
+METRIC
+COMPS
+3 /
+EQLDIMS
+1 /
+PROPS
+ZMFVD
+ 2000.0 0.30001 0.29999 0.39998 /
+END
+)");
+
+ const auto tmgr = Opm::TableManager{ deck };
+ const auto& zmfvd = tmgr.getZmfvdTables();
+ const auto& table = zmfvd.getTable(0);
+
+ double sum = 0.0;
+ for (int c = 0; c < 3; ++c) {
+ sum += table.getMoleFractionColumn(c)[0];
+ }
+ BOOST_CHECK_CLOSE(sum, 1.0, 1.0e-12);
+}
+
+BOOST_AUTO_TEST_CASE(ZmfvdTable_NonFiniteMoleFractionIsRejected) {
+ // NaN is an acceptable floating-point token to the parser. It must not
+ // reach the normalization: a NaN sum makes every tolerance comparison
+ // false, so an unguarded helper would divide the row through by it and
+ // store a composition of NaNs.
+ const auto deck = Opm::Parser{}.parseString(R"(
+RUNSPEC
+METRIC
+COMPS
+2 /
+EQLDIMS
+1 /
+PROPS
+ZMFVD
+ 100.0 NaN 0.5 /
+END
+)");
+
+ BOOST_CHECK_THROW(Opm::TableManager{ deck }, Opm::OpmInputError);
+}
+
BOOST_AUTO_TEST_CASE(CompvdTable_MoleFractionsMustSumToOne) {
const auto deck = Opm::Parser{}.parseString(R"(
RUNSPEC