78 using typename FlowProblemType::Scalar;
79 using typename FlowProblemType::Simulator;
80 using typename FlowProblemType::GridView;
81 using typename FlowProblemType::FluidSystem;
82 using typename FlowProblemType::Vanguard;
83 using typename FlowProblemType::GlobalEqVector;
84 using typename FlowProblemType::EqVector;
85 using FlowProblemType::dim;
86 using FlowProblemType::dimWorld;
87 using FlowProblemType::numEq;
88 using FlowProblemType::numPhases;
89 using FlowProblemType::numComponents;
92 using FlowProblemType::enableConvectiveMixing;
93 using FlowProblemType::enableBrine;
94 using FlowProblemType::enableDiffusion;
95 using FlowProblemType::enableDispersion;
96 using FlowProblemType::enableEnergy;
97 using FlowProblemType::enableExperiments;
98 using FlowProblemType::enableExtbo;
99 using FlowProblemType::enableFoam;
100 using FlowProblemType::enableMICP;
101 using FlowProblemType::enablePolymer;
102 using FlowProblemType::enablePolymerMolarWeight;
103 using FlowProblemType::enableSaltPrecipitation;
104 using FlowProblemType::enableSolvent;
105 using FlowProblemType::enableTemperature;
106 using FlowProblemType::enableThermalFluxBoundaries;
108 using FlowProblemType::gasPhaseIdx;
109 using FlowProblemType::oilPhaseIdx;
110 using FlowProblemType::waterPhaseIdx;
112 using FlowProblemType::waterCompIdx;
113 using FlowProblemType::oilCompIdx;
114 using FlowProblemType::gasCompIdx;
117 using typename FlowProblemType::RateVector;
118 using typename FlowProblemType::PrimaryVariables;
119 using typename FlowProblemType::Indices;
120 using typename FlowProblemType::IntensiveQuantities;
121 using typename FlowProblemType::ElementContext;
123 using typename FlowProblemType::MaterialLaw;
124 using typename FlowProblemType::DimMatrix;
137 using InitialFluidState =
typename EquilInitializer<TypeTag>::ScalarFluidState;
155 EclWriterType::registerParameters();
157 DamarisWriterType::registerParameters();
167 , thresholdPressures_(simulator)
168 , mixControls_(simulator.vanguard().schedule())
169 , actionHandler_(simulator.vanguard().eclState(),
170 simulator.vanguard().schedule(),
171 simulator.vanguard().actionState(),
172 simulator.vanguard().summaryState(),
174 simulator.vanguard().grid().comm())
179 const auto& vanguard = simulator.vanguard();
183 enableSaltPrecipitation>(vanguard.eclState());
186 DiffusionModule::initFromState(vanguard.eclState());
187 DispersionModule::initFromState(vanguard.eclState());
210 eclWriter_ = std::make_unique<EclWriterType>(simulator);
211 enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
215 enableDamarisOutput_ = Parameters::Get<Parameters::EnableDamarisOutput>();
226 auto& simulator = this->simulator();
228 const auto& schedule = simulator.vanguard().schedule();
232 actionHandler_.evalUDQAssignments(
episodeIdx, simulator.vanguard().udqState());
236 if (
oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
239 FluidSystem::setVapPars(0.0, 0.0);
243 ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), simulator.vanguard().schedule(),
episodeIdx, moduleParams_.convectiveMixingModuleParam);
254 FlowProblemType::finishInit();
255 auto& simulator = this->simulator();
260 this->transmissibilities_.finishInit([&
vg = this->simulator().vanguard()](
const unsigned int it) {
261 return vg.gridIdxToEquilGridIdx(it);
270 if (enableEclOutput_) {
271 if (simulator.vanguard().grid().comm().size() > 1) {
272 if (simulator.vanguard().grid().comm().rank() == 0)
273 eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
276 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
279 std::function<
unsigned int(
unsigned int)>
equilGridToGrid = [&simulator](
unsigned int i) {
280 return simulator.vanguard().gridEquilIdxToGridIdx(i);
284 simulator.vanguard().releaseGlobalTransmissibilities();
286 const auto& eclState = simulator.vanguard().eclState();
287 const auto& schedule = simulator.vanguard().schedule();
290 simulator.setStartTime(schedule.getStartTime());
291 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
297 simulator.setEpisodeIndex(-1);
298 simulator.setEpisodeLength(0.0);
303 this->gravity_ = 0.0;
304 if (Parameters::Get<Parameters::EnableGravity>())
305 this->gravity_[dim - 1] = 9.80665;
306 if (!eclState.getInitConfig().hasGravity())
307 this->gravity_[dim - 1] = 0.0;
309 if (this->enableTuning_) {
312 const auto&
tuning = schedule[0].tuning();
313 this->initialTimeStepSize_ =
tuning.TSINIT.has_value() ?
tuning.TSINIT.value() : -1.0;
314 this->maxTimeStepAfterWellEvent_ =
tuning.TMAXWC;
317 this->initFluidSystem_();
321 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
322 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
323 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
326 this->readRockParameters_(simulator.vanguard().cellCenterDepths(),
327 [&simulator](
const unsigned idx)
329 std::array<int,dim> coords;
330 simulator.vanguard().cartesianCoordinate(idx, coords);
331 for (auto& c : coords) {
336 this->readMaterialParameters_();
337 this->readThermalParameters_();
340 if (enableEclOutput_) {
341 eclWriter_->writeInit();
346 const auto&
initconfig = eclState.getInitConfig();
347 this->tracerModel_.init(
initconfig.restartRequested());
349 readEclRestartSolution_();
351 readInitialCondition_();
353 this->tracerModel_.prepareTracerBatches();
355 this->updatePffDofData_();
358 const auto& vanguard = this->simulator().vanguard();
359 const auto& gridView = vanguard.gridView();
361 this->polymer_.maxAdsorption.resize(
numElements, 0.0);
364 this->readBoundaryConditions_();
367 computeAndSetEqWeights_();
369 if (this->enableDriftCompensation_) {
370 this->drift_.resize(this->model().numGridDof());
374 if (this->enableVtkOutput_ && eclState.getIOConfig().initOnly()) {
375 simulator.setTimeStepSize(0.0);
383 simulator.startNextEpisode(schedule.seconds(1));
384 simulator.setEpisodeIndex(0);
385 simulator.setTimeStepIndex(0);
390 this->mixControls_.init(this->model().numGridDof(),
391 this->episodeIndex(),
392 eclState.runspec().tabdims().getNumPVTTables());
400 FlowProblemType::endTimeStep();
403 this->eclWriter()->mutableOutputModule().invalidateLocalData();
405 const bool isSubStep = !this->simulator().episodeWillBeOver();
408 const auto& grid = this->simulator().vanguard().gridView().grid();
410 using GridType = std::remove_cv_t<std::remove_reference_t<
decltype(grid)>>;
411 constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
412 if (!
isCpGrid || (grid.maxLevel() == 0)) {
413 this->eclWriter_->evalSummaryState(
isSubStep);
420 auto& simulator = this->simulator();
424 .applyActions(
episodeIdx, simulator.time() + simulator.timeStepSize(),
425 [
this](
const bool global)
427 using TransUpdateQuantities = typename Vanguard::TransmissibilityType::TransUpdateQuantities;
428 this->transmissibilities_
429 .update(global, TransUpdateQuantities::All, [&vg = this->simulator().vanguard()]
430 (const unsigned int i)
432 return vg.gridIdxToEquilGridIdx(i);
438 if constexpr (enableMICP) {
439 auto& model = this->model();
440 const auto& residual = model.linearizer().residual();
444 MICPModule::checkCloggingMICP(model, phi,
globalDofIdx);
454 OPM_TIMEBLOCK(endEpisode);
455 const int episodeIdx = this->episodeIndex();
467 .evalUDQAssignments(episodeIdx, this->simulator().vanguard().udqState());
469 FlowProblemType::endEpisode();
473 if (enableEclOutput_){
474 eclWriter_->writeReports(timer);
485 FlowProblemType::writeOutput(verbose);
487 bool isSubStep = !this->simulator().episodeWillBeOver();
489 data::Solution localCellData = {};
493 if (enableDamarisOutput_) {
494 damarisWriter_->writeOutput(localCellData, isSubStep) ;
497 if (enableEclOutput_){
498 eclWriter_->writeOutput(std::move(localCellData), isSubStep);
502 void finalizeOutput()
504 OPM_TIMEBLOCK(finalizeOutput);
516 FlowProblemType::initialSolutionApplied();
521 this->thresholdPressures_.finishInit();
523 if (this->simulator().episodeIndex() == 0) {
524 eclWriter_->writeInitialFIPReport();
528 void addToSourceDense(RateVector& rate,
529 unsigned globalDofIdx,
530 unsigned timeIdx)
const override
532 this->aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
535 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
536 std::array<int,3> ijk;
537 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
539 if (source.hasSource(ijk)) {
540 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
541 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
542 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
543 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
545 for (
unsigned i = 0; i < phidx_map.size(); ++i) {
546 const auto phaseIdx = phidx_map[i];
547 const auto sourceComp = sc_map[i];
548 const auto compIdx = cidx_map[i];
549 if (!FluidSystem::phaseIsActive(phaseIdx)) {
552 Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
553 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
554 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
556 rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate;
559 if constexpr (enableSolvent) {
560 Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
561 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
562 const auto& solventPvt = SolventModule::solventPvt();
563 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
565 rate[Indices::contiSolventEqIdx] += mass_rate;
567 if constexpr (enablePolymer) {
568 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
570 if constexpr (enableEnergy) {
571 for (
unsigned i = 0; i < phidx_map.size(); ++i) {
572 const auto phaseIdx = phidx_map[i];
573 if (!FluidSystem::phaseIsActive(phaseIdx)) {
576 const auto sourceComp = sc_map[i];
577 if (source.hasHrate({ijk, sourceComp})) {
578 rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
580 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, 0);
581 auto fs = intQuants.fluidState();
583 if (source.hasTemperature({ijk, sourceComp})) {
584 Scalar temperature = source.temperature({ijk, sourceComp});
585 fs.setTemperature(temperature);
587 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
588 Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx);
589 Scalar energy_rate = getValue(h)*mass_rate;
590 rate[Indices::contiEnergyEqIdx] += energy_rate;
598 if (this->enableDriftCompensation_) {
599 const auto& simulator = this->simulator();
600 const auto& model = this->model();
605 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
606 Scalar poro = this->porosity(globalDofIdx, timeIdx);
607 Scalar dt = simulator.timeStepSize();
608 EqVector dofDriftRate = this->drift_[globalDofIdx];
609 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
612 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
613 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
614 if (cnv > maxCompensation) {
615 dofDriftRate[eqIdx] *= maxCompensation/cnv;
619 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
620 rate[eqIdx] -= dofDriftRate[eqIdx];
629 template <
class LhsEval>
632 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
633 if (!enableSaltPrecipitation)
636 const auto& fs = intQuants.fluidState();
637 unsigned tableIdx = fs.pvtRegionIndex();
638 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
639 porosityFactor = min(porosityFactor, 1.0);
640 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
641 return permfactTable.eval(porosityFactor,
true);
645 const InitialFluidState& initialFluidState(
unsigned globalDofIdx)
const
646 {
return initialFluidStates_[globalDofIdx]; }
648 std::vector<InitialFluidState>& initialFluidStates()
649 {
return initialFluidStates_; }
651 const std::vector<InitialFluidState>& initialFluidStates()
const
652 {
return initialFluidStates_; }
654 const EclipseIO& eclIO()
const
655 {
return eclWriter_->eclIO(); }
657 void setSubStepReport(
const SimulatorReportSingle& report)
658 {
return eclWriter_->setSubStepReport(report); }
660 void setSimulationReport(
const SimulatorReport& report)
661 {
return eclWriter_->setSimulationReport(report); }
663 InitialFluidState boundaryFluidState(
unsigned globalDofIdx,
const int directionId)
const
665 OPM_TIMEBLOCK_LOCAL(boundaryFluidState);
666 const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop;
667 if (bcprop.size() > 0) {
668 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
672 if (this->bcindex_(dir)[globalDofIdx] == 0)
673 return initialFluidStates_[globalDofIdx];
675 const auto& bc = bcprop[this->bcindex_(dir)[globalDofIdx]];
676 if (bc.bctype == BCType::DIRICHLET )
678 InitialFluidState fluidState;
679 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
680 fluidState.setPvtRegionIndex(pvtRegionIdx);
682 switch (bc.component) {
683 case BCComponent::OIL:
684 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
685 throw std::logic_error(
"oil is not active and you're trying to add oil BC");
687 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
689 case BCComponent::GAS:
690 if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
691 throw std::logic_error(
"gas is not active and you're trying to add gas BC");
693 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
695 case BCComponent::WATER:
696 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
697 throw std::logic_error(
"water is not active and you're trying to add water BC");
699 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
701 case BCComponent::SOLVENT:
702 case BCComponent::POLYMER:
703 case BCComponent::NONE:
704 throw std::logic_error(
"you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
706 fluidState.setTotalSaturation(1.0);
707 double pressure = initialFluidStates_[globalDofIdx].pressure(this->refPressurePhaseIdx_());
708 const auto pressure_input = bc.pressure;
709 if (pressure_input) {
710 pressure = *pressure_input;
713 std::array<Scalar, numPhases> pc = {0};
714 const auto& matParams = this->materialLawParams(globalDofIdx);
715 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
716 Valgrind::CheckDefined(pressure);
717 Valgrind::CheckDefined(pc);
718 for (
unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
719 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
721 fluidState.setPc(phaseIdx, pc[phaseIdx]);
722 if (Indices::oilEnabled)
723 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
724 else if (Indices::gasEnabled)
725 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
726 else if (Indices::waterEnabled)
728 fluidState.setPressure(phaseIdx, pressure);
731 double temperature = initialFluidStates_[globalDofIdx].temperature(0);
732 const auto temperature_input = bc.temperature;
733 if(temperature_input)
734 temperature = *temperature_input;
735 fluidState.setTemperature(temperature);
737 if (FluidSystem::enableDissolvedGas()) {
738 fluidState.setRs(0.0);
739 fluidState.setRv(0.0);
741 if (FluidSystem::enableDissolvedGasInWater()) {
742 fluidState.setRsw(0.0);
744 if (FluidSystem::enableVaporizedWater())
745 fluidState.setRvw(0.0);
747 for (
unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
748 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
750 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
751 fluidState.setInvB(phaseIdx, b);
753 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
754 fluidState.setDensity(phaseIdx, rho);
756 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
757 fluidState.setEnthalpy(phaseIdx, h);
760 fluidState.checkDefined();
764 return initialFluidStates_[globalDofIdx];
768 const std::unique_ptr<EclWriterType>& eclWriter()
const
773 void setConvData(
const std::vector<std::vector<int>>& data)
775 eclWriter_->mutableOutputModule().setCnvData(data);
784 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
785 this->episodeIndex(),
786 this->pvtRegionIndex(globalDofIdx));
795 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
796 this->episodeIndex(),
797 this->pvtRegionIndex(globalDofIdx));
810 int episodeIdx = this->episodeIndex();
811 return !this->mixControls_.drsdtActive(episodeIdx) &&
812 !this->mixControls_.drvdtActive(episodeIdx) &&
813 this->rockCompPoroMultWc_.empty() &&
814 this->rockCompPoroMult_.empty();
823 template <
class Context>
824 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
826 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
828 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
829 values.assignNaive(initialFluidStates_[globalDofIdx]);
831 SolventModule::assignPrimaryVars(values,
832 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
833 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
835 if constexpr (enablePolymer)
836 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
838 if constexpr (enablePolymerMolarWeight)
839 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
841 if constexpr (enableBrine) {
842 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
843 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
846 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
850 if constexpr (enableMICP){
851 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
852 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
853 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
854 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
855 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
858 values.checkDefined();
862 Scalar drsdtcon(
unsigned elemIdx,
int episodeIdx)
const
864 return this->mixControls_.drsdtcon(elemIdx, episodeIdx,
865 this->pvtRegionIndex(elemIdx));
873 template <
class Context>
875 const Context& context,
877 unsigned timeIdx)
const
879 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
880 if (!context.intersection(spaceIdx).boundary())
883 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
891 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
892 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
893 values.setThermalFlow(context, spaceIdx, timeIdx, this->initialFluidStates_[globalDofIdx] );
896 if (this->nonTrivialBoundaryConditions()) {
897 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
898 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
899 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
900 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
901 const auto [type, massrate] = this->boundaryCondition(globalDofIdx, indexInInside);
902 if (type == BCType::THERMAL)
903 values.setThermalFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
904 else if (type == BCType::FREE || type == BCType::DIRICHLET)
905 values.setFreeFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
906 else if (type == BCType::RATE)
907 values.setMassRate(massrate, pvtRegionIdx);
916 {
return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
919 {
return thresholdPressures_; }
921 FlowThresholdPressure<TypeTag>& thresholdPressure()
922 {
return thresholdPressures_; }
924 const ModuleParams& moduleParams()
const
926 return moduleParams_;
929 template<
class Serializer>
930 void serializeOp(Serializer& serializer)
932 serializer(
static_cast<FlowProblemType&
>(*
this));
933 serializer(mixControls_);
934 serializer(*eclWriter_);
938 void updateExplicitQuantities_(
int episodeIdx,
int timeStepSize,
const bool first_step_after_restart)
override
940 this->updateExplicitQuantities_(first_step_after_restart);
942 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
943 updateMaxPolymerAdsorption_();
945 mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize);
948 void updateMaxPolymerAdsorption_()
951 this->updateProperty_(
"FlowProblemBlackoil::updateMaxPolymerAdsorption_() failed:",
952 [
this](
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
954 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
958 bool updateMaxPolymerAdsorption_(
unsigned compressedDofIdx,
const IntensiveQuantities& iq)
960 const Scalar pa = scalarValue(iq.polymerAdsorption());
961 auto& mpa = this->polymer_.maxAdsorption;
962 if (mpa[compressedDofIdx] < pa) {
963 mpa[compressedDofIdx] = pa;
970 void computeAndSetEqWeights_()
972 std::vector<Scalar> sumInvB(numPhases, 0.0);
973 const auto& gridView = this->gridView();
974 ElementContext elemCtx(this->simulator());
975 for(
const auto& elem: elements(gridView, Dune::Partitions::interior)) {
976 elemCtx.updatePrimaryStencil(elem);
977 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
978 const auto& dofFluidState = this->initialFluidStates_[elemIdx];
979 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
980 if (!FluidSystem::phaseIsActive(phaseIdx))
983 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
987 std::size_t numDof = this->model().numGridDof();
988 const auto& comm = this->simulator().vanguard().grid().comm();
989 comm.sum(sumInvB.data(),sumInvB.size());
990 Scalar numTotalDof = comm.sum(numDof);
992 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
993 if (!FluidSystem::phaseIsActive(phaseIdx))
996 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
997 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
998 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
999 this->model().setEqWeight(activeSolventCompIdx, avgB);
1004 bool updateCompositionChangeLimits_()
1006 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1009 int episodeIdx = this->episodeIndex();
1010 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1011 this->mixControls_.drsdtActive(episodeIdx),
1012 this->mixControls_.drvdtActive(episodeIdx)};
1013 if (!active[0] && !active[1] && !active[2]) {
1017 this->updateProperty_(
"FlowProblemBlackoil::updateCompositionChangeLimits_()) failed:",
1018 [
this,episodeIdx,active](
unsigned compressedDofIdx,
1019 const IntensiveQuantities& iq)
1021 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1022 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1023 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1024 this->mixControls_.update(compressedDofIdx,
1027 this->gravity_[dim - 1],
1028 perm[dim - 1][dim - 1],
1038 void readEclRestartSolution_()
1041 if(this->simulator().vanguard().grid().maxLevel() > 0) {
1042 throw std::invalid_argument(
"Refined grids are not yet supported for restart ");
1046 auto& simulator = this->simulator();
1047 const auto& schedule = simulator.vanguard().schedule();
1048 const auto& eclState = simulator.vanguard().eclState();
1049 const auto& initconfig = eclState.getInitConfig();
1050 const int restart_step = initconfig.getRestartStep();
1052 simulator.setTime(schedule.seconds(restart_step));
1054 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
1055 schedule.stepLength(restart_step));
1056 simulator.setEpisodeIndex(restart_step);
1058 this->eclWriter_->beginRestart();
1060 Scalar dt = std::min(this->eclWriter_->restartTimeStepSize(), simulator.episodeLength());
1061 simulator.setTimeStepSize(dt);
1063 std::size_t numElems = this->model().numGridDof();
1064 this->initialFluidStates_.resize(numElems);
1065 if constexpr (enableSolvent) {
1066 this->solventSaturation_.resize(numElems, 0.0);
1067 this->solventRsw_.resize(numElems, 0.0);
1070 if constexpr (enablePolymer)
1071 this->polymer_.concentration.resize(numElems, 0.0);
1073 if constexpr (enablePolymerMolarWeight) {
1074 const std::string msg {
"Support of the RESTART for polymer molecular weight "
1075 "is not implemented yet. The polymer weight value will be "
1076 "zero when RESTART begins"};
1077 OpmLog::warning(
"NO_POLYMW_RESTART", msg);
1078 this->polymer_.moleWeight.resize(numElems, 0.0);
1081 if constexpr (enableMICP) {
1082 this->micp_.resize(numElems);
1086 this->mixControls_.init(numElems, restart_step, eclState.runspec().tabdims().getNumPVTTables());
1088 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1089 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1090 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
1091 this->eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
1092 this->eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
1101 auto ssol = enableSolvent
1102 ? this->eclWriter_->outputModule().getSolventSaturation(elemIdx)
1105 this->processRestartSaturations_(elemFluidState, ssol);
1107 if constexpr (enableSolvent) {
1108 this->solventSaturation_[elemIdx] = ssol;
1109 this->solventRsw_[elemIdx] = this->eclWriter_->outputModule().getSolventRsw(elemIdx);
1114 bool isThermal = eclState.getSimulationConfig().isThermal();
1115 bool needTemperature = (eclState.runspec().co2Storage() || eclState.runspec().h2Storage());
1116 if (!isThermal && needTemperature) {
1117 const auto& fp = simulator.vanguard().eclState().fieldProps();
1118 elemFluidState.setTemperature(fp.get_double(
"TEMPI")[elemIdx]);
1121 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
1123 if constexpr (enablePolymer)
1124 this->polymer_.concentration[elemIdx] = this->eclWriter_->outputModule().getPolymerConcentration(elemIdx);
1125 if constexpr (enableMICP){
1126 this->micp_.microbialConcentration[elemIdx] = this->eclWriter_->outputModule().getMicrobialConcentration(elemIdx);
1127 this->micp_.oxygenConcentration[elemIdx] = this->eclWriter_->outputModule().getOxygenConcentration(elemIdx);
1128 this->micp_.ureaConcentration[elemIdx] = this->eclWriter_->outputModule().getUreaConcentration(elemIdx);
1129 this->micp_.biofilmConcentration[elemIdx] = this->eclWriter_->outputModule().getBiofilmConcentration(elemIdx);
1130 this->micp_.calciteConcentration[elemIdx] = this->eclWriter_->outputModule().getCalciteConcentration(elemIdx);
1135 const int episodeIdx = this->episodeIndex();
1136 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
1141 auto& sol = this->model().solution(0);
1142 const auto& gridView = this->gridView();
1143 ElementContext elemCtx(simulator);
1144 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1145 elemCtx.updatePrimaryStencil(elem);
1146 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
1147 this->initial(sol[elemIdx], elemCtx, 0, 0);
1155 this->model().syncOverlap();
1157 this->eclWriter_->endRestart();
1160 void readEquilInitialCondition_()
override
1162 const auto& simulator = this->simulator();
1165 EquilInitializer<TypeTag> equilInitializer(simulator, *(this->materialLawManager_));
1167 std::size_t numElems = this->model().numGridDof();
1168 this->initialFluidStates_.resize(numElems);
1169 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1170 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1171 elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
1175 void readExplicitInitialCondition_()
override
1177 const auto& simulator = this->simulator();
1178 const auto& vanguard = simulator.vanguard();
1179 const auto& eclState = vanguard.eclState();
1180 const auto& fp = eclState.fieldProps();
1181 bool has_swat = fp.has_double(
"SWAT");
1182 bool has_sgas = fp.has_double(
"SGAS");
1183 bool has_rs = fp.has_double(
"RS");
1184 bool has_rv = fp.has_double(
"RV");
1185 bool has_rvw = fp.has_double(
"RVW");
1186 bool has_pressure = fp.has_double(
"PRESSURE");
1187 bool has_salt = fp.has_double(
"SALT");
1188 bool has_saltp = fp.has_double(
"SALTP");
1191 if (Indices::numPhases > 1) {
1192 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
1193 throw std::runtime_error(
"The ECL input file requires the presence of the SWAT keyword if "
1194 "the water phase is active");
1195 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
1196 throw std::runtime_error(
"The ECL input file requires the presence of the SGAS keyword if "
1197 "the gas phase is active");
1200 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
1201 "keyword if the model is initialized explicitly");
1202 if (FluidSystem::enableDissolvedGas() && !has_rs)
1203 throw std::runtime_error(
"The ECL input file requires the RS keyword to be present if"
1204 " dissolved gas is enabled");
1205 if (FluidSystem::enableVaporizedOil() && !has_rv)
1206 throw std::runtime_error(
"The ECL input file requires the RV keyword to be present if"
1207 " vaporized oil is enabled");
1208 if (FluidSystem::enableVaporizedWater() && !has_rvw)
1209 throw std::runtime_error(
"The ECL input file requires the RVW keyword to be present if"
1210 " vaporized water is enabled");
1211 if (enableBrine && !has_salt)
1212 throw std::runtime_error(
"The ECL input file requires the SALT keyword to be present if"
1213 " brine is enabled and the model is initialized explicitly");
1214 if (enableSaltPrecipitation && !has_saltp)
1215 throw std::runtime_error(
"The ECL input file requires the SALTP keyword to be present if"
1216 " salt precipitation is enabled and the model is initialized explicitly");
1218 std::size_t numDof = this->model().numGridDof();
1220 initialFluidStates_.resize(numDof);
1222 std::vector<double> waterSaturationData;
1223 std::vector<double> gasSaturationData;
1224 std::vector<double> pressureData;
1225 std::vector<double> rsData;
1226 std::vector<double> rvData;
1227 std::vector<double> rvwData;
1228 std::vector<double> tempiData;
1229 std::vector<double> saltData;
1230 std::vector<double> saltpData;
1232 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
1233 waterSaturationData = fp.get_double(
"SWAT");
1235 waterSaturationData.resize(numDof);
1237 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
1238 gasSaturationData = fp.get_double(
"SGAS");
1240 gasSaturationData.resize(numDof);
1242 pressureData = fp.get_double(
"PRESSURE");
1243 if (FluidSystem::enableDissolvedGas())
1244 rsData = fp.get_double(
"RS");
1246 if (FluidSystem::enableVaporizedOil())
1247 rvData = fp.get_double(
"RV");
1249 if (FluidSystem::enableVaporizedWater())
1250 rvwData = fp.get_double(
"RVW");
1253 tempiData = fp.get_double(
"TEMPI");
1256 if constexpr (enableBrine)
1257 saltData = fp.get_double(
"SALT");
1260 if constexpr (enableSaltPrecipitation)
1261 saltpData = fp.get_double(
"SALTP");
1264 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1265 auto& dofFluidState = initialFluidStates_[dofIdx];
1267 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
1272 Scalar temperatureLoc = tempiData[dofIdx];
1273 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
1274 temperatureLoc = FluidSystem::surfaceTemperature;
1275 dofFluidState.setTemperature(temperatureLoc);
1280 if constexpr (enableBrine)
1281 dofFluidState.setSaltConcentration(saltData[dofIdx]);
1286 if constexpr (enableSaltPrecipitation)
1287 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
1292 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1293 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
1294 waterSaturationData[dofIdx]);
1296 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
1297 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
1298 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1300 - waterSaturationData[dofIdx]);
1303 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1304 gasSaturationData[dofIdx]);
1306 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
1307 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
1309 - waterSaturationData[dofIdx]
1310 - gasSaturationData[dofIdx]);
1315 Scalar pressure = pressureData[dofIdx];
1319 std::array<Scalar, numPhases> pc = {0};
1320 const auto& matParams = this->materialLawParams(dofIdx);
1321 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
1322 Valgrind::CheckDefined(pressure);
1323 Valgrind::CheckDefined(pc);
1324 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1325 if (!FluidSystem::phaseIsActive(phaseIdx))
1328 if (Indices::oilEnabled)
1329 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1330 else if (Indices::gasEnabled)
1331 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1332 else if (Indices::waterEnabled)
1334 dofFluidState.setPressure(phaseIdx, pressure);
1337 if (FluidSystem::enableDissolvedGas())
1338 dofFluidState.setRs(rsData[dofIdx]);
1339 else if (Indices::gasEnabled && Indices::oilEnabled)
1340 dofFluidState.setRs(0.0);
1342 if (FluidSystem::enableVaporizedOil())
1343 dofFluidState.setRv(rvData[dofIdx]);
1344 else if (Indices::gasEnabled && Indices::oilEnabled)
1345 dofFluidState.setRv(0.0);
1347 if (FluidSystem::enableVaporizedWater())
1348 dofFluidState.setRvw(rvwData[dofIdx]);
1353 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1354 if (!FluidSystem::phaseIsActive(phaseIdx))
1357 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1358 dofFluidState.setInvB(phaseIdx, b);
1360 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1361 dofFluidState.setDensity(phaseIdx, rho);
1368 void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
1372 const Scalar smallSaturationTolerance = 1.e-6;
1373 Scalar sumSaturation = 0.0;
1374 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1375 if (FluidSystem::phaseIsActive(phaseIdx)) {
1376 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
1377 elemFluidState.setSaturation(phaseIdx, 0.0);
1379 sumSaturation += elemFluidState.saturation(phaseIdx);
1383 if constexpr (enableSolvent) {
1384 if (solventSaturation < smallSaturationTolerance)
1385 solventSaturation = 0.0;
1387 sumSaturation += solventSaturation;
1390 assert(sumSaturation > 0.0);
1392 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1393 if (FluidSystem::phaseIsActive(phaseIdx)) {
1394 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
1395 elemFluidState.setSaturation(phaseIdx, saturation);
1398 if constexpr (enableSolvent) {
1399 solventSaturation = solventSaturation / sumSaturation;
1403 void readInitialCondition_()
override
1405 FlowProblemType::readInitialCondition_();
1407 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
1408 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
1411 enablePolymerMolarWeight,
1416 void handleSolventBC(
const BCProp::BCFace& bc, RateVector& rate)
const override
1418 if constexpr (!enableSolvent)
1419 throw std::logic_error(
"solvent is disabled and you're trying to add solvent to BC");
1421 rate[Indices::solventSaturationIdx] = bc.rate;
1424 void handlePolymerBC(
const BCProp::BCFace& bc, RateVector& rate)
const override
1426 if constexpr (!enablePolymer)
1427 throw std::logic_error(
"polymer is disabled and you're trying to add polymer to BC");
1429 rate[Indices::polymerConcentrationIdx] = bc.rate;
1432 void updateExplicitQuantities_(
const bool first_step_after_restart)
1434 OPM_TIMEBLOCK(updateExplicitQuantities);
1435 const bool invalidateFromMaxWaterSat = this->updateMaxWaterSaturation_();
1436 const bool invalidateFromMinPressure = this->updateMinPressure_();
1439 const bool invalidateFromHyst = this->updateHysteresis_();
1440 const bool invalidateFromMaxOilSat = this->updateMaxOilSaturation_();
1443 const bool invalidateDRDT = !first_step_after_restart && this->updateCompositionChangeLimits_();
1446 const bool invalidateIntensiveQuantities
1447 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat || invalidateDRDT;
1448 if (invalidateIntensiveQuantities) {
1449 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1450 this->model().invalidateAndUpdateIntensiveQuantities(0);
1453 this->updateRockCompTransMultVal_();
1456 FlowThresholdPressure<TypeTag> thresholdPressures_;
1458 std::vector<InitialFluidState> initialFluidStates_;
1460 bool enableEclOutput_;
1461 std::unique_ptr<EclWriterType> eclWriter_;
1463 bool enableDamarisOutput_ = false ;
1464 std::unique_ptr<DamarisWriterType> damarisWriter_;
1466 MixingRateControls<FluidSystem> mixControls_;
1468 ActionHandler<Scalar> actionHandler_;
1470 ModuleParams moduleParams_;