Merge pull request #5596 from bska/enable-new-satfunc-consistency-checks

Enable New Saturation Function Consistency Checks

opm/simulators/flow/FlowProblem.hpp

@@ -37,8 +37,8 @@
 #include <opm/common/utility/TimeService.hpp>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/input/eclipse/Parser/ParserKeywords/E.hpp>
 #include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/input/eclipse/Units/Units.hpp>
 
 #include <opm/material/common/ConditionalStorage.hpp>
 #include <opm/material/common/Valgrind.hpp>
@@ -72,8 +72,10 @@
 #include <opm/utility/CopyablePtr.hpp>
 
 #include <algorithm>
+#include <cstddef>
 #include <functional>
 #include <set>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
@@ -224,38 +226,41 @@ public:
         , pffDofData_(simulator.gridView(), this->elementMapper())
         , tracerModel_(simulator)
     {
-        const auto& vanguard = simulator.vanguard();
-
-        enableDriftCompensation_ = Parameters::Get<Parameters::EnableDriftCompensation>();
-
-        enableVtkOutput_ = Parameters::Get<Parameters::EnableVtkOutput>();
-
+        this->enableDriftCompensation_ = Parameters::Get<Parameters::EnableDriftCompensation>();
+        this->enableVtkOutput_ = Parameters::Get<Parameters::EnableVtkOutput>();
         this->enableTuning_ = Parameters::Get<Parameters::EnableTuning>();
+
         this->initialTimeStepSize_ = Parameters::Get<Parameters::InitialTimeStepSize<Scalar>>();
-        this->maxTimeStepAfterWellEvent_ = Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60;
+        this->maxTimeStepAfterWellEvent_ = unit::convert::from
+            (Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>(), unit::day);
 
-        // The value N for this parameter is defined in the following order of presedence:
+        // The value N for this parameter is defined in the following order of precedence:
+        //
         // 1. Command line value (--num-pressure-points-equil=N)
-        // 2. EQLDIMS item 2
-        // Default value is defined in opm-common/src/opm/input/eclipse/share/keywords/000_Eclipse100/E/EQLDIMS
-        if (Parameters::IsSet<Parameters::NumPressurePointsEquil>())
-        {
-            this->numPressurePointsEquil_ = Parameters::Get<Parameters::NumPressurePointsEquil>();
-        } else {
-            this->numPressurePointsEquil_ = simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP();
+        //
+        // 2. EQLDIMS item 2.  Default value from
+        //    opm-common/opm/input/eclipse/share/keywords/000_Eclipse100/E/EQLDIMS
+
+        this->numPressurePointsEquil_ = Parameters::IsSet<Parameters::NumPressurePointsEquil>()
+            ? Parameters::Get<Parameters::NumPressurePointsEquil>()
+            : simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP();
+
+        this->explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
+
+        if (! Parameters::Get<Parameters::CheckSatfuncConsistency>()) {
+            // User did not enable the "new" saturation function consistency
+            // check module.  Run the original checker instead.  This is a
+            // temporary measure.
+            RelpermDiagnostics relpermDiagnostics{};
+            relpermDiagnostics.diagnosis(simulator.vanguard().eclState(),
+                                         simulator.vanguard().cartesianIndexMapper());
         }
-
-        explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
-
-
-        RelpermDiagnostics relpermDiagnostics;
-        relpermDiagnostics.diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
     }
 
     virtual ~FlowProblem() = default;
 
     void prefetch(const Element& elem) const
-    { pffDofData_.prefetch(elem); }
+    { this->pffDofData_.prefetch(elem); }
 
     /*!
      * \brief This method restores the complete state of the problem and its sub-objects

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -43,22 +43,28 @@
 
 #include <opm/output/eclipse/EclipseIO.hpp>
 
+#include <opm/input/eclipse/Units/Units.hpp>
+
+#include <opm/simulators/flow/ActionHandler.hpp>
 #include <opm/simulators/flow/FlowProblem.hpp>
 #include <opm/simulators/flow/FlowProblemBlackoilProperties.hpp>
 #include <opm/simulators/flow/FlowThresholdPressure.hpp>
+#include <opm/simulators/flow/MixingRateControls.hpp>
 #include <opm/simulators/flow/VtkTracerModule.hpp>
 
-#include <opm/simulators/flow/MixingRateControls.hpp>
+#include <opm/simulators/utils/satfunc/SatfuncConsistencyCheckManager.hpp>
 
-#include <opm/simulators/flow/ActionHandler.hpp>
 #if HAVE_DAMARIS
 #include <opm/simulators/flow/DamarisWriter.hpp>
 #endif
 
 #include <algorithm>
+#include <cstddef>
 #include <functional>
 #include <set>
+#include <stdexcept>
 #include <string>
+#include <string_view>
 #include <vector>
 
 namespace Opm {
@@ -209,6 +215,7 @@ public:
         // create the ECL writer
         eclWriter_ = std::make_unique<EclWriterType>(simulator);
         enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
+
 #if HAVE_DAMARIS
         // create Damaris writer
         damarisWriter_ = std::make_unique<DamarisWriterType>(simulator);
@@ -224,23 +231,27 @@ public:
         FlowProblemType::beginEpisode();
 
         auto& simulator = this->simulator();
-        int episodeIdx = simulator.episodeIndex();
+
+        const int episodeIdx = simulator.episodeIndex();
         const auto& schedule = simulator.vanguard().schedule();
 
         // Evaluate UDQ assign statements to make sure the settings are
         // available as UDA controls for the current report step.
-        actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
+        this->actionHandler_
+            .evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
 
         if (episodeIdx >= 0) {
             const auto& oilVap = schedule[episodeIdx].oilvap();
             if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
                 FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2());
-            } else {
+            }
+            else {
                 FluidSystem::setVapPars(0.0, 0.0);
             }
         }
 
-        ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), simulator.vanguard().schedule(), episodeIdx, moduleParams_.convectiveMixingModuleParam);
+        ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), schedule, episodeIdx,
+                                             this->moduleParams_.convectiveMixingModuleParam);
     }
 
     /*!
@@ -248,18 +259,21 @@ public:
      */
     void finishInit()
     {
-        // TODO: there should be room to remove duplication for this function,
-        // but there is relatively complicated logic in the function calls in this function
-        // some refactoring is needed for this function
+        // TODO: there should be room to remove duplication for this
+        // function, but there is relatively complicated logic in the
+        // function calls here.  Some refactoring is needed.
         FlowProblemType::finishInit();
+
         auto& simulator = this->simulator();
 
         auto finishTransmissibilities = [updated = false, this]() mutable
         {
             if (updated) { return; }
+
             this->transmissibilities_.finishInit([&vg = this->simulator().vanguard()](const unsigned int it) {
                 return vg.gridIdxToEquilGridIdx(it);
             });
+
             updated = true;
         };
 
@@ -279,7 +293,8 @@ public:
             std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
                 return simulator.vanguard().gridEquilIdxToGridIdx(i);
             };
-            eclWriter_->extractOutputTransAndNNC(equilGridToGrid);
+
+            this->eclWriter_->extractOutputTransAndNNC(equilGridToGrid);
         }
         simulator.vanguard().releaseGlobalTransmissibilities();
 
