Test the PSAT phase selection and reduce its failure count over ranks

Move the choice between the cell pressure, the bubble point and the dew
point into CompositionalContainer::cellSaturationPressure(), which takes
plain values instead of a fluid state, and cover its four outcomes in
the container test: both phases present, oil only, gas only, and a
supercritical mixture with no saturation pressure. The output module
only unpacks the fluid state and writes the result.

Count the cells without a saturation pressure in validateLocalData(),
which every rank reaches after the per-cell pass, and reduce it over
the grid communicator before logging the total once from rank 0. The
previous report at output time was the I/O rank's count alone.

Say at both sites what zero stands for: the solver reports that it
could not determine a saturation pressure, which a supercritical
mixture also does; it does not assert that none exists.

Author: GitPaean

opm/simulators/flow/CompositionalContainer.cpp

@@ -23,11 +23,15 @@
 #include <config.h>
 #include <opm/simulators/flow/CompositionalContainer.hpp>
 
+#include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
+
+#include <opm/material/constraintsolvers/SaturationPressure.hpp>
 #include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
 
 #include <opm/output/data/Solution.hpp>
 
 #include <algorithm>
+#include <optional>
 #include <tuple>
 
 #include <fmt/format.h>
@@ -237,6 +241,36 @@ outputRestart(data::Solution& sol,
     this->allocated_ = false;
 }
 
+template<class FluidSystem>
+std::optional<typename CompositionalContainer<FluidSystem>::Scalar>
+CompositionalContainer<FluidSystem>::
+cellSaturationPressure(const Scalar oilSaturation,
+                       const Scalar gasSaturation,
+                       const Scalar oilPressure,
+                       const std::array<Scalar, numComponents>& moleFractions,
+                       const Scalar temperature,
+                       const CompositionalConfig::EOSType eosType)
+{
+    if (oilSaturation > 0.0 && gasSaturation > 0.0) {
+        return oilPressure;
+    }
+
+    // The zero-component instantiation exists only to register a parameter
+    // and has no equation of state to solve.
+    if constexpr (numComponents > 0) {
+        using SatP = SaturationPressure<Scalar, FluidSystem>;
+        typename SatP::CompVec incipient;
+        Scalar psat = 0.0;
+        const bool found = (gasSaturation <= 0.0)
+            ? SatP::bubblePressure(moleFractions, temperature, eosType, psat, incipient)
+            : SatP::dewPressure(moleFractions, temperature, eosType, psat, incipient);
+        if (found) {
+            return psat;
+        }
+    }
+    return std::nullopt;
+}
+
 #define INSTANTIATE_COMP_THREEPHASE(NUM) \
     template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM, true>; \
     template class CompositionalContainer<FS##NUM<double>>;

opm/simulators/flow/CompositionalContainer.hpp

@@ -26,9 +26,12 @@
 #ifndef OPM_COMPOSITIONAL_CONTAINER_HPP
 #define OPM_COMPOSITIONAL_CONTAINER_HPP
 
+#include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
+
 #include <array>
 #include <functional>
 #include <map>
+#include <optional>
 #include <string>
 #include <vector>
 
@@ -72,6 +75,20 @@ class CompositionalContainer
     void assignSaturationPressure(const unsigned globalDofIdx,
                                   const Scalar psat);
 
+    /// The saturation pressure of a cell, or nothing when the solver could
+    /// not determine one.  A cell holding both hydrocarbon phases is at its
+    /// saturation pressure; a single-phase cell has the bubble- (dew-) point
+    /// pressure of its total composition.  Takes plain values rather than a
+    /// fluid state so the phase selection and the solver wiring can be
+    /// exercised without a simulator.
+    static std::optional<Scalar>
+    cellSaturationPressure(const Scalar oilSaturation,
+                           const Scalar gasSaturation,
+                           const Scalar oilPressure,
+                           const std::array<Scalar, numComponents>& moleFractions,
+                           const Scalar temperature,
+                           const CompositionalConfig::EOSType eosType);
+
     void assignVaporFraction(const unsigned globalDofIdx,
                              const Scalar vmf);
 

