Write compositional saturation pressure to restart output

Handle PSAT requests in compositional restart output. Use the cell
pressure when both hydrocarbon phases are present, solve the bubble or
dew pressure for exact single-phase flash labels, and write zero when no
boundary can be resolved.

Allocate PSAT only for restart output and calculate it in the OpenMP
cell loop. Make output validation virtual and run the MPI reductions
after rank-wide exception handling, preventing a local failure from
leaving ranks in different collectives. Reduce the failure count across
ranks and log it once with clear singular and plural wording.

Keep the phase decision independent of simulator state and test restart
allocation, phase labels near round-off boundaries, solver results, and
the failure contract.

Author: GitPaean

CMakeLists_files.cmake

@@ -468,6 +468,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_aqantrc_flow_keyword.cpp
   tests/test_blackoil_amg.cpp
   tests/test_blackoilprimaryvariables.cpp
+  tests/test_compositionalcontainer.cpp
   tests/test_compwell_equations.cpp
   tests/test_compwell_jacobian.cpp
   tests/test_convergenceoutputconfiguration.cpp

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>
@@ -37,7 +41,8 @@ namespace Opm {
 template<class FluidSystem>
 void CompositionalContainer<FluidSystem>::
 allocate(const unsigned bufferSize,
-         std::map<std::string, int>& rstKeywords)
+         std::map<std::string, int>& rstKeywords,
+         const bool isRestartOutput)
 {
     if (auto& zmf = rstKeywords["ZMF"]; zmf > 0) {
         this->allocated_ = true;
@@ -75,6 +80,17 @@ allocate(const unsigned bufferSize,
         gasPressure_.resize(bufferSize, 0.0);
     }
 
+    // Summary-only substeps must not retain a PSAT buffer from a previous
+    // preparation pass: its presence enables a nonlinear solve in every cell.
+    saturationPressure_.clear();
+    if (auto& psat = rstKeywords["PSAT"]; psat > 0) {
+        psat = 0;
+        if (isRestartOutput) {
+            this->allocated_ = true;
+            saturationPressure_.resize(bufferSize, 0.0);
+        }
+    }
+
     if (auto& vmf = rstKeywords["VMF"]; vmf > 0) {
         this->allocated_ = true;
         vmf = 0;
@@ -140,6 +156,16 @@ assignPhasePressures(const unsigned globalDofIdx,
     }
 }
 
+template<class FluidSystem>
+void CompositionalContainer<FluidSystem>::
+assignSaturationPressure(const unsigned globalDofIdx,
+                         const Scalar psat)
+{
+    if (!saturationPressure_.empty()) {
+        saturationPressure_[globalDofIdx] = psat;
+    }
+}
+
 template<class FluidSystem>
 void CompositionalContainer<FluidSystem>::
 assignVaporFraction(const unsigned globalDofIdx,
@@ -205,6 +231,7 @@ outputRestart(data::Solution& sol,
 
     entries.emplace_back("POIL", UnitSystem::measure::pressure, oilPressure_);
     entries.emplace_back("PGAS", UnitSystem::measure::pressure, gasPressure_);
+    entries.emplace_back("PSAT", UnitSystem::measure::pressure, saturationPressure_);
     entries.emplace_back("VMF", UnitSystem::measure::identity, vaporFraction_);
 
     std::ranges::for_each(entries,
@@ -214,6 +241,38 @@ outputRestart(data::Solution& sol,
     this->allocated_ = false;
 }
 
+template<class FluidSystem>
+auto CompositionalContainer<FluidSystem>::
+cellSaturationPressure(const Scalar liquidFraction,
+                       const Scalar oilPressure,
+                       const std::array<Scalar, numComponents>& moleFractions,
+                       const Scalar temperature,
+                       const CompositionalConfig::EOSType eosType) -> std::optional<Scalar>
+{
+    // Compare the flash's exact single-phase labels: a two-phase Rachford-Rice
+    // result can round slightly outside (0, 1).
+    const bool liquidOnly = (liquidFraction == 1.0);
+    const bool vapourOnly = (liquidFraction == 0.0);
+    if (!liquidOnly && !vapourOnly) {
+        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 Solver = SaturationPressure<Scalar, FluidSystem>;
+        typename Solver::CompVec incipient;
+        Scalar psat = 0.0;
+        const bool found = liquidOnly
+            ? Solver::bubblePressure(moleFractions, temperature, eosType, psat, incipient)
+            : Solver::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>
 
@@ -51,7 +54,8 @@ class CompositionalContainer
 
 public:
     void allocate(const unsigned bufferSize,
-                  std::map<std::string, int>& rstKeywords);
+                  std::map<std::string, int>& rstKeywords,
+                  const bool isRestartOutput);
 
     using AssignFunction = std::function<Scalar(const unsigned)>;
 
@@ -68,6 +72,24 @@ class CompositionalContainer
                               const Scalar oilPressure,
                               const Scalar gasPressure);
 
+    void assignSaturationPressure(const unsigned globalDofIdx,
+                                  const Scalar psat);
+
+    /// Return the cell pressure for two hydrocarbon phases, or the bubble/dew
+    /// pressure of the total composition for a single phase. Return std::nullopt
+    /// if the solver cannot determine a saturation pressure.
+    ///
+    /// \p liquidFraction is the flash liquid fraction L: exactly one denotes
+    /// liquid only, exactly zero vapour only, and other values denote two phases.
+    /// Use L because computed phase saturations can contain round-off residuals.
+    /// Plain-value arguments allow testing phase selection without a simulator.
+    [[nodiscard]] static std::optional<Scalar>
+    cellSaturationPressure(const Scalar liquidFraction,
+                           const Scalar oilPressure,
+                           const std::array<Scalar, numComponents>& moleFractions,
+                           const Scalar temperature,
+                           const CompositionalConfig::EOSType eosType);
+
     void assignVaporFraction(const unsigned globalDofIdx,
                              const Scalar vmf);
 
@@ -89,6 +111,9 @@ class CompositionalContainer
     bool vaporFractionAllocated() const
     { return !vaporFraction_.empty(); }
 
+    bool saturationPressureAllocated() const
+    { return !saturationPressure_.empty(); }
+
     bool allocated() const
     { return allocated_; }
 
@@ -101,6 +126,8 @@ class CompositionalContainer
     // phase pressures (POIL, PGAS)
     ScalarBuffer oilPressure_;
     ScalarBuffer gasPressure_;
+    // saturation pressure (PSAT)
+    ScalarBuffer saturationPressure_;
     // vapour mole fraction of the total mixture (VMF)
     ScalarBuffer vaporFraction_;
 };

opm/simulators/flow/EclWriter.hpp

@@ -863,8 +863,6 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
                 this->outputModule_->processElementBlockData(elemCtx);
             }
             this->outputModule_->clearExtractors();
-
-            this->outputModule_->accumulateDensityParallel();
         }
 
         {
@@ -881,10 +879,12 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
             }
         }
 
-        this->outputModule_->validateLocalData();
-
         OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ",
                                    this->simulator_.vanguard().grid().comm());