@@ -301,10 +316,12 @@ public:
         // disables gravity, else the standard value of the gravity constant at sea level
         // on earth is used
         this->gravity_ = 0.0;
-        if (Parameters::Get<Parameters::EnableGravity>())
-            this->gravity_[dim - 1] = 9.80665;
-        if (!eclState.getInitConfig().hasGravity())
-            this->gravity_[dim - 1] = 0.0;
+        if (Parameters::Get<Parameters::EnableGravity>() &&
+            eclState.getInitConfig().hasGravity())
+        {
+            // unit::gravity is 9.80665 m^2/s--i.e., standard measure at Tellus equator.
+            this->gravity_[dim - 1] = unit::gravity;
+        }
 
         if (this->enableTuning_) {
             // if support for the TUNING keyword is enabled, we get the initial time
@@ -316,8 +333,6 @@ public:
 
         this->initFluidSystem_();
 
-
-
         if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
             FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
             this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
@@ -333,22 +348,25 @@ public:
                                       }
                                       return coords;
                                   });
+
         this->readMaterialParameters_();
         this->readThermalParameters_();
 
         // write the static output files (EGRID, INIT)
         if (enableEclOutput_) {
-            eclWriter_->writeInit();
+            this->eclWriter_->writeInit();
         }
 
         finishTransmissibilities();
 
         const auto& initconfig = eclState.getInitConfig();
         this->tracerModel_.init(initconfig.restartRequested());
