Loading...
Searching...
No Matches
FlowProblem.hpp
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
Copyright 2023 INRIA
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef OPM_FLOW_PROBLEM_HPP
#define OPM_FLOW_PROBLEM_HPP
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <opm/common/utility/TimeService.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.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>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/material/thermal/EclThermalLawManager.hpp>
#include <opm/models/common/directionalmobility.hh>
#include <opm/models/utils/pffgridvector.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/simulators/flow/CpGridVanguard.hpp>
#include <opm/simulators/flow/DummyGradientCalculator.hpp>
#include <opm/simulators/flow/EclWriter.hpp>
#include <opm/simulators/flow/EquilInitializer.hpp>
#include <opm/simulators/flow/FlowGenericProblem.hpp>
// TODO: maybe we can name it FlowProblemProperties.hpp
#include <opm/simulators/flow/FlowBaseProblemProperties.hpp>
#include <opm/simulators/flow/FlowUtils.hpp>
#include <opm/simulators/flow/TracerModel.hpp>
#include <opm/simulators/flow/Transmissibility.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#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>
namespace Opm {
template <class TypeTag>
class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
, public FlowGenericProblem<GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::FluidSystem>>
{
protected:
using BaseType = FlowGenericProblem<GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::FluidSystem>>;
using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
using Implementation = GetPropType<TypeTag, Properties::Problem>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Stencil = GetPropType<TypeTag, Properties::Stencil>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
using EqVector = GetPropType<TypeTag, Properties::EqVector>;
using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
// Grid and world dimension
enum { dim = GridView::dimension };
enum { dimWorld = GridView::dimensionworld };
// copy some indices for convenience
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
enum { enableMICP = Indices::enableMICP };
enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
// TODO: later, gasCompIdx, oilCompIdx and waterCompIdx should go to the FlowProblemBlackoil in the future
// we do not want them in the compositional setting
enum { gasCompIdx = FluidSystem::gasCompIdx };
enum { oilCompIdx = FluidSystem::oilCompIdx };
enum { waterCompIdx = FluidSystem::waterCompIdx };
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Element = typename GridView::template Codim<0>::Entity;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using EclMaterialLawManager = typename GetProp<TypeTag, Properties::MaterialLaw>::EclMaterialLawManager;
using EclThermalLawManager = typename GetProp<TypeTag, Properties::SolidEnergyLaw>::EclThermalLawManager;
using MaterialLawParams = typename EclMaterialLawManager::MaterialLawParams;
using SolidEnergyLawParams = typename EclThermalLawManager::SolidEnergyLawParams;
using ThermalConductionLawParams = typename EclThermalLawManager::ThermalConductionLawParams;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using WellModel = GetPropType<TypeTag, Properties::WellModel>;
using AquiferModel = GetPropType<TypeTag, Properties::AquiferModel>;
using Toolbox = MathToolbox<Evaluation>;
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using TracerModel = GetPropType<TypeTag, Properties::TracerModel>;
using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag>>;
public:
using BaseType::briefDescription;
using BaseType::helpPreamble;
using BaseType::shouldWriteOutput;
using BaseType::shouldWriteRestartFile;
using BaseType::rockCompressibility;
using BaseType::porosity;
static void registerParameters()
{
ParentType::registerParameters();
registerFlowProblemParameters<Scalar>();
}
static int handlePositionalParameter(std::function<void(const std::string&,
const std::string&)> addKey,
std::set<std::string>& seenParams,
std::string& errorMsg,
int,
const char** argv,
int paramIdx,
int)
{
return detail::eclPositionalParameter(addKey,
seenParams,
errorMsg,
argv,
paramIdx);
}
explicit FlowProblem(Simulator& simulator)
: ParentType(simulator)
, BaseType(simulator.vanguard().eclState(),
simulator.vanguard().schedule(),
simulator.vanguard().gridView())
, transmissibilities_(simulator.vanguard().eclState(),
simulator.vanguard().gridView(),
simulator.vanguard().cartesianIndexMapper(),
simulator.vanguard().grid(),
simulator.vanguard().cellCentroids(),
enableEnergy,
enableDiffusion,
enableDispersion)
, wellModel_(simulator)
, aquiferModel_(simulator)
, pffDofData_(simulator.gridView(), this->elementMapper())
, tracerModel_(simulator)
{
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().levelCartesianIndexMapper());
}
}
virtual ~FlowProblem() = default;
void prefetch(const Element& elem) const
{ this->pffDofData_.prefetch(elem); }
template <class Restarter>
void deserialize(Restarter& res)
{
// reload the current episode/report step from the deck
this->beginEpisode();
// deserialize the wells
wellModel_.deserialize(res);
// deserialize the aquifer
aquiferModel_.deserialize(res);
}
template <class Restarter>
void serialize(Restarter& res)
{
wellModel_.serialize(res);
aquiferModel_.serialize(res);
}
int episodeIndex() const
{
return std::max(this->simulator().episodeIndex(), 0);
}
virtual void beginEpisode()
{
OPM_TIMEBLOCK(beginEpisode);
// Proceed to the next report step
auto& simulator = this->simulator();
int episodeIdx = simulator.episodeIndex();
auto& eclState = simulator.vanguard().eclState();
const auto& schedule = simulator.vanguard().schedule();
const auto& events = schedule[episodeIdx].events();
if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
// bring the contents of the keywords to the current state of the SCHEDULE
// section.
//
// TODO (?): make grid topology changes possible (depending on what exactly
// has changed, the grid may need be re-created which has some serious
// implications on e.g., the solution of the simulation.)
const auto& miniDeck = schedule[episodeIdx].geo_keywords();
const auto& cc = simulator.vanguard().grid().comm();
eclState.apply_schedule_keywords( miniDeck );
eclBroadcast(cc, eclState.getTransMult() );
// Re-ordering in case of ALUGrid
std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridEquilIdxToGridIdx(i);
};
// re-compute all quantities which may possibly be affected.
using TransUpdateQuantities = typename Vanguard::TransmissibilityType::TransUpdateQuantities;
transmissibilities_.update(true, TransUpdateQuantities::All, equilGridToGrid);
this->referencePorosity_[1] = this->referencePorosity_[0];
updateReferencePorosity_();
updatePffDofData_();
this->model().linearizer().updateDiscretizationParameters();
}
bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
// set up the wells for the next episode.
wellModel_.beginEpisode();
// set up the aquifers for the next episode.
aquiferModel_.beginEpisode();
// set the size of the initial time step of the episode
Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
// negative value of initialTimeStepSize_ indicates no active limit from TSINIT or NEXTSTEP
if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
// allow the size of the initial time step to be set via an external parameter
// if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
dt = std::min(dt, this->initialTimeStepSize_);
simulator.setTimeStepSize(dt);
}
virtual void beginTimeStep()
{
OPM_TIMEBLOCK(beginTimeStep);
const int episodeIdx = this->episodeIndex();
const int timeStepSize = this->simulator().timeStepSize();
this->beginTimeStep_(enableExperiments,
episodeIdx,
this->simulator().timeStepIndex(),
this->simulator().startTime(),
this->simulator().time(),
timeStepSize,
this->simulator().endTime());
// update maximum water saturation and minimum pressure
// used when ROCKCOMP is activated
// Do not update max RS first step after a restart
this->updateExplicitQuantities_(episodeIdx, timeStepSize, first_step_ && (episodeIdx > 0));
first_step_ = false;
if (nonTrivialBoundaryConditions()) {
this->model().linearizer().updateBoundaryConditionData();
}
wellModel_.beginTimeStep();
aquiferModel_.beginTimeStep();
tracerModel_.beginTimeStep();
}
void beginIteration()
{
OPM_TIMEBLOCK(beginIteration);
wellModel_.beginIteration();
aquiferModel_.beginIteration();
}
void endIteration()
{
OPM_TIMEBLOCK(endIteration);
wellModel_.endIteration();
aquiferModel_.endIteration();
}
virtual void endTimeStep()
{
OPM_TIMEBLOCK(endTimeStep);
#ifndef NDEBUG
if constexpr (getPropValue<TypeTag, Properties::EnableDebuggingChecks>()) {
// in debug mode, we don't care about performance, so we check
// if the model does the right thing (i.e., the mass change
// inside the whole reservoir must be equivalent to the fluxes
// over the grid's boundaries plus the source rates specified by
// the problem).
const int rank = this->simulator().gridView().comm().rank();
if (rank == 0) {
std::cout << "checking conservativeness of solution\n";
}
this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
if (rank == 0) {
std::cout << "solution is sufficiently conservative\n";
}
}
#endif // NDEBUG
auto& simulator = this->simulator();
simulator.setTimeStepIndex(simulator.timeStepIndex()+1);
this->wellModel_.endTimeStep();
this->aquiferModel_.endTimeStep();
this->tracerModel_.endTimeStep();
// Compute flux for output
this->model().linearizer().updateFlowsInfo();
if (this->enableDriftCompensation_) {
OPM_TIMEBLOCK(driftCompansation);
const auto& residual = this->model().linearizer().residual();
for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(globalDofIdx);
this->drift_[sfcdofIdx] = residual[sfcdofIdx] * simulator.timeStepSize();
if constexpr (getPropValue<TypeTag, Properties::UseVolumetricResidual>()) {
this->drift_[sfcdofIdx] *= this->model().dofTotalVolume(sfcdofIdx);
}
}
}
}
virtual void endEpisode()
{
const int episodeIdx = this->episodeIndex();
this->wellModel_.endEpisode();
this->aquiferModel_.endEpisode();
const auto& schedule = this->simulator().vanguard().schedule();
// End simulation when completed.
if (episodeIdx + 1 >= static_cast<int>(schedule.size()) - 1) {
this->simulator().setFinished(true);
return;
}
// Otherwise, start next episode (report step).
this->simulator().startNextEpisode(schedule.stepLength(episodeIdx + 1));
}
virtual void writeOutput(bool verbose)
{
OPM_TIMEBLOCK(problemWriteOutput);
if (Parameters::Get<Parameters::EnableWriteAllSolutions>() ||
this->episodeWillBeOver())
{
// Create VTK output as needed.
ParentType::writeOutput(verbose);
}
}
template <class Context>
const DimMatrix& intrinsicPermeability(const Context& context,
unsigned spaceIdx,
unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return transmissibilities_.permeability(globalSpaceIdx);
}
const DimMatrix& intrinsicPermeability(unsigned globalElemIdx) const
{ return transmissibilities_.permeability(globalElemIdx); }
template <class Context>
Scalar transmissibility(const Context& context,
[[maybe_unused]] unsigned fromDofLocalIdx,
unsigned toDofLocalIdx) const
{
assert(fromDofLocalIdx == 0);
return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
}
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
{
return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
}
template <class Context>
Scalar diffusivity(const Context& context,
[[maybe_unused]] unsigned fromDofLocalIdx,
unsigned toDofLocalIdx) const
{
assert(fromDofLocalIdx == 0);
return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
}
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
}
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
}
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx,
const unsigned boundaryFaceIdx) const
{
return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
}
template <class Context>
Scalar transmissibilityBoundary(const Context& elemCtx,
unsigned boundaryFaceIdx) const
{
unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
}
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx,
const unsigned boundaryFaceIdx) const
{
return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
}
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn,
const unsigned globalSpaceIdxOut) const
{
return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
}
template <class Context>
Scalar thermalHalfTransmissibilityIn(const Context& context,
unsigned faceIdx,
unsigned timeIdx) const
{
const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
unsigned toDofLocalIdx = face.exteriorIndex();
return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
}
template <class Context>
Scalar thermalHalfTransmissibilityOut(const Context& context,
unsigned faceIdx,
unsigned timeIdx) const
{
const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
unsigned toDofLocalIdx = face.exteriorIndex();
return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
}
template <class Context>
Scalar thermalHalfTransmissibilityBoundary(const Context& elemCtx,
unsigned boundaryFaceIdx) const
{
unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
}
const typename Vanguard::TransmissibilityType& eclTransmissibilities() const
{ return transmissibilities_; }
const TracerModel& tracerModel() const
{ return tracerModel_; }
TracerModel& tracerModel()
{ return tracerModel_; }
template <class Context>
Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return this->porosity(globalSpaceIdx, timeIdx);
}
template <class Context>
Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return this->dofCenterDepth(globalSpaceIdx);
}
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
{
return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
}
template <class Context>
Scalar rockCompressibility(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return this->rockCompressibility(globalSpaceIdx);
}
template <class Context>
Scalar rockReferencePressure(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return rockReferencePressure(globalSpaceIdx);
}
Scalar rockReferencePressure(unsigned globalSpaceIdx) const
{
const auto& rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
if (rock_config.store()) {
return asImp_().initialFluidState(globalSpaceIdx).pressure(refPressurePhaseIdx_());
}
else {
if (this->rockParams_.empty())
return 1e5;
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty()) {
tableIdx = this->rockTableIdx_[globalSpaceIdx];
}
return this->rockParams_[tableIdx].referencePressure;
}
}
template <class Context>
const MaterialLawParams& materialLawParams(const Context& context,
unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return this->materialLawParams(globalSpaceIdx);
}
const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const
{
return materialLawManager_->materialLawParams(globalDofIdx);
}
const MaterialLawParams& materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const
{
return materialLawManager_->materialLawParams(globalDofIdx, facedir);
}
template <class Context>
const SolidEnergyLawParams&
solidEnergyLawParams(const Context& context,
unsigned spaceIdx,
unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
}
template <class Context>
const ThermalConductionLawParams &
thermalConductionLawParams(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
}
std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
{ return materialLawManager_; }
template <class FluidState, class ...Args>
void updateRelperms(
std::array<Evaluation,numPhases> &mobility,
DirectionalMobilityPtr &dirMob,
FluidState &fluidState,
unsigned globalSpaceIdx) const
{
using ContainerT = std::array<Evaluation, numPhases>;
OPM_TIMEBLOCK_LOCAL(updateRelperms, Subsystem::SatProps);
{
// calculate relative permeabilities. note that we store the result into the
// mobility_ class attribute. the division by the phase viscosity happens later.
const auto& materialParams = materialLawParams(globalSpaceIdx);
MaterialLaw::template relativePermeabilities<ContainerT, FluidState, Args...>(mobility, materialParams, fluidState);
Valgrind::CheckDefined(mobility);
}
if (materialLawManager_->hasDirectionalRelperms()
|| materialLawManager_->hasDirectionalImbnum())
{
using Dir = FaceDir::DirEnum;
constexpr int ndim = 3;
dirMob = std::make_unique<DirectionalMobility<TypeTag>>();
Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
for (int i = 0; i<ndim; i++) {
const auto& materialParams = materialLawParams(globalSpaceIdx, facedirs[i]);
auto& mob_array = dirMob->getArray(i);
MaterialLaw::template relativePermeabilities<ContainerT, FluidState, Args...>(mob_array, materialParams, fluidState);
}
}
}
std::shared_ptr<EclMaterialLawManager> materialLawManager()
{ return materialLawManager_; }
using BaseType::pvtRegionIndex;
template <class Context>
unsigned pvtRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{ return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
using BaseType::satnumRegionIndex;
template <class Context>
unsigned satnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{ return this->satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
using BaseType::miscnumRegionIndex;
template <class Context>
unsigned miscnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{ return this->miscnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
using BaseType::plmixnumRegionIndex;
template <class Context>
unsigned plmixnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{ return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
// TODO: polymer related might need to go to the blackoil side
using BaseType::maxPolymerAdsorption;
template <class Context>
Scalar maxPolymerAdsorption(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{ return this->maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); }
std::string name() const
{ return this->simulator().vanguard().caseName(); }
template <class Context>
Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
// use the initial temperature of the DOF if temperature is not a primary
// variable
unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0);
}
Scalar temperature(unsigned globalDofIdx, unsigned /*timeIdx*/) const
{
// use the initial temperature of the DOF if temperature is not a primary
// variable
return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0);
}
const SolidEnergyLawParams&
solidEnergyLawParams(unsigned globalSpaceIdx,
unsigned /*timeIdx*/) const
{
return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
}
const ThermalConductionLawParams &
thermalConductionLawParams(unsigned globalSpaceIdx,
unsigned /*timeIdx*/)const
{
return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
}
Scalar maxOilSaturation(unsigned globalDofIdx) const
{
if (!this->vapparsActive(this->episodeIndex()))
return 0.0;
return this->maxOilSaturation_[globalDofIdx];
}
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
{
if (!this->vapparsActive(this->episodeIndex()))
return;
this->maxOilSaturation_[globalDofIdx] = value;
}
virtual void initialSolutionApplied()
{
// Calculate all intensive quantities.
this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/0);
// We also need the intensive quantities for timeIdx == 1
// corresponding to the start of the current timestep, if we
// do not use the storage cache, or if we cannot recycle the
// first iteration storage.
if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/1);
}
// initialize the wells. Note that this needs to be done after initializing the
// intrinsic permeabilities and the after applying the initial solution because
// the well model uses these...
wellModel_.init();
aquiferModel_.initialSolutionApplied();
const bool invalidateFromHyst = updateHysteresis_();
if (invalidateFromHyst) {
OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
}
}
template <class Context>
void source(RateVector& rate,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx) const
{
const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
source(rate, globalDofIdx, timeIdx);
}
void source(RateVector& rate,
unsigned globalDofIdx,
unsigned timeIdx) const
{
OPM_TIMEBLOCK_LOCAL(eclProblemSource, Subsystem::Assembly);
rate = 0.0;
// Add well contribution to source here.
wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
// convert the source term from the total mass rate of the
// cell to the one per unit of volume as used by the model.
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
Valgrind::CheckDefined(rate[eqIdx]);
assert(isfinite(rate[eqIdx]));
}
// Add non-well sources.
addToSourceDense(rate, globalDofIdx, timeIdx);
}
virtual void addToSourceDense(RateVector& rate,
unsigned globalDofIdx,
unsigned timeIdx) const = 0;
const WellModel& wellModel() const
{ return wellModel_; }
WellModel& wellModel()
{ return wellModel_; }
const AquiferModel& aquiferModel() const
{ return aquiferModel_; }
AquiferModel& mutableAquiferModel()
{ return aquiferModel_; }
bool nonTrivialBoundaryConditions() const
{ return nonTrivialBoundaryConditions_; }
Scalar nextTimeStepSize() const
{
OPM_TIMEBLOCK(nexTimeStepSize);
// allow external code to do the timestepping
if (this->nextTimeStepSize_ > 0.0)
return this->nextTimeStepSize_;
const auto& simulator = this->simulator();
int episodeIdx = simulator.episodeIndex();
// for the initial episode, we use a fixed time step size
if (episodeIdx < 0)
return this->initialTimeStepSize_;
// ask the newton method for a suggestion. This suggestion will be based on how
// well the previous time step converged. After that, apply the runtime time
// stepping constraints.
const auto& newtonMethod = this->model().newtonMethod();
return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
}
template <class LhsEval>
LhsEval rockCompPoroMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
{
OPM_TIMEBLOCK_LOCAL(rockCompPoroMultiplier, Subsystem::PvtProps);
if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
return 1.0;
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty())
tableIdx = this->rockTableIdx_[elementIdx];
const auto& fs = intQuants.fluidState();
LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
const auto& rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
if (!this->minRefPressure_.empty())
// The pore space change is irreversible
effectivePressure =
min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
this->minRefPressure_[elementIdx]);
if (!this->overburdenPressure_.empty())
effectivePressure -= this->overburdenPressure_[elementIdx];
if (rock_config.store()) {
effectivePressure -= asImp_().initialFluidState(elementIdx).pressure(refPressurePhaseIdx_());
}
if (!this->rockCompPoroMult_.empty()) {
return this->rockCompPoroMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
}
// water compaction
assert(!this->rockCompPoroMultWc_.empty());
LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
}
template <class LhsEval>
LhsEval rockCompTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
{
const bool implicit = !this->explicitRockCompaction_;
return implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<LhsEval>(intQuants, elementIdx)
: this->simulator().problem().getRockCompTransMultVal(elementIdx);
}
template <class LhsEval>
LhsEval wellTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
{
OPM_TIMEBLOCK_LOCAL(wellTransMultiplier, Subsystem::Wells);
const bool implicit = !this->explicitRockCompaction_;
double trans_mult = implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<double>(intQuants, elementIdx)
: this->simulator().problem().getRockCompTransMultVal(elementIdx);
trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants, elementIdx);
return trans_mult;
}
std::pair<BCType, RateVector> boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) const
{
OPM_TIMEBLOCK_LOCAL(boundaryCondition, Subsystem::Assembly);
if (!nonTrivialBoundaryConditions_) {
return { BCType::NONE, RateVector(0.0) };
}
FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
const auto& schedule = this->simulator().vanguard().schedule();
if (bcindex_(dir)[globalSpaceIdx] == 0) {
return { BCType::NONE, RateVector(0.0) };
}
if (schedule[this->episodeIndex()].bcprop.size() == 0) {
return { BCType::NONE, RateVector(0.0) };
}
const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
if (bc.bctype!=BCType::RATE) {
return { bc.bctype, RateVector(0.0) };
}
RateVector rate = 0.0;
switch (bc.component) {
case BCComponent::OIL:
rate[FluidSystem::canonicalToActiveCompIdx(oilCompIdx)] = bc.rate;
break;
case BCComponent::GAS:
rate[FluidSystem::canonicalToActiveCompIdx(gasCompIdx)] = bc.rate;
break;
case BCComponent::WATER:
rate[FluidSystem::canonicalToActiveCompIdx(waterCompIdx)] = bc.rate;
break;
case BCComponent::SOLVENT:
this->handleSolventBC(bc, rate);
break;
case BCComponent::POLYMER:
this->handlePolymerBC(bc, rate);
break;
case BCComponent::MICR:
this->handleMicrBC(bc, rate);
break;
case BCComponent::OXYG:
this->handleOxygBC(bc, rate);
break;
case BCComponent::UREA:
this->handleUreaBC(bc, rate);
break;
case BCComponent::NONE:
throw std::logic_error("you need to specify the component when RATE type is set in BC");
break;
}
//TODO add support for enthalpy rate
return {bc.bctype, rate};
}
template<class Serializer>
void serializeOp(Serializer& serializer)
{
serializer(static_cast<BaseType&>(*this));
serializer(drift_);
serializer(wellModel_);
serializer(aquiferModel_);
serializer(tracerModel_);
serializer(*materialLawManager_);
}
private:
Implementation& asImp_()
{ return *static_cast<Implementation *>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation *>(this); }
protected:
template<class UpdateFunc>
void updateProperty_(const std::string& failureMsg,
UpdateFunc func)
{
OPM_TIMEBLOCK(updateProperty);
const auto& model = this->simulator().model();
const auto& primaryVars = model.solution(/*timeIdx*/0);
const auto& vanguard = this->simulator().vanguard();
std::size_t numGridDof = primaryVars.size();
OPM_BEGIN_PARALLEL_TRY_CATCH();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
const auto& iq = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/ 0);
func(dofIdx, iq);
}
OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
}
bool updateMaxOilSaturation_()
{
OPM_TIMEBLOCK(updateMaxOilSaturation);
int episodeIdx = this->episodeIndex();
// we use VAPPARS
if (this->vapparsActive(episodeIdx)) {
this->updateProperty_("FlowProblem::updateMaxOilSaturation_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
this->updateMaxOilSaturation_(compressedDofIdx,iq);
});
return true;
}
return false;
}
bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation, Subsystem::SatProps);
const auto& fs = iq.fluidState();
const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
auto& mos = this->maxOilSaturation_;
if(mos[compressedDofIdx] < So){
mos[compressedDofIdx] = So;
return true;
}else{
return false;
}
}
bool updateMaxWaterSaturation_()
{
OPM_TIMEBLOCK(updateMaxWaterSaturation);
// water compaction is activated in ROCKCOMP
if (this->maxWaterSaturation_.empty())
return false;
this->maxWaterSaturation_[/*timeIdx=*/1] = this->maxWaterSaturation_[/*timeIdx=*/0];
this->updateProperty_("FlowProblem::updateMaxWaterSaturation_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
this->updateMaxWaterSaturation_(compressedDofIdx,iq);
});
return true;
}
bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation, Subsystem::SatProps);
const auto& fs = iq.fluidState();
const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
auto& mow = this->maxWaterSaturation_;
if(mow[compressedDofIdx]< Sw){
mow[compressedDofIdx] = Sw;
return true;
}else{
return false;
}
}
bool updateMinPressure_()
{
OPM_TIMEBLOCK(updateMinPressure);
// IRREVERS option is used in ROCKCOMP
if (this->minRefPressure_.empty())
return false;
this->updateProperty_("FlowProblem::updateMinPressure_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
this->updateMinPressure_(compressedDofIdx,iq);
});
return true;
}
bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
OPM_TIMEBLOCK_LOCAL(updateMinPressure, Subsystem::PvtProps);
const auto& fs = iq.fluidState();
const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
auto& min_pressures = this->minRefPressure_;
if(min_pressures[compressedDofIdx]> min_pressure){
min_pressures[compressedDofIdx] = min_pressure;
return true;
}else{
return false;
}
}
// \brief Function to assign field properties of type double, on the leaf grid view.
//
// For CpGrid with local grid refinement, the field property of a cell on the leaf
// is inherited from its parent or equivalent (when has no parent) cell on level zero.
std::function<std::vector<double>(const FieldPropsManager&, const std::string&)>
fieldPropDoubleOnLeafAssigner_()
{
const auto& lookup = this->lookUpData_;
return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString)
{
return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
};
}
// \brief Function to assign field properties of type int, unsigned int, ..., on the leaf grid view.
//
// For CpGrid with local grid refinement, the field property of a cell on the leaf
// is inherited from its parent or equivalent (when has no parent) cell on level zero.
template<typename IntType>
std::function<std::vector<IntType>(const FieldPropsManager&, const std::string&, bool)>
fieldPropIntTypeOnLeafAssigner_()
{
const auto& lookup = this->lookUpData_;
return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString, bool needsTranslation)
{
return lookup.template assignFieldPropsIntOnLeaf<IntType>(fieldPropManager, propString, needsTranslation);
};
}
void readMaterialParameters_()
{
OPM_TIMEBLOCK(readMaterialParameters);
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
// the PVT and saturation region numbers
OPM_BEGIN_PARALLEL_TRY_CATCH();
this->updatePvtnum_();
this->updateSatnum_();
// the MISC region numbers (solvent model)
this->updateMiscnum_();
// the PLMIX region numbers (polymer model)
this->updatePlmixnum_();
OPM_END_PARALLEL_TRY_CATCH("Invalid region numbers: ", vanguard.gridView().comm());
// porosity
updateReferencePorosity_();
this->referencePorosity_[1] = this->referencePorosity_[0];
// fluid-matrix interactions (saturation functions; relperm/capillary pressure)
materialLawManager_ = std::make_shared<EclMaterialLawManager>();
materialLawManager_->initFromState(eclState);
materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
this-> template fieldPropIntTypeOnLeafAssigner_<int>(),
this-> lookupIdxOnLevelZeroAssigner_());
}
void readThermalParameters_()
{
if constexpr (enableEnergy)
{
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
// fluid-matrix interactions (saturation functions; relperm/capillary pressure)
thermalLawManager_ = std::make_shared<EclThermalLawManager>();
thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
this-> fieldPropDoubleOnLeafAssigner_(),
this-> template fieldPropIntTypeOnLeafAssigner_<unsigned int>());
}
}
void updateReferencePorosity_()
{
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
std::size_t numDof = this->model().numGridDof();
this->referencePorosity_[/*timeIdx=*/0].resize(numDof);
const auto& fp = eclState.fieldProps();
const std::vector<double> porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV");
for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(dofIdx);
Scalar poreVolume = porvData[dofIdx];
// we define the porosity as the accumulated pore volume divided by the
// geometric volume of the element. Note that -- in pathetic cases -- it can
// be larger than 1.0!
Scalar dofVolume = simulator.model().dofTotalVolume(sfcdofIdx);
assert(dofVolume > 0.0);
this->referencePorosity_[/*timeIdx=*/0][sfcdofIdx] = poreVolume/dofVolume;
}
}
virtual void readInitialCondition_()
{
// TODO: whether we should move this to FlowProblemBlackoil
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
if (eclState.getInitConfig().hasEquil())
readEquilInitialCondition_();
else
readExplicitInitialCondition_();
//initialize min/max values
std::size_t numElems = this->model().numGridDof();
for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
const auto& fs = asImp_().initialFluidStates()[elemIdx];
if (!this->maxWaterSaturation_.empty() && waterPhaseIdx > -1)
this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
if (!this->maxOilSaturation_.empty() && oilPhaseIdx > -1)
this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
if (!this->minRefPressure_.empty() && refPressurePhaseIdx_() > -1)
this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
}
}
virtual void readEquilInitialCondition_() = 0;
virtual void readExplicitInitialCondition_() = 0;
// update the hysteresis parameters of the material laws for the whole grid
bool updateHysteresis_()
{
if (!materialLawManager_->enableHysteresis())
return false;
// we need to update the hysteresis data for _all_ elements (i.e., not just the
// interior ones) to avoid desynchronization of the processes in the parallel case!
this->updateProperty_("FlowProblem::updateHysteresis_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
});
return true;
}
bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
OPM_TIMEBLOCK_LOCAL(updateHysteresis_, Subsystem::SatProps);
materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
//TODO change materials to give a bool
return true;
}
Scalar getRockCompTransMultVal(std::size_t dofIdx) const
{
if (this->rockCompTransMultVal_.empty())
return 1.0;
return this->rockCompTransMultVal_[dofIdx];
}
protected:
struct PffDofData_
{
ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
ConditionalStorage<enableDiffusion, Scalar> diffusivity;
ConditionalStorage<enableDispersion, Scalar> dispersivity;
Scalar transmissibility;
};
// update the prefetch friendly data object
void updatePffDofData_()
{
const auto& distFn =
[this](PffDofData_& dofData,
const Stencil& stencil,
unsigned localDofIdx)
-> void
{
const auto& elementMapper = this->model().elementMapper();
unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
if (localDofIdx != 0) {
unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(/*dofIdx=*/0));
dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
if constexpr (enableEnergy) {
*dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
*dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
}
if constexpr (enableDiffusion)
*dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
if (enableDispersion)
dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
}
};
pffDofData_.update(distFn);
}
virtual void updateExplicitQuantities_(int episodeIdx, int timeStepSize, bool first_step_after_restart) = 0;
void readBoundaryConditions_()
{
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
if (bcconfig.size() > 0) {
nonTrivialBoundaryConditions_ = true;
std::size_t numCartDof = vanguard.cartesianSize();
unsigned numElems = vanguard.gridView().size(/*codim=*/0);
std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
bcindex_.resize(numElems, 0);
auto loopAndApply = [&cartesianToCompressedElemIdx,
&vanguard](const auto& bcface,
auto apply)
{
for (int i = bcface.i1; i <= bcface.i2; ++i) {
for (int j = bcface.j1; j <= bcface.j2; ++j) {
for (int k = bcface.k1; k <= bcface.k2; ++k) {
std::array<int, 3> tmp = {i,j,k};
auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
if (elemIdx >= 0)
apply(elemIdx);
}
}
}
};
for (const auto& bcface : bcconfig) {
std::vector<int>& data = bcindex_(bcface.dir);
const int index = bcface.index;
loopAndApply(bcface,
[&data,index](int elemIdx)
{ data[elemIdx] = index; });
}
}
}
// this method applies the runtime constraints specified via the deck and/or command
// line parameters for the size of the next time step.
Scalar limitNextTimeStepSize_(Scalar dtNext) const
{
if constexpr (enableExperiments) {
const auto& simulator = this->simulator();
const auto& schedule = simulator.vanguard().schedule();
int episodeIdx = simulator.episodeIndex();
// first thing in the morning, limit the time step size to the maximum size
Scalar maxTimeStepSize = Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60;
int reportStepIdx = std::max(episodeIdx, 0);
if (this->enableTuning_) {
const auto& tuning = schedule[reportStepIdx].tuning();
maxTimeStepSize = tuning.TSMAXZ;
}
dtNext = std::min(dtNext, maxTimeStepSize);
Scalar remainingEpisodeTime =
simulator.episodeStartTime() + simulator.episodeLength()
- (simulator.startTime() + simulator.time());
assert(remainingEpisodeTime >= 0.0);
// if we would have a small amount of time left over in the current episode, make
// two equal time steps instead of a big and a small one
if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
// note: limiting to the maximum time step size here is probably not strictly
// necessary, but it should not hurt and is more fool-proof
dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
if (simulator.episodeStarts()) {
// if a well event occurred, respect the limit for the maximum time step after
// that, too
const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
bool wellEventOccured =
events.hasEvent(ScheduleEvents::NEW_WELL)
|| events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
|| events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
|| events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
}
}
return dtNext;
}
int refPressurePhaseIdx_() const {
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
return oilPhaseIdx;
}
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
return gasPhaseIdx;
}
else {
return waterPhaseIdx;
}
}
void updateRockCompTransMultVal_()
{
const auto& model = this->simulator().model();
std::size_t numGridDof = this->model().numGridDof();
this->rockCompTransMultVal_.resize(numGridDof, 1.0);
for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
const auto& iq = *model.cachedIntensiveQuantities(elementIdx, /*timeIdx=*/ 0);
Scalar trans_mult = computeRockCompTransMultiplier_<Scalar>(iq, elementIdx);
this->rockCompTransMultVal_[elementIdx] = trans_mult;
}
}
template <class LhsEval>
LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities& intQuants, unsigned elementIdx) const
{
OPM_TIMEBLOCK_LOCAL(computeRockCompTransMultiplier, Subsystem::PvtProps);
if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
return 1.0;
unsigned tableIdx = 0;
if (!this->rockTableIdx_.empty())
tableIdx = this->rockTableIdx_[elementIdx];
const auto& fs = intQuants.fluidState();
LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
const auto& rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
if (!this->minRefPressure_.empty())
// The pore space change is irreversible
effectivePressure =
min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
this->minRefPressure_[elementIdx]);
if (!this->overburdenPressure_.empty())
effectivePressure -= this->overburdenPressure_[elementIdx];
if (rock_config.store()) {
effectivePressure -= asImp_().initialFluidState(elementIdx).pressure(refPressurePhaseIdx_());
}
if (!this->rockCompTransMult_.empty())
return this->rockCompTransMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
// water compaction
assert(!this->rockCompTransMultWc_.empty());
LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
}
typename Vanguard::TransmissibilityType transmissibilities_;
std::shared_ptr<EclMaterialLawManager> materialLawManager_;
std::shared_ptr<EclThermalLawManager> thermalLawManager_;
GlobalEqVector drift_;
WellModel wellModel_;
AquiferModel aquiferModel_;
PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
TracerModel tracerModel_;
template<class T>
struct BCData
{
std::array<std::vector<T>,6> data;
void resize(std::size_t size, T defVal)
{
for (auto& d : data)
d.resize(size, defVal);
}
const std::vector<T>& operator()(FaceDir::DirEnum dir) const
{
if (dir == FaceDir::DirEnum::Unknown)
throw std::runtime_error("Tried to access BC data for the 'Unknown' direction");
int idx = 0;
int div = static_cast<int>(dir);
while ((div /= 2) >= 1)
++idx;
assert(idx >= 0 && idx <= 5);
return data[idx];
}
std::vector<T>& operator()(FaceDir::DirEnum dir)
{
return const_cast<std::vector<T>&>(std::as_const(*this)(dir));
}
};
virtual void handleSolventBC(const BCProp::BCFace&, RateVector&) const = 0;
virtual void handlePolymerBC(const BCProp::BCFace&, RateVector&) const = 0;
virtual void handleMicrBC(const BCProp::BCFace&, RateVector&) const = 0;
virtual void handleOxygBC(const BCProp::BCFace&, RateVector&) const = 0;
virtual void handleUreaBC(const BCProp::BCFace&, RateVector&) const = 0;
BCData<int> bcindex_;
bool nonTrivialBoundaryConditions_ = false;
bool first_step_ = true;
virtual bool episodeWillBeOver() const
{
return this->simulator().episodeWillBeOver();
}
};
} // namespace Opm
#endif // OPM_FLOW_PROBLEM_HPP
CpGridVanguard.hpp
Helper class for grid instantiation of ECL file-format using problems.
DummyGradientCalculator.hpp
This is a "dummy" gradient calculator which does not do anything.
EclWriter.hpp
Collects necessary output values and pass it to opm-common's ECL output.
EquilInitializer.hpp
Computes the initial condition based on the EQUIL keyword from ECL.
FlowBaseProblemProperties.hpp
FlowGenericProblem.hpp
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
TracerModel.hpp
A class which handles tracers as specified in by ECL.
Transmissibility.hpp
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::maxPolymerAdsorption
Scalar maxPolymerAdsorption(unsigned elemIdx) const
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::pvtRegionIndex
unsigned pvtRegionIndex(unsigned elemIdx) const
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::briefDescription
static std::string briefDescription()
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::porosity
Scalar porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::rockCompressibility
Scalar rockCompressibility(unsigned globalSpaceIdx) const
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::miscnumRegionIndex
unsigned miscnumRegionIndex(unsigned elemIdx) const
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::satnumRegionIndex
unsigned satnumRegionIndex(unsigned elemIdx) const
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::plmixnumRegionIndex
unsigned plmixnumRegionIndex(unsigned elemIdx) const
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::shouldWriteOutput
bool shouldWriteOutput() const
Definition FlowGenericProblem.hpp:277
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::helpPreamble
static std::string helpPreamble(int, const char **argv)
Opm::FlowGenericProblem< GetPropType< TypeTag, Properties::GridView >, GetPropType< TypeTag, Properties::FluidSystem > >::shouldWriteRestartFile
bool shouldWriteRestartFile() const
Definition FlowGenericProblem.hpp:286
Opm::FlowProblem
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:94
Opm::FlowProblem::episodeWillBeOver
virtual bool episodeWillBeOver() const
Whether or not the current episode will end at the end of the current time step.
Definition FlowProblem.hpp:1733
Opm::FlowProblem::wellModel
const WellModel & wellModel() const
Returns a reference to the ECL well manager used by the problem.
Definition FlowProblem.hpp:1012
Opm::FlowProblem::transmissibility
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
Direct access to the transmissibility between two elements.
Definition FlowProblem.hpp:527
Opm::FlowProblem::handlePositionalParameter
static int handlePositionalParameter(std::function< void(const std::string &, const std::string &)> addKey, std::set< std::string > &seenParams, std::string &errorMsg, int, const char **argv, int paramIdx, int)
Handles positional command line parameters.
Definition FlowProblem.hpp:197
Opm::FlowProblem::writeOutput
virtual void writeOutput(bool verbose)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:479
Opm::FlowProblem::intrinsicPermeability
const DimMatrix & intrinsicPermeability(unsigned globalElemIdx) const
This method returns the intrinsic permeability tensor given a global element index.
Definition FlowProblem.hpp:509
Opm::FlowProblem::pvtRegionIndex
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:826
Opm::FlowProblem::porosity
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:659
Opm::FlowProblem::beginIteration
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition FlowProblem.hpp:384
Opm::FlowProblem::satnumRegionIndex
unsigned satnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:834
Opm::FlowProblem::thermalHalfTransmissibilityOut
Scalar thermalHalfTransmissibilityOut(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:617
Opm::FlowProblem::transmissibilityBoundary
Scalar transmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:574
Opm::FlowProblem::thermalHalfTransmissibility
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn, const unsigned globalSpaceIdxOut) const
Definition FlowProblem.hpp:594
Opm::FlowProblem::rockCompressibility
Scalar rockCompressibility(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:693
Opm::FlowProblem::maxOilSaturation
Scalar maxOilSaturation(unsigned globalDofIdx) const
Returns an element's historic maximum oil phase saturation that was observed during the simulation.
Definition FlowProblem.hpp:910
Opm::FlowProblem::name
std::string name() const
The problem name.
Definition FlowProblem.hpp:865
Opm::FlowProblem::rockCompTransMultiplier
LhsEval rockCompTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1104
Opm::FlowProblem::endIteration
void endIteration()
Called by the simulator after each Newton-Raphson iteration.
Definition FlowProblem.hpp:394
Opm::FlowProblem::FlowProblem
FlowProblem(Simulator &simulator)
Definition FlowProblem.hpp:216
Opm::FlowProblem::eclTransmissibilities
const Vanguard::TransmissibilityType & eclTransmissibilities() const
Return a reference to the object that handles the "raw" transmissibilities.
Definition FlowProblem.hpp:640
Opm::FlowProblem::nextTimeStepSize
Scalar nextTimeStepSize() const
Propose the size of the next time step to the simulator.
Definition FlowProblem.hpp:1033
Opm::FlowProblem::computeRockCompTransMultiplier_
LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1638
Opm::FlowProblem::thermalConductionLawParams
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:769
Opm::FlowProblem::transmissibility
Scalar transmissibility(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:516
Opm::FlowProblem::beginEpisode
virtual void beginEpisode()
Called by the simulator before an episode begins.
Definition FlowProblem.hpp:295
Opm::FlowProblem::beginTimeStep
virtual void beginTimeStep()
Called by the simulator before each time integration.
Definition FlowProblem.hpp:352
Opm::FlowProblem::dofCenterDepth
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
Direct indexed acces to the depth of an degree of freedom [m].
Definition FlowProblem.hpp:684
Opm::FlowProblem::wellTransMultiplier
LhsEval wellTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Return the well transmissibility multiplier due to rock changes.
Definition FlowProblem.hpp:1116
Opm::FlowProblem::materialLawManager
std::shared_ptr< const EclMaterialLawManager > materialLawManager() const
Returns the ECL material law manager.
Definition FlowProblem.hpp:781
Opm::FlowProblem::plmixnumRegionIndex
unsigned plmixnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:850
Opm::FlowProblem::temperature
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:872
Opm::FlowProblem::thermalHalfTransmissibilityBoundary
Scalar thermalHalfTransmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:630
Opm::FlowProblem::miscnumRegionIndex
unsigned miscnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:842
Opm::FlowProblem::intrinsicPermeability
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:495
Opm::FlowProblem::maxPolymerAdsorption
Scalar maxPolymerAdsorption(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the max polymer adsorption value.
Definition FlowProblem.hpp:859
Opm::FlowProblem::thermalHalfTransmissibilityIn
Scalar thermalHalfTransmissibilityIn(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:604
Opm::FlowProblem::setMaxOilSaturation
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
Sets an element's maximum oil phase saturation observed during the simulation.
Definition FlowProblem.hpp:927
Opm::FlowProblem::solidEnergyLawParams
const SolidEnergyLawParams & solidEnergyLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the parameters for the energy storage law of the rock.
Definition FlowProblem.hpp:756
Opm::FlowProblem::dispersivity
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the dispersivity for a face i.e.
Definition FlowProblem.hpp:554
Opm::FlowProblem::registerParameters
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblem.hpp:181
Opm::FlowProblem::endEpisode
virtual void endEpisode()
Called by the simulator after the end of an episode.
Definition FlowProblem.hpp:456
Opm::FlowProblem::source
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition FlowProblem.hpp:971
Opm::FlowProblem::endTimeStep
virtual void endTimeStep()
Called by the simulator after each time integration.
Definition FlowProblem.hpp:404
Opm::FlowProblem::initialSolutionApplied
virtual void initialSolutionApplied()
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition FlowProblem.hpp:938
Opm::FlowProblem::materialLawParams
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:734
Opm::FlowProblem::thermalTransmissibilityBoundary
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:561
Opm::FlowProblem::diffusivity
Scalar diffusivity(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:536
Opm::FlowProblem::rockReferencePressure
Scalar rockReferencePressure(unsigned globalSpaceIdx) const
Definition FlowProblem.hpp:712
Opm::FlowProblem::deserialize
void deserialize(Restarter &res)
This method restores the complete state of the problem and its sub-objects from disk.
Definition FlowProblem.hpp:261
Opm::FlowProblem::materialLawManager
std::shared_ptr< EclMaterialLawManager > materialLawManager()
Definition FlowProblem.hpp:818
Opm::FlowProblem::dofCenterDepth
Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the depth of an degree of freedom [m].
Definition FlowProblem.hpp:672
Opm::FlowProblem::rockCompPoroMultiplier
LhsEval rockCompPoroMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the porosity multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1060
Opm::FlowProblem::transmissibilityBoundary
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:584
Opm::FlowProblem::serialize
void serialize(Restarter &res)
This method writes the complete state of the problem and its subobjects to disk.
Definition FlowProblem.hpp:280
Opm::FlowProblem::rockReferencePressure
Scalar rockReferencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:703
Opm::FlowProblem::diffusivity
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the transmissibility for a face i.e.
Definition FlowProblem.hpp:547
Opm::PffGridVector
A random-access container which stores data attached to a grid's degrees of freedom in a prefetch fri...
Definition pffgridvector.hh:50
Opm::RelpermDiagnostics
This class is intend to be a relperm diagnostics, to detect wrong input of relperm table and endpoint...
Definition RelpermDiagnostics.hpp:51
Opm::RelpermDiagnostics::diagnosis
void diagnosis(const EclipseState &eclState, const LevelCartesianIndexMapper &levelCartesianIndexMapper)
This function is used to diagnosis relperm in eclipse data file.
Definition RelpermDiagnostics.cpp:830
Opm::TracerModel
A class which handles tracers as specified in by ECL.
Definition TracerModel.hpp:75
directionalmobility.hh
This file contains definitions related to directional mobilities.
ecfvdiscretization.hh
The base class for the element-centered finite-volume discretization scheme.
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
Opm::prefetch
void prefetch(const T &val, unsigned n=1)
Template function which emits prefetch instructions for a range of memory.
Definition prefetch.hh:39
Opm::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:233
Opm::getPropValue
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:240
Opm::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition propertysystem.hh:224
Opm::Parameters::Get
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition parametersystem.hpp:187
pffgridvector.hh
Opm::FlowProblem::BCData
Definition FlowProblem.hpp:1690
Opm::FlowProblem::PffDofData_
Definition FlowProblem.hpp:1479