Enable New Saturation Function Consistency Checks

The new parameter CheckSatfuncConsistency, command line option
--check-satfunc-consistency, allows users to opt into running the
checks.  The option currently defaults to 'false' to reflect the
somewhat experimental nature of the new facility.

The new parameter --num-satfunc-consistency-sample-points allows the
user to select the maximum number of reported failures for each
individual consistency check.  By default, the simulator will report
at most five failures for each check.

We check the unscaled curves if the run does not activate the
end-point scaling option and the scaled curves otherwise.  At
present we're limited to reversible and non-directional saturation
functions for the drainage process, but those restrictions will be
lifted in due time.

opm/simulators/flow/FlowProblem.hpp

@@ -247,9 +247,14 @@ public:
 
         this->explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
 
-        RelpermDiagnostics relpermDiagnostics{};
-        relpermDiagnostics.diagnosis(simulator.vanguard().eclState(),
-                                     simulator.vanguard().cartesianIndexMapper());
+        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());
+        }
     }
 
     virtual ~FlowProblem() = default;

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -52,6 +52,7 @@
 #include <opm/simulators/flow/MixingRateControls.hpp>
 #include <opm/simulators/flow/VtkTracerModule.hpp>
 
+#include <opm/simulators/utils/satfunc/SatfuncConsistencyCheckManager.hpp>
 
 #if HAVE_DAMARIS
 #include <opm/simulators/flow/DamarisWriter.hpp>
@@ -63,6 +64,7 @@
 #include <set>
 #include <stdexcept>
 #include <string>
+#include <string_view>
 #include <vector>
 
 namespace Opm {
@@ -401,6 +403,24 @@ 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
@@ -1467,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

@@ -119,6 +119,7 @@ initSimulator(const char *filename)
     const char* argv[] = {
         "test_equil",
         filenameArg.c_str(),
+        "--check-satfunc-consistency=false",
     };
 
     Opm::setupParameters_<TypeTag>(/*argc=*/sizeof(argv)/sizeof(argv[0]), argv, /*registerParams=*/false);