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>;
184 using InitialFluidState =
typename EquilInitializer<TypeTag>::ScalarFluidState;
186 using Toolbox = MathToolbox<Evaluation>;
187 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
195 using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
211 ParentType::registerParameters();
212 EclWriterType::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();
260 return eclPositionalParameter(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>();
319 this->initialTimeStepSize_ = Parameters::get<TypeTag, Properties::InitialTimeStepSize>();
320 this->maxTimeStepAfterWellEvent_ = Parameters::get<TypeTag, Properties::TimeStepAfterEventInDays>() * 24 * 60 * 60;
326 if (Parameters::isSet<TypeTag, Properties::NumPressurePointsEquil>())
328 this->numPressurePointsEquil_ = Parameters::get<TypeTag, Properties::NumPressurePointsEquil>();
330 this->numPressurePointsEquil_ = simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP();
334 relpermDiagnostics.
diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
342 ParentType::finishInit();
344 auto& simulator = this->simulator();
345 const auto& eclState = simulator.vanguard().eclState();
346 const auto& schedule = simulator.vanguard().schedule();
349 simulator.setStartTime(schedule.getStartTime());
350 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
356 simulator.setEpisodeIndex(-1);
357 simulator.setEpisodeLength(0.0);
362 this->gravity_ = 0.0;
363 if (Parameters::get<TypeTag, Properties::EnableGravity>())
364 this->gravity_[dim - 1] = 9.80665;
365 if (!eclState.getInitConfig().hasGravity())
366 this->gravity_[dim - 1] = 0.0;
368 if (this->enableTuning_) {
371 const auto& tuning = schedule[0].tuning();
372 this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
373 this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
376 this->initFluidSystem_();
379 this->mixControls_.init(this->model().numGridDof(),
380 this->episodeIndex(),
381 eclState.runspec().tabdims().getNumPVTTables());
383 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
384 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
385 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
388 this->readRockParameters_(simulator.vanguard().cellCenterDepths(),
389 [&simulator](
const unsigned idx)
391 std::array<int,dim> coords;
392 simulator.vanguard().cartesianCoordinate(idx, coords);
393 for (auto& c : coords) {
398 readMaterialParameters_();
399 readThermalParameters_();
402 std::function<
unsigned int(
unsigned int)> gridToEquilGrid = [&simulator](
unsigned int i) {
403 return simulator.vanguard().gridIdxToEquilGridIdx(i);
405 transmissibilities_.finishInit(gridToEquilGrid);
407 const auto& initconfig = eclState.getInitConfig();
408 tracerModel_.init(initconfig.restartRequested());
409 if (initconfig.restartRequested())
410 readEclRestartSolution_();
412 readInitialCondition_();
414 tracerModel_.prepareTracerBatches();
418 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
419 const auto& vanguard = this->simulator().vanguard();
420 const auto& gridView = vanguard.gridView();
421 int numElements = gridView.size(0);
422 this->polymer_.maxAdsorption.resize(numElements, 0.0);
425 readBoundaryConditions_();
428 computeAndSetEqWeights_();
430 if (enableDriftCompensation_) {
431 drift_.resize(this->model().numGridDof());
436 if (enableEclOutput_) {
437 if (simulator.vanguard().grid().comm().size() > 1) {
438 if (simulator.vanguard().grid().comm().rank() == 0)
439 eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
441 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
444 std::function<
unsigned int(
unsigned int)> equilGridToGrid = [&simulator](
unsigned int i) {
445 return simulator.vanguard().gridEquilIdxToGridIdx(i);
447 eclWriter_->writeInit(equilGridToGrid);
450 simulator.vanguard().releaseGlobalTransmissibilities();
455 if (!initconfig.restartRequested()) {
456 simulator.startNextEpisode(schedule.seconds(1));
457 simulator.setEpisodeIndex(0);
461 void prefetch(
const Element& elem)
const
462 { pffDofData_.prefetch(elem); }
475 template <
class Restarter>
482 wellModel_.deserialize(res);
485 aquiferModel_.deserialize(res);
494 template <
class Restarter>
497 wellModel_.serialize(res);
499 aquiferModel_.serialize(res);
502 int episodeIndex()
const
504 return std::max(this->simulator().episodeIndex(), 0);
512 OPM_TIMEBLOCK(beginEpisode);
514 auto& simulator = this->simulator();
515 int episodeIdx = simulator.episodeIndex();
516 auto& eclState = simulator.vanguard().eclState();
517 const auto& schedule = simulator.vanguard().schedule();
518 const auto& events = schedule[episodeIdx].events();
520 if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
527 const auto& miniDeck = schedule[episodeIdx].geo_keywords();
528 const auto& cc = simulator.vanguard().grid().comm();
529 eclState.apply_schedule_keywords( miniDeck );
530 eclBroadcast(cc, eclState.getTransMult() );
533 std::function<
unsigned int(
unsigned int)> equilGridToGrid = [&simulator](
unsigned int i) {
534 return simulator.vanguard().gridEquilIdxToGridIdx(i);
538 transmissibilities_.update(
true, equilGridToGrid);
539 this->referencePorosity_[1] = this->referencePorosity_[0];
540 updateReferencePorosity_();
542 this->model().linearizer().updateDiscretizationParameters();
545 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
548 wellModel_.beginEpisode();
551 aquiferModel_.beginEpisode();
554 Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
556 if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
559 dt = std::min(dt, this->initialTimeStepSize_);
560 simulator.setTimeStepSize(dt);
564 actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
566 if (episodeIdx >= 0) {
567 const auto& oilVap = schedule[episodeIdx].oilvap();
568 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
569 FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2());
571 FluidSystem::setVapPars(0.0, 0.0);
581 OPM_TIMEBLOCK(beginTimeStep);
582 int episodeIdx = this->episodeIndex();
584 this->beginTimeStep_(enableExperiments,
586 this->simulator().timeStepIndex(),
587 this->simulator().startTime(),
588 this->simulator().time(),
589 this->simulator().timeStepSize(),
590 this->simulator().endTime());
594 asImp_().updateExplicitQuantities_();
596 if (nonTrivialBoundaryConditions()) {
597 this->model().linearizer().updateBoundaryConditionData();
600 wellModel_.beginTimeStep();
601 aquiferModel_.beginTimeStep();
602 tracerModel_.beginTimeStep();
611 OPM_TIMEBLOCK(beginIteration);
612 wellModel_.beginIteration();
613 aquiferModel_.beginIteration();
621 OPM_TIMEBLOCK(endIteration);
622 wellModel_.endIteration();
623 aquiferModel_.endIteration();
631 OPM_TIMEBLOCK(endTimeStep);
633 if constexpr (getPropValue<TypeTag, Properties::EnableDebuggingChecks>()) {
638 int rank = this->simulator().gridView().comm().rank();
640 std::cout <<
"checking conservativeness of solution\n";
641 this->model().checkConservativeness(-1,
true);
643 std::cout <<
"solution is sufficiently conservative\n";
647 auto& simulator = this->simulator();
648 wellModel_.endTimeStep();
649 aquiferModel_.endTimeStep();
650 tracerModel_.endTimeStep();
654 this->model().linearizer().updateFlowsInfo();
657 asImp_().updateCompositionChangeLimits_();
659 OPM_TIMEBLOCK(driftCompansation);
660 if (enableDriftCompensation_) {
661 const auto& residual = this->model().linearizer().residual();
662 for (
unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
663 drift_[globalDofIdx] = residual[globalDofIdx];
664 drift_[globalDofIdx] *= simulator.timeStepSize();
665 if constexpr (getPropValue<TypeTag, Properties::UseVolumetricResidual>())
666 drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
670 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
671 !this->simulator().episodeWillBeOver();
673 const auto& grid = this->simulator().vanguard().gridView().grid();
674 using GridType = std::remove_cv_t<
typename std::remove_reference<
decltype(grid)>::type>;
675 bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
676 if ( !isCpGrid || (this->simulator().vanguard().gridView().grid().maxLevel()==0)) {
677 eclWriter_->evalSummaryState(isSubStep);
680 int episodeIdx = this->episodeIndex();
683 std::function<
unsigned int(
unsigned int)> gridToEquilGrid = [&simulator](
unsigned int i) {
684 return simulator.vanguard().gridIdxToEquilGridIdx(i);
687 std::function<void(
bool)> transUp =
688 [
this,gridToEquilGrid](
bool global) {
689 this->transmissibilities_.update(global,gridToEquilGrid);
692 OPM_TIMEBLOCK(applyActions);
693 actionHandler_.applyActions(episodeIdx,
694 simulator.time() + simulator.timeStepSize(),
698 if constexpr (enableMICP){
699 auto& model = this->model();
700 const auto& residual = this->model().linearizer().residual();
701 for (
unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
702 auto& phi = this->referencePorosity_[1][globalDofIdx];
703 MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
713 OPM_TIMEBLOCK(endEpisode);
714 auto& simulator = this->simulator();
715 auto& schedule = simulator.vanguard().schedule();
717 wellModel_.endEpisode();
718 aquiferModel_.endEpisode();
720 int episodeIdx = this->episodeIndex();
722 if (episodeIdx + 1 >=
static_cast<int>(schedule.size() - 1)) {
723 simulator.setFinished(
true);
728 simulator.startNextEpisode(schedule.stepLength(episodeIdx + 1));
737 OPM_TIMEBLOCK(problemWriteOutput);
740 if (Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() ||
741 this->simulator().episodeWillBeOver()) {
742 ParentType::writeOutput(verbose);
745 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
746 !this->simulator().episodeWillBeOver();
748 data::Solution localCellData = {};
752 if (enableDamarisOutput_) {
753 damarisWriter_->writeOutput(localCellData, isSubStep) ;
756 if (enableEclOutput_){
757 eclWriter_->writeOutput(std::move(localCellData), timer, isSubStep);
761 void finalizeOutput() {
762 OPM_TIMEBLOCK(finalizeOutput);
772 template <
class Context>
775 unsigned timeIdx)
const
777 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
778 return transmissibilities_.permeability(globalSpaceIdx);
788 {
return transmissibilities_.permeability(globalElemIdx); }
793 template <
class Context>
795 [[maybe_unused]]
unsigned fromDofLocalIdx,
796 unsigned toDofLocalIdx)
const
798 assert(fromDofLocalIdx == 0);
799 return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
807 return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
813 template <
class Context>
815 [[maybe_unused]]
unsigned fromDofLocalIdx,
816 unsigned toDofLocalIdx)
const
818 assert(fromDofLocalIdx == 0);
819 return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
825 Scalar
diffusivity(
const unsigned globalCellIn,
const unsigned globalCellOut)
const{
826 return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
832 Scalar
dispersivity(
const unsigned globalCellIn,
const unsigned globalCellOut)
const{
833 return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
840 const unsigned boundaryFaceIdx)
const
842 return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
851 template <
class Context>
853 unsigned boundaryFaceIdx)
const
855 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
856 return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
863 const unsigned boundaryFaceIdx)
const
865 return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
873 const unsigned globalSpaceIdxOut)
const
875 return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
881 template <
class Context>
884 unsigned timeIdx)
const
886 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
887 unsigned toDofLocalIdx = face.exteriorIndex();
888 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
894 template <
class Context>
897 unsigned timeIdx)
const
899 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
900 unsigned toDofLocalIdx = face.exteriorIndex();
901 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
907 template <
class Context>
909 unsigned boundaryFaceIdx)
const
911 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
912 return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
919 {
return transmissibilities_; }
925 {
return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
928 {
return thresholdPressures_; }
930 FlowThresholdPressure<TypeTag>& thresholdPressure()
931 {
return thresholdPressures_; }
933 const TracerModel& tracerModel()
const
934 {
return tracerModel_; }
936 TracerModel& tracerModel()
937 {
return tracerModel_; }
947 template <
class Context>
948 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
950 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
951 return this->porosity(globalSpaceIdx, timeIdx);
960 template <
class Context>
961 Scalar
dofCenterDepth(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
963 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
964 return this->dofCenterDepth(globalSpaceIdx);
975 return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
981 template <
class Context>
984 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
985 return this->rockCompressibility(globalSpaceIdx);
991 template <
class Context>
994 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
995 return this->rockReferencePressure(globalSpaceIdx);
1001 template <
class Context>
1003 unsigned spaceIdx,
unsigned timeIdx)
const
1005 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1006 return this->materialLawParams(globalSpaceIdx);
1009 const MaterialLawParams& materialLawParams(
unsigned globalDofIdx)
const
1011 return materialLawManager_->materialLawParams(globalDofIdx);
1014 const MaterialLawParams& materialLawParams(
unsigned globalDofIdx, FaceDir::DirEnum facedir)
const
1016 return materialLawManager_->materialLawParams(globalDofIdx, facedir);
1022 template <
class Context>
1023 const SolidEnergyLawParams&
1026 unsigned timeIdx)
const
1028 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1029 return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1035 template <
class Context>
1036 const ThermalConductionLawParams &
1039 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1040 return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1050 {
return materialLawManager_; }
1052 template <
class Flu
idState>
1053 void updateRelperms(
1054 std::array<Evaluation,numPhases> &mobility,
1055 DirectionalMobilityPtr &dirMob,
1056 FluidState &fluidState,
1057 unsigned globalSpaceIdx)
const
1059 OPM_TIMEBLOCK_LOCAL(updateRelperms);
1063 const auto& materialParams = materialLawParams(globalSpaceIdx);
1064 MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
1065 Valgrind::CheckDefined(mobility);
1067 if (materialLawManager_->hasDirectionalRelperms()
1068 || materialLawManager_->hasDirectionalImbnum())
1070 using Dir = FaceDir::DirEnum;
1071 constexpr int ndim = 3;
1072 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
1073 Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
1074 for (
int i = 0; i<ndim; i++) {
1075 const auto& materialParams = materialLawParams(globalSpaceIdx, facedirs[i]);
1076 auto& mob_array = dirMob->getArray(i);
1077 MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
1086 {
return materialLawManager_; }
1088 using BaseType::pvtRegionIndex;
1092 template <
class Context>
1093 unsigned pvtRegionIndex(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
1094 {
return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1096 using BaseType::satnumRegionIndex;
1100 template <
class Context>
1102 {
return this->satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1104 using BaseType::miscnumRegionIndex;
1108 template <
class Context>
1110 {
return this->miscnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1112 using BaseType::plmixnumRegionIndex;
1116 template <
class Context>
1118 {
return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1120 using BaseType::maxPolymerAdsorption;
1124 template <
class Context>
1126 {
return this->maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1132 {
return this->simulator().vanguard().caseName(); }
1137 template <
class Context>
1138 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
1142 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1143 return initialFluidStates_[globalDofIdx].temperature(0);
1147 Scalar temperature(
unsigned globalDofIdx,
unsigned )
const
1151 return initialFluidStates_[globalDofIdx].temperature(0);
1154 const SolidEnergyLawParams&
1155 solidEnergyLawParams(
unsigned globalSpaceIdx,
1158 return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1160 const ThermalConductionLawParams &
1161 thermalConductionLawParams(
unsigned globalSpaceIdx,
1164 return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1172 template <
class Context>
1174 const Context& context,
1176 unsigned timeIdx)
const
1178 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
1179 if (!context.intersection(spaceIdx).boundary())
1182 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
1190 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1191 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1192 values.setThermalFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
1195 if (nonTrivialBoundaryConditions()) {
1196 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
1197 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1198 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1199 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
1200 const auto [type, massrate] = boundaryCondition(globalDofIdx, indexInInside);
1201 if (type == BCType::THERMAL)
1202 values.setThermalFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1203 else if (type == BCType::FREE || type == BCType::DIRICHLET)
1204 values.setFreeFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1205 else if (type == BCType::RATE)
1206 values.setMassRate(massrate, pvtRegionIdx);
1221 if (!this->vapparsActive(this->episodeIndex()))
1224 return this->maxOilSaturation_[globalDofIdx];
1238 if (!this->vapparsActive(this->episodeIndex()))
1241 this->maxOilSaturation_[globalDofIdx] = value;
1250 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
1251 this->episodeIndex(),
1252 this->pvtRegionIndex(globalDofIdx));
1261 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
1262 this->episodeIndex(),
1263 this->pvtRegionIndex(globalDofIdx));
1276 int episodeIdx = this->episodeIndex();
1277 return !this->mixControls_.drsdtActive(episodeIdx) &&
1278 !this->mixControls_.drvdtActive(episodeIdx) &&
1279 this->rockCompPoroMultWc_.empty() &&
1280 this->rockCompPoroMult_.empty();
1289 template <
class Context>
1290 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
1292 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1294 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
1295 values.assignNaive(initialFluidStates_[globalDofIdx]);
1297 SolventModule::assignPrimaryVars(values,
1298 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
1299 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
1301 if constexpr (enablePolymer)
1302 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
1304 if constexpr (enablePolymerMolarWeight)
1305 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
1307 if constexpr (enableBrine) {
1308 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
1309 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
1312 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
1316 if constexpr (enableMICP){
1317 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
1318 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
1319 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
1320 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
1321 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
1324 values.checkDefined();
1333 this->model().invalidateAndUpdateIntensiveQuantities(0);
1339 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
1340 this->model().invalidateAndUpdateIntensiveQuantities(1);
1351 thresholdPressures_.finishInit();
1353 updateCompositionChangeLimits_();
1355 aquiferModel_.initialSolutionApplied();
1357 if (this->simulator().episodeIndex() == 0) {
1358 eclWriter_->writeInitialFIPReport();
1367 template <
class Context>
1369 const Context& context,
1371 unsigned timeIdx)
const
1373 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1374 source(rate, globalDofIdx, timeIdx);
1377 void source(RateVector& rate,
1378 unsigned globalDofIdx,
1379 unsigned timeIdx)
const
1381 OPM_TIMEBLOCK_LOCAL(eclProblemSource);
1385 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
1389 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1390 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
1392 Valgrind::CheckDefined(rate[eqIdx]);
1393 assert(isfinite(rate[eqIdx]));
1397 addToSourceDense(rate, globalDofIdx, timeIdx);
1400 void addToSourceDense(RateVector& rate,
1401 unsigned globalDofIdx,
1402 unsigned timeIdx)
const
1404 aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
1407 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
1408 std::array<int,3> ijk;
1409 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
1411 if (source.hasSource(ijk)) {
1412 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1413 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
1414 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
1415 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
1417 for (
unsigned i = 0; i < phidx_map.size(); ++i) {
1418 const auto phaseIdx = phidx_map[i];
1419 const auto sourceComp = sc_map[i];
1420 const auto compIdx = cidx_map[i];
1421 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1424 Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1425 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1426 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
1428 rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate;
1431 if constexpr (enableSolvent) {
1432 Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
1433 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1434 const auto& solventPvt = SolventModule::solventPvt();
1435 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
1437 rate[Indices::contiSolventEqIdx] += mass_rate;
1439 if constexpr (enablePolymer) {
1440 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
1442 if constexpr (enableEnergy) {
1443 for (
unsigned i = 0; i < phidx_map.size(); ++i) {
1444 const auto phaseIdx = phidx_map[i];
1445 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1448 const auto sourceComp = sc_map[i];
1449 if (source.hasHrate({ijk, sourceComp})) {
1450 rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1452 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, 0);
1453 auto fs = intQuants.fluidState();
1455 if (source.hasTemperature({ijk, sourceComp})) {
1456 Scalar temperature = source.temperature({ijk, sourceComp});
1457 fs.setTemperature(temperature);
1459 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
1460 Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx);
1461 Scalar energy_rate = getValue(h)*mass_rate;
1462 rate[Indices::contiEnergyEqIdx] += energy_rate;
1471 const bool compensateDrift = wellModel_.wellsActive();
1472 if (enableDriftCompensation_ && compensateDrift) {
1473 const auto& simulator = this->simulator();
1474 const auto& model = this->model();
1479 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
1480 Scalar poro = this->porosity(globalDofIdx, timeIdx);
1481 Scalar dt = simulator.timeStepSize();
1482 EqVector dofDriftRate = drift_[globalDofIdx];
1483 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
1486 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1487 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
1488 if (cnv > maxCompensation) {
1489 dofDriftRate[eqIdx] *= maxCompensation/cnv;
1493 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
1494 rate[eqIdx] -= dofDriftRate[eqIdx];
1504 {
return wellModel_; }
1506 WellModel& wellModel()
1507 {
return wellModel_; }
1509 const AquiferModel& aquiferModel()
const
1510 {
return aquiferModel_; }
1512 AquiferModel& mutableAquiferModel()
1513 {
return aquiferModel_; }
1516 const InitialFluidState& initialFluidState(
unsigned globalDofIdx)
const
1517 {
return initialFluidStates_[globalDofIdx]; }
1519 const EclipseIO& eclIO()
const
1520 {
return eclWriter_->eclIO(); }
1522 void setSubStepReport(
const SimulatorReportSingle& report)
1523 {
return eclWriter_->setSubStepReport(report); }
1525 void setSimulationReport(
const SimulatorReport& report)
1526 {
return eclWriter_->setSimulationReport(report); }
1528 bool nonTrivialBoundaryConditions()
const
1529 {
return nonTrivialBoundaryConditions_; }
1531 const InitialFluidState boundaryFluidState(
unsigned globalDofIdx,
const int directionId)
const
1533 OPM_TIMEBLOCK_LOCAL(boundaryFluidState);
1534 const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop;
1535 if (bcprop.size() > 0) {
1536 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1540 if (bcindex_(dir)[globalDofIdx] == 0)
1541 return initialFluidStates_[globalDofIdx];
1543 const auto& bc = bcprop[bcindex_(dir)[globalDofIdx]];
1544 if (bc.bctype == BCType::DIRICHLET )
1546 InitialFluidState fluidState;
1547 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1548 fluidState.setPvtRegionIndex(pvtRegionIdx);
1550 switch (bc.component) {
1551 case BCComponent::OIL:
1552 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
1553 throw std::logic_error(
"oil is not active and you're trying to add oil BC");
1555 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
1557 case BCComponent::GAS:
1558 if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1559 throw std::logic_error(
"gas is not active and you're trying to add gas BC");
1561 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
1563 case BCComponent::WATER:
1564 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1565 throw std::logic_error(
"water is not active and you're trying to add water BC");
1567 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
1569 case BCComponent::SOLVENT:
1570 case BCComponent::POLYMER:
1571 case BCComponent::NONE:
1572 throw std::logic_error(
"you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
1575 double pressure = initialFluidStates_[globalDofIdx].pressure(refPressurePhaseIdx_());
1576 const auto pressure_input = bc.pressure;
1577 if (pressure_input) {
1578 pressure = *pressure_input;
1581 std::array<Scalar, numPhases> pc = {0};
1582 const auto& matParams = materialLawParams(globalDofIdx);
1583 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
1584 Valgrind::CheckDefined(pressure);
1585 Valgrind::CheckDefined(pc);
1586 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1587 if (!FluidSystem::phaseIsActive(phaseIdx))
1590 if (Indices::oilEnabled)
1591 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1592 else if (Indices::gasEnabled)
1593 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1594 else if (Indices::waterEnabled)
1596 fluidState.setPressure(phaseIdx, pressure);
1599 double temperature = initialFluidStates_[globalDofIdx].temperature(0);
1600 const auto temperature_input = bc.temperature;
1601 if(temperature_input)
1602 temperature = *temperature_input;
1603 fluidState.setTemperature(temperature);
1605 if (FluidSystem::enableDissolvedGas()) {
1606 fluidState.setRs(0.0);
1607 fluidState.setRv(0.0);
1609 if (FluidSystem::enableDissolvedGasInWater()) {
1610 fluidState.setRsw(0.0);
1612 if (FluidSystem::enableVaporizedWater())
1613 fluidState.setRvw(0.0);
1615 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1616 if (!FluidSystem::phaseIsActive(phaseIdx))
1619 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
1620 fluidState.setInvB(phaseIdx, b);
1622 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
1623 fluidState.setDensity(phaseIdx, rho);
1625 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
1626 fluidState.setEnthalpy(phaseIdx, h);
1629 fluidState.checkDefined();
1633 return initialFluidStates_[globalDofIdx];
1644 OPM_TIMEBLOCK(nexTimeStepSize);
1646 if (this->nextTimeStepSize_ > 0.0)
1647 return this->nextTimeStepSize_;
1649 const auto& simulator = this->simulator();
1650 int episodeIdx = simulator.episodeIndex();
1654 return this->initialTimeStepSize_;
1659 const auto& newtonMethod = this->model().newtonMethod();
1660 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1668 template <
class LhsEval>
1671 OPM_TIMEBLOCK_LOCAL(rockCompPoroMultiplier);
1672 if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
1675 unsigned tableIdx = 0;
1676 if (!this->rockTableIdx_.empty())
1677 tableIdx = this->rockTableIdx_[elementIdx];
1679 const auto& fs = intQuants.fluidState();
1680 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1681 if (!this->minRefPressure_.empty())
1684 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1685 this->minRefPressure_[elementIdx]);
1687 if (!this->overburdenPressure_.empty())
1688 effectivePressure -= this->overburdenPressure_[elementIdx];
1691 if (!this->rockCompPoroMult_.empty()) {
1692 return this->rockCompPoroMult_[tableIdx].eval(effectivePressure,
true);
1696 assert(!this->rockCompPoroMultWc_.empty());
1697 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1698 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
1700 return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax,
true);
1708 template <
class LhsEval>
1711 bool implicit = !Parameters::get<TypeTag, Properties::ExplicitRockCompaction>();
1712 return implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<LhsEval>(intQuants, elementIdx)
1713 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1721 template <
class LhsEval>
1724 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
1725 if (!enableSaltPrecipitation)
1728 const auto& fs = intQuants.fluidState();
1729 unsigned tableIdx = fs.pvtRegionIndex();
1730 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
1731 porosityFactor = min(porosityFactor, 1.0);
1732 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
1733 return permfactTable.eval(porosityFactor,
true);
1739 template <
class LhsEval>
1742 OPM_TIMEBLOCK_LOCAL(wellTransMultiplier);
1744 bool implicit = !Parameters::get<TypeTag, Properties::ExplicitRockCompaction>();
1745 double trans_mult = implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<double>(intQuants, elementIdx)
1746 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1747 trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants);
1752 std::pair<BCType, RateVector> boundaryCondition(
const unsigned int globalSpaceIdx,
const int directionId)
const
1754 OPM_TIMEBLOCK_LOCAL(boundaryCondition);
1755 if (!nonTrivialBoundaryConditions_) {
1756 return { BCType::NONE, RateVector(0.0) };
1758 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1759 const auto& schedule = this->simulator().vanguard().schedule();
1760 if (bcindex_(dir)[globalSpaceIdx] == 0) {
1761 return { BCType::NONE, RateVector(0.0) };
1763 if (schedule[this->episodeIndex()].bcprop.size() == 0) {
1764 return { BCType::NONE, RateVector(0.0) };
1766 const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
1767 if (bc.bctype!=BCType::RATE) {
1768 return { bc.bctype, RateVector(0.0) };
1771 RateVector rate = 0.0;
1772 switch (bc.component) {
1773 case BCComponent::OIL:
1774 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = bc.rate;
1776 case BCComponent::GAS:
1777 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = bc.rate;
1779 case BCComponent::WATER:
1780 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate;
1782 case BCComponent::SOLVENT:
1783 if constexpr (!enableSolvent)
1784 throw std::logic_error(
"solvent is disabled and you're trying to add solvent to BC");
1786 rate[Indices::solventSaturationIdx] = bc.rate;
1788 case BCComponent::POLYMER:
1789 if constexpr (!enablePolymer)
1790 throw std::logic_error(
"polymer is disabled and you're trying to add polymer to BC");
1792 rate[Indices::polymerConcentrationIdx] = bc.rate;
1794 case BCComponent::NONE:
1795 throw std::logic_error(
"you need to specify the component when RATE type is set in BC");
1799 return {bc.bctype, rate};
1802 const std::unique_ptr<EclWriterType>& eclWriter()
const
1807 void setConvData(
const std::vector<std::vector<int>>& data)
1809 eclWriter_->mutableOutputModule().setCnvData(data);
1812 template<
class Serializer>
1813 void serializeOp(Serializer& serializer)
1815 serializer(
static_cast<BaseType&
>(*
this));
1817 serializer(wellModel_);
1818 serializer(aquiferModel_);
1819 serializer(tracerModel_);
1820 serializer(*materialLawManager_);
1821 serializer(*eclWriter_);
1824 Implementation& asImp_()
1825 {
return *
static_cast<Implementation *
>(
this); }
1827 void updateExplicitQuantities_()
1829 OPM_TIMEBLOCK(updateExplicitQuantities);
1830 const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
1831 const bool invalidateFromMinPressure = updateMinPressure_();
1834 const bool invalidateFromHyst = updateHysteresis_();
1835 const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
1838 bool invalidateIntensiveQuantities
1839 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
1840 if (invalidateIntensiveQuantities) {
1841 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1842 this->model().invalidateAndUpdateIntensiveQuantities(0);
1845 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
1846 updateMaxPolymerAdsorption_();
1848 updateRockCompTransMultVal_();
1851 template<
class UpdateFunc>
1852 void updateProperty_(
const std::string& failureMsg,
1855 OPM_TIMEBLOCK(updateProperty);
1856 const auto& model = this->simulator().model();
1857 const auto& primaryVars = model.solution(0);
1858 const auto& vanguard = this->simulator().vanguard();
1859 std::size_t numGridDof = primaryVars.size();
1860 OPM_BEGIN_PARALLEL_TRY_CATCH();
1862#pragma omp parallel for
1864 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1865 const auto& iq = *model.cachedIntensiveQuantities(dofIdx, 0);
1868 OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
1872 void updateCompositionChangeLimits_()
1874 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1877 int episodeIdx = this->episodeIndex();
1878 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1879 this->mixControls_.drsdtActive(episodeIdx),
1880 this->mixControls_.drvdtActive(episodeIdx)};
1881 if (!active[0] && !active[1] && !active[2]) {
1885 this->updateProperty_(
"FlowProblem::updateCompositionChangeLimits_()) failed:",
1886 [
this,episodeIdx,active](
unsigned compressedDofIdx,
1887 const IntensiveQuantities& iq)
1889 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1890 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1891 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1892 this->mixControls_.update(compressedDofIdx,
1895 this->gravity_[dim - 1],
1896 perm[dim - 1][dim - 1],
1904 bool updateMaxOilSaturation_()
1906 OPM_TIMEBLOCK(updateMaxOilSaturation);
1907 int episodeIdx = this->episodeIndex();
1910 if (this->vapparsActive(episodeIdx)) {
1911 this->updateProperty_(
"FlowProblem::updateMaxOilSaturation_() failed:",
1912 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1914 this->updateMaxOilSaturation_(compressedDofIdx,iq);
1922 bool updateMaxOilSaturation_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1924 OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
1925 const auto& fs = iq.fluidState();
1926 const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1927 auto& mos = this->maxOilSaturation_;
1928 if(mos[compressedDofIdx] < So){
1929 mos[compressedDofIdx] = So;
1936 bool updateMaxWaterSaturation_()
1938 OPM_TIMEBLOCK(updateMaxWaterSaturation);
1940 if (this->maxWaterSaturation_.empty())
1943 this->maxWaterSaturation_[1] = this->maxWaterSaturation_[0];
1944 this->updateProperty_(
"FlowProblem::updateMaxWaterSaturation_() failed:",
1945 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1947 this->updateMaxWaterSaturation_(compressedDofIdx,iq);
1953 bool updateMaxWaterSaturation_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1955 OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation);
1956 const auto& fs = iq.fluidState();
1957 const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
1958 auto& mow = this->maxWaterSaturation_;
1959 if(mow[compressedDofIdx]< Sw){
1960 mow[compressedDofIdx] = Sw;
1967 bool updateMinPressure_()
1969 OPM_TIMEBLOCK(updateMinPressure);
1971 if (this->minRefPressure_.empty())
1974 this->updateProperty_(
"FlowProblem::updateMinPressure_() failed:",
1975 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1977 this->updateMinPressure_(compressedDofIdx,iq);
1982 bool updateMinPressure_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq){
1983 OPM_TIMEBLOCK_LOCAL(updateMinPressure);
1984 const auto& fs = iq.fluidState();
1985 const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
1986 auto& min_pressures = this->minRefPressure_;
1987 if(min_pressures[compressedDofIdx]> min_pressure){
1988 min_pressures[compressedDofIdx] = min_pressure;
1999 std::function<std::vector<double>(
const FieldPropsManager&,
const std::string&)>
2000 fieldPropDoubleOnLeafAssigner_()
2002 const auto& lookup = this->lookUpData_;
2003 return [&lookup](
const FieldPropsManager& fieldPropManager,
const std::string& propString)
2005 return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
2013 template<
typename IntType>
2014 std::function<std::vector<IntType>(
const FieldPropsManager&,
const std::string&,
bool)>
2015 fieldPropIntTypeOnLeafAssigner_()
2017 const auto& lookup = this->lookUpData_;
2018 return [&lookup](
const FieldPropsManager& fieldPropManager,
const std::string& propString,
bool needsTranslation)
2020 return lookup.template assignFieldPropsIntOnLeaf<IntType>(fieldPropManager, propString, needsTranslation);
2024 void readMaterialParameters_()
2026 OPM_TIMEBLOCK(readMaterialParameters);
2027 const auto& simulator = this->simulator();
2028 const auto& vanguard = simulator.vanguard();
2029 const auto& eclState = vanguard.eclState();
2032 OPM_BEGIN_PARALLEL_TRY_CATCH();
2033 this->updatePvtnum_();
2034 this->updateSatnum_();
2037 this->updateMiscnum_();
2039 this->updatePlmixnum_();
2042 this->updateKrnum_();
2043 OPM_END_PARALLEL_TRY_CATCH(
"Invalid region numbers: ", vanguard.gridView().comm());
2046 updateReferencePorosity_();
2047 this->referencePorosity_[1] = this->referencePorosity_[0];
2052 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
2053 materialLawManager_->initFromState(eclState);
2054 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2055 this->
template fieldPropIntTypeOnLeafAssigner_<int>(),
2056 this-> lookupIdxOnLevelZeroAssigner_());
2060 void readThermalParameters_()
2062 if constexpr (enableEnergy)
2064 const auto& simulator = this->simulator();
2065 const auto& vanguard = simulator.vanguard();
2066 const auto& eclState = vanguard.eclState();
2069 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
2070 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2071 this-> fieldPropDoubleOnLeafAssigner_(),
2072 this->
template fieldPropIntTypeOnLeafAssigner_<unsigned int>());
2076 void updateReferencePorosity_()
2078 const auto& simulator = this->simulator();
2079 const auto& vanguard = simulator.vanguard();
2080 const auto& eclState = vanguard.eclState();
2082 std::size_t numDof = this->model().numGridDof();
2084 this->referencePorosity_[0].resize(numDof);
2086 const auto& fp = eclState.fieldProps();
2087 const std::vector<double> porvData =
this -> fieldPropDoubleOnLeafAssigner_()(fp,
"PORV");
2088 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2089 Scalar poreVolume = porvData[dofIdx];
2094 Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx);
2095 assert(dofVolume > 0.0);
2096 this->referencePorosity_[0][dofIdx] = poreVolume/dofVolume;
2100 void readInitialCondition_()
2102 const auto& simulator = this->simulator();
2103 const auto& vanguard = simulator.vanguard();
2104 const auto& eclState = vanguard.eclState();
2106 if (eclState.getInitConfig().hasEquil())
2107 readEquilInitialCondition_();
2109 readExplicitInitialCondition_();
2111 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
2112 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
2115 enablePolymerMolarWeight,
2119 std::size_t numElems = this->model().numGridDof();
2120 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2121 const auto& fs = initialFluidStates_[elemIdx];
2122 if (!this->maxWaterSaturation_.empty())
2123 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
2124 if (!this->maxOilSaturation_.empty())
2125 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
2126 if (!this->minRefPressure_.empty())
2127 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
2133 void readEquilInitialCondition_()
2135 const auto& simulator = this->simulator();
2138 EquilInitializer<TypeTag> equilInitializer(simulator, *materialLawManager_);
2140 std::size_t numElems = this->model().numGridDof();
2141 initialFluidStates_.resize(numElems);
2142 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2143 auto& elemFluidState = initialFluidStates_[elemIdx];
2144 elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
2148 void readEclRestartSolution_()
2151 if(this->simulator().vanguard().grid().maxLevel() > 0) {
2152 throw std::invalid_argument(
"Refined grids are not yet supported for restart ");
2156 auto& simulator = this->simulator();
2157 const auto& schedule = simulator.vanguard().schedule();
2158 const auto& eclState = simulator.vanguard().eclState();
2159 const auto& initconfig = eclState.getInitConfig();
2161 int restart_step = initconfig.getRestartStep();
2163 simulator.setTime(schedule.seconds(restart_step));
2165 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
2166 schedule.stepLength(restart_step));
2167 simulator.setEpisodeIndex(restart_step);
2169 eclWriter_->beginRestart();
2171 Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength());
2172 simulator.setTimeStepSize(dt);
2174 std::size_t numElems = this->model().numGridDof();
2175 initialFluidStates_.resize(numElems);
2176 if constexpr (enableSolvent) {
2177 this->solventSaturation_.resize(numElems, 0.0);
2178 this->solventRsw_.resize(numElems, 0.0);
2181 if constexpr (enablePolymer)
2182 this->polymer_.concentration.resize(numElems, 0.0);
2184 if constexpr (enablePolymerMolarWeight) {
2185 const std::string msg {
"Support of the RESTART for polymer molecular weight "
2186 "is not implemented yet. The polymer weight value will be "
2187 "zero when RESTART begins"};
2188 OpmLog::warning(
"NO_POLYMW_RESTART", msg);
2189 this->polymer_.moleWeight.resize(numElems, 0.0);
2192 if constexpr (enableMICP) {
2193 this->micp_.resize(numElems);
2196 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2197 auto& elemFluidState = initialFluidStates_[elemIdx];
2198 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
2199 eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
2200 eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
2209 auto ssol = enableSolvent
2210 ? eclWriter_->outputModule().getSolventSaturation(elemIdx)
2213 processRestartSaturations_(elemFluidState, ssol);
2215 if constexpr (enableSolvent) {
2216 this->solventSaturation_[elemIdx] = ssol;
2217 this->solventRsw_[elemIdx] = eclWriter_->outputModule().getSolventRsw(elemIdx);
2221 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
2223 if constexpr (enablePolymer)
2224 this->polymer_.concentration[elemIdx] = eclWriter_->outputModule().getPolymerConcentration(elemIdx);
2225 if constexpr (enableMICP){
2226 this->micp_.microbialConcentration[elemIdx] = eclWriter_->outputModule().getMicrobialConcentration(elemIdx);
2227 this->micp_.oxygenConcentration[elemIdx] = eclWriter_->outputModule().getOxygenConcentration(elemIdx);
2228 this->micp_.ureaConcentration[elemIdx] = eclWriter_->outputModule().getUreaConcentration(elemIdx);
2229 this->micp_.biofilmConcentration[elemIdx] = eclWriter_->outputModule().getBiofilmConcentration(elemIdx);
2230 this->micp_.calciteConcentration[elemIdx] = eclWriter_->outputModule().getCalciteConcentration(elemIdx);
2235 const int episodeIdx = this->episodeIndex();
2236 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
2241 auto& sol = this->model().solution(0);
2242 const auto& gridView = this->gridView();
2243 ElementContext elemCtx(simulator);
2244 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2245 elemCtx.updatePrimaryStencil(elem);
2246 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
2247 initial(sol[elemIdx], elemCtx, 0, 0);
2255 this->model().syncOverlap();
2257 eclWriter_->endRestart();
2260 void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
2264 const Scalar smallSaturationTolerance = 1.e-6;
2265 Scalar sumSaturation = 0.0;
2266 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2267 if (FluidSystem::phaseIsActive(phaseIdx)) {
2268 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
2269 elemFluidState.setSaturation(phaseIdx, 0.0);
2271 sumSaturation += elemFluidState.saturation(phaseIdx);
2275 if constexpr (enableSolvent) {
2276 if (solventSaturation < smallSaturationTolerance)
2277 solventSaturation = 0.0;
2279 sumSaturation += solventSaturation;
2282 assert(sumSaturation > 0.0);
2284 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2285 if (FluidSystem::phaseIsActive(phaseIdx)) {
2286 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
2287 elemFluidState.setSaturation(phaseIdx, saturation);
2290 if constexpr (enableSolvent) {
2291 solventSaturation = solventSaturation / sumSaturation;
2295 void readExplicitInitialCondition_()
2297 const auto& simulator = this->simulator();
2298 const auto& vanguard = simulator.vanguard();
2299 const auto& eclState = vanguard.eclState();
2300 const auto& fp = eclState.fieldProps();
2301 bool has_swat = fp.has_double(
"SWAT");
2302 bool has_sgas = fp.has_double(
"SGAS");
2303 bool has_rs = fp.has_double(
"RS");
2304 bool has_rv = fp.has_double(
"RV");
2305 bool has_rvw = fp.has_double(
"RVW");
2306 bool has_pressure = fp.has_double(
"PRESSURE");
2307 bool has_salt = fp.has_double(
"SALT");
2308 bool has_saltp = fp.has_double(
"SALTP");
2311 if (Indices::numPhases > 1) {
2312 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
2313 throw std::runtime_error(
"The ECL input file requires the presence of the SWAT keyword if "
2314 "the water phase is active");
2315 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
2316 throw std::runtime_error(
"The ECL input file requires the presence of the SGAS keyword if "
2317 "the gas phase is active");
2320 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
2321 "keyword if the model is initialized explicitly");
2322 if (FluidSystem::enableDissolvedGas() && !has_rs)
2323 throw std::runtime_error(
"The ECL input file requires the RS keyword to be present if"
2324 " dissolved gas is enabled");
2325 if (FluidSystem::enableVaporizedOil() && !has_rv)
2326 throw std::runtime_error(
"The ECL input file requires the RV keyword to be present if"
2327 " vaporized oil is enabled");
2328 if (FluidSystem::enableVaporizedWater() && !has_rvw)
2329 throw std::runtime_error(
"The ECL input file requires the RVW keyword to be present if"
2330 " vaporized water is enabled");
2331 if (enableBrine && !has_salt)
2332 throw std::runtime_error(
"The ECL input file requires the SALT keyword to be present if"
2333 " brine is enabled and the model is initialized explicitly");
2334 if (enableSaltPrecipitation && !has_saltp)
2335 throw std::runtime_error(
"The ECL input file requires the SALTP keyword to be present if"
2336 " salt precipitation is enabled and the model is initialized explicitly");
2338 std::size_t numDof = this->model().numGridDof();
2340 initialFluidStates_.resize(numDof);
2342 std::vector<double> waterSaturationData;
2343 std::vector<double> gasSaturationData;
2344 std::vector<double> pressureData;
2345 std::vector<double> rsData;
2346 std::vector<double> rvData;
2347 std::vector<double> rvwData;
2348 std::vector<double> tempiData;
2349 std::vector<double> saltData;
2350 std::vector<double> saltpData;
2352 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
2353 waterSaturationData = fp.get_double(
"SWAT");
2355 waterSaturationData.resize(numDof);
2357 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
2358 gasSaturationData = fp.get_double(
"SGAS");
2360 gasSaturationData.resize(numDof);
2362 pressureData = fp.get_double(
"PRESSURE");
2363 if (FluidSystem::enableDissolvedGas())
2364 rsData = fp.get_double(
"RS");
2366 if (FluidSystem::enableVaporizedOil())
2367 rvData = fp.get_double(
"RV");
2369 if (FluidSystem::enableVaporizedWater())
2370 rvwData = fp.get_double(
"RVW");
2373 tempiData = fp.get_double(
"TEMPI");
2376 if constexpr (enableBrine)
2377 saltData = fp.get_double(
"SALT");
2380 if constexpr (enableSaltPrecipitation)
2381 saltpData = fp.get_double(
"SALTP");
2384 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2385 auto& dofFluidState = initialFluidStates_[dofIdx];
2387 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
2392 Scalar temperatureLoc = tempiData[dofIdx];
2393 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
2394 temperatureLoc = FluidSystem::surfaceTemperature;
2395 dofFluidState.setTemperature(temperatureLoc);
2400 if constexpr (enableBrine)
2401 dofFluidState.setSaltConcentration(saltData[dofIdx]);
2406 if constexpr (enableSaltPrecipitation)
2407 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
2412 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
2413 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
2414 waterSaturationData[dofIdx]);
2416 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
2417 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
2418 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2420 - waterSaturationData[dofIdx]);
2423 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2424 gasSaturationData[dofIdx]);
2426 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
2427 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
2429 - waterSaturationData[dofIdx]
2430 - gasSaturationData[dofIdx]);
2435 Scalar pressure = pressureData[dofIdx];
2439 std::array<Scalar, numPhases> pc = {0};
2440 const auto& matParams = materialLawParams(dofIdx);
2441 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
2442 Valgrind::CheckDefined(pressure);
2443 Valgrind::CheckDefined(pc);
2444 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2445 if (!FluidSystem::phaseIsActive(phaseIdx))
2448 if (Indices::oilEnabled)
2449 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
2450 else if (Indices::gasEnabled)
2451 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
2452 else if (Indices::waterEnabled)
2454 dofFluidState.setPressure(phaseIdx, pressure);
2457 if (FluidSystem::enableDissolvedGas())
2458 dofFluidState.setRs(rsData[dofIdx]);
2459 else if (Indices::gasEnabled && Indices::oilEnabled)
2460 dofFluidState.setRs(0.0);
2462 if (FluidSystem::enableVaporizedOil())
2463 dofFluidState.setRv(rvData[dofIdx]);
2464 else if (Indices::gasEnabled && Indices::oilEnabled)
2465 dofFluidState.setRv(0.0);
2467 if (FluidSystem::enableVaporizedWater())
2468 dofFluidState.setRvw(rvwData[dofIdx]);
2473 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2474 if (!FluidSystem::phaseIsActive(phaseIdx))
2477 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2478 dofFluidState.setInvB(phaseIdx, b);
2480 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2481 dofFluidState.setDensity(phaseIdx, rho);
2488 bool updateHysteresis_()
2490 if (!materialLawManager_->enableHysteresis())
2495 this->updateProperty_(
"FlowProblem::updateHysteresis_() failed:",
2496 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2498 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2504 bool updateHysteresis_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2506 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
2507 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2512 void updateMaxPolymerAdsorption_()
2515 this->updateProperty_(
"FlowProblem::updateMaxPolymerAdsorption_() failed:",
2516 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2518 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
2522 bool updateMaxPolymerAdsorption_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
2524 const Scalar pa = scalarValue(iq.polymerAdsorption());
2525 auto& mpa = this->polymer_.maxAdsorption;
2526 if (mpa[compressedDofIdx] < pa) {
2527 mpa[compressedDofIdx] = pa;
2534 Scalar getRockCompTransMultVal(std::size_t dofIdx)
const
2536 if (this->rockCompTransMultVal_.empty())
2539 return this->rockCompTransMultVal_[dofIdx];
2546 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
2547 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
2548 ConditionalStorage<enableDiffusion, Scalar> diffusivity;
2549 ConditionalStorage<enableDispersion, Scalar> dispersivity;
2550 Scalar transmissibility;
2554 void updatePffDofData_()
2556 const auto& distFn =
2557 [
this](PffDofData_& dofData,
2558 const Stencil& stencil,
2559 unsigned localDofIdx)
2562 const auto& elementMapper = this->model().elementMapper();
2564 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
2565 if (localDofIdx != 0) {
2566 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(0));
2567 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
2569 if constexpr (enableEnergy) {
2570 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
2571 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
2573 if constexpr (enableDiffusion)
2574 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
2575 if (enableDispersion)
2576 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
2580 pffDofData_.update(distFn);
2583 void readBoundaryConditions_()
2585 const auto& simulator = this->simulator();
2586 const auto& vanguard = simulator.vanguard();
2587 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
2588 if (bcconfig.size() > 0) {
2589 nonTrivialBoundaryConditions_ =
true;
2591 std::size_t numCartDof = vanguard.cartesianSize();
2592 unsigned numElems = vanguard.gridView().size(0);
2593 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
2595 for (
unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
2596 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
2598 bcindex_.resize(numElems, 0);
2599 auto loopAndApply = [&cartesianToCompressedElemIdx,
2600 &vanguard](
const auto& bcface,
2603 for (
int i = bcface.i1; i <= bcface.i2; ++i) {
2604 for (
int j = bcface.j1; j <= bcface.j2; ++j) {
2605 for (
int k = bcface.k1; k <= bcface.k2; ++k) {
2606 std::array<int, 3> tmp = {i,j,k};
2607 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
2614 for (
const auto& bcface : bcconfig) {
2615 std::vector<int>& data = bcindex_(bcface.dir);
2616 const int index = bcface.index;
2617 loopAndApply(bcface,
2618 [&data,index](
int elemIdx)
2619 { data[elemIdx] = index; });
2626 Scalar limitNextTimeStepSize_(Scalar dtNext)
const
2628 if constexpr (enableExperiments) {
2629 const auto& simulator = this->simulator();
2630 const auto& schedule = simulator.vanguard().schedule();
2631 int episodeIdx = simulator.episodeIndex();
2634 Scalar maxTimeStepSize = Parameters::get<TypeTag, Properties::SolverMaxTimeStepInDays>() * 24 * 60 * 60;
2635 int reportStepIdx = std::max(episodeIdx, 0);
2636 if (this->enableTuning_) {
2637 const auto& tuning = schedule[reportStepIdx].tuning();
2638 maxTimeStepSize = tuning.TSMAXZ;
2641 dtNext = std::min(dtNext, maxTimeStepSize);
2643 Scalar remainingEpisodeTime =
2644 simulator.episodeStartTime() + simulator.episodeLength()
2645 - (simulator.startTime() + simulator.time());
2646 assert(remainingEpisodeTime >= 0.0);
2650 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
2653 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
2655 if (simulator.episodeStarts()) {
2658 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
2659 bool wellEventOccured =
2660 events.hasEvent(ScheduleEvents::NEW_WELL)
2661 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
2662 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
2663 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
2664 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
2665 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
2672 void computeAndSetEqWeights_()
2674 std::vector<Scalar> sumInvB(numPhases, 0.0);
2675 const auto& gridView = this->gridView();
2676 ElementContext elemCtx(this->simulator());
2677 for(
const auto& elem: elements(gridView, Dune::Partitions::interior)) {
2678 elemCtx.updatePrimaryStencil(elem);
2679 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
2680 const auto& dofFluidState = initialFluidStates_[elemIdx];
2681 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2682 if (!FluidSystem::phaseIsActive(phaseIdx))
2685 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
2689 std::size_t numDof = this->model().numGridDof();
2690 const auto& comm = this->simulator().vanguard().grid().comm();
2691 comm.sum(sumInvB.data(),sumInvB.size());
2692 Scalar numTotalDof = comm.sum(numDof);
2694 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2695 if (!FluidSystem::phaseIsActive(phaseIdx))
2698 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
2699 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
2700 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
2701 this->model().setEqWeight(activeSolventCompIdx, avgB);
2705 int refPressurePhaseIdx_()
const {
2706 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2709 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2713 return waterPhaseIdx;
2717 void updateRockCompTransMultVal_()
2719 const auto& model = this->simulator().model();
2720 std::size_t numGridDof = this->model().numGridDof();
2721 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
2722 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
2723 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, 0);
2724 Scalar trans_mult = computeRockCompTransMultiplier_<Scalar>(iq, elementIdx);
2725 this->rockCompTransMultVal_[elementIdx] = trans_mult;
2734 template <
class LhsEval>
2735 LhsEval computeRockCompTransMultiplier_(
const IntensiveQuantities& intQuants,
unsigned elementIdx)
const
2737 OPM_TIMEBLOCK_LOCAL(computeRockCompTransMultiplier);
2738 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
2741 unsigned tableIdx = 0;
2742 if (!this->rockTableIdx_.empty())
2743 tableIdx = this->rockTableIdx_[elementIdx];
2745 const auto& fs = intQuants.fluidState();
2746 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
2748 if (!this->minRefPressure_.empty())
2751 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
2752 this->minRefPressure_[elementIdx]);
2754 if (!this->overburdenPressure_.empty())
2755 effectivePressure -= this->overburdenPressure_[elementIdx];
2757 if (!this->rockCompTransMult_.empty())
2758 return this->rockCompTransMult_[tableIdx].eval(effectivePressure,
true);
2761 assert(!this->rockCompTransMultWc_.empty());
2762 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
2763 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
2765 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax,
true);
2768 typename Vanguard::TransmissibilityType transmissibilities_;
2770 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
2771 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
2773 FlowThresholdPressure<TypeTag> thresholdPressures_;
2775 std::vector<InitialFluidState> initialFluidStates_;
2777 bool enableDriftCompensation_;
2778 GlobalEqVector drift_;
2780 WellModel wellModel_;
2781 AquiferModel aquiferModel_;
2783 bool enableEclOutput_;
2784 std::unique_ptr<EclWriterType> eclWriter_;
2787 bool enableDamarisOutput_ = false ;
2788 std::unique_ptr<DamarisWriterType> damarisWriter_;
2791 PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
2792 TracerModel tracerModel_;
2794 ActionHandler actionHandler_;
2799 std::array<std::vector<T>,6> data;
2801 void resize(std::size_t size, T defVal)
2803 for (
auto& d : data)
2804 d.resize(size, defVal);
2807 const std::vector<T>& operator()(FaceDir::DirEnum dir)
const
2809 if (dir == FaceDir::DirEnum::Unknown)
2810 throw std::runtime_error(
"Tried to access BC data for the 'Unknown' direction");
2812 int div =
static_cast<int>(dir);
2813 while ((div /= 2) >= 1)
2815 assert(idx >= 0 && idx <= 5);
2819 std::vector<T>& operator()(FaceDir::DirEnum dir)
2821 return const_cast<std::vector<T>&
>(std::as_const(*
this)(dir));
2825 BCData<int> bcindex_;
2826 bool nonTrivialBoundaryConditions_ =
false;