Enable New Saturation Function Consistency Checks

Currently opt-in by the new parameter CheckSatfuncConsistency
(command line option --check-satfunc-consistency), this commit hooks
up the SatfuncConsistencyCheckManager<> pass to the FlowProblem.  I
hope to make the parameter be "opt-out" at some point, but for now
the new pass breaks too many regression tests.

opm/simulators/flow/FlowProblem.hpp

@@ -71,8 +71,10 @@
 #include <opm/utility/CopyablePtr.hpp>
 
 #include <algorithm>
+#include <cstddef>
 #include <functional>
 #include <set>
+#include <stdexcept>
 #include <string>
 #include <vector>
 

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -43,22 +43,26 @@
 
 #include <opm/output/eclipse/EclipseIO.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 +213,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 +229,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 +291,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();
 
@@ -316,8 +329,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,12 +344,13 @@ 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();
@@ -357,14 +369,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 +397,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 +1470,64 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
  this->updateRockCompTransMultVal_();
  }
 
+ bool satfuncConsistencyRequirementsNotMet() const
+ {
+ 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;
+
+ if (sfuncConsistencyChecks.anyFailedChecks()) {
+ if (isIoRank) {
+ OpmLog::warning("Benign Saturation Function "
+ "End-point Consistency Failures");
+ }
+
+ sfuncConsistencyChecks.reportFailures
+ (ViolationLevel::Standard,
+ [](std::string_view record)
+ { OpmLog::info(std::string { record }); });
+
+ this->simulator().vanguard().grid().comm().barrier();
+ }
+
+ if (sfuncConsistencyChecks.anyFailedCriticalChecks()) {
+ if (isIoRank) {
+ OpmLog::error("Saturation Function "
+ "End-point Consistency Failures");
+ }
+
+ sfuncConsistencyChecks.reportFailures
+ (ViolationLevel::Critical,
+ [](std::string_view record)
+ { OpmLog::info(std::string { record }); });
+
+ this->simulator().vanguard().grid().comm().barrier();
+
+ return true;
+ }
+
+ return false;
+ }
+
  FlowThresholdPressure<TypeTag> thresholdPressures_;
 
  std::vector<InitialFluidState> initialFluidStates_;

opm/simulators/flow/FlowProblemParameters.cpp

@@ -66,6 +66,9 @@ 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");
+
  // 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,9 @@ 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; };
+
 // Parameterize equilibration accuracy
 struct NumPressurePointsEquil
 { static constexpr int value = ParserKeywords::EQLDIMS::DEPTH_NODES_P::defaultValue; };