opm/simulators/flow/OutputCompositionalModule.hpp

@@ -42,7 +42,6 @@
 #include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 
 #include <opm/material/common/Valgrind.hpp>
-#include <opm/material/constraintsolvers/SaturationPressure.hpp>
 
 #include <opm/models/blackoil/blackoilproperties.hh>
 #include <opm/models/common/multiphasebaseproperties.hh>
@@ -202,15 +201,6 @@ class OutputCompositionalModule : public GenericOutputModule<GetPropType<TypeTag
 
     void assignToSolution(data::Solution& sol) override
     {
-        if (this->numFailedSaturationPressures_ > 0) {
-            // The solver reports both an absent root and a failed search as
-            // zero. This count covers the current local output snapshot.
-            OpmLog::info(fmt::format("No saturation pressure for {} cells; "
-                                     "PSAT is reported as zero there.",
-                                     this->numFailedSaturationPressures_));
-            this->numFailedSaturationPressures_ = 0;
-        }
-
         this->compC_.outputRestart(sol, this->saturation_[oilPhaseIdx]);
         BaseType::assignToSolution(sol);
 
@@ -738,6 +728,22 @@ class OutputCompositionalModule : public GenericOutputModule<GetPropType<TypeTag
         this->assignSaturationPressure_(globalDofIdx, intQuants.fluidState());
     }
 
+    /// Runs on every rank once the per-cell pass is complete, so this is
+    /// where the count of cells without a saturation pressure is reduced.
+    /// Reporting it at output time would report the I/O rank's cells alone:
+    /// the output is assembled there, while the count is gathered per rank.
+    void validateLocalData()
+    {
+        const auto& comm = this->simulator_.gridView().comm();
+        const auto total = comm.sum(this->numFailedSaturationPressures_);
+        this->numFailedSaturationPressures_ = 0;
+        if (total > 0 && comm.rank() == 0) {
+            OpmLog::info(fmt::format("No saturation pressure could be determined for {} "
+                                     "cells; PSAT is written as zero there.", total));
+        }
+        BaseType::validateLocalData();
+    }
+
 protected:
     //! \brief Allocate compositional relative-permeability buffers.
     void allocFormulationBuffers(std::map<std::string, int>& rstKeywords,
@@ -812,37 +818,29 @@ class OutputCompositionalModule : public GenericOutputModule<GetPropType<TypeTag
             return;
         }
 