-        if (initconfig.restartRequested())
-            readEclRestartSolution_();
-        else
-            readInitialCondition_();
+        if (initconfig.restartRequested()) {
+            this->readEclRestartSolution_();
+        }
+        else {
+            this->readInitialCondition_();
+        }
 
         this->tracerModel_.prepareTracerBatches();
 
@@ -357,14 +375,14 @@ public:
         if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
             const auto& vanguard = this->simulator().vanguard();
             const auto& gridView = vanguard.gridView();
-            int numElements = gridView.size(/*codim=*/0);
+            const int numElements = gridView.size(/*codim=*/0);
             this->polymer_.maxAdsorption.resize(numElements, 0.0);
         }
 
         this->readBoundaryConditions_();
 
         // compute and set eq weights based on initial b values
-        computeAndSetEqWeights_();
+        this->computeAndSetEqWeights_();
 
         if (this->enableDriftCompensation_) {
             this->drift_.resize(this->model().numGridDof());
@@ -385,11 +403,30 @@ public:
             simulator.setTimeStepIndex(0);
         }
 
+        if (Parameters::Get<Parameters::CheckSatfuncConsistency>() &&
+            ! this->satfuncConsistencyRequirementsMet())
+        {
+            // User requested that saturation functions be checked for
+            // consistency and essential/critical requirements are not met.
+            // Abort simulation run.
+            //
+            // Note: We need synchronisation here lest ranks other than the
+            // I/O rank throw exceptions too early thereby risking an
+            // incomplete failure report being shown to the user.
+            this->simulator().vanguard().grid().comm().barrier();
+
+            throw std::domain_error {
+                "Saturation function end-points do not "
+                "meet requisite consistency conditions"
+            };
+        }
+
         // TODO: move to the end for later refactoring of the function finishInit()
