30#ifndef OPM_FLOW_PROBLEM_HPP
31#define OPM_FLOW_PROBLEM_HPP
33#include <dune/common/version.hh>
34#include <dune/common/fvector.hh>
35#include <dune/common/fmatrix.hh>
37#include <opm/common/utility/TimeService.hpp>
41#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
42#include <opm/input/eclipse/Parser/ParserKeywords/E.hpp>
43#include <opm/input/eclipse/Schedule/Schedule.hpp>
45#include <opm/material/common/ConditionalStorage.hpp>
46#include <opm/material/common/Valgrind.hpp>
47#include <opm/material/densead/Evaluation.hpp>
48#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
49#include <opm/material/fluidstates/CompositionalFluidState.hpp>
50#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
51#include <opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp>
52#include <opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp>
53#include <opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp>
54#include <opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp>
55#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp>
56#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp>
57#include <opm/material/thermal/EclThermalLawManager.hpp>
59#include <opm/models/common/directionalmobility.hh>
60#include <opm/models/utils/pffgridvector.hh>
61#include <opm/models/blackoil/blackoilmodel.hh>
62#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
64#include <opm/output/eclipse/EclipseIO.hpp>
86#include <opm/utility/CopyablePtr.hpp>
88#include <opm/common/OpmLog/OpmLog.hpp>
108template <
class TypeTag>
109class FlowProblem :
public GetPropType<TypeTag, Properties::BaseProblem>
111 GetPropType<TypeTag, Properties::FluidSystem>>
114 GetPropType<TypeTag, Properties::FluidSystem>>;
115 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
116 using Implementation = GetPropType<TypeTag, Properties::Problem>;
118 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
119 using GridView = GetPropType<TypeTag, Properties::GridView>;
120 using Stencil = GetPropType<TypeTag, Properties::Stencil>;
121 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
122 using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
123 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
124 using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
127 enum { dim = GridView::dimension };
128 enum { dimWorld = GridView::dimensionworld };
131 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
132 enum { numPhases = FluidSystem::numPhases };
133 enum { numComponents = FluidSystem::numComponents };
134 enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
135 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
136 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
137 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
138 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
139 enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
140 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
141 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
142 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
143 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
144 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
145 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
146 enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
147 enum { enableApiTracking = getPropValue<TypeTag, Properties::EnableApiTracking>() };
148 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
149 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
150 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
151 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
152 enum { gasCompIdx = FluidSystem::gasCompIdx };
153 enum { oilCompIdx = FluidSystem::oilCompIdx };
154 enum { waterCompIdx = FluidSystem::waterCompIdx };
156 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
157 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
158 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
159 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
160 using Element =
typename GridView::template Codim<0>::Entity;
161 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
162 using EclMaterialLawManager =
typename GetProp<TypeTag, Properties::MaterialLaw>::EclMaterialLawManager;
163 using EclThermalLawManager =
typename GetProp<TypeTag, Properties::SolidEnergyLaw>::EclThermalLawManager;
164 using MaterialLawParams =
typename EclMaterialLawManager::MaterialLawParams;
165 using SolidEnergyLawParams =
typename EclThermalLawManager::SolidEnergyLawParams;
166 using ThermalConductionLawParams =
typename EclThermalLawManager::ThermalConductionLawParams;
167 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
168 using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
169 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
170 using Indices = GetPropType<TypeTag, Properties::Indices>;
171 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
172 using WellModel = GetPropType<TypeTag, Properties::WellModel>;
173 using AquiferModel = GetPropType<TypeTag, Properties::AquiferModel>;
175 using SolventModule = BlackOilSolventModule<TypeTag>;
176 using PolymerModule = BlackOilPolymerModule<TypeTag>;
177 using FoamModule = BlackOilFoamModule<TypeTag>;
178 using BrineModule = BlackOilBrineModule<TypeTag>;
179 using ExtboModule = BlackOilExtboModule<TypeTag>;
180 using MICPModule = BlackOilMICPModule<TypeTag>;
181 using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
182 using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
186 using Toolbox = MathToolbox<Evaluation>;
187 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
195 using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
211 ParentType::registerParameters();
214 DamarisWriterType::registerParameters();
219 Parameters::registerParam<TypeTag, Properties::EnableWriteAllSolutions>
220 (
"Write all solutions to disk instead of only the ones for the "
222 Parameters::registerParam<TypeTag, Properties::EnableEclOutput>
223 (
"Write binary output which is compatible with the commercial "
224 "Eclipse simulator");
226 Parameters::registerParam<TypeTag, Properties::EnableDamarisOutput>
227 (
"Write a specific variable using Damaris in a separate core");
229 Parameters::registerParam<TypeTag, Properties::EclOutputDoublePrecision>
230 (
"Tell the output writer to use double precision. Useful for 'perfect' restarts");
231 Parameters::registerParam<TypeTag, Properties::RestartWritingInterval>
232 (
"The frequencies of which time steps are serialized to disk");
233 Parameters::registerParam<TypeTag, Properties::EnableDriftCompensation>
234 (
"Enable partial compensation of systematic mass losses via "
235 "the source term of the next time step");
236 Parameters::registerParam<TypeTag, Properties::OutputMode>
237 (
"Specify which messages are going to be printed. "
238 "Valid values are: none, log, all (default)");
239 Parameters::registerParam<TypeTag, Properties::NumPressurePointsEquil>
240 (
"Number of pressure points (in each direction) in tables used for equilibration");
241 Parameters::hideParam<TypeTag, Properties::NumPressurePointsEquil>();
242 Parameters::registerParam<TypeTag, Properties::ExplicitRockCompaction>
243 (
"Use pressure from end of the last time step when evaluating rock compaction");
244 Parameters::hideParam<TypeTag, Properties::ExplicitRockCompaction>();
252 std::string& errorMsg,
258 using ParamsMeta = GetProp<TypeTag, Properties::ParameterMetaData>;
259 Dune::ParameterTree& tree = ParamsMeta::tree();
271 : ParentType(simulator)
272 ,
BaseType(simulator.vanguard().eclState(),
273 simulator.vanguard().schedule(),
274 simulator.vanguard().gridView())
275 , transmissibilities_(simulator.vanguard().eclState(),
276 simulator.vanguard().gridView(),
277 simulator.vanguard().cartesianIndexMapper(),
278 simulator.vanguard().grid(),
279 simulator.vanguard().cellCentroids(),
283 , thresholdPressures_(simulator)
284 , wellModel_(simulator)
285 , aquiferModel_(simulator)
286 , pffDofData_(simulator.gridView(), this->elementMapper())
287 , tracerModel_(simulator)
288 , actionHandler_(simulator.vanguard().eclState(),
289 simulator.vanguard().schedule(),
290 simulator.vanguard().actionState(),
291 simulator.vanguard().summaryState(),
293 simulator.vanguard().grid().comm())
297 const auto& vanguard = simulator.vanguard();
298 SolventModule::initFromState(vanguard.eclState(), vanguard.schedule());
299 PolymerModule::initFromState(vanguard.eclState());
300 FoamModule::initFromState(vanguard.eclState());
301 BrineModule::initFromState(vanguard.eclState());
302 ExtboModule::initFromState(vanguard.eclState());
303 MICPModule::initFromState(vanguard.eclState());
304 DispersionModule::initFromState(vanguard.eclState());
305 DiffusionModule::initFromState(vanguard.eclState());
308 eclWriter_ = std::make_unique<EclWriterType>(simulator);
311 damarisWriter_ = std::make_unique<DamarisWriterType>(simulator);
312 enableDamarisOutput_ = Parameters::get<TypeTag, Properties::EnableDamarisOutput>();
314 enableDriftCompensation_ = Parameters::get<TypeTag, Properties::EnableDriftCompensation>();
316 enableEclOutput_ = Parameters::get<TypeTag, Properties::EnableEclOutput>();
318 this->
enableTuning_ = Parameters::get<TypeTag, Properties::EnableTuning>();
326 if (Parameters::isSet<TypeTag, Properties::NumPressurePointsEquil>())
330 this->
numPressurePointsEquil_ = simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP();
333 explicitRockCompaction_ = Parameters::get<TypeTag, Properties::ExplicitRockCompaction>();
337 relpermDiagnostics.
diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
345 ParentType::finishInit();
347 auto& simulator = this->simulator();
348 const auto& eclState = simulator.vanguard().eclState();
349 const auto& schedule = simulator.vanguard().schedule();
352 simulator.setStartTime(schedule.getStartTime());
353 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
359 simulator.setEpisodeIndex(-1);
360 simulator.setEpisodeLength(0.0);
365 this->gravity_ = 0.0;
366 if (Parameters::get<TypeTag, Properties::EnableGravity>())
367 this->gravity_[dim - 1] = 9.80665;
368 if (!eclState.getInitConfig().hasGravity())
369 this->gravity_[dim - 1] = 0.0;
374 const auto& tuning = schedule[0].tuning();
384 eclState.runspec().tabdims().getNumPVTTables());
386 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
387 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
392 [&simulator](
const unsigned idx)
394 std::array<int,dim> coords;
395 simulator.vanguard().cartesianCoordinate(idx, coords);
396 for (auto& c : coords) {
405 std::function<
unsigned int(
unsigned int)> gridToEquilGrid = [&simulator](
unsigned int i) {
406 return simulator.vanguard().gridIdxToEquilGridIdx(i);
408 transmissibilities_.finishInit(gridToEquilGrid);
410 const auto& initconfig = eclState.getInitConfig();
411 tracerModel_.init(initconfig.restartRequested());
412 if (initconfig.restartRequested())
417 tracerModel_.prepareTracerBatches();
421 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
422 const auto& vanguard = this->simulator().vanguard();
423 const auto& gridView = vanguard.gridView();
424 int numElements = gridView.size(0);
425 this->
polymer_.maxAdsorption.resize(numElements, 0.0);
428 readBoundaryConditions_();
431 computeAndSetEqWeights_();
433 if (enableDriftCompensation_) {
434 drift_.resize(this->model().numGridDof());
439 if (enableEclOutput_) {
440 if (simulator.vanguard().grid().comm().size() > 1) {
441 if (simulator.vanguard().grid().comm().rank() == 0)
442 eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
444 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
447 std::function<
unsigned int(
unsigned int)> equilGridToGrid = [&simulator](
unsigned int i) {
448 return simulator.vanguard().gridEquilIdxToGridIdx(i);
450 eclWriter_->writeInit(equilGridToGrid);
453 simulator.vanguard().releaseGlobalTransmissibilities();
458 if (!initconfig.restartRequested()) {
459 simulator.startNextEpisode(schedule.seconds(1));
460 simulator.setEpisodeIndex(0);
465 { pffDofData_.prefetch(elem); }
478 template <
class Restarter>
485 wellModel_.deserialize(res);
488 aquiferModel_.deserialize(res);
497 template <
class Restarter>
500 wellModel_.serialize(res);
502 aquiferModel_.serialize(res);
507 return std::max(this->simulator().episodeIndex(), 0);
515 OPM_TIMEBLOCK(beginEpisode);
517 auto& simulator = this->simulator();
518 int episodeIdx = simulator.episodeIndex();
519 auto& eclState = simulator.vanguard().eclState();
520 const auto& schedule = simulator.vanguard().schedule();
521 const auto& events = schedule[episodeIdx].events();
523 if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
530 const auto& miniDeck = schedule[episodeIdx].geo_keywords();
531 const auto& cc = simulator.vanguard().grid().comm();
532 eclState.apply_schedule_keywords( miniDeck );
536 std::function<
unsigned int(
unsigned int)> equilGridToGrid = [&simulator](
unsigned int i) {
537 return simulator.vanguard().gridEquilIdxToGridIdx(i);
541 transmissibilities_.update(
true, equilGridToGrid);
542 this->referencePorosity_[1] = this->referencePorosity_[0];
543 updateReferencePorosity_();
545 this->model().linearizer().updateDiscretizationParameters();
548 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
551 wellModel_.beginEpisode();
554 aquiferModel_.beginEpisode();
557 Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
559 if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
562 dt = std::min(dt, this->initialTimeStepSize_);
563 simulator.setTimeStepSize(dt);
567 actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
569 if (episodeIdx >= 0) {
570 const auto& oilVap = schedule[episodeIdx].oilvap();
571 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
572 FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2());
574 FluidSystem::setVapPars(0.0, 0.0);
584 OPM_TIMEBLOCK(beginTimeStep);
585 int episodeIdx = this->episodeIndex();
587 this->beginTimeStep_(enableExperiments,
589 this->simulator().timeStepIndex(),
590 this->simulator().startTime(),
591 this->simulator().time(),
592 this->simulator().timeStepSize(),
593 this->simulator().endTime());
597 asImp_().updateExplicitQuantities_();
599 if (nonTrivialBoundaryConditions()) {
600 this->model().linearizer().updateBoundaryConditionData();
603 wellModel_.beginTimeStep();
604 aquiferModel_.beginTimeStep();
605 tracerModel_.beginTimeStep();
614 OPM_TIMEBLOCK(beginIteration);
615 wellModel_.beginIteration();
616 aquiferModel_.beginIteration();
625 wellModel_.endIteration();
626 aquiferModel_.endIteration();
634 OPM_TIMEBLOCK(endTimeStep);
636 if constexpr (getPropValue<TypeTag, Properties::EnableDebuggingChecks>()) {
641 int rank = this->simulator().gridView().comm().rank();
643 std::cout <<
"checking conservativeness of solution\n";
644 this->model().checkConservativeness(-1,
true);
646 std::cout <<
"solution is sufficiently conservative\n";
650 auto& simulator = this->simulator();
651 wellModel_.endTimeStep();
652 aquiferModel_.endTimeStep();
653 tracerModel_.endTimeStep();
657 this->model().linearizer().updateFlowsInfo();
660 asImp_().updateCompositionChangeLimits_();
662 OPM_TIMEBLOCK(driftCompansation);
663 if (enableDriftCompensation_) {
664 const auto& residual = this->model().linearizer().residual();
665 for (
unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
666 drift_[globalDofIdx] = residual[globalDofIdx];
667 drift_[globalDofIdx] *= simulator.timeStepSize();
668 if constexpr (getPropValue<TypeTag, Properties::UseVolumetricResidual>())
669 drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
673 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
674 !this->simulator().episodeWillBeOver();
676 const auto& grid = this->simulator().vanguard().gridView().grid();
677 using GridType = std::remove_cv_t<
typename std::remove_reference<
decltype(grid)>::type>;
678 bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
679 if ( !isCpGrid || (this->simulator().vanguard().gridView().grid().maxLevel()==0)) {
680 eclWriter_->evalSummaryState(isSubStep);
683 int episodeIdx = this->episodeIndex();
686 std::function<
unsigned int(
unsigned int)> gridToEquilGrid = [&simulator](
unsigned int i) {
687 return simulator.vanguard().gridIdxToEquilGridIdx(i);
690 std::function<void(
bool)> transUp =
691 [
this,gridToEquilGrid](
bool global) {
692 this->transmissibilities_.update(global,gridToEquilGrid);
695 OPM_TIMEBLOCK(applyActions);
696 actionHandler_.applyActions(episodeIdx,
697 simulator.time() + simulator.timeStepSize(),
701 if constexpr (enableMICP){
702 auto& model = this->model();
703 const auto& residual = this->model().linearizer().residual();
704 for (
unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
705 auto& phi = this->referencePorosity_[1][globalDofIdx];
706 MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
716 OPM_TIMEBLOCK(endEpisode);
717 auto& simulator = this->simulator();
718 auto& schedule = simulator.vanguard().schedule();
720 wellModel_.endEpisode();
721 aquiferModel_.endEpisode();
723 int episodeIdx = this->episodeIndex();
725 if (episodeIdx + 1 >=
static_cast<int>(schedule.size() - 1)) {
726 simulator.setFinished(
true);
731 simulator.startNextEpisode(schedule.stepLength(episodeIdx + 1));
740 OPM_TIMEBLOCK(problemWriteOutput);
743 if (Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() ||
744 this->simulator().episodeWillBeOver()) {
745 ParentType::writeOutput(verbose);
748 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
749 !this->simulator().episodeWillBeOver();
751 data::Solution localCellData = {};
755 if (enableDamarisOutput_) {
756 damarisWriter_->writeOutput(localCellData, isSubStep) ;
759 if (enableEclOutput_){
760 eclWriter_->writeOutput(std::move(localCellData), timer, isSubStep);
765 OPM_TIMEBLOCK(finalizeOutput);
775 template <
class Context>
778 unsigned timeIdx)
const
780 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
781 return transmissibilities_.permeability(globalSpaceIdx);
791 {
return transmissibilities_.permeability(globalElemIdx); }
796 template <
class Context>
798 [[maybe_unused]]
unsigned fromDofLocalIdx,
799 unsigned toDofLocalIdx)
const
801 assert(fromDofLocalIdx == 0);
802 return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
810 return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
816 template <
class Context>
818 [[maybe_unused]]
unsigned fromDofLocalIdx,
819 unsigned toDofLocalIdx)
const
821 assert(fromDofLocalIdx == 0);
822 return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
828 Scalar
diffusivity(
const unsigned globalCellIn,
const unsigned globalCellOut)
const{
829 return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
835 Scalar
dispersivity(
const unsigned globalCellIn,
const unsigned globalCellOut)
const{
836 return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
843 const unsigned boundaryFaceIdx)
const
845 return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
854 template <
class Context>
856 unsigned boundaryFaceIdx)
const
858 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
859 return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
866 const unsigned boundaryFaceIdx)
const
868 return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
876 const unsigned globalSpaceIdxOut)
const
878 return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
884 template <
class Context>
887 unsigned timeIdx)
const
889 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
890 unsigned toDofLocalIdx = face.exteriorIndex();
891 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
897 template <
class Context>
900 unsigned timeIdx)
const
902 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
903 unsigned toDofLocalIdx = face.exteriorIndex();
904 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
910 template <
class Context>
912 unsigned boundaryFaceIdx)
const
914 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
915 return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
922 {
return transmissibilities_; }
928 {
return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
931 {
return thresholdPressures_; }
934 {
return thresholdPressures_; }
937 {
return tracerModel_; }
940 {
return tracerModel_; }
950 template <
class Context>
951 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
953 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
954 return this->porosity(globalSpaceIdx, timeIdx);
963 template <
class Context>
964 Scalar
dofCenterDepth(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
966 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
967 return this->dofCenterDepth(globalSpaceIdx);
978 return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
984 template <
class Context>
987 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
988 return this->rockCompressibility(globalSpaceIdx);
994 template <
class Context>
997 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
998 return this->rockReferencePressure(globalSpaceIdx);
1004 template <
class Context>
1006 unsigned spaceIdx,
unsigned timeIdx)
const
1008 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1009 return this->materialLawParams(globalSpaceIdx);
1014 return materialLawManager_->materialLawParams(globalDofIdx);
1019 return materialLawManager_->materialLawParams(globalDofIdx, facedir);
1025 template <
class Context>
1026 const SolidEnergyLawParams&
1029 unsigned timeIdx)
const
1031 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1032 return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1038 template <
class Context>
1039 const ThermalConductionLawParams &
1042 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1043 return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1053 {
return materialLawManager_; }
1055 template <
class Flu
idState>
1057 std::array<Evaluation,numPhases> &mobility,
1058 DirectionalMobilityPtr &dirMob,
1059 FluidState &fluidState,
1060 unsigned globalSpaceIdx)
const
1062 OPM_TIMEBLOCK_LOCAL(updateRelperms);
1066 const auto& materialParams = materialLawParams(globalSpaceIdx);
1067 MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
1068 Valgrind::CheckDefined(mobility);
1070 if (materialLawManager_->hasDirectionalRelperms()
1071 || materialLawManager_->hasDirectionalImbnum())
1073 using Dir = FaceDir::DirEnum;
1074 constexpr int ndim = 3;
1075 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
1076 Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
1077 for (
int i = 0; i<ndim; i++) {
1078 const auto& materialParams = materialLawParams(globalSpaceIdx, facedirs[i]);
1079 auto& mob_array = dirMob->getArray(i);
1080 MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
1089 {
return materialLawManager_; }
1091 using BaseType::pvtRegionIndex;
1095 template <
class Context>
1096 unsigned pvtRegionIndex(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
1097 {
return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1099 using BaseType::satnumRegionIndex;
1103 template <
class Context>
1105 {
return this->satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1107 using BaseType::miscnumRegionIndex;
1111 template <
class Context>
1113 {
return this->miscnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1115 using BaseType::plmixnumRegionIndex;
1119 template <
class Context>
1121 {
return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1123 using BaseType::maxPolymerAdsorption;
1127 template <
class Context>
1129 {
return this->maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1135 {
return this->simulator().vanguard().caseName(); }
1140 template <
class Context>
1141 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
1145 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1146 return initialFluidStates_[globalDofIdx].temperature(0);
1154 return initialFluidStates_[globalDofIdx].temperature(0);
1157 const SolidEnergyLawParams&
1161 return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1163 const ThermalConductionLawParams &
1167 return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1175 template <
class Context>
1177 const Context& context,
1179 unsigned timeIdx)
const
1181 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
1182 if (!context.intersection(spaceIdx).boundary())
1185 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
1193 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1194 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1195 values.setThermalFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
1198 if (nonTrivialBoundaryConditions()) {
1199 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
1200 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1201 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1202 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
1203 const auto [type, massrate] = boundaryCondition(globalDofIdx, indexInInside);
1204 if (type == BCType::THERMAL)
1205 values.setThermalFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1206 else if (type == BCType::FREE || type == BCType::DIRICHLET)
1207 values.setFreeFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1208 else if (type == BCType::RATE)
1209 values.setMassRate(massrate, pvtRegionIdx);
1224 if (!this->vapparsActive(this->episodeIndex()))
1227 return this->maxOilSaturation_[globalDofIdx];
1241 if (!this->vapparsActive(this->episodeIndex()))
1244 this->maxOilSaturation_[globalDofIdx] = value;
1253 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
1254 this->episodeIndex(),
1255 this->pvtRegionIndex(globalDofIdx));
1264 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
1265 this->episodeIndex(),
1266 this->pvtRegionIndex(globalDofIdx));
1279 int episodeIdx = this->episodeIndex();
1280 return !this->mixControls_.drsdtActive(episodeIdx) &&
1281 !this->mixControls_.drvdtActive(episodeIdx) &&
1282 this->rockCompPoroMultWc_.empty() &&
1283 this->rockCompPoroMult_.empty();
1292 template <
class Context>
1293 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
1295 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1297 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
1298 values.assignNaive(initialFluidStates_[globalDofIdx]);
1300 SolventModule::assignPrimaryVars(values,
1301 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
1302 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
1304 if constexpr (enablePolymer)
1305 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
1307 if constexpr (enablePolymerMolarWeight)
1308 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
1310 if constexpr (enableBrine) {
1311 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
1312 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
1315 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
1319 if constexpr (enableMICP){
1320 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
1321 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
1322 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
1323 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
1324 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
1327 values.checkDefined();
1336 this->model().invalidateAndUpdateIntensiveQuantities(0);
1342 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
1343 this->model().invalidateAndUpdateIntensiveQuantities(1);
1354 thresholdPressures_.finishInit();
1356 updateCompositionChangeLimits_();
1358 aquiferModel_.initialSolutionApplied();
1360 if (this->simulator().episodeIndex() == 0) {
1361 eclWriter_->writeInitialFIPReport();
1370 template <
class Context>
1372 const Context& context,
1374 unsigned timeIdx)
const
1376 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1377 source(rate, globalDofIdx, timeIdx);
1381 unsigned globalDofIdx,
1382 unsigned timeIdx)
const
1384 OPM_TIMEBLOCK_LOCAL(eclProblemSource);
1388 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
1392 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1393 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
1395 Valgrind::CheckDefined(rate[eqIdx]);
1396 assert(isfinite(rate[eqIdx]));
1400 addToSourceDense(rate, globalDofIdx, timeIdx);
1404 unsigned globalDofIdx,
1405 unsigned timeIdx)
const
1407 aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
1410 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
1411 std::array<int,3> ijk;
1412 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
1414 if (source.hasSource(ijk)) {
1415 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1416 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
1417 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
1418 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
1420 for (
unsigned i = 0; i < phidx_map.size(); ++i) {
1421 const auto phaseIdx = phidx_map[i];
1422 const auto sourceComp = sc_map[i];
1423 const auto compIdx = cidx_map[i];
1424 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1427 Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1428 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1429 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
1431 rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate;
1434 if constexpr (enableSolvent) {
1435 Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
1436 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1437 const auto& solventPvt = SolventModule::solventPvt();
1438 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
1440 rate[Indices::contiSolventEqIdx] += mass_rate;
1442 if constexpr (enablePolymer) {
1443 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
1445 if constexpr (enableEnergy) {
1446 for (
unsigned i = 0; i < phidx_map.size(); ++i) {
1447 const auto phaseIdx = phidx_map[i];
1448 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1451 const auto sourceComp = sc_map[i];
1452 if (source.hasHrate({ijk, sourceComp})) {
1453 rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1455 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, 0);
1456 auto fs = intQuants.fluidState();
1458 if (source.hasTemperature({ijk, sourceComp})) {
1459 Scalar temperature = source.temperature({ijk, sourceComp});
1460 fs.setTemperature(temperature);
1462 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
1463 Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx);
1464 Scalar energy_rate = getValue(h)*mass_rate;
1465 rate[Indices::contiEnergyEqIdx] += energy_rate;
1474 const bool compensateDrift = wellModel_.wellsActive();
1475 if (enableDriftCompensation_ && compensateDrift) {
1476 const auto& simulator = this->simulator();
1477 const auto& model = this->model();
1482 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
1483 Scalar poro = this->porosity(globalDofIdx, timeIdx);
1484 Scalar dt = simulator.timeStepSize();
1485 EqVector dofDriftRate = drift_[globalDofIdx];
1486 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
1489 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1490 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
1491 if (cnv > maxCompensation) {
1492 dofDriftRate[eqIdx] *= maxCompensation/cnv;
1496 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
1497 rate[eqIdx] -= dofDriftRate[eqIdx];
1507 {
return wellModel_; }
1510 {
return wellModel_; }
1513 {
return aquiferModel_; }
1516 {
return aquiferModel_; }
1520 {
return initialFluidStates_[globalDofIdx]; }
1523 {
return eclWriter_->eclIO(); }
1526 {
return eclWriter_->setSubStepReport(report); }
1529 {
return eclWriter_->setSimulationReport(report); }
1532 {
return nonTrivialBoundaryConditions_; }
1536 OPM_TIMEBLOCK_LOCAL(boundaryFluidState);
1537 const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop;
1538 if (bcprop.size() > 0) {
1539 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1543 if (bcindex_(dir)[globalDofIdx] == 0)
1544 return initialFluidStates_[globalDofIdx];
1546 const auto& bc = bcprop[bcindex_(dir)[globalDofIdx]];
1547 if (bc.bctype == BCType::DIRICHLET )
1549 InitialFluidState fluidState;
1550 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1551 fluidState.setPvtRegionIndex(pvtRegionIdx);
1553 switch (bc.component) {
1554 case BCComponent::OIL:
1555 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
1556 throw std::logic_error(
"oil is not active and you're trying to add oil BC");
1558 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
1560 case BCComponent::GAS:
1561 if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1562 throw std::logic_error(
"gas is not active and you're trying to add gas BC");
1564 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
1566 case BCComponent::WATER:
1567 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1568 throw std::logic_error(
"water is not active and you're trying to add water BC");
1570 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
1572 case BCComponent::SOLVENT:
1573 case BCComponent::POLYMER:
1575 throw std::logic_error(
"you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
1578 double pressure = initialFluidStates_[globalDofIdx].pressure(refPressurePhaseIdx_());
1579 const auto pressure_input = bc.pressure;
1580 if (pressure_input) {
1581 pressure = *pressure_input;
1584 std::array<Scalar, numPhases> pc = {0};
1585 const auto& matParams = materialLawParams(globalDofIdx);
1586 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
1587 Valgrind::CheckDefined(pressure);
1588 Valgrind::CheckDefined(pc);
1589 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1590 if (!FluidSystem::phaseIsActive(phaseIdx))
1593 if (Indices::oilEnabled)
1594 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1595 else if (Indices::gasEnabled)
1596 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1597 else if (Indices::waterEnabled)
1599 fluidState.setPressure(phaseIdx, pressure);
1602 double temperature = initialFluidStates_[globalDofIdx].temperature(0);
1603 const auto temperature_input = bc.temperature;
1604 if(temperature_input)
1605 temperature = *temperature_input;
1606 fluidState.setTemperature(temperature);
1608 if (FluidSystem::enableDissolvedGas()) {
1609 fluidState.setRs(0.0);
1610 fluidState.setRv(0.0);
1612 if (FluidSystem::enableDissolvedGasInWater()) {
1613 fluidState.setRsw(0.0);
1615 if (FluidSystem::enableVaporizedWater())
1616 fluidState.setRvw(0.0);
1618 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1619 if (!FluidSystem::phaseIsActive(phaseIdx))
1622 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
1623 fluidState.setInvB(phaseIdx, b);
1625 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
1626 fluidState.setDensity(phaseIdx, rho);
1628 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
1629 fluidState.setEnthalpy(phaseIdx, h);
1632 fluidState.checkDefined();
1636 return initialFluidStates_[globalDofIdx];
1647 OPM_TIMEBLOCK(nexTimeStepSize);
1649 if (this->nextTimeStepSize_ > 0.0)
1650 return this->nextTimeStepSize_;
1652 const auto& simulator = this->simulator();
1653 int episodeIdx = simulator.episodeIndex();
1657 return this->initialTimeStepSize_;
1662 const auto& newtonMethod = this->model().newtonMethod();
1663 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1671 template <
class LhsEval>
1674 OPM_TIMEBLOCK_LOCAL(rockCompPoroMultiplier);
1675 if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
1678 unsigned tableIdx = 0;
1679 if (!this->rockTableIdx_.empty())
1680 tableIdx = this->rockTableIdx_[elementIdx];
1682 const auto& fs = intQuants.fluidState();
1683 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1684 if (!this->minRefPressure_.empty())
1687 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1688 this->minRefPressure_[elementIdx]);
1690 if (!this->overburdenPressure_.empty())
1691 effectivePressure -= this->overburdenPressure_[elementIdx];
1694 if (!this->rockCompPoroMult_.empty()) {
1695 return this->rockCompPoroMult_[tableIdx].eval(effectivePressure,
true);
1699 assert(!this->rockCompPoroMultWc_.empty());
1700 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1701 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
1703 return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax,
true);
1711 template <
class LhsEval>
1714 const bool implicit = !this->explicitRockCompaction_;
1715 return implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<LhsEval>(intQuants, elementIdx)
1716 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1724 template <
class LhsEval>
1727 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
1728 if (!enableSaltPrecipitation)
1731 const auto& fs = intQuants.fluidState();
1732 unsigned tableIdx = fs.pvtRegionIndex();
1733 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
1734 porosityFactor = min(porosityFactor, 1.0);
1735 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
1736 return permfactTable.eval(porosityFactor,
true);
1742 template <
class LhsEval>
1745 OPM_TIMEBLOCK_LOCAL(wellTransMultiplier);
1747 const bool implicit = !this->explicitRockCompaction_;
1748 double trans_mult = implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<double>(intQuants, elementIdx)
1749 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1750 trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants);
1755 std::pair<BCType, RateVector>
boundaryCondition(
const unsigned int globalSpaceIdx,
const int directionId)
const
1757 OPM_TIMEBLOCK_LOCAL(boundaryCondition);
1758 if (!nonTrivialBoundaryConditions_) {
1761 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1762 const auto& schedule = this->simulator().vanguard().schedule();
1763 if (bcindex_(dir)[globalSpaceIdx] == 0) {
1766 if (schedule[this->episodeIndex()].bcprop.size() == 0) {
1769 const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
1770 if (bc.bctype!=BCType::RATE) {
1771 return { bc.bctype, RateVector(0.0) };
1774 RateVector rate = 0.0;
1775 switch (bc.component) {
1776 case BCComponent::OIL:
1777 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = bc.rate;
1779 case BCComponent::GAS:
1780 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = bc.rate;
1782 case BCComponent::WATER:
1783 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate;
1785 case BCComponent::SOLVENT:
1786 if constexpr (!enableSolvent)
1787 throw std::logic_error(
"solvent is disabled and you're trying to add solvent to BC");
1789 rate[Indices::solventSaturationIdx] = bc.rate;
1791 case BCComponent::POLYMER:
1792 if constexpr (!enablePolymer)
1793 throw std::logic_error(
"polymer is disabled and you're trying to add polymer to BC");
1795 rate[Indices::polymerConcentrationIdx] = bc.rate;
1798 throw std::logic_error(
"you need to specify the component when RATE type is set in BC");
1802 return {bc.bctype, rate};
1812 eclWriter_->mutableOutputModule().setCnvData(data);
1815 template<
class Serializer>
1818 serializer(
static_cast<BaseType&
>(*
this));
1820 serializer(wellModel_);
1821 serializer(aquiferModel_);
1822 serializer(tracerModel_);
1823 serializer(*materialLawManager_);
1824 serializer(*eclWriter_);
1827 Implementation& asImp_()
1828 {
return *
static_cast<Implementation *
>(
this); }
1832 OPM_TIMEBLOCK(updateExplicitQuantities);
1833 const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
1834 const bool invalidateFromMinPressure = updateMinPressure_();
1837 const bool invalidateFromHyst = updateHysteresis_();
1838 const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
1841 bool invalidateIntensiveQuantities
1842 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
1843 if (invalidateIntensiveQuantities) {
1844 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1845 this->model().invalidateAndUpdateIntensiveQuantities(0);
1848 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
1849 updateMaxPolymerAdsorption_();
1851 updateRockCompTransMultVal_();
1854 template<
class UpdateFunc>
1858 OPM_TIMEBLOCK(updateProperty);
1859 const auto& model = this->simulator().model();
1860 const auto& primaryVars = model.solution(0);
1861 const auto& vanguard = this->simulator().vanguard();
1862 std::size_t numGridDof = primaryVars.size();
1865#pragma omp parallel for
1867 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1868 const auto& iq = *model.cachedIntensiveQuantities(dofIdx, 0);
1877 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1880 int episodeIdx = this->episodeIndex();
1881 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1882 this->mixControls_.drsdtActive(episodeIdx),
1883 this->mixControls_.drvdtActive(episodeIdx)};
1884 if (!active[0] && !active[1] && !active[2]) {
1888 this->updateProperty_(
"FlowProblem::updateCompositionChangeLimits_()) failed:",
1889 [
this,episodeIdx,active](
unsigned compressedDofIdx,
1890 const IntensiveQuantities& iq)
1892 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1893 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1894 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1895 this->mixControls_.update(compressedDofIdx,
1898 this->gravity_[dim - 1],
1899 perm[dim - 1][dim - 1],
1909 OPM_TIMEBLOCK(updateMaxOilSaturation);
1910 int episodeIdx = this->episodeIndex();
1913 if (this->vapparsActive(episodeIdx)) {
1914 this->updateProperty_(
"FlowProblem::updateMaxOilSaturation_() failed:",
1915 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1917 this->updateMaxOilSaturation_(compressedDofIdx,iq);
1927 OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
1928 const auto& fs = iq.fluidState();
1929 const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1930 auto& mos = this->maxOilSaturation_;
1931 if(mos[compressedDofIdx] < So){
1932 mos[compressedDofIdx] = So;
1941 OPM_TIMEBLOCK(updateMaxWaterSaturation);
1943 if (this->maxWaterSaturation_.empty())
1946 this->maxWaterSaturation_[1] = this->maxWaterSaturation_[0];
1947 this->updateProperty_(
"FlowProblem::updateMaxWaterSaturation_() failed:",
1948 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1950 this->updateMaxWaterSaturation_(compressedDofIdx,iq);
1958 OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation);
1959 const auto& fs = iq.fluidState();
1960 const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
1961 auto& mow = this->maxWaterSaturation_;
1962 if(mow[compressedDofIdx]< Sw){
1963 mow[compressedDofIdx] = Sw;
1972 OPM_TIMEBLOCK(updateMinPressure);
1974 if (this->minRefPressure_.empty())
1977 this->updateProperty_(
"FlowProblem::updateMinPressure_() failed:",
1978 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1980 this->updateMinPressure_(compressedDofIdx,iq);
1986 OPM_TIMEBLOCK_LOCAL(updateMinPressure);
1987 const auto& fs = iq.fluidState();
1988 const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
1989 auto& min_pressures = this->minRefPressure_;
1990 if(min_pressures[compressedDofIdx]> min_pressure){
1991 min_pressures[compressedDofIdx] = min_pressure;
2002 std::function<std::vector<double>(
const FieldPropsManager&,
const std::string&)>
2005 const auto& lookup = this->lookUpData_;
2006 return [&lookup](
const FieldPropsManager& fieldPropManager,
const std::string& propString)
2008 return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
2016 template<
typename IntType>
2017 std::function<std::vector<IntType>(
const FieldPropsManager&,
const std::string&,
bool)>
2020 const auto& lookup = this->lookUpData_;
2021 return [&lookup](
const FieldPropsManager& fieldPropManager,
const std::string& propString,
bool needsTranslation)
2023 return lookup.template assignFieldPropsIntOnLeaf<IntType>(fieldPropManager, propString, needsTranslation);
2029 OPM_TIMEBLOCK(readMaterialParameters);
2030 const auto& simulator = this->simulator();
2031 const auto& vanguard = simulator.vanguard();
2032 const auto& eclState = vanguard.eclState();
2036 this->updatePvtnum_();
2037 this->updateSatnum_();
2040 this->updateMiscnum_();
2042 this->updatePlmixnum_();
2047 updateReferencePorosity_();
2048 this->referencePorosity_[1] = this->referencePorosity_[0];
2053 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
2054 materialLawManager_->initFromState(eclState);
2055 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2056 this->
template fieldPropIntTypeOnLeafAssigner_<int>(),
2057 this-> lookupIdxOnLevelZeroAssigner_());
2063 if constexpr (enableEnergy)
2065 const auto& simulator = this->simulator();
2066 const auto& vanguard = simulator.vanguard();
2067 const auto& eclState = vanguard.eclState();
2070 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
2071 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2072 this-> fieldPropDoubleOnLeafAssigner_(),
2073 this->
template fieldPropIntTypeOnLeafAssigner_<unsigned int>());
2079 const auto& simulator = this->simulator();
2080 const auto& vanguard = simulator.vanguard();
2081 const auto& eclState = vanguard.eclState();
2083 std::size_t numDof = this->model().numGridDof();
2085 this->referencePorosity_[0].resize(numDof);
2087 const auto& fp = eclState.fieldProps();
2088 const std::vector<double> porvData =
this -> fieldPropDoubleOnLeafAssigner_()(fp,
"PORV");
2089 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2090 Scalar poreVolume = porvData[dofIdx];
2095 Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx);
2096 assert(dofVolume > 0.0);
2097 this->referencePorosity_[0][dofIdx] = poreVolume/dofVolume;
2103 const auto& simulator = this->simulator();
2104 const auto& vanguard = simulator.vanguard();
2105 const auto& eclState = vanguard.eclState();
2107 if (eclState.getInitConfig().hasEquil())
2108 readEquilInitialCondition_();
2110 readExplicitInitialCondition_();
2112 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
2113 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
2116 enablePolymerMolarWeight,
2120 std::size_t numElems = this->model().numGridDof();
2121 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2122 const auto& fs = initialFluidStates_[elemIdx];
2123 if (!this->maxWaterSaturation_.empty())
2124 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
2125 if (!this->maxOilSaturation_.empty())
2126 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
2127 if (!this->minRefPressure_.empty())
2128 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
2136 const auto& simulator = this->simulator();
2141 std::size_t numElems = this->model().numGridDof();
2142 initialFluidStates_.resize(numElems);
2143 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2144 auto& elemFluidState = initialFluidStates_[elemIdx];
2152 if(this->simulator().vanguard().grid().maxLevel() > 0) {
2153 throw std::invalid_argument(
"Refined grids are not yet supported for restart ");
2157 auto& simulator = this->simulator();
2158 const auto& schedule = simulator.vanguard().schedule();
2159 const auto& eclState = simulator.vanguard().eclState();
2160 const auto& initconfig = eclState.getInitConfig();
2162 int restart_step = initconfig.getRestartStep();
2164 simulator.setTime(schedule.seconds(restart_step));
2166 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
2167 schedule.stepLength(restart_step));
2168 simulator.setEpisodeIndex(restart_step);
2170 eclWriter_->beginRestart();
2172 Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength());
2173 simulator.setTimeStepSize(dt);
2175 std::size_t numElems = this->model().numGridDof();
2176 initialFluidStates_.resize(numElems);
2177 if constexpr (enableSolvent) {
2178 this->solventSaturation_.resize(numElems, 0.0);
2179 this->solventRsw_.resize(numElems, 0.0);
2182 if constexpr (enablePolymer)
2183 this->polymer_.concentration.resize(numElems, 0.0);
2185 if constexpr (enablePolymerMolarWeight) {
2186 const std::string msg {
"Support of the RESTART for polymer molecular weight "
2187 "is not implemented yet. The polymer weight value will be "
2188 "zero when RESTART begins"};
2189 OpmLog::warning(
"NO_POLYMW_RESTART", msg);
2190 this->polymer_.moleWeight.resize(numElems, 0.0);
2193 if constexpr (enableMICP) {
2194 this->micp_.resize(numElems);
2197 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2198 auto& elemFluidState = initialFluidStates_[elemIdx];
2199 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
2200 eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
2201 eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
2210 auto ssol = enableSolvent
2211 ? eclWriter_->outputModule().getSolventSaturation(elemIdx)
2214 processRestartSaturations_(elemFluidState, ssol);
2216 if constexpr (enableSolvent) {
2217 this->solventSaturation_[elemIdx] = ssol;
2218 this->solventRsw_[elemIdx] = eclWriter_->outputModule().getSolventRsw(elemIdx);
2223 bool isThermal = eclState.getSimulationConfig().isThermal();
2224 bool needTemperature = (eclState.runspec().co2Storage() || eclState.runspec().h2Storage());
2225 if (!isThermal && needTemperature) {
2226 const auto& fp = simulator.vanguard().eclState().fieldProps();
2227 elemFluidState.setTemperature(fp.get_double(
"TEMPI")[elemIdx]);
2230 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
2232 if constexpr (enablePolymer)
2233 this->polymer_.concentration[elemIdx] = eclWriter_->outputModule().getPolymerConcentration(elemIdx);
2234 if constexpr (enableMICP){
2235 this->micp_.microbialConcentration[elemIdx] = eclWriter_->outputModule().getMicrobialConcentration(elemIdx);
2236 this->micp_.oxygenConcentration[elemIdx] = eclWriter_->outputModule().getOxygenConcentration(elemIdx);
2237 this->micp_.ureaConcentration[elemIdx] = eclWriter_->outputModule().getUreaConcentration(elemIdx);
2238 this->micp_.biofilmConcentration[elemIdx] = eclWriter_->outputModule().getBiofilmConcentration(elemIdx);
2239 this->micp_.calciteConcentration[elemIdx] = eclWriter_->outputModule().getCalciteConcentration(elemIdx);
2244 const int episodeIdx = this->episodeIndex();
2245 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
2250 auto& sol = this->model().solution(0);
2251 const auto& gridView = this->gridView();
2252 ElementContext elemCtx(simulator);
2253 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2254 elemCtx.updatePrimaryStencil(elem);
2255 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
2256 initial(sol[elemIdx], elemCtx, 0, 0);
2264 this->model().syncOverlap();
2266 eclWriter_->endRestart();
2273 const Scalar smallSaturationTolerance = 1.e-6;
2274 Scalar sumSaturation = 0.0;
2275 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2276 if (FluidSystem::phaseIsActive(phaseIdx)) {
2277 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
2278 elemFluidState.setSaturation(phaseIdx, 0.0);
2280 sumSaturation += elemFluidState.saturation(phaseIdx);
2284 if constexpr (enableSolvent) {
2285 if (solventSaturation < smallSaturationTolerance)
2286 solventSaturation = 0.0;
2288 sumSaturation += solventSaturation;
2291 assert(sumSaturation > 0.0);
2293 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2294 if (FluidSystem::phaseIsActive(phaseIdx)) {
2295 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
2296 elemFluidState.setSaturation(phaseIdx, saturation);
2299 if constexpr (enableSolvent) {
2300 solventSaturation = solventSaturation / sumSaturation;
2306 const auto& simulator = this->simulator();
2307 const auto& vanguard = simulator.vanguard();
2308 const auto& eclState = vanguard.eclState();
2309 const auto& fp = eclState.fieldProps();
2310 bool has_swat = fp.has_double(
"SWAT");
2311 bool has_sgas = fp.has_double(
"SGAS");
2312 bool has_rs = fp.has_double(
"RS");
2313 bool has_rv = fp.has_double(
"RV");
2314 bool has_rvw = fp.has_double(
"RVW");
2315 bool has_pressure = fp.has_double(
"PRESSURE");
2316 bool has_salt = fp.has_double(
"SALT");
2317 bool has_saltp = fp.has_double(
"SALTP");
2320 if (Indices::numPhases > 1) {
2321 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
2322 throw std::runtime_error(
"The ECL input file requires the presence of the SWAT keyword if "
2323 "the water phase is active");
2324 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
2325 throw std::runtime_error(
"The ECL input file requires the presence of the SGAS keyword if "
2326 "the gas phase is active");
2329 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
2330 "keyword if the model is initialized explicitly");
2331 if (FluidSystem::enableDissolvedGas() && !has_rs)
2332 throw std::runtime_error(
"The ECL input file requires the RS keyword to be present if"
2333 " dissolved gas is enabled");
2334 if (FluidSystem::enableVaporizedOil() && !has_rv)
2335 throw std::runtime_error(
"The ECL input file requires the RV keyword to be present if"
2336 " vaporized oil is enabled");
2337 if (FluidSystem::enableVaporizedWater() && !has_rvw)
2338 throw std::runtime_error(
"The ECL input file requires the RVW keyword to be present if"
2339 " vaporized water is enabled");
2340 if (enableBrine && !has_salt)
2341 throw std::runtime_error(
"The ECL input file requires the SALT keyword to be present if"
2342 " brine is enabled and the model is initialized explicitly");
2343 if (enableSaltPrecipitation && !has_saltp)
2344 throw std::runtime_error(
"The ECL input file requires the SALTP keyword to be present if"
2345 " salt precipitation is enabled and the model is initialized explicitly");
2347 std::size_t numDof = this->model().numGridDof();
2349 initialFluidStates_.resize(numDof);
2351 std::vector<double> waterSaturationData;
2352 std::vector<double> gasSaturationData;
2353 std::vector<double> pressureData;
2354 std::vector<double> rsData;
2355 std::vector<double> rvData;
2356 std::vector<double> rvwData;
2357 std::vector<double> tempiData;
2358 std::vector<double> saltData;
2359 std::vector<double> saltpData;
2361 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
2362 waterSaturationData = fp.get_double(
"SWAT");
2364 waterSaturationData.resize(numDof);
2366 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
2367 gasSaturationData = fp.get_double(
"SGAS");
2369 gasSaturationData.resize(numDof);
2371 pressureData = fp.get_double(
"PRESSURE");
2372 if (FluidSystem::enableDissolvedGas())
2373 rsData = fp.get_double(
"RS");
2375 if (FluidSystem::enableVaporizedOil())
2376 rvData = fp.get_double(
"RV");
2378 if (FluidSystem::enableVaporizedWater())
2379 rvwData = fp.get_double(
"RVW");
2382 tempiData = fp.get_double(
"TEMPI");
2385 if constexpr (enableBrine)
2386 saltData = fp.get_double(
"SALT");
2389 if constexpr (enableSaltPrecipitation)
2390 saltpData = fp.get_double(
"SALTP");
2393 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2394 auto& dofFluidState = initialFluidStates_[dofIdx];
2396 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
2401 Scalar temperatureLoc = tempiData[dofIdx];
2402 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
2403 temperatureLoc = FluidSystem::surfaceTemperature;
2404 dofFluidState.setTemperature(temperatureLoc);
2409 if constexpr (enableBrine)
2410 dofFluidState.setSaltConcentration(saltData[dofIdx]);
2415 if constexpr (enableSaltPrecipitation)
2416 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
2421 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
2422 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
2423 waterSaturationData[dofIdx]);
2425 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
2426 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
2427 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2429 - waterSaturationData[dofIdx]);
2432 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2433 gasSaturationData[dofIdx]);
2435 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
2436 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
2438 - waterSaturationData[dofIdx]
2439 - gasSaturationData[dofIdx]);
2444 Scalar pressure = pressureData[dofIdx];
2448 std::array<Scalar, numPhases> pc = {0};
2449 const auto& matParams = materialLawParams(dofIdx);
2450 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
2451 Valgrind::CheckDefined(pressure);
2452 Valgrind::CheckDefined(pc);
2453 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2454 if (!FluidSystem::phaseIsActive(phaseIdx))
2457 if (Indices::oilEnabled)
2458 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
2459 else if (Indices::gasEnabled)
2460 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
2461 else if (Indices::waterEnabled)
2463 dofFluidState.setPressure(phaseIdx, pressure);
2466 if (FluidSystem::enableDissolvedGas())
2467 dofFluidState.setRs(rsData[dofIdx]);
2468 else if (Indices::gasEnabled && Indices::oilEnabled)
2469 dofFluidState.setRs(0.0);
2471 if (FluidSystem::enableVaporizedOil())
2472 dofFluidState.setRv(rvData[dofIdx]);
2473 else if (Indices::gasEnabled && Indices::oilEnabled)
2474 dofFluidState.setRv(0.0);
2476 if (FluidSystem::enableVaporizedWater())
2477 dofFluidState.setRvw(rvwData[dofIdx]);
2482 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2483 if (!FluidSystem::phaseIsActive(phaseIdx))
2486 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2487 dofFluidState.setInvB(phaseIdx, b);
2489 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2490 dofFluidState.setDensity(phaseIdx, rho);
2499 if (!materialLawManager_->enableHysteresis())
2504 this->updateProperty_(
"FlowProblem::updateHysteresis_() failed:",
2505 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2507 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2515 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
2516 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2524 this->updateProperty_(
"FlowProblem::updateMaxPolymerAdsorption_() failed:",
2525 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2527 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
2533 const Scalar pa = scalarValue(iq.polymerAdsorption());
2534 auto& mpa = this->polymer_.maxAdsorption;
2535 if (mpa[compressedDofIdx] < pa) {
2536 mpa[compressedDofIdx] = pa;
2545 if (this->rockCompTransMultVal_.empty())
2548 return this->rockCompTransMultVal_[dofIdx];
2555 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
2556 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
2557 ConditionalStorage<enableDiffusion, Scalar> diffusivity;
2558 ConditionalStorage<enableDispersion, Scalar> dispersivity;
2559 Scalar transmissibility;
2563 void updatePffDofData_()
2565 const auto& distFn =
2566 [
this](PffDofData_& dofData,
2567 const Stencil& stencil,
2568 unsigned localDofIdx)
2571 const auto& elementMapper = this->model().elementMapper();
2573 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
2574 if (localDofIdx != 0) {
2575 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(0));
2576 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
2578 if constexpr (enableEnergy) {
2579 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
2580 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
2582 if constexpr (enableDiffusion)
2583 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
2584 if (enableDispersion)
2585 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
2589 pffDofData_.update(distFn);
2592 void readBoundaryConditions_()
2594 const auto& simulator = this->simulator();
2595 const auto& vanguard = simulator.vanguard();
2596 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
2597 if (bcconfig.size() > 0) {
2598 nonTrivialBoundaryConditions_ =
true;
2600 std::size_t numCartDof = vanguard.cartesianSize();
2601 unsigned numElems = vanguard.gridView().size(0);
2602 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
2604 for (
unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
2605 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
2607 bcindex_.resize(numElems, 0);
2608 auto loopAndApply = [&cartesianToCompressedElemIdx,
2609 &vanguard](
const auto& bcface,
2612 for (
int i = bcface.i1; i <= bcface.i2; ++i) {
2613 for (
int j = bcface.j1; j <= bcface.j2; ++j) {
2614 for (
int k = bcface.k1; k <= bcface.k2; ++k) {
2615 std::array<int, 3> tmp = {i,j,k};
2616 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
2623 for (
const auto& bcface : bcconfig) {
2624 std::vector<int>& data = bcindex_(bcface.dir);
2625 const int index = bcface.index;
2626 loopAndApply(bcface,
2627 [&data,index](
int elemIdx)
2628 { data[elemIdx] = index; });
2635 Scalar limitNextTimeStepSize_(Scalar dtNext)
const
2637 if constexpr (enableExperiments) {
2638 const auto& simulator = this->simulator();
2639 const auto& schedule = simulator.vanguard().schedule();
2640 int episodeIdx = simulator.episodeIndex();
2643 Scalar maxTimeStepSize = Parameters::get<TypeTag, Properties::SolverMaxTimeStepInDays>() * 24 * 60 * 60;
2644 int reportStepIdx = std::max(episodeIdx, 0);
2645 if (this->enableTuning_) {
2646 const auto& tuning = schedule[reportStepIdx].tuning();
2647 maxTimeStepSize = tuning.TSMAXZ;
2650 dtNext = std::min(dtNext, maxTimeStepSize);
2652 Scalar remainingEpisodeTime =
2653 simulator.episodeStartTime() + simulator.episodeLength()
2654 - (simulator.startTime() + simulator.time());
2655 assert(remainingEpisodeTime >= 0.0);
2659 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
2662 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
2664 if (simulator.episodeStarts()) {
2667 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
2668 bool wellEventOccured =
2669 events.hasEvent(ScheduleEvents::NEW_WELL)
2670 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
2671 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
2672 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
2673 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
2674 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
2681 void computeAndSetEqWeights_()
2683 std::vector<Scalar> sumInvB(numPhases, 0.0);
2684 const auto& gridView = this->gridView();
2685 ElementContext elemCtx(this->simulator());
2686 for(
const auto& elem: elements(gridView, Dune::Partitions::interior)) {
2687 elemCtx.updatePrimaryStencil(elem);
2688 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
2689 const auto& dofFluidState = initialFluidStates_[elemIdx];
2690 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2691 if (!FluidSystem::phaseIsActive(phaseIdx))
2694 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
2698 std::size_t numDof = this->model().numGridDof();
2699 const auto& comm = this->simulator().vanguard().grid().comm();
2700 comm.sum(sumInvB.data(),sumInvB.size());
2701 Scalar numTotalDof = comm.sum(numDof);
2703 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2704 if (!FluidSystem::phaseIsActive(phaseIdx))
2707 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
2708 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
2709 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
2710 this->model().setEqWeight(activeSolventCompIdx, avgB);
2714 int refPressurePhaseIdx_()
const {
2715 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2718 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2722 return waterPhaseIdx;
2726 void updateRockCompTransMultVal_()
2728 const auto& model = this->simulator().model();
2729 std::size_t numGridDof = this->model().numGridDof();
2730 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
2731 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
2732 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, 0);
2733 Scalar trans_mult = computeRockCompTransMultiplier_<Scalar>(iq, elementIdx);
2734 this->rockCompTransMultVal_[elementIdx] = trans_mult;
2743 template <
class LhsEval>
2744 LhsEval computeRockCompTransMultiplier_(
const IntensiveQuantities& intQuants,
unsigned elementIdx)
const
2746 OPM_TIMEBLOCK_LOCAL(computeRockCompTransMultiplier);
2747 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
2750 unsigned tableIdx = 0;
2751 if (!this->rockTableIdx_.empty())
2752 tableIdx = this->rockTableIdx_[elementIdx];
2754 const auto& fs = intQuants.fluidState();
2755 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
2757 if (!this->minRefPressure_.empty())
2760 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
2761 this->minRefPressure_[elementIdx]);
2763 if (!this->overburdenPressure_.empty())
2764 effectivePressure -= this->overburdenPressure_[elementIdx];
2766 if (!this->rockCompTransMult_.empty())
2767 return this->rockCompTransMult_[tableIdx].eval(effectivePressure,
true);
2770 assert(!this->rockCompTransMultWc_.empty());
2771 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
2772 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
2774 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax,
true);
2777 typename Vanguard::TransmissibilityType transmissibilities_;
2779 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
2780 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
2782 FlowThresholdPressure<TypeTag> thresholdPressures_;
2784 std::vector<InitialFluidState> initialFluidStates_;
2786 bool enableDriftCompensation_;
2787 GlobalEqVector drift_;
2789 WellModel wellModel_;
2790 AquiferModel aquiferModel_;
2792 bool enableEclOutput_;
2793 std::unique_ptr<EclWriterType> eclWriter_;
2796 bool enableDamarisOutput_ = false ;
2797 std::unique_ptr<DamarisWriterType> damarisWriter_;
2800 PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
2801 TracerModel tracerModel_;
2803 ActionHandler actionHandler_;
2808 std::array<std::vector<T>,6> data;
2810 void resize(std::size_t size, T defVal)
2812 for (
auto& d : data)
2813 d.resize(size, defVal);
2816 const std::vector<T>& operator()(FaceDir::DirEnum dir)
const
2818 if (dir == FaceDir::DirEnum::Unknown)
2819 throw std::runtime_error(
"Tried to access BC data for the 'Unknown' direction");
2821 int div =
static_cast<int>(dir);
2822 while ((div /= 2) >= 1)
2824 assert(idx >= 0 && idx <= 5);
2828 std::vector<T>& operator()(FaceDir::DirEnum dir)
2830 return const_cast<std::vector<T>&
>(std::as_const(*
this)(dir));
2834 BCData<int> bcindex_;
2835 bool nonTrivialBoundaryConditions_ =
false;
2836 bool explicitRockCompaction_ =
false;
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:172
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:138
This file contains the flux module which is used for flow problems.
Collects necessary output values and pass them to Damaris server processes.
Definition: DamarisWriter.hpp:88
Collects necessary output values and pass it to opm-output.
Definition: EclWriter.hpp:100
static void registerParameters()
Definition: EclWriter.hpp:123
Computes the initial condition based on the EQUIL keyword from ECL.
Definition: EquilInitializer.hpp:58
const ScalarFluidState & initialFluidState(unsigned elemIdx) const
Return the initial thermodynamic state which should be used as the initial condition.
Definition: EquilInitializer.hpp:195
BlackOilFluidState< Scalar, FluidSystem, enableTemperature, enableEnergy, enableDissolution, enableVapwat, enableBrine, enableSaltPrecipitation, has_disgas_in_water, Indices::numPhases > ScalarFluidState
Definition: EquilInitializer.hpp:98
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowGenericProblem.hpp:70
PolymerSolutionContainer< Scalar > polymer_
Definition: FlowGenericProblem.hpp:358
Scalar rockReferencePressure(unsigned globalSpaceIdx) const
Definition: FlowGenericProblem_impl.hpp:324
static std::string briefDescription()
Definition: FlowGenericProblem_impl.hpp:133
Scalar porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
Direct indexed access to the porosity.
Definition: FlowGenericProblem_impl.hpp:339
Scalar rockCompressibility(unsigned globalSpaceIdx) const
Definition: FlowGenericProblem_impl.hpp:309
Scalar initialTimeStepSize_
Definition: FlowGenericProblem.hpp:371
bool enableTuning_
Definition: FlowGenericProblem.hpp:370
MixingRateControls< GetPropType< TypeTag, Properties::FluidSystem > > mixControls_
Definition: FlowGenericProblem.hpp:367
void readRockParameters_(const std::vector< Scalar > &cellCenterDepths, std::function< std::array< int, 3 >(const unsigned)> ijkIndex)
Definition: FlowGenericProblem_impl.hpp:140
std::vector< Scalar > maxOilSaturation_
Definition: FlowGenericProblem.hpp:359
int numPressurePointsEquil_
Definition: FlowGenericProblem.hpp:375
void initFluidSystem_()
Definition: FlowGenericProblem_impl.hpp:485
Scalar maxTimeStepAfterWellEvent_
Definition: FlowGenericProblem.hpp:372
bool shouldWriteOutput() const
Always returns true. The ecl output writer takes care of the rest.
Definition: FlowGenericProblem.hpp:301
static std::string helpPreamble(int, const char **argv)
Definition: FlowGenericProblem_impl.hpp:118
bool shouldWriteRestartFile() const
Returns true if an eWoms restart file should be written to disk.
Definition: FlowGenericProblem.hpp:310
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowProblem.hpp:112
const WellModel & wellModel() const
Returns a reference to the ECL well manager used by the problem.
Definition: FlowProblem.hpp:1506
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
Direct access to the transmissibility between two elements.
Definition: FlowProblem.hpp:808
Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
Definition: FlowProblem.hpp:927
void readEquilInitialCondition_()
Definition: FlowProblem.hpp:2134
void readExplicitInitialCondition_()
Definition: FlowProblem.hpp:2304
bool nonTrivialBoundaryConditions() const
Definition: FlowProblem.hpp:1531
std::function< std::vector< IntType >(const FieldPropsManager &, const std::string &, bool)> fieldPropIntTypeOnLeafAssigner_()
Definition: FlowProblem.hpp:2018
const DimMatrix & intrinsicPermeability(unsigned globalElemIdx) const
This method returns the intrinsic permeability tensor given a global element index.
Definition: FlowProblem.hpp:790
void endEpisode()
Called by the simulator after the end of an episode.
Definition: FlowProblem.hpp:714
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition: FlowProblem.hpp:1096
Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
Returns the maximum value of the oil vaporization factor at the current time for a given degree of fr...
Definition: FlowProblem.hpp:1262
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:951
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition: FlowProblem.hpp:612
const std::unique_ptr< EclWriterType > & eclWriter() const
Definition: FlowProblem.hpp:1805
void setSubStepReport(const SimulatorReportSingle &report)
Definition: FlowProblem.hpp:1525
bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities &iq)
Definition: FlowProblem.hpp:2513
unsigned satnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition: FlowProblem.hpp:1104
Scalar thermalHalfTransmissibilityOut(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:898
const ThermalConductionLawParams & thermalConductionLawParams(unsigned globalSpaceIdx, unsigned) const
Definition: FlowProblem.hpp:1164
bool updateMinPressure_()
Definition: FlowProblem.hpp:1970
std::function< std::vector< double >(const FieldPropsManager &, const std::string &)> fieldPropDoubleOnLeafAssigner_()
Definition: FlowProblem.hpp:2003
Scalar transmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition: FlowProblem.hpp:855
void updateReferencePorosity_()
Definition: FlowProblem.hpp:2077
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn, const unsigned globalSpaceIdxOut) const
Definition: FlowProblem.hpp:875
void endTimeStep()
Called by the simulator after each time integration.
Definition: FlowProblem.hpp:632
Scalar rockCompressibility(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:985
bool updateMaxWaterSaturation_()
Definition: FlowProblem.hpp:1939
bool recycleFirstIterationStorage() const
Return if the storage term of the first iteration is identical to the storage term for the solution o...
Definition: FlowProblem.hpp:1277
const InitialFluidState & initialFluidState(unsigned globalDofIdx) const
Definition: FlowProblem.hpp:1519
Scalar maxOilSaturation(unsigned globalDofIdx) const
Returns an element's historic maximum oil phase saturation that was observed during the simulation.
Definition: FlowProblem.hpp:1222
bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities &iq)
Definition: FlowProblem.hpp:1985
std::string name() const
Definition: FlowProblem.hpp:1134
int episodeIndex() const
Definition: FlowProblem.hpp:505
void readInitialCondition_()
Definition: FlowProblem.hpp:2101
LhsEval rockCompTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition: FlowProblem.hpp:1712
void endIteration()
Called by the simulator after each Newton-Raphson iteration.
Definition: FlowProblem.hpp:622
FlowProblem(Simulator &simulator)
Definition: FlowProblem.hpp:270
const Vanguard::TransmissibilityType & eclTransmissibilities() const
Return a reference to the object that handles the "raw" transmissibilities.
Definition: FlowProblem.hpp:921
void addToSourceDense(RateVector &rate, unsigned globalDofIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1403
void source(RateVector &rate, unsigned globalDofIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1380
Scalar nextTimeStepSize() const
Propose the size of the next time step to the simulator.
Definition: FlowProblem.hpp:1645
std::pair< BCType, RateVector > boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) const
Definition: FlowProblem.hpp:1755
void writeOutput(const SimulatorTimer &timer, bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition: FlowProblem.hpp:738
void updateExplicitQuantities_()
Definition: FlowProblem.hpp:1830
void updateRelperms(std::array< Evaluation, numPhases > &mobility, DirectionalMobilityPtr &dirMob, FluidState &fluidState, unsigned globalSpaceIdx) const
Definition: FlowProblem.hpp:1056
const FlowThresholdPressure< TypeTag > & thresholdPressure() const
Definition: FlowProblem.hpp:930
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1040
Scalar transmissibility(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition: FlowProblem.hpp:797
const TracerModel & tracerModel() const
Definition: FlowProblem.hpp:936
const SolidEnergyLawParams & solidEnergyLawParams(unsigned globalSpaceIdx, unsigned) const
Definition: FlowProblem.hpp:1158
Scalar getRockCompTransMultVal(std::size_t dofIdx) const
Definition: FlowProblem.hpp:2543
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
Direct indexed acces to the depth of an degree of freedom [m].
Definition: FlowProblem.hpp:976
LhsEval wellTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Return the well transmissibility multiplier due to rock changues.
Definition: FlowProblem.hpp:1743
const MaterialLawParams & materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const
Definition: FlowProblem.hpp:1017
std::shared_ptr< const EclMaterialLawManager > materialLawManager() const
Returns the ECL material law manager.
Definition: FlowProblem.hpp:1052
unsigned plmixnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition: FlowProblem.hpp:1120
const MaterialLawParams & materialLawParams(unsigned globalDofIdx) const
Definition: FlowProblem.hpp:1012
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1141
Scalar thermalHalfTransmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition: FlowProblem.hpp:911
unsigned miscnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition: FlowProblem.hpp:1112
Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
Returns the maximum value of the gas dissolution factor at the current time for a given degree of fre...
Definition: FlowProblem.hpp:1251
FlowThresholdPressure< TypeTag > & thresholdPressure()
Definition: FlowProblem.hpp:933
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:776
bool updateHysteresis_()
Definition: FlowProblem.hpp:2497
void readThermalParameters_()
Definition: FlowProblem.hpp:2061
void serializeOp(Serializer &serializer)
Definition: FlowProblem.hpp:1816
void finalizeOutput()
Definition: FlowProblem.hpp:764
Scalar maxPolymerAdsorption(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the max polymer adsorption value.
Definition: FlowProblem.hpp:1128
void beginTimeStep()
Called by the simulator before each time integration.
Definition: FlowProblem.hpp:582
Scalar thermalHalfTransmissibilityIn(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:885
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
Sets an element's maximum oil phase saturation observed during the simulation.
Definition: FlowProblem.hpp:1239
void setSimulationReport(const SimulatorReport &report)
Definition: FlowProblem.hpp:1528
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:1027
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const
Definition: FlowProblem.hpp:835
AquiferModel & mutableAquiferModel()
Definition: FlowProblem.hpp:1515
Scalar temperature(unsigned globalDofIdx, unsigned) const
Definition: FlowProblem.hpp:1150
void finishInit()
Definition: FlowProblem.hpp:343
bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities &iq)
Definition: FlowProblem.hpp:2531
void readEclRestartSolution_()
Definition: FlowProblem.hpp:2149
WellModel & wellModel()
Definition: FlowProblem.hpp:1509
void updateProperty_(const std::string &failureMsg, UpdateFunc func)
Definition: FlowProblem.hpp:1855
void setConvData(const std::vector< std::vector< int > > &data)
Definition: FlowProblem.hpp:1810
const EclipseIO & eclIO() const
Definition: FlowProblem.hpp:1522
void initialSolutionApplied()
Definition: FlowProblem.hpp:1333
bool updateMaxOilSaturation_()
Definition: FlowProblem.hpp:1907
static void registerParameters()
Definition: FlowProblem.hpp:209
static int handlePositionalParameter(std::set< std::string > &seenParams, std::string &errorMsg, int, const char **argv, int paramIdx, int)
Definition: FlowProblem.hpp:251
void updateMaxPolymerAdsorption_()
Definition: FlowProblem.hpp:2521
void updateCompositionChangeLimits_()
Definition: FlowProblem.hpp:1875
void beginEpisode()
Called by the simulator before an episode begins.
Definition: FlowProblem.hpp:513
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1371
const InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const
Definition: FlowProblem.hpp:1534
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1005
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition: FlowProblem.hpp:842
Scalar diffusivity(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition: FlowProblem.hpp:817
void deserialize(Restarter &res)
This method restores the complete state of the problem and its sub-objects from disk.
Definition: FlowProblem.hpp:479
void processRestartSaturations_(InitialFluidState &elemFluidState, Scalar &solventSaturation)
Definition: FlowProblem.hpp:2269
std::shared_ptr< EclMaterialLawManager > materialLawManager()
Definition: FlowProblem.hpp:1088
bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities &iq)
Definition: FlowProblem.hpp:1956
Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the depth of an degree of freedom [m].
Definition: FlowProblem.hpp:964
bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities &iq)
Definition: FlowProblem.hpp:1925
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1293
const AquiferModel & aquiferModel() const
Definition: FlowProblem.hpp:1512
LhsEval permFactTransMultiplier(const IntensiveQuantities &intQuants) const
Calculate the transmissibility multiplier due to porosity reduction.
Definition: FlowProblem.hpp:1725
LhsEval rockCompPoroMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the porosity multiplier due to water induced rock compaction.
Definition: FlowProblem.hpp:1672
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition: FlowProblem.hpp:865
TracerModel & tracerModel()
Definition: FlowProblem.hpp:939
void readMaterialParameters_()
Definition: FlowProblem.hpp:2027
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:1176
void serialize(Restarter &res)
This method writes the complete state of the problem and its subobjects to disk.
Definition: FlowProblem.hpp:498
Scalar rockReferencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:995
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const
Definition: FlowProblem.hpp:828
void prefetch(const Element &elem) const
Definition: FlowProblem.hpp:464
This class calculates the threshold pressure for grid faces according to the Eclipse Reference Manual...
Definition: FlowThresholdPressure.hpp:59
void init(std::size_t numDof, int episodeIdx, const unsigned ntpvt)
Definition: RelpermDiagnostics.hpp:50
void diagnosis(const EclipseState &eclState, const CartesianIndexMapper &cartesianIndexMapper)
Definition: SimulatorTimer.hpp:39
A class which handles tracers as specified in by ECL.
Definition: TracerModel.hpp:67
VTK output module for the tracer model's parameters.
Definition: VtkTracerModule.hpp:70
static void registerParameters()
Register all run-time parameters for the tracer VTK output module.
Definition: VtkTracerModule.hpp:96
int endIteration(int rank)
@ NONE
Definition: DeferredLogger.hpp:46
Definition: BlackoilPhases.hpp:27
void eclBroadcast(Parallel::Communication, T &)
int eclPositionalParameter(Dune::ParameterTree &tree, std::set< std::string > &seenParams, std::string &errorMsg, const char **argv, int paramIdx)
Definition: FlowGenericProblem_impl.hpp:48
Definition: SimulatorReport.hpp:100
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:34