-        const Scalar sOil = getValue(fs.saturation(oilPhaseIdx));
-        const Scalar sGas = getValue(fs.saturation(gasPhaseIdx));
-
-        Scalar psat = 0.0;
-        if (sOil > 0.0 && sGas > 0.0) {
-            psat = getValue(fs.pressure(oilPhaseIdx));
+        std::array<Scalar, numComponents> z;
+        for (int c = 0; c < numComponents; ++c) {
+            z[c] = getValue(fs.moleFraction(c));
         }
-        else {
-            using SatP = SaturationPressure<Scalar, FluidSystem>;
-            typename SatP::CompVec z;
-            typename SatP::CompVec incipient;
-            for (int c = 0; c < numComponents; ++c) {
-                z[c] = getValue(fs.moleFraction(c));
-            }
-            const Scalar temp = getValue(fs.temperature(oilPhaseIdx));
-            const bool converged = (sGas <= 0.0)
-                ? SatP::bubblePressure(z, temp, eosType_, psat, incipient)
-                : SatP::dewPressure(z, temp, eosType_, psat, incipient);
-            if (!converged) {
-                // A supercritical mixture genuinely has no saturation pressure,
-                // and so does a search that gave up; both report zero.  Count
-                // them so that a systematically failing solve is not silent.
-                psat = 0.0;
+        const auto psat = CompositionalContainer<FluidSystem>::cellSaturationPressure(
+            getValue(fs.saturation(oilPhaseIdx)),
+            getValue(fs.saturation(gasPhaseIdx)),
+            getValue(fs.pressure(oilPhaseIdx)),
+            z,
+            getValue(fs.temperature(oilPhaseIdx)),
+            eosType_);
+        if (!psat) {
+            // The solver never asserts that no saturation pressure exists; it
+            // reports that it could not determine one, which a supercritical
+            // mixture, where none exists, also does.  Zero stands for that,
+            // and the cells are counted so a field-wide total is visible.
 #ifdef _OPENMP
 #pragma omp atomic
 #endif
-                ++this->numFailedSaturationPressures_;
-            }
+            ++this->numFailedSaturationPressures_;
         }
 
-        this->compC_.assignSaturationPressure(globalDofIdx, psat);
+        this->compC_.assignSaturationPressure(globalDofIdx, psat.value_or(Scalar{0}));
     }
 
     const Simulator& simulator_;

tests/test_compositionalcontainer.cpp

@@ -24,17 +24,37 @@
 
 #include <opm/simulators/flow/CompositionalContainer.hpp>
 
+#include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
+
 #include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
 #include <opm/output/data/Solution.hpp>
 
+#include <array>
 #include <map>
+#include <optional>
 #include <string>
 #include <vector>
 
 namespace {
 using FluidSystem = Opm::GenericOilGasWaterFluidSystem<double, 3, false>;
 using Container = Opm::CompositionalContainer<FluidSystem>;
-}
+
+// The three-component CO2/methane/decane fluid of the saturation-pressure
+// solver's own tests, so their reference values can be reused here.
+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});
+    }
+};
+} // anonymous namespace
+
+BOOST_GLOBAL_FIXTURE(Fixture);
 
 BOOST_AUTO_TEST_CASE(SaturationPressureRequiresRestartOutput)
 {
@@ -85,3 +105,47 @@ BOOST_AUTO_TEST_CASE(DisabledSaturationPressureClearsPreviousRequest)
     container.allocate(2, keywords, /*isRestartOutput=*/true);
     BOOST_CHECK(!container.saturationPressureAllocated());
 }
+
+BOOST_AUTO_TEST_CASE(CellSaturationPressureSelectsThePhaseAndTheBranch)
+{
+    // 100 degC with Peng-Robinson, as in the solver's tests.
+    constexpr double temperature = 373.15;
+    constexpr auto eos = Opm::CompositionalConfig::EOSType::PR;
+    using CompVec = std::array<double, 3>;
+
+    // Both phases present: the cell is at its saturation pressure and the
+    // solver is not consulted.
+    {
+        const auto psat = Container::cellSaturationPressure(
+            0.3, 0.7, 75.0e5, CompVec{0.0, 0.5, 0.5}, temperature, eos);
+        BOOST_REQUIRE(psat.has_value());
+        BOOST_CHECK_CLOSE(*psat, 75.0e5, 1.0e-12);
+    }
+
+    // Oil only: the bubble point of the total composition.  The reference
+    // simulator puts this liquid at 160.5601 bar.
+    {
+        const auto psat = Container::cellSaturationPressure(
+            1.0, 0.0, 150.0e5, CompVec{0.0, 0.5, 0.5}, temperature, eos);
+        BOOST_REQUIRE(psat.has_value());
+        BOOST_CHECK_CLOSE(*psat / 1.0e5, 160.56010, 1.0e-3);
+    }
+
+    // Gas only: the retrograde dew point of that liquid's equilibrium vapour,
+    // which recovers the same pressure.  The vapour is known to single
+    // precision, hence the tolerance.
+    {
+        const auto psat = Container::cellSaturationPressure(
+            0.0, 1.0, 150.0e5, CompVec{0.0, 0.987784, 0.012216}, temperature, eos);
+        BOOST_REQUIRE(psat.has_value());
+        BOOST_CHECK_CLOSE(*psat / 1.0e5, 160.56010, 1.0e-2);
+    }
+
+    // Pure methane far above its critical temperature has no saturation
+    // pressure, and the cell pressure must not stand in for one.
+    {
+        const auto psat = Container::cellSaturationPressure(
+            0.0, 1.0, 150.0e5, CompVec{0.0, 1.0, 0.0}, temperature, eos);
+        BOOST_CHECK(!psat.has_value());
+    }
+}