-	// deal with DRSDT
+        //
+        // deal with DRSDT
         this->mixControls_.init(this->model().numGridDof(),
                                 this->episodeIndex(),
-                                eclState.runspec().tabdims().getNumPVTTables());   
+                                eclState.runspec().tabdims().getNumPVTTables());
     }
 
     /*!
@@ -1450,6 +1487,78 @@ protected:
         this->updateRockCompTransMultVal_();
     }
 
+    bool satfuncConsistencyRequirementsMet() const
+    {
+        if (const auto nph = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
+            + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
+            + FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
+            nph < 2)
+        {
+            // Single phase runs don't need saturation functions and there's
+            // nothing to do here.  Return 'true' to tell caller that the
+            // consistency requirements are Met.
+            return true;
+        }
+
+        const auto numSamplePoints = static_cast<std::size_t>
+            (Parameters::Get<Parameters::NumSatfuncConsistencySamplePoints>());
+
+        auto sfuncConsistencyChecks =
+            Satfunc::PhaseChecks::SatfuncConsistencyCheckManager<Scalar> {
+            numSamplePoints, this->simulator().vanguard().eclState(),
+            [&cmap = this->simulator().vanguard().cartesianIndexMapper()](const int elemIdx)
+            { return cmap.cartesianIndex(elemIdx); }
+        };
+
+        const auto ioRank = 0;
+        const auto isIoRank = this->simulator().vanguard()
+            .grid().comm().rank() == ioRank;
+
+        sfuncConsistencyChecks.collectFailuresTo(ioRank)
+            .run(this->simulator().vanguard().grid().leafGridView(),
+                 [&vg   = this->simulator().vanguard(),
+                  &emap = this->simulator().model().elementMapper()]
+                 (const auto& elem)
+                 { return vg.gridIdxToEquilGridIdx(emap.index(elem)); });
+
+        using ViolationLevel = typename Satfunc::PhaseChecks::
+            SatfuncConsistencyCheckManager<Scalar>::ViolationLevel;
+
+        auto reportFailures = [&sfuncConsistencyChecks, this]
+            (const ViolationLevel level)
+        {
+            sfuncConsistencyChecks.reportFailures
+                (level, [](std::string_view record)
+                { OpmLog::info(std::string { record }); });
+        };
+
+        if (sfuncConsistencyChecks.anyFailedStandardChecks()) {
+            if (isIoRank) {
+                OpmLog::warning("Saturation Function "
+                                "End-point Consistency Problems");
+
+                reportFailures(ViolationLevel::Standard);
+            }
+        }
+
+        if (sfuncConsistencyChecks.anyFailedCriticalChecks()) {
+            if (isIoRank) {
+                OpmLog::error("Saturation Function "
+                              "End-point Consistency Failures");
+
+                reportFailures(ViolationLevel::Critical);
+            }
+
+            // There are "critical" check failures.  Report that consistency
+            // requirements are not Met.
+            return false;
+        }
+
+        // If we get here then there are no critical failures.  Report
+        // Met = true, i.e., that the consistency requirements ARE met.
+        return true;
+    }
+
     FlowThresholdPressure<TypeTag> thresholdPressures_;
 
     std::vector<InitialFluidState> initialFluidStates_;

opm/simulators/flow/FlowProblemParameters.cpp

@@ -66,6 +66,12 @@ void registerFlowProblemParameters()
         ("Use pressure from end of the last time step when evaluating rock compaction");
     Parameters::Hide<Parameters::ExplicitRockCompaction>(); // Users will typically not need to modify this parameter..
 
+    Parameters::Register<Parameters::CheckSatfuncConsistency>
+        ("Whether or not to check saturation function consistency requirements");
+
+    Parameters::Register<Parameters::NumSatfuncConsistencySamplePoints>
+        ("Maximum number of reported failures for each individual saturation function consistency check");
+
     // By default, stop it after the universe will probably have stopped
     // to exist. (the ECL problem will finish the simulation explicitly
     // after it simulated the last episode specified in the deck.)

opm/simulators/flow/FlowProblemParameters.hpp

@@ -39,6 +39,13 @@ struct EnableDriftCompensation { static constexpr bool value = false; };
 // implicit or explicit pressure in rock compaction
 struct ExplicitRockCompaction { static constexpr bool value = false; };
 
+// Whether or not to check saturation function consistency requirements.
+struct CheckSatfuncConsistency { static constexpr bool value = false; };
+
+// Maximum number of reported failures for each saturation function
+// consistency check.
+struct NumSatfuncConsistencySamplePoints { static constexpr int value = 5; };
+
 // Parameterize equilibration accuracy
 struct NumPressurePointsEquil
 { static constexpr int value = ParserKeywords::EQLDIMS::DEPTH_NODES_P::defaultValue; };

opm/simulators/utils/satfunc/SatfuncConsistencyCheckManager.hpp

@@ -152,8 +152,8 @@ namespace Opm::Satfunc::PhaseChecks {
         /// Reports only those conditions/checks for which there is at least
         /// one violation.
         ///
-        /// In a parallel run it is only safe to call this function on the
-        /// root process defined by collectFailuresTo().
+        /// In a parallel run this function does nothing on ranks other than
+        /// the root process defined by collectFailuresTo().
         ///
         /// \param[in] level Report's severity level.
         ///

tests/test_equil.cpp

@@ -20,10 +20,20 @@
   module for the precise wording of the license and the list of
   copyright holders.
 */
+
 #include "config.h"
 
 #define BOOST_TEST_MODULE Equil
 
+#include <boost/test/unit_test.hpp>
+
+#include <boost/version.hpp>
+#if (BOOST_VERSION / 100000 == 1) && ((BOOST_VERSION / 100) % 1000 < 71)
+#include <boost/test/floating_point_comparison.hpp>
+#else
+#include <boost/test/tools/floating_point_comparison.hpp>
+#endif
+
 #include <opm/grid/UnstructuredGrid.h>
 #include <opm/grid/GridManager.hpp>
 #include <opm/grid/cpgrid/GridHelpers.hpp>
@@ -58,42 +68,37 @@
 #include <string>
 #include <vector>
 
