91 GetPropType<TypeTag, Properties::FluidSystem>>
108 enum { dim = GridView::dimension };
109 enum { dimWorld = GridView::dimensionworld };
113 enum { numPhases = FluidSystem::numPhases };
114 enum { numComponents = FluidSystem::numComponents };
132 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
133 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
134 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
138 enum { gasCompIdx = FluidSystem::gasCompIdx };
139 enum { oilCompIdx = FluidSystem::oilCompIdx };
140 enum { waterCompIdx = FluidSystem::waterCompIdx };
145 using Element =
typename GridView::template Codim<0>::Entity;
149 using MaterialLawParams =
typename EclMaterialLawManager::MaterialLawParams;
150 using SolidEnergyLawParams =
typename EclThermalLawManager::SolidEnergyLawParams;
151 using ThermalConductionLawParams =
typename EclThermalLawManager::ThermalConductionLawParams;
160 using Toolbox = MathToolbox<Evaluation>;
161 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
165 using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
181 ParentType::registerParameters();
183 registerFlowProblemParameters<Scalar>();
191 const std::string&)> addKey,
192 std::set<std::string>& seenParams,
193 std::string& errorMsg,
199 return detail::eclPositionalParameter(addKey,
210 : ParentType(simulator)
211 ,
BaseType(simulator.vanguard().eclState(),
212 simulator.vanguard().schedule(),
213 simulator.vanguard().gridView())
214 , transmissibilities_(simulator.vanguard().eclState(),
215 simulator.vanguard().gridView(),
216 simulator.vanguard().cartesianIndexMapper(),
217 simulator.vanguard().grid(),
218 simulator.vanguard().cellCentroids(),
222 , wellModel_(simulator)
223 , aquiferModel_(simulator)
224 , pffDofData_(simulator.gridView(), this->elementMapper())
225 , tracerModel_(simulator)
227 const auto& vanguard = simulator.vanguard();
229 enableDriftCompensation_ = Parameters::Get<Parameters::EnableDriftCompensation>();
231 enableVtkOutput_ = Parameters::Get<Parameters::EnableVtkOutput>();
233 this->enableTuning_ = Parameters::Get<Parameters::EnableTuning>();
234 this->initialTimeStepSize_ = Parameters::Get<Parameters::InitialTimeStepSize<Scalar>>();
235 this->maxTimeStepAfterWellEvent_ = Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60;
241 if (Parameters::IsSet<Parameters::NumPressurePointsEquil>())
243 this->numPressurePointsEquil_ = Parameters::Get<Parameters::NumPressurePointsEquil>();
245 this->numPressurePointsEquil_ = simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP();
248 explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
252 relpermDiagnostics.
diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
257 void prefetch(
const Element& elem)
const
258 { pffDofData_.prefetch(elem); }
271 template <
class Restarter>
278 wellModel_.deserialize(res);
281 aquiferModel_.deserialize(res);
290 template <
class Restarter>
293 wellModel_.serialize(res);
295 aquiferModel_.serialize(res);
298 int episodeIndex()
const
300 return std::max(this->simulator().episodeIndex(), 0);
310 auto& simulator = this->simulator();
311 int episodeIdx = simulator.episodeIndex();
312 auto& eclState = simulator.vanguard().eclState();
313 const auto& schedule = simulator.vanguard().schedule();
314 const auto& events = schedule[episodeIdx].events();
316 if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
323 const auto& miniDeck = schedule[episodeIdx].geo_keywords();
324 const auto& cc = simulator.vanguard().grid().comm();
325 eclState.apply_schedule_keywords( miniDeck );
326 eclBroadcast(cc, eclState.getTransMult() );
329 std::function<
unsigned int(
unsigned int)> equilGridToGrid = [&simulator](
unsigned int i) {
330 return simulator.vanguard().gridEquilIdxToGridIdx(i);
334 using TransUpdateQuantities =
typename Vanguard::TransmissibilityType::TransUpdateQuantities;
335 transmissibilities_.update(
true, TransUpdateQuantities::All, equilGridToGrid);
336 this->referencePorosity_[1] = this->referencePorosity_[0];
337 updateReferencePorosity_();
339 this->model().linearizer().updateDiscretizationParameters();
342 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
345 wellModel_.beginEpisode();
348 aquiferModel_.beginEpisode();
351 Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
353 if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
356 dt = std::min(dt, this->initialTimeStepSize_);
357 simulator.setTimeStepSize(dt);
366 const int episodeIdx = this->episodeIndex();
367 const int timeStepSize = this->simulator().timeStepSize();
369 this->beginTimeStep_(enableExperiments,
371 this->simulator().timeStepIndex(),
372 this->simulator().startTime(),
373 this->simulator().time(),
375 this->simulator().endTime());
380 this->updateExplicitQuantities_(episodeIdx, timeStepSize, first_step_ && (episodeIdx > 0));
383 if (nonTrivialBoundaryConditions()) {
384 this->model().linearizer().updateBoundaryConditionData();
387 wellModel_.beginTimeStep();
388 aquiferModel_.beginTimeStep();
389 tracerModel_.beginTimeStep();
399 wellModel_.beginIteration();
400 aquiferModel_.beginIteration();
409 wellModel_.endIteration();
410 aquiferModel_.endIteration();
427 const int rank = this->simulator().gridView().comm().rank();
429 std::cout <<
"checking conservativeness of solution\n";
432 this->model().checkConservativeness(-1,
true);
434 std::cout <<
"solution is sufficiently conservative\n";
439 auto& simulator = this->simulator();
440 simulator.setTimeStepIndex(simulator.timeStepIndex()+1);
442 this->wellModel_.endTimeStep();
443 this->aquiferModel_.endTimeStep();
444 this->tracerModel_.endTimeStep();
447 this->model().linearizer().updateFlowsInfo();
449 if (this->enableDriftCompensation_) {
450 OPM_TIMEBLOCK(driftCompansation);
452 const auto& residual = this->model().linearizer().residual();
454 for (
unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
455 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(globalDofIdx);
456 this->drift_[sfcdofIdx] = residual[sfcdofIdx] * simulator.timeStepSize();
459 this->drift_[sfcdofIdx] *= this->model().dofTotalVolume(sfcdofIdx);
470 const int episodeIdx = this->episodeIndex();
472 this->wellModel_.endEpisode();
473 this->aquiferModel_.endEpisode();
475 const auto& schedule = this->simulator().vanguard().schedule();
478 if (episodeIdx + 1 >=
static_cast<int>(schedule.size()) - 1) {
479 this->simulator().setFinished(
true);
484 this->simulator().startNextEpisode(schedule.stepLength(episodeIdx + 1));
493 OPM_TIMEBLOCK(problemWriteOutput);
496 if (Parameters::Get<Parameters::EnableWriteAllSolutions>() ||
497 this->simulator().episodeWillBeOver()) {
498 ParentType::writeOutput(verbose);
505 template <
class Context>
508 unsigned timeIdx)
const
510 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
511 return transmissibilities_.permeability(globalSpaceIdx);
521 {
return transmissibilities_.permeability(globalElemIdx); }
526 template <
class Context>
528 [[maybe_unused]]
unsigned fromDofLocalIdx,
529 unsigned toDofLocalIdx)
const
531 assert(fromDofLocalIdx == 0);
532 return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
540 return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
546 template <
class Context>
548 [[maybe_unused]]
unsigned fromDofLocalIdx,
549 unsigned toDofLocalIdx)
const
551 assert(fromDofLocalIdx == 0);
552 return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
558 Scalar
diffusivity(
const unsigned globalCellIn,
const unsigned globalCellOut)
const{
559 return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
565 Scalar
dispersivity(
const unsigned globalCellIn,
const unsigned globalCellOut)
const{
566 return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
573 const unsigned boundaryFaceIdx)
const
575 return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
584 template <
class Context>
586 unsigned boundaryFaceIdx)
const
588 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
589 return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
596 const unsigned boundaryFaceIdx)
const
598 return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
606 const unsigned globalSpaceIdxOut)
const
608 return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
614 template <
class Context>
617 unsigned timeIdx)
const
619 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
620 unsigned toDofLocalIdx = face.exteriorIndex();
621 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
627 template <
class Context>
630 unsigned timeIdx)
const
632 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
633 unsigned toDofLocalIdx = face.exteriorIndex();
634 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
640 template <
class Context>
642 unsigned boundaryFaceIdx)
const
644 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
645 return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
652 {
return transmissibilities_; }
656 {
return tracerModel_; }
658 TracerModel& tracerModel()
659 {
return tracerModel_; }
669 template <
class Context>
670 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
672 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
673 return this->
porosity(globalSpaceIdx, timeIdx);
682 template <
class Context>
683 Scalar
dofCenterDepth(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
685 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
697 return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
703 template <
class Context>
706 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
713 template <
class Context>
716 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
723 template <
class Context>
725 unsigned spaceIdx,
unsigned timeIdx)
const
727 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
733 return materialLawManager_->materialLawParams(globalDofIdx);
736 const MaterialLawParams&
materialLawParams(
unsigned globalDofIdx, FaceDir::DirEnum facedir)
const
738 return materialLawManager_->materialLawParams(globalDofIdx, facedir);
744 template <
class Context>
745 const SolidEnergyLawParams&
748 unsigned timeIdx)
const
750 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
751 return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
757 template <
class Context>
758 const ThermalConductionLawParams &
761 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
762 return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
772 {
return materialLawManager_; }
774 template <
class Flu
idState>
776 std::array<Evaluation,numPhases> &mobility,
777 DirectionalMobilityPtr &dirMob,
778 FluidState &fluidState,
779 unsigned globalSpaceIdx)
const
781 OPM_TIMEBLOCK_LOCAL(updateRelperms);
786 MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
787 Valgrind::CheckDefined(mobility);
789 if (materialLawManager_->hasDirectionalRelperms()
790 || materialLawManager_->hasDirectionalImbnum())
792 using Dir = FaceDir::DirEnum;
793 constexpr int ndim = 3;
794 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
795 Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
796 for (
int i = 0; i<ndim; i++) {
798 auto& mob_array = dirMob->getArray(i);
799 MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
808 {
return materialLawManager_; }
814 template <
class Context>
815 unsigned pvtRegionIndex(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
816 {
return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
822 template <
class Context>
830 template <
class Context>
838 template <
class Context>
847 template <
class Context>
855 {
return this->simulator().vanguard().caseName(); }
860 template <
class Context>
861 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
865 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
866 return asImp_().initialFluidState(globalDofIdx).temperature(0);
870 Scalar
temperature(
unsigned globalDofIdx,
unsigned )
const
874 return asImp_().initialFluidState(globalDofIdx).temperature(0);
877 const SolidEnergyLawParams&
881 return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
883 const ThermalConductionLawParams &
887 return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
901 if (!this->vapparsActive(this->episodeIndex()))
904 return this->maxOilSaturation_[globalDofIdx];
918 if (!this->vapparsActive(this->episodeIndex()))
921 this->maxOilSaturation_[globalDofIdx] = value;
930 this->model().invalidateAndUpdateIntensiveQuantities(0);
936 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
937 this->model().invalidateAndUpdateIntensiveQuantities(1);
945 aquiferModel_.initialSolutionApplied();
947 const bool invalidateFromHyst = updateHysteresis_();
948 if (invalidateFromHyst) {
949 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
950 this->model().invalidateAndUpdateIntensiveQuantities(0);
959 template <
class Context>
961 const Context& context,
963 unsigned timeIdx)
const
965 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
966 source(rate, globalDofIdx, timeIdx);
969 void source(RateVector& rate,
970 unsigned globalDofIdx,
971 unsigned timeIdx)
const
973 OPM_TIMEBLOCK_LOCAL(eclProblemSource);
977 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
981 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
982 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
984 Valgrind::CheckDefined(rate[eqIdx]);
985 assert(isfinite(rate[eqIdx]));
989 addToSourceDense(rate, globalDofIdx, timeIdx);
992 virtual void addToSourceDense(RateVector& rate,
993 unsigned globalDofIdx,
994 unsigned timeIdx)
const = 0;
1002 {
return wellModel_; }
1005 {
return wellModel_; }
1007 const AquiferModel& aquiferModel()
const
1008 {
return aquiferModel_; }
1010 AquiferModel& mutableAquiferModel()
1011 {
return aquiferModel_; }
1013 bool nonTrivialBoundaryConditions()
const
1014 {
return nonTrivialBoundaryConditions_; }
1024 OPM_TIMEBLOCK(nexTimeStepSize);
1026 if (this->nextTimeStepSize_ > 0.0)
1027 return this->nextTimeStepSize_;
1029 const auto& simulator = this->simulator();
1030 int episodeIdx = simulator.episodeIndex();
1034 return this->initialTimeStepSize_;
1039 const auto& newtonMethod = this->model().newtonMethod();
1040 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1048 template <
class LhsEval>
1052 if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
1055 unsigned tableIdx = 0;
1056 if (!this->rockTableIdx_.empty())
1057 tableIdx = this->rockTableIdx_[elementIdx];
1059 const auto& fs = intQuants.fluidState();
1060 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1061 if (!this->minRefPressure_.empty())
1064 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1065 this->minRefPressure_[elementIdx]);
1067 if (!this->overburdenPressure_.empty())
1068 effectivePressure -= this->overburdenPressure_[elementIdx];
1071 if (!this->rockCompPoroMult_.empty()) {
1072 return this->rockCompPoroMult_[tableIdx].eval(effectivePressure,
true);
1076 assert(!this->rockCompPoroMultWc_.empty());
1077 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1078 LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
1080 return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax,
true);
1088 template <
class LhsEval>
1091 const bool implicit = !this->explicitRockCompaction_;
1093 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1100 template <
class LhsEval>
1105 const bool implicit = !this->explicitRockCompaction_;
1107 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1108 trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants);
1113 std::pair<BCType, RateVector> boundaryCondition(
const unsigned int globalSpaceIdx,
const int directionId)
const
1115 OPM_TIMEBLOCK_LOCAL(boundaryCondition);
1116 if (!nonTrivialBoundaryConditions_) {
1117 return { BCType::NONE, RateVector(0.0) };
1119 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1120 const auto& schedule = this->simulator().vanguard().schedule();
1121 if (bcindex_(dir)[globalSpaceIdx] == 0) {
1122 return { BCType::NONE, RateVector(0.0) };
1124 if (schedule[this->episodeIndex()].bcprop.size() == 0) {
1125 return { BCType::NONE, RateVector(0.0) };
1127 const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
1128 if (bc.bctype!=BCType::RATE) {
1129 return { bc.bctype, RateVector(0.0) };
1132 RateVector rate = 0.0;
1133 switch (bc.component) {
1134 case BCComponent::OIL:
1135 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = bc.rate;
1137 case BCComponent::GAS:
1138 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = bc.rate;
1140 case BCComponent::WATER:
1141 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate;
1143 case BCComponent::SOLVENT:
1144 this->handleSolventBC(bc, rate);
1146 case BCComponent::POLYMER:
1147 this->handlePolymerBC(bc, rate);
1149 case BCComponent::NONE:
1150 throw std::logic_error(
"you need to specify the component when RATE type is set in BC");
1154 return {bc.bctype, rate};
1158 template<
class Serializer>
1159 void serializeOp(Serializer& serializer)
1161 serializer(
static_cast<BaseType&
>(*
this));
1163 serializer(wellModel_);
1164 serializer(aquiferModel_);
1165 serializer(tracerModel_);
1166 serializer(*materialLawManager_);
1170 Implementation& asImp_()
1171 {
return *
static_cast<Implementation *
>(
this); }
1173 const Implementation& asImp_()
const
1174 {
return *
static_cast<const Implementation *
>(
this); }
1177 template<
class UpdateFunc>
1178 void updateProperty_(
const std::string& failureMsg,
1181 OPM_TIMEBLOCK(updateProperty);
1182 const auto& model = this->simulator().model();
1183 const auto& primaryVars = model.solution(0);
1184 const auto& vanguard = this->simulator().vanguard();
1185 std::size_t numGridDof = primaryVars.size();
1186 OPM_BEGIN_PARALLEL_TRY_CATCH();
1188#pragma omp parallel for
1190 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1191 const auto& iq = *model.cachedIntensiveQuantities(dofIdx, 0);
1194 OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
1197 bool updateMaxOilSaturation_()
1199 OPM_TIMEBLOCK(updateMaxOilSaturation);
1200 int episodeIdx = this->episodeIndex();
1203 if (this->vapparsActive(episodeIdx)) {
1204 this->updateProperty_(
"FlowProblem::updateMaxOilSaturation_() failed:",
1205 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1207 this->updateMaxOilSaturation_(compressedDofIdx,iq);
1215 bool updateMaxOilSaturation_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1217 OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
1218 const auto& fs = iq.fluidState();
1219 const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1220 auto& mos = this->maxOilSaturation_;
1221 if(mos[compressedDofIdx] < So){
1222 mos[compressedDofIdx] = So;
1229 bool updateMaxWaterSaturation_()
1231 OPM_TIMEBLOCK(updateMaxWaterSaturation);
1233 if (this->maxWaterSaturation_.empty())
1236 this->maxWaterSaturation_[1] = this->maxWaterSaturation_[0];
1237 this->updateProperty_(
"FlowProblem::updateMaxWaterSaturation_() failed:",
1238 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1240 this->updateMaxWaterSaturation_(compressedDofIdx,iq);
1246 bool updateMaxWaterSaturation_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1248 OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation);
1249 const auto& fs = iq.fluidState();
1250 const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
1251 auto& mow = this->maxWaterSaturation_;
1252 if(mow[compressedDofIdx]< Sw){
1253 mow[compressedDofIdx] = Sw;
1260 bool updateMinPressure_()
1262 OPM_TIMEBLOCK(updateMinPressure);
1264 if (this->minRefPressure_.empty())
1267 this->updateProperty_(
"FlowProblem::updateMinPressure_() failed:",
1268 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1270 this->updateMinPressure_(compressedDofIdx,iq);
1275 bool updateMinPressure_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq){
1276 OPM_TIMEBLOCK_LOCAL(updateMinPressure);
1277 const auto& fs = iq.fluidState();
1278 const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
1279 auto& min_pressures = this->minRefPressure_;
1280 if(min_pressures[compressedDofIdx]> min_pressure){
1281 min_pressures[compressedDofIdx] = min_pressure;
1292 std::function<std::vector<double>(
const FieldPropsManager&,
const std::string&)>
1293 fieldPropDoubleOnLeafAssigner_()
1295 const auto& lookup = this->lookUpData_;
1296 return [&lookup](
const FieldPropsManager& fieldPropManager,
const std::string& propString)
1298 return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
1306 template<
typename IntType>
1307 std::function<std::vector<IntType>(
const FieldPropsManager&,
const std::string&,
bool)>
1308 fieldPropIntTypeOnLeafAssigner_()
1310 const auto& lookup = this->lookUpData_;
1311 return [&lookup](
const FieldPropsManager& fieldPropManager,
const std::string& propString,
bool needsTranslation)
1313 return lookup.template assignFieldPropsIntOnLeaf<IntType>(fieldPropManager, propString, needsTranslation);
1317 void readMaterialParameters_()
1319 OPM_TIMEBLOCK(readMaterialParameters);
1320 const auto& simulator = this->simulator();
1321 const auto& vanguard = simulator.vanguard();
1322 const auto& eclState = vanguard.eclState();
1325 OPM_BEGIN_PARALLEL_TRY_CATCH();
1326 this->updatePvtnum_();
1327 this->updateSatnum_();
1330 this->updateMiscnum_();
1332 this->updatePlmixnum_();
1334 OPM_END_PARALLEL_TRY_CATCH(
"Invalid region numbers: ", vanguard.gridView().comm());
1337 updateReferencePorosity_();
1338 this->referencePorosity_[1] = this->referencePorosity_[0];
1343 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
1344 materialLawManager_->initFromState(eclState);
1345 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1346 this->
template fieldPropIntTypeOnLeafAssigner_<int>(),
1347 this-> lookupIdxOnLevelZeroAssigner_());
1351 void readThermalParameters_()
1353 if constexpr (enableEnergy)
1355 const auto& simulator = this->simulator();
1356 const auto& vanguard = simulator.vanguard();
1357 const auto& eclState = vanguard.eclState();
1360 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
1361 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1362 this-> fieldPropDoubleOnLeafAssigner_(),
1363 this->
template fieldPropIntTypeOnLeafAssigner_<unsigned int>());
1367 void updateReferencePorosity_()
1369 const auto& simulator = this->simulator();
1370 const auto& vanguard = simulator.vanguard();
1371 const auto& eclState = vanguard.eclState();
1373 std::size_t numDof = this->model().numGridDof();
1375 this->referencePorosity_[0].resize(numDof);
1377 const auto& fp = eclState.fieldProps();
1378 const std::vector<double> porvData =
this -> fieldPropDoubleOnLeafAssigner_()(fp,
"PORV");
1379 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1380 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(dofIdx);
1381 Scalar poreVolume = porvData[dofIdx];
1386 Scalar dofVolume = simulator.model().dofTotalVolume(sfcdofIdx);
1387 assert(dofVolume > 0.0);
1388 this->referencePorosity_[0][sfcdofIdx] = poreVolume/dofVolume;
1392 virtual void readInitialCondition_()
1395 const auto& simulator = this->simulator();
1396 const auto& vanguard = simulator.vanguard();
1397 const auto& eclState = vanguard.eclState();
1399 if (eclState.getInitConfig().hasEquil())
1400 readEquilInitialCondition_();
1402 readExplicitInitialCondition_();
1405 std::size_t numElems = this->model().numGridDof();
1406 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1407 const auto& fs = asImp_().initialFluidStates()[elemIdx];
1408 if (!this->maxWaterSaturation_.empty() && waterPhaseIdx > -1)
1409 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
1410 if (!this->maxOilSaturation_.empty() && oilPhaseIdx > -1)
1411 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
1412 if (!this->minRefPressure_.empty() && refPressurePhaseIdx_() > -1)
1413 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
1417 virtual void readEquilInitialCondition_() = 0;
1418 virtual void readExplicitInitialCondition_() = 0;
1421 bool updateHysteresis_()
1423 if (!materialLawManager_->enableHysteresis())
1428 this->updateProperty_(
"FlowProblem::updateHysteresis_() failed:",
1429 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1431 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
1437 bool updateHysteresis_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
1439 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
1440 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
1445 Scalar getRockCompTransMultVal(std::size_t dofIdx)
const
1447 if (this->rockCompTransMultVal_.empty())
1450 return this->rockCompTransMultVal_[dofIdx];
1456 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
1457 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
1458 ConditionalStorage<enableDiffusion, Scalar> diffusivity;
1459 ConditionalStorage<enableDispersion, Scalar> dispersivity;
1460 Scalar transmissibility;
1464 void updatePffDofData_()
1466 const auto& distFn =
1468 const Stencil& stencil,
1469 unsigned localDofIdx)
1472 const auto& elementMapper = this->model().elementMapper();
1474 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
1475 if (localDofIdx != 0) {
1476 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(0));
1477 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
1479 if constexpr (enableEnergy) {
1480 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
1481 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
1483 if constexpr (enableDiffusion)
1484 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
1485 if (enableDispersion)
1486 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
1490 pffDofData_.update(distFn);
1493 virtual void updateExplicitQuantities_(
int episodeIdx,
int timeStepSize,
bool first_step_after_restart) = 0;
1495 void readBoundaryConditions_()
1497 const auto& simulator = this->simulator();
1498 const auto& vanguard = simulator.vanguard();
1499 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
1500 if (bcconfig.size() > 0) {
1501 nonTrivialBoundaryConditions_ =
true;
1503 std::size_t numCartDof = vanguard.cartesianSize();
1504 unsigned numElems = vanguard.gridView().size(0);
1505 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
1507 for (
unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
1508 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
1510 bcindex_.resize(numElems, 0);
1511 auto loopAndApply = [&cartesianToCompressedElemIdx,
1512 &vanguard](
const auto& bcface,
1515 for (
int i = bcface.i1; i <= bcface.i2; ++i) {
1516 for (
int j = bcface.j1; j <= bcface.j2; ++j) {
1517 for (
int k = bcface.k1; k <= bcface.k2; ++k) {
1518 std::array<int, 3> tmp = {i,j,k};
1519 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
1526 for (
const auto& bcface : bcconfig) {
1527 std::vector<int>& data = bcindex_(bcface.dir);
1528 const int index = bcface.index;
1529 loopAndApply(bcface,
1530 [&data,index](
int elemIdx)
1531 { data[elemIdx] = index; });
1538 Scalar limitNextTimeStepSize_(Scalar dtNext)
const
1540 if constexpr (enableExperiments) {
1541 const auto& simulator = this->simulator();
1542 const auto& schedule = simulator.vanguard().schedule();
1543 int episodeIdx = simulator.episodeIndex();
1546 Scalar maxTimeStepSize = Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60;
1547 int reportStepIdx = std::max(episodeIdx, 0);
1548 if (this->enableTuning_) {
1549 const auto& tuning = schedule[reportStepIdx].tuning();
1550 maxTimeStepSize = tuning.TSMAXZ;
1553 dtNext = std::min(dtNext, maxTimeStepSize);
1555 Scalar remainingEpisodeTime =
1556 simulator.episodeStartTime() + simulator.episodeLength()
1557 - (simulator.startTime() + simulator.time());
1558 assert(remainingEpisodeTime >= 0.0);
1562 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
1565 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
1567 if (simulator.episodeStarts()) {
1570 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
1571 bool wellEventOccured =
1572 events.hasEvent(ScheduleEvents::NEW_WELL)
1573 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
1574 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
1575 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
1576 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
1577 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
1584 int refPressurePhaseIdx_()
const {
1585 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1588 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1592 return waterPhaseIdx;
1596 void updateRockCompTransMultVal_()
1598 const auto& model = this->simulator().model();
1599 std::size_t numGridDof = this->model().numGridDof();
1600 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
1601 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
1602 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, 0);
1604 this->rockCompTransMultVal_[elementIdx] = trans_mult;
1613 template <
class LhsEval>
1616 OPM_TIMEBLOCK_LOCAL(computeRockCompTransMultiplier);
1617 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
1620 unsigned tableIdx = 0;
1621 if (!this->rockTableIdx_.empty())
1622 tableIdx = this->rockTableIdx_[elementIdx];
1624 const auto& fs = intQuants.fluidState();
1625 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1627 if (!this->minRefPressure_.empty())
1630 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1631 this->minRefPressure_[elementIdx]);
1633 if (!this->overburdenPressure_.empty())
1634 effectivePressure -= this->overburdenPressure_[elementIdx];
1636 if (!this->rockCompTransMult_.empty())
1637 return this->rockCompTransMult_[tableIdx].eval(effectivePressure,
true);
1640 assert(!this->rockCompTransMultWc_.empty());
1641 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1642 LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
1644 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax,
true);
1647 typename Vanguard::TransmissibilityType transmissibilities_;
1649 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
1650 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
1652 bool enableDriftCompensation_;
1653 GlobalEqVector drift_;
1655 WellModel wellModel_;
1656 AquiferModel aquiferModel_;
1658 bool enableVtkOutput_;
1667 std::array<std::vector<T>,6> data;
1669 void resize(std::size_t size, T defVal)
1671 for (
auto& d : data)
1672 d.resize(size, defVal);
1675 const std::vector<T>& operator()(FaceDir::DirEnum dir)
const
1677 if (dir == FaceDir::DirEnum::Unknown)
1678 throw std::runtime_error(
"Tried to access BC data for the 'Unknown' direction");
1680 int div =
static_cast<int>(dir);
1681 while ((div /= 2) >= 1)
1683 assert(idx >= 0 && idx <= 5);
1687 std::vector<T>& operator()(FaceDir::DirEnum dir)
1689 return const_cast<std::vector<T>&
>(std::as_const(*
this)(dir));
1693 virtual void handleSolventBC(
const BCProp::BCFace&, RateVector&)
const = 0;
1695 virtual void handlePolymerBC(
const BCProp::BCFace&, RateVector&)
const = 0;
1698 bool nonTrivialBoundaryConditions_ =
false;
1699 bool explicitRockCompaction_ =
false;
1700 bool first_step_ =
true;