Enable New Saturation Function Consistency Checks

CMakeLists_files.cmake

@@ -221,8 +221,11 @@ if (HAVE_ECL_INPUT)
     opm/simulators/utils/satfunc/GasPhaseConsistencyChecks.cpp
     opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.cpp
     opm/simulators/utils/satfunc/PhaseCheckBase.cpp
+    opm/simulators/utils/satfunc/SatfuncConsistencyCheckManager.cpp
     opm/simulators/utils/satfunc/SatfuncConsistencyChecks.cpp
+    opm/simulators/utils/satfunc/ScaledSatfuncCheckPoint.cpp
     opm/simulators/utils/satfunc/ThreePointHorizontalConsistencyChecks.cpp
+    opm/simulators/utils/satfunc/UnscaledSatfuncCheckPoint.cpp
     opm/simulators/utils/satfunc/WaterPhaseConsistencyChecks.cpp
   )
 endif()
@@ -1039,8 +1042,12 @@ if (HAVE_ECL_INPUT)
     opm/simulators/utils/satfunc/GasPhaseConsistencyChecks.hpp
     opm/simulators/utils/satfunc/OilPhaseConsistencyChecks.hpp
     opm/simulators/utils/satfunc/PhaseCheckBase.hpp
+    opm/simulators/utils/satfunc/SatfuncCheckPointInterface.hpp
+    opm/simulators/utils/satfunc/SatfuncConsistencyCheckManager.hpp
     opm/simulators/utils/satfunc/SatfuncConsistencyChecks.hpp
+    opm/simulators/utils/satfunc/ScaledSatfuncCheckPoint.hpp
     opm/simulators/utils/satfunc/ThreePointHorizontalConsistencyChecks.hpp
+    opm/simulators/utils/satfunc/UnscaledSatfuncCheckPoint.hpp
     opm/simulators/utils/satfunc/WaterPhaseConsistencyChecks.hpp
   )
 endif()

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>
@@ -66,13 +66,14 @@
 
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/utils/ParallelSerialization.hpp>
-#include <opm/simulators/utils/satfunc/RelpermDiagnostics.hpp>
 
 #include <opm/utility/CopyablePtr.hpp>
 
 #include <algorithm>
+#include <cstddef>
 #include <functional>
 #include <set>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
@@ -223,38 +224,32 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
         , 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:
+        //
         // 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();
-        }
-
-        explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
+        //
+        // 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();
 
-        RelpermDiagnostics relpermDiagnostics;
-        relpermDiagnostics.diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
+        this->explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
     }
 
     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 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         // 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,42 +231,49 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         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);
     }
 
     /*!
      * \copydoc FvBaseProblem::finishInit
      */
     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 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
             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 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         // 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 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
 
         this->initFluidSystem_();
 
-
-
         if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
             FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
             this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
@@ -333,22 +348,25 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
                                       }
                                       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 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         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,26 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
             simulator.setTimeStepIndex(0);
         }
 
+        if (Parameters::Get<Parameters::CheckSatfuncConsistency>() &&
+            this->satfuncConsistencyRequirementsNotMet())
+        {
+            if (simulator.vanguard().grid().comm().rank() == 0) {
+                throw std::domain_error {
+                    "Saturation function end-points do not "
+                    "meet requisite consistency conditions"
+                };
+            }
+            else {
+                throw std::domain_error {""};
+            }
+        }
+
         // 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());
     }
 
     /*!
@@ -1443,6 +1476,80 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         this->updateRockCompTransMultVal_();
     }
 
+    bool satfuncConsistencyRequirementsNotMet() 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 'false' to tell caller that the
+            // consistency requirements ARE met, i.e., NotMet = false, in
+            // this case.
+            return false;
+        }
+
+        const auto numSamplePoints = std::size_t {5};
+
+        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 }); });
+
+            this->simulator().vanguard().grid().comm().barrier();
+        };
+
+        if (sfuncConsistencyChecks.anyFailedChecks()) {
+            if (isIoRank) {
+                OpmLog::warning("Benign Saturation Function "
+                                "End-point Consistency Failures");
+            }
+
+            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 NotMet.
+            return true;
+        }
+
+        // If we get here then there are no critical failures.  Report
+        // NotMet = false, i.e., that the consistency requirements ARE met.
+        return false;
+    }
+
     FlowThresholdPressure<TypeTag> thresholdPressures_;
 
     std::vector<InitialFluidState> initialFluidStates_;