-#include <boost/test/unit_test.hpp>
-#include <boost/version.hpp>
-#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 71
-#include <boost/test/floating_point_comparison.hpp>
-#else
-#include <boost/test/tools/floating_point_comparison.hpp>
-#endif
-
-
 namespace Opm::Properties {
 
 namespace TTag {
 
-
 struct TestEquilTypeTag {
     using InheritsFrom = std::tuple<FlowBaseProblemBlackoil,
                                     BlackOilModel>;
 };
+
 struct TestEquilVapwatTypeTag {
     using InheritsFrom = std::tuple<FlowBaseProblemBlackoil,
                                     BlackOilModel>;
 };
-}
+
+} // namespace TTag
 
 template<class TypeTag>
 struct WellModel<TypeTag, TTag::TestEquilTypeTag> {
     using type = BlackoilWellModel<TypeTag>;
 };
+
 template<class TypeTag>
 struct EnableVapwat<TypeTag, TTag::TestEquilTypeTag> {
     static constexpr bool value = true;
 };
+
 template<class TypeTag>
 struct WellModel<TypeTag, TTag::TestEquilVapwatTypeTag> {
     using type = BlackoilWellModel<TypeTag>;
 };
+
 template<class TypeTag>
 struct EnableVapwat<TypeTag, TTag::TestEquilVapwatTypeTag> {
     static constexpr bool value = true;
@@ -101,18 +106,20 @@ struct EnableVapwat<TypeTag, TTag::TestEquilVapwatTypeTag> {
 
 } // namespace Opm::Properties
 
+namespace {
+
 template <class TypeTag>
 std::unique_ptr<Opm::GetPropType<TypeTag, Opm::Properties::Simulator>>
 initSimulator(const char *filename)
 {
     using Simulator = Opm::GetPropType<TypeTag, Opm::Properties::Simulator>;
 
-    std::string filenameArg = "--ecl-deck-file-name=";
-    filenameArg += filename;
+    const auto filenameArg = std::string {"--ecl-deck-file-name="} + filename;
 
     const char* argv[] = {
         "test_equil",
-        filenameArg.c_str()
+        filenameArg.c_str(),
+        "--check-satfunc-consistency=false",
     };
 
     Opm::setupParameters_<TypeTag>(/*argc=*/sizeof(argv)/sizeof(argv[0]), argv, /*registerParams=*/false);
@@ -123,7 +130,8 @@ initSimulator(const char *filename)
 }
 
 template <class GridView>
-static std::vector<std::pair<double,double>> cellVerticalExtent(const GridView& gridView)
+std::vector<std::pair<double,double>>
+cellVerticalExtent(const GridView& gridView)
 {
     using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
     ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
@@ -142,7 +150,7 @@ static std::vector<std::pair<double,double>> cellVerticalExtent(const GridView&
 }
 
 template <class TypeTag>
-static void initDefaultFluidSystem()
+void initDefaultFluidSystem()
 {
     using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
 
@@ -207,11 +215,12 @@ static void initDefaultFluidSystem()
     FluidSystem::initEnd();
 }
 
-static Opm::EquilRecord mkEquilRecord( double datd, double datp,
-                                       double zwoc, double pcow_woc,
-                                       double zgoc, double pcgo_goc )
+Opm::EquilRecord
+mkEquilRecord(const double datd, const double datp,
+              const double zwoc, const double pcow_woc,
+              const double zgoc, const double pcgo_goc)
 {
-    return Opm::EquilRecord( datd, datp, zwoc, pcow_woc, zgoc, pcgo_goc, true, true, 0, true);
+    return { datd, datp, zwoc, pcow_woc, zgoc, pcgo_goc, true, true, 0, true};
 }
 
 template <typename Simulator>
@@ -220,8 +229,6 @@ double centerDepth(const Simulator& sim, const std::size_t cell)
     return Opm::UgGridHelpers::cellCenterDepth(sim.vanguard().grid(), cell);
 }
 
-namespace {
-
 struct EquilFixture {
     EquilFixture() {
         int argc = boost::unit_test::framework::master_test_suite().argc;
@@ -253,7 +260,7 @@ struct EquilFixture {
                                                                         CartesianIndexMapper>;
 };
 
-}
+} // Anonymous namespace
 
 BOOST_GLOBAL_FIXTURE(EquilFixture);
 
@@ -845,7 +852,6 @@ BOOST_AUTO_TEST_CASE(DeckWithCO2STORE)
     }
 }
 
-
 BOOST_AUTO_TEST_CASE(DeckWithWetGas)
 {
     using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;