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_();
2045 this->updateKrnum_();
2049 updateReferencePorosity_();
2050 this->referencePorosity_[1] = this->referencePorosity_[0];
2055 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
2056 materialLawManager_->initFromState(eclState);
2057 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2058 this->
template fieldPropIntTypeOnLeafAssigner_<int>(),
2059 this-> lookupIdxOnLevelZeroAssigner_());
2065 if constexpr (enableEnergy)
2067 const auto& simulator = this->simulator();
2068 const auto& vanguard = simulator.vanguard();
2069 const auto& eclState = vanguard.eclState();
2072 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
2073 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2074 this-> fieldPropDoubleOnLeafAssigner_(),
2075 this->
template fieldPropIntTypeOnLeafAssigner_<unsigned int>());
2081 const auto& simulator = this->simulator();
2082 const auto& vanguard = simulator.vanguard();
2083 const auto& eclState = vanguard.eclState();
2085 std::size_t numDof = this->model().numGridDof();
2087 this->referencePorosity_[0].resize(numDof);
2089 const auto& fp = eclState.fieldProps();
2090 const std::vector<double> porvData =
this -> fieldPropDoubleOnLeafAssigner_()(fp,
"PORV");
2091 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2092 Scalar poreVolume = porvData[dofIdx];
2097 Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx);
2098 assert(dofVolume > 0.0);
2099 this->referencePorosity_[0][dofIdx] = poreVolume/dofVolume;
2105 const auto& simulator = this->simulator();
2106 const auto& vanguard = simulator.vanguard();
2107 const auto& eclState = vanguard.eclState();
2109 if (eclState.getInitConfig().hasEquil())
2110 readEquilInitialCondition_();
2112 readExplicitInitialCondition_();
2114 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
2115 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
2118 enablePolymerMolarWeight,
2122 std::size_t numElems = this->model().numGridDof();
2123 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2124 const auto& fs = initialFluidStates_[elemIdx];
2125 if (!this->maxWaterSaturation_.empty())
2126 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
2127 if (!this->maxOilSaturation_.empty())
2128 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
2129 if (!this->minRefPressure_.empty())
2130 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
2138 const auto& simulator = this->simulator();
2143 std::size_t numElems = this->model().numGridDof();
2144 initialFluidStates_.resize(numElems);
2145 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2146 auto& elemFluidState = initialFluidStates_[elemIdx];
2154 if(this->simulator().vanguard().grid().maxLevel() > 0) {
2155 throw std::invalid_argument(
"Refined grids are not yet supported for restart ");
2159 auto& simulator = this->simulator();
2160 const auto& schedule = simulator.vanguard().schedule();
2161 const auto& eclState = simulator.vanguard().eclState();
2162 const auto& initconfig = eclState.getInitConfig();
2164 int restart_step = initconfig.getRestartStep();
2166 simulator.setTime(schedule.seconds(restart_step));
2168 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
2169 schedule.stepLength(restart_step));
2170 simulator.setEpisodeIndex(restart_step);
2172 eclWriter_->beginRestart();
2174 Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength());
2175 simulator.setTimeStepSize(dt);
2177 std::size_t numElems = this->model().numGridDof();
2178 initialFluidStates_.resize(numElems);
2179 if constexpr (enableSolvent) {
2180 this->solventSaturation_.resize(numElems, 0.0);
2181 this->solventRsw_.resize(numElems, 0.0);
2184 if constexpr (enablePolymer)
2185 this->polymer_.concentration.resize(numElems, 0.0);
2187 if constexpr (enablePolymerMolarWeight) {
2188 const std::string msg {
"Support of the RESTART for polymer molecular weight "
2189 "is not implemented yet. The polymer weight value will be "
2190 "zero when RESTART begins"};
2191 OpmLog::warning(
"NO_POLYMW_RESTART", msg);
2192 this->polymer_.moleWeight.resize(numElems, 0.0);
2195 if constexpr (enableMICP) {
2196 this->micp_.resize(numElems);
2199 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2200 auto& elemFluidState = initialFluidStates_[elemIdx];
2201 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
2202 eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
2203 eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
2212 auto ssol = enableSolvent
2213 ? eclWriter_->outputModule().getSolventSaturation(elemIdx)
2216 processRestartSaturations_(elemFluidState, ssol);
2218 if constexpr (enableSolvent) {
2219 this->solventSaturation_[elemIdx] = ssol;
2220 this->solventRsw_[elemIdx] = eclWriter_->outputModule().getSolventRsw(elemIdx);
2224 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
2226 if constexpr (enablePolymer)
2227 this->polymer_.concentration[elemIdx] = eclWriter_->outputModule().getPolymerConcentration(elemIdx);
2228 if constexpr (enableMICP){
2229 this->micp_.microbialConcentration[elemIdx] = eclWriter_->outputModule().getMicrobialConcentration(elemIdx);
2230 this->micp_.oxygenConcentration[elemIdx] = eclWriter_->outputModule().getOxygenConcentration(elemIdx);
2231 this->micp_.ureaConcentration[elemIdx] = eclWriter_->outputModule().getUreaConcentration(elemIdx);
2232 this->micp_.biofilmConcentration[elemIdx] = eclWriter_->outputModule().getBiofilmConcentration(elemIdx);
2233 this->micp_.calciteConcentration[elemIdx] = eclWriter_->outputModule().getCalciteConcentration(elemIdx);
2238 const int episodeIdx = this->episodeIndex();
2239 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
2244 auto& sol = this->model().solution(0);
2245 const auto& gridView = this->gridView();
2246 ElementContext elemCtx(simulator);
2247 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2248 elemCtx.updatePrimaryStencil(elem);
2249 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
2250 initial(sol[elemIdx], elemCtx, 0, 0);
2258 this->model().syncOverlap();
2260 eclWriter_->endRestart();
2267 const Scalar smallSaturationTolerance = 1.e-6;
2268 Scalar sumSaturation = 0.0;
2269 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2270 if (FluidSystem::phaseIsActive(phaseIdx)) {
2271 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
2272 elemFluidState.setSaturation(phaseIdx, 0.0);
2274 sumSaturation += elemFluidState.saturation(phaseIdx);
2278 if constexpr (enableSolvent) {
2279 if (solventSaturation < smallSaturationTolerance)
2280 solventSaturation = 0.0;
2282 sumSaturation += solventSaturation;
2285 assert(sumSaturation > 0.0);
2287 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2288 if (FluidSystem::phaseIsActive(phaseIdx)) {
2289 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
2290 elemFluidState.setSaturation(phaseIdx, saturation);
2293 if constexpr (enableSolvent) {
2294 solventSaturation = solventSaturation / sumSaturation;
2300 const auto& simulator = this->simulator();
2301 const auto& vanguard = simulator.vanguard();
2302 const auto& eclState = vanguard.eclState();
2303 const auto& fp = eclState.fieldProps();
2304 bool has_swat = fp.has_double(
"SWAT");
2305 bool has_sgas = fp.has_double(
"SGAS");
2306 bool has_rs = fp.has_double(
"RS");
2307 bool has_rv = fp.has_double(
"RV");
2308 bool has_rvw = fp.has_double(
"RVW");
2309 bool has_pressure = fp.has_double(
"PRESSURE");
2310 bool has_salt = fp.has_double(
"SALT");
2311 bool has_saltp = fp.has_double(
"SALTP");
2314 if (Indices::numPhases > 1) {
2315 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
2316 throw std::runtime_error(
"The ECL input file requires the presence of the SWAT keyword if "
2317 "the water phase is active");
2318 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
2319 throw std::runtime_error(
"The ECL input file requires the presence of the SGAS keyword if "
2320 "the gas phase is active");
2323 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
2324 "keyword if the model is initialized explicitly");
2325 if (FluidSystem::enableDissolvedGas() && !has_rs)
2326 throw std::runtime_error(
"The ECL input file requires the RS keyword to be present if"
2327 " dissolved gas is enabled");
2328 if (FluidSystem::enableVaporizedOil() && !has_rv)
2329 throw std::runtime_error(
"The ECL input file requires the RV keyword to be present if"
2330 " vaporized oil is enabled");
2331 if (FluidSystem::enableVaporizedWater() && !has_rvw)
2332 throw std::runtime_error(
"The ECL input file requires the RVW keyword to be present if"
2333 " vaporized water is enabled");
2334 if (enableBrine && !has_salt)
2335 throw std::runtime_error(
"The ECL input file requires the SALT keyword to be present if"
2336 " brine is enabled and the model is initialized explicitly");
2337 if (enableSaltPrecipitation && !has_saltp)
2338 throw std::runtime_error(
"The ECL input file requires the SALTP keyword to be present if"
2339 " salt precipitation is enabled and the model is initialized explicitly");
2341 std::size_t numDof = this->model().numGridDof();
2343 initialFluidStates_.resize(numDof);
2345 std::vector<double> waterSaturationData;
2346 std::vector<double> gasSaturationData;
2347 std::vector<double> pressureData;
2348 std::vector<double> rsData;
2349 std::vector<double> rvData;
2350 std::vector<double> rvwData;
2351 std::vector<double> tempiData;
2352 std::vector<double> saltData;
2353 std::vector<double> saltpData;
2355 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
2356 waterSaturationData = fp.get_double(
"SWAT");
2358 waterSaturationData.resize(numDof);
2360 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
2361 gasSaturationData = fp.get_double(
"SGAS");
2363 gasSaturationData.resize(numDof);
2365 pressureData = fp.get_double(
"PRESSURE");
2366 if (FluidSystem::enableDissolvedGas())
2367 rsData = fp.get_double(
"RS");
2369 if (FluidSystem::enableVaporizedOil())
2370 rvData = fp.get_double(
"RV");
2372 if (FluidSystem::enableVaporizedWater())
2373 rvwData = fp.get_double(
"RVW");
2376 tempiData = fp.get_double(
"TEMPI");
2379 if constexpr (enableBrine)
2380 saltData = fp.get_double(
"SALT");
2383 if constexpr (enableSaltPrecipitation)
2384 saltpData = fp.get_double(
"SALTP");
2387 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2388 auto& dofFluidState = initialFluidStates_[dofIdx];
2390 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
2395 Scalar temperatureLoc = tempiData[dofIdx];
2396 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
2397 temperatureLoc = FluidSystem::surfaceTemperature;
2398 dofFluidState.setTemperature(temperatureLoc);
2403 if constexpr (enableBrine)
2404 dofFluidState.setSaltConcentration(saltData[dofIdx]);
2409 if constexpr (enableSaltPrecipitation)
2410 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
2415 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
2416 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
2417 waterSaturationData[dofIdx]);
2419 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
2420 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
2421 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2423 - waterSaturationData[dofIdx]);
2426 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2427 gasSaturationData[dofIdx]);
2429 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
2430 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
2432 - waterSaturationData[dofIdx]
2433 - gasSaturationData[dofIdx]);
2438 Scalar pressure = pressureData[dofIdx];
2442 std::array<Scalar, numPhases> pc = {0};
2443 const auto& matParams = materialLawParams(dofIdx);
2444 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
2445 Valgrind::CheckDefined(pressure);
2446 Valgrind::CheckDefined(pc);
2447 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2448 if (!FluidSystem::phaseIsActive(phaseIdx))
2451 if (Indices::oilEnabled)
2452 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
2453 else if (Indices::gasEnabled)
2454 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
2455 else if (Indices::waterEnabled)
2457 dofFluidState.setPressure(phaseIdx, pressure);
2460 if (FluidSystem::enableDissolvedGas())
2461 dofFluidState.setRs(rsData[dofIdx]);
2462 else if (Indices::gasEnabled && Indices::oilEnabled)
2463 dofFluidState.setRs(0.0);
2465 if (FluidSystem::enableVaporizedOil())
2466 dofFluidState.setRv(rvData[dofIdx]);
2467 else if (Indices::gasEnabled && Indices::oilEnabled)
2468 dofFluidState.setRv(0.0);
2470 if (FluidSystem::enableVaporizedWater())
2471 dofFluidState.setRvw(rvwData[dofIdx]);
2476 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2477 if (!FluidSystem::phaseIsActive(phaseIdx))
2480 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2481 dofFluidState.setInvB(phaseIdx, b);
2483 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2484 dofFluidState.setDensity(phaseIdx, rho);
2493 if (!materialLawManager_->enableHysteresis())
2498 this->updateProperty_(
"FlowProblem::updateHysteresis_() failed:",
2499 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2501 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2509 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
2510 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2518 this->updateProperty_(
"FlowProblem::updateMaxPolymerAdsorption_() failed:",
2519 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2521 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
2527 const Scalar pa = scalarValue(iq.polymerAdsorption());
2528 auto& mpa = this->polymer_.maxAdsorption;
2529 if (mpa[compressedDofIdx] < pa) {
2530 mpa[compressedDofIdx] = pa;
2539 if (this->rockCompTransMultVal_.empty())
2542 return this->rockCompTransMultVal_[dofIdx];
2549 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
2550 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
2551 ConditionalStorage<enableDiffusion, Scalar> diffusivity;
2552 ConditionalStorage<enableDispersion, Scalar> dispersivity;
2553 Scalar transmissibility;
2557 void updatePffDofData_()
2559 const auto& distFn =
2560 [
this](PffDofData_& dofData,
2561 const Stencil& stencil,
2562 unsigned localDofIdx)
2565 const auto& elementMapper = this->model().elementMapper();
2567 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
2568 if (localDofIdx != 0) {
2569 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(0));
2570 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
2572 if constexpr (enableEnergy) {
2573 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
2574 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
2576 if constexpr (enableDiffusion)
2577 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
2578 if (enableDispersion)
2579 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
2583 pffDofData_.update(distFn);
2586 void readBoundaryConditions_()
2588 const auto& simulator = this->simulator();
2589 const auto& vanguard = simulator.vanguard();
2590 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
2591 if (bcconfig.size() > 0) {
2592 nonTrivialBoundaryConditions_ =
true;
2594 std::size_t numCartDof = vanguard.cartesianSize();
2595 unsigned numElems = vanguard.gridView().size(0);
2596 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
2598 for (
unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
2599 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
2601 bcindex_.resize(numElems, 0);
2602 auto loopAndApply = [&cartesianToCompressedElemIdx,
2603 &vanguard](
const auto& bcface,
2606 for (
int i = bcface.i1; i <= bcface.i2; ++i) {
2607 for (
int j = bcface.j1; j <= bcface.j2; ++j) {
2608 for (
int k = bcface.k1; k <= bcface.k2; ++k) {
2609 std::array<int, 3> tmp = {i,j,k};
2610 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
2617 for (
const auto& bcface : bcconfig) {
2618 std::vector<int>& data = bcindex_(bcface.dir);
2619 const int index = bcface.index;
2620 loopAndApply(bcface,
2621 [&data,index](
int elemIdx)
2622 { data[elemIdx] = index; });
2629 Scalar limitNextTimeStepSize_(Scalar dtNext)
const
2631 if constexpr (enableExperiments) {
2632 const auto& simulator = this->simulator();
2633 const auto& schedule = simulator.vanguard().schedule();
2634 int episodeIdx = simulator.episodeIndex();
2637 Scalar maxTimeStepSize = Parameters::get<TypeTag, Properties::SolverMaxTimeStepInDays>() * 24 * 60 * 60;
2638 int reportStepIdx = std::max(episodeIdx, 0);
2639 if (this->enableTuning_) {
2640 const auto& tuning = schedule[reportStepIdx].tuning();
2641 maxTimeStepSize = tuning.TSMAXZ;
2644 dtNext = std::min(dtNext, maxTimeStepSize);
2646 Scalar remainingEpisodeTime =
2647 simulator.episodeStartTime() + simulator.episodeLength()
2648 - (simulator.startTime() + simulator.time());
2649 assert(remainingEpisodeTime >= 0.0);
2653 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
2656 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
2658 if (simulator.episodeStarts()) {
2661 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
2662 bool wellEventOccured =
2663 events.hasEvent(ScheduleEvents::NEW_WELL)
2664 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
2665 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
2666 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
2667 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
2668 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
2675 void computeAndSetEqWeights_()
2677 std::vector<Scalar> sumInvB(numPhases, 0.0);
2678 const auto& gridView = this->gridView();
2679 ElementContext elemCtx(this->simulator());
2680 for(
const auto& elem: elements(gridView, Dune::Partitions::interior)) {
2681 elemCtx.updatePrimaryStencil(elem);
2682 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
2683 const auto& dofFluidState = initialFluidStates_[elemIdx];
2684 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2685 if (!FluidSystem::phaseIsActive(phaseIdx))
2688 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
2692 std::size_t numDof = this->model().numGridDof();
2693 const auto& comm = this->simulator().vanguard().grid().comm();
2694 comm.sum(sumInvB.data(),sumInvB.size());
2695 Scalar numTotalDof = comm.sum(numDof);
2697 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2698 if (!FluidSystem::phaseIsActive(phaseIdx))
2701 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
2702 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
2703 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
2704 this->model().setEqWeight(activeSolventCompIdx, avgB);
2708 int refPressurePhaseIdx_()
const {
2709 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2712 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2716 return waterPhaseIdx;
2720 void updateRockCompTransMultVal_()
2722 const auto& model = this->simulator().model();
2723 std::size_t numGridDof = this->model().numGridDof();
2724 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
2725 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
2726 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, 0);
2727 Scalar trans_mult = computeRockCompTransMultiplier_<Scalar>(iq, elementIdx);
2728 this->rockCompTransMultVal_[elementIdx] = trans_mult;
2737 template <
class LhsEval>
2738 LhsEval computeRockCompTransMultiplier_(
const IntensiveQuantities& intQuants,
unsigned elementIdx)
const
2740 OPM_TIMEBLOCK_LOCAL(computeRockCompTransMultiplier);
2741 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
2744 unsigned tableIdx = 0;
2745 if (!this->rockTableIdx_.empty())
2746 tableIdx = this->rockTableIdx_[elementIdx];
2748 const auto& fs = intQuants.fluidState();
2749 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
2751 if (!this->minRefPressure_.empty())
2754 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
2755 this->minRefPressure_[elementIdx]);
2757 if (!this->overburdenPressure_.empty())
2758 effectivePressure -= this->overburdenPressure_[elementIdx];
2760 if (!this->rockCompTransMult_.empty())
2761 return this->rockCompTransMult_[tableIdx].eval(effectivePressure,
true);
2764 assert(!this->rockCompTransMultWc_.empty());
2765 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
2766 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
2768 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax,
true);
2771 typename Vanguard::TransmissibilityType transmissibilities_;
2773 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
2774 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
2776 FlowThresholdPressure<TypeTag> thresholdPressures_;
2778 std::vector<InitialFluidState> initialFluidStates_;
2780 bool enableDriftCompensation_;
2781 GlobalEqVector drift_;
2783 WellModel wellModel_;
2784 AquiferModel aquiferModel_;
2786 bool enableEclOutput_;
2787 std::unique_ptr<EclWriterType> eclWriter_;
2790 bool enableDamarisOutput_ = false ;
2791 std::unique_ptr<DamarisWriterType> damarisWriter_;
2794 PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
2795 TracerModel tracerModel_;
2797 ActionHandler actionHandler_;
2802 std::array<std::vector<T>,6> data;
2804 void resize(std::size_t size, T defVal)
2806 for (
auto& d : data)
2807 d.resize(size, defVal);
2810 const std::vector<T>& operator()(FaceDir::DirEnum dir)
const
2812 if (dir == FaceDir::DirEnum::Unknown)
2813 throw std::runtime_error(
"Tried to access BC data for the 'Unknown' direction");
2815 int div =
static_cast<int>(dir);
2816 while ((div /= 2) >= 1)
2818 assert(idx >= 0 && idx <= 5);
2822 std::vector<T>& operator()(FaceDir::DirEnum dir)
2824 return const_cast<std::vector<T>&
>(std::as_const(*
this)(dir));
2828 BCData<int> bcindex_;
2829 bool nonTrivialBoundaryConditions_ =
false;
2830 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:366
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:379
bool enableTuning_
Definition: FlowGenericProblem.hpp:378
MixingRateControls< GetPropType< TypeTag, Properties::FluidSystem > > mixControls_
Definition: FlowGenericProblem.hpp:375
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:367
int numPressurePointsEquil_
Definition: FlowGenericProblem.hpp:383
void initFluidSystem_()
Definition: FlowGenericProblem_impl.hpp:505
Scalar maxTimeStepAfterWellEvent_
Definition: FlowGenericProblem.hpp:380
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:2136
void readExplicitInitialCondition_()
Definition: FlowProblem.hpp:2298
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:2507
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:2079
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:2103
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:2537
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:2491
void readThermalParameters_()
Definition: FlowProblem.hpp:2063
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:2525
void readEclRestartSolution_()
Definition: FlowProblem.hpp:2151
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:2515
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:2263
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