Skip to content

Commit

Permalink
Enable New Saturation Function Consistency Checks
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@bska
bska committed Sep 20, 2024
1 parent a7157e0 commit d7d435b
Show file tree
Hide file tree
Showing 5 changed files with 185 additions and 71 deletions.
39 changes: 17 additions & 22 deletions opm/simulators/flow/FlowProblem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -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>

Expand Down Expand Up @@ -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
Expand Down
157 changes: 132 additions & 25 deletions opm/simulators/flow/FlowProblemBlackoil.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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);
Expand All @@ -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;
};

Expand All @@ -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();

Expand All @@ -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
Expand All @@ -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);
Expand All @@ -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();

Expand All @@ -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());
Expand All @@ -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());
}

/*!
Expand Down Expand Up @@ -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_;
Expand Down
3 changes: 3 additions & 0 deletions opm/simulators/flow/FlowProblemParameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.)
Expand Down
3 changes: 3 additions & 0 deletions opm/simulators/flow/FlowProblemParameters.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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; };
Expand Down
Loading

0 comments on commit d7d435b

Please sign in to comment.