+
+        // Complete rank-wide exception handling before entering output collectives.
+        this->outputModule_->accumulateDensityParallel();
+        this->outputModule_->validateLocalData();
     }
 
     void captureLocalFluxData()

opm/simulators/flow/GenericOutputModule.hpp

@@ -286,7 +286,10 @@ class GenericOutputModule {
         local_data_valid_ = false;
     }
 
-    void validateLocalData(){
+    /// Mark per-cell output data valid. All ranks must call this method,
+    /// since overrides may reduce per-rank data.
+    virtual void validateLocalData()
+    {
         local_data_valid_ = true;
     }
 

opm/simulators/flow/OutputCompositionalModule.hpp

@@ -36,6 +36,7 @@
 #include <opm/common/TimingMacros.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
 
+#include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
 #include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 
 #include <opm/material/common/Valgrind.hpp>
@@ -52,6 +53,7 @@
 #include <opm/simulators/flow/OutputExtractor.hpp>
 
 #include <algorithm>
+#include <array>
 #include <cstddef>
 #include <fstream>
 #include <memory>
@@ -62,6 +64,7 @@
 #include <utility>
 #include <vector>
 
+#include <fmt/format.h>
 
 namespace Opm {
 
@@ -127,6 +130,7 @@ class OutputCompositionalModule : public GenericOutputModule<GetPropType<TypeTag
                    getPropValue<TypeTag, Properties::EnableBioeffects>(),
                    getPropValue<TypeTag, Properties::EnableGeochemistry>())
         , simulator_(simulator)
+        , eosType_(simulator.vanguard().eclState().compositionalConfig().eosType(0))
     {
         for (auto& region_pair : this->regions_) {
             this->createLocalRegion_(region_pair.second);
@@ -184,7 +188,10 @@ class OutputCompositionalModule : public GenericOutputModule<GetPropType<TypeTag
         }
 
         auto rstKeywords = this->schedule_.rst_keywords(reportStepNum);
-        this->compC_.allocate(bufferSize, rstKeywords);
+        const bool isRestartOutput =
+            isRestart || (!substep && this->schedule_.write_rst_file(reportStepNum));
+        this->compC_.allocate(bufferSize, rstKeywords, isRestartOutput);
+        this->numFailedSaturationPressures_ = 0;
 
         this->doAllocBuffers(bufferSize, reportStepNum, substep, log, isRestart,
                              /* hysteresisConfig = */ nullptr,
@@ -714,6 +721,26 @@ class OutputCompositionalModule : public GenericOutputModule<GetPropType<TypeTag
                                               intQuants,
                                               totVolume,
                                               referencePorosity);
+
+        // Run the nonlinear PSAT solve in the caller's OpenMP loop.
+        // The assignment skips cells when no restart buffer is allocated.
+        this->assignSaturationPressure_(globalDofIdx, intQuants.fluidState());
+    }
+
+    /// Reduce PSAT failures across all ranks before marking the output data valid.
+    void validateLocalData() override
+    {
+        const auto& comm = this->simulator_.gridView().comm();
+        const auto totalFailures = comm.sum(this->numFailedSaturationPressures_);
+        this->numFailedSaturationPressures_ = 0;
+        if (totalFailures > 0 && comm.rank() == 0) {
+            const auto* const cell = totalFailures == 1 ? "cell" : "cells";
+            OpmLog::info(fmt::format("Could not determine saturation pressure in {} {}; "
+                                     "PSAT is written as zero for every affected cell.",
+                                     totalFailures,
+                                     cell));
+        }
+        BaseType::validateLocalData();
     }
 
 protected:
@@ -779,8 +806,41 @@ class OutputCompositionalModule : public GenericOutputModule<GetPropType<TypeTag
         }
     }
 
+    /// Store the cell's saturation pressure, using zero for an unsuccessful solve.
+    /// Concurrent calls must use distinct cell indices; the failure count is atomic.
+    template<class FluidState>
+    void assignSaturationPressure_(const unsigned globalDofIdx, const FluidState& fluidState)
+    {
+        if (!this->compC_.saturationPressureAllocated()) {
+            return;
+        }
+
+        std::array<Scalar, numComponents> moleFractions;
+        for (int c = 0; c < numComponents; ++c) {
+            moleFractions[c] = getValue(fluidState.moleFraction(c));
+        }
+        const auto psat = CompositionalContainer<FluidSystem>::cellSaturationPressure(
+            getValue(fluidState.L()),
+            getValue(fluidState.pressure(oilPhaseIdx)),
+            moleFractions,
+            getValue(fluidState.temperature(oilPhaseIdx)),
+            this->eosType_);
+        if (!psat) {
+            // Failure includes supercritical mixtures; it does not establish
+            // whether a saturation pressure exists.
+#ifdef _OPENMP
+#pragma omp atomic
+#endif
+            ++this->numFailedSaturationPressures_;
+        }
+
+        this->compC_.assignSaturationPressure(globalDofIdx, psat.value_or(Scalar{0}));
+    }
+
     const Simulator& simulator_;
     CompositionalContainer<FluidSystem> compC_;
+    CompositionalConfig::EOSType eosType_;
+    std::size_t numFailedSaturationPressures_{};
     std::vector<typename Extractor::Entry> extractors_;
     typename BlockExtractor::ExecMap blockExtractors_;
 };