From 08dc4e11c9e5e9329459d35af59b50ba83b787a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 10 Sep 2024 18:18:11 +0200 Subject: [PATCH] Enable New Saturation Function Consistency Checks The new parameter CheckSatfuncConsistency, command line option --check-satfunc-consistency, allows users to opt out of running the checks, e.g., because the curves are not of primary interest in the use case. 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 | 39 ++--- opm/simulators/flow/FlowProblemBlackoil.hpp | 157 +++++++++++++++--- opm/simulators/flow/FlowProblemParameters.cpp | 3 + opm/simulators/flow/FlowProblemParameters.hpp | 3 + tests/test_equil.cpp | 54 +++--- 5 files changed, 185 insertions(+), 71 deletions(-) diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index a9d359bb7c8..af694033926 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -37,8 +37,8 @@ #include #include -#include #include +#include #include #include @@ -66,13 +66,14 @@ #include #include -#include #include #include +#include #include #include +#include #include #include @@ -223,38 +224,32 @@ class FlowProblem : public GetPropType , pffDofData_(simulator.gridView(), this->elementMapper()) , tracerModel_(simulator) { - const auto& vanguard = simulator.vanguard(); - - enableDriftCompensation_ = Parameters::Get(); - - enableVtkOutput_ = Parameters::Get(); - + this->enableDriftCompensation_ = Parameters::Get(); + this->enableVtkOutput_ = Parameters::Get(); this->enableTuning_ = Parameters::Get(); + this->initialTimeStepSize_ = Parameters::Get>(); - this->maxTimeStepAfterWellEvent_ = Parameters::Get>() * 24 * 60 * 60; + this->maxTimeStepAfterWellEvent_ = unit::convert::from + (Parameters::Get>(), 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()) - { - this->numPressurePointsEquil_ = Parameters::Get(); - } else { - this->numPressurePointsEquil_ = simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP(); - } - - explicitRockCompaction_ = Parameters::Get(); + // + // 2. EQLDIMS item 2. Default value from + // opm-common/opm/input/eclipse/share/keywords/000_Eclipse100/E/EQLDIMS + this->numPressurePointsEquil_ = Parameters::IsSet() + ? Parameters::Get() + : simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP(); - RelpermDiagnostics relpermDiagnostics; - relpermDiagnostics.diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper()); + this->explicitRockCompaction_ = Parameters::Get(); } 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 diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp index 63fe0153a37..77f65e834b5 100644 --- a/opm/simulators/flow/FlowProblemBlackoil.hpp +++ b/opm/simulators/flow/FlowProblemBlackoil.hpp @@ -43,22 +43,28 @@ #include +#include + +#include #include #include #include +#include #include -#include +#include -#include #if HAVE_DAMARIS #include #endif #include +#include #include #include +#include #include +#include #include namespace Opm { @@ -209,6 +215,7 @@ class FlowProblemBlackoil : public FlowProblem // create the ECL writer eclWriter_ = std::make_unique(simulator); enableEclOutput_ = Parameters::Get(); + #if HAVE_DAMARIS // create Damaris writer damarisWriter_ = std::make_unique(simulator); @@ -224,23 +231,27 @@ class FlowProblemBlackoil : public FlowProblem 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 @@ class FlowProblemBlackoil : public FlowProblem */ 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 std::function 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 // disables gravity, else the standard value of the gravity constant at sea level // on earth is used this->gravity_ = 0.0; - if (Parameters::Get()) - this->gravity_[dim - 1] = 9.80665; - if (!eclState.getInitConfig().hasGravity()) - this->gravity_[dim - 1] = 0.0; + if (Parameters::Get() && + 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 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 } 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 if constexpr (getPropValue()) { 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 simulator.setTimeStepIndex(0); } + if (Parameters::Get() && + 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()); } /*! @@ -1450,6 +1483,80 @@ class FlowProblemBlackoil : public FlowProblem 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 { + 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::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 thresholdPressures_; std::vector initialFluidStates_; diff --git a/opm/simulators/flow/FlowProblemParameters.cpp b/opm/simulators/flow/FlowProblemParameters.cpp index 92c1e9a8828..3a04f844e7e 100644 --- a/opm/simulators/flow/FlowProblemParameters.cpp +++ b/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(); // Users will typically not need to modify this parameter.. + Parameters::Register + ("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.) diff --git a/opm/simulators/flow/FlowProblemParameters.hpp b/opm/simulators/flow/FlowProblemParameters.hpp index 3ac0bf6f39b..3442a9366c1 100644 --- a/opm/simulators/flow/FlowProblemParameters.hpp +++ b/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 = true; }; + // Parameterize equilibration accuracy struct NumPressurePointsEquil { static constexpr int value = ParserKeywords::EQLDIMS::DEPTH_NODES_P::defaultValue; }; diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 75c2284c93e..cfcc22c60cc 100644 --- a/tests/test_equil.cpp +++ b/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 + +#include +#if (BOOST_VERSION / 100000 == 1) && ((BOOST_VERSION / 100) % 1000 < 71) +#include +#else +#include +#endif + #include #include #include @@ -58,42 +68,37 @@ #include #include -#include -#include -#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 71 -#include -#else -#include -#endif - - namespace Opm::Properties { namespace TTag { - struct TestEquilTypeTag { using InheritsFrom = std::tuple; }; + struct TestEquilVapwatTypeTag { using InheritsFrom = std::tuple; }; -} + +} // namespace TTag template struct WellModel { using type = BlackoilWellModel; }; + template struct EnableVapwat { static constexpr bool value = true; }; + template struct WellModel { using type = BlackoilWellModel; }; + template struct EnableVapwat { static constexpr bool value = true; @@ -101,18 +106,20 @@ struct EnableVapwat { } // namespace Opm::Properties +namespace { + template std::unique_ptr> initSimulator(const char *filename) { using Simulator = Opm::GetPropType; - 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_(/*argc=*/sizeof(argv)/sizeof(argv[0]), argv, /*registerParams=*/false); @@ -123,7 +130,8 @@ initSimulator(const char *filename) } template -static std::vector> cellVerticalExtent(const GridView& gridView) +std::vector> +cellVerticalExtent(const GridView& gridView) { using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); @@ -142,7 +150,7 @@ static std::vector> cellVerticalExtent(const GridView& } template -static void initDefaultFluidSystem() +void initDefaultFluidSystem() { using FluidSystem = Opm::GetPropType; @@ -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 @@ -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;