[test] compositional equilibration combination

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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.hpp>
+
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/utility/OpmInputError.hpp>
+
+#include <cmath>
+#include <limits>
+#include <numeric>
+
+#include <fmt/format.h>
+
+namespace Opm {
+
+double moleFractionTolerance()
+{
+    return 1.0e-4;
+}
+
+double exactSumSlack(const std::size_t numValues)
+{
+    return 2.0 * numValues * std::numeric_limits<double>::epsilon();
+}
+
+void normalizeMoleFractions(std::vector<double>& 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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_NORMALIZE_MOLE_FRACTIONS_HPP
+#define OPM_NORMALIZE_MOLE_FRACTIONS_HPP
+
+#include <cstddef>
+#include <string>
+#include <vector>
+
+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<double>& fractions,
+                            const std::string& what,
+                            const KeywordLocation& location);
+
+} // namespace Opm
+
+#endif // OPM_NORMALIZE_MOLE_FRACTIONS_HPP

opm/input/eclipse/EclipseState/InitConfig/Equil.cpp

@@ -78,7 +78,10 @@ namespace Opm {
         if (compositional) {
             comp_init_type = record.getItem<ParserKeywords::EQUIL::COMP_INIT_TYPE>().get<int>(0);
             if (comp_init_type == 2 || comp_init_type == 3) {
-                set_to_saturation_pressure = record.getItem<ParserKeywords::EQUIL::COMP_NOT_SET_SAT_PRESSURE>().get<int>(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<ParserKeywords::EQUIL::COMP_NOT_SET_SAT_PRESSURE>();
+                set_to_saturation_pressure = !item.hasValue(0) || (item.get<int>(0) != 1);
             }
         }
     }

opm/input/eclipse/EclipseState/Tables/Tables.cpp

@@ -95,6 +95,8 @@
 #include <opm/input/eclipse/Units/UnitSystem.hpp>
 #include <opm/input/eclipse/Units/Units.hpp>
 
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.hpp>
 #include <opm/common/utility/OpmInputError.hpp>
 
 #include <opm/input/eclipse/Deck/Deck.hpp>
@@ -114,6 +116,7 @@
 #include <cstddef>
 #include <functional>
 #include <initializer_list>
+#include <limits>
 #include <numeric>
 #include <stdexcept>
 #include <string>
@@ -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<double> 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<double> 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<double>(compIdx);
-            moles[c] = mole_fraction;
-            getColumn(1 + c).addValue(mole_fraction, tableName);
+            moles[c] = item.get<double>(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<double>();
 
     const std::string tableName{"COMPVD"};
+    std::vector<double> 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<double> 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.

opm/input/eclipse/Schedule/Well/WellKeywordHandlers.cpp

@@ -36,6 +36,7 @@
 #include <opm/input/eclipse/Schedule/Well/WListManager.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPDP.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPEXP.hpp>
+#include <opm/input/eclipse/EclipseState/Compositional/NormalizeMoleFractions.hpp>
 #include <opm/input/eclipse/Schedule/Well/Well.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellEconProductionLimits.hpp>
@@ -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<ParserKeywords::WELLSTRE::STREAM>().getTrimmedString(0);
-        const auto& composition = record.getItem<ParserKeywords::WELLSTRE::COMPOSITIONS>().getSIDoubleData();
+        auto composition = record.getItem<ParserKeywords::WELLSTRE::COMPOSITIONS>().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<double>::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<std::vector<double>>(composition);
         inj_streams.update(stream_name, std::move(composition_ptr));
     }