30#ifndef OPM_FLOW_PROBLEM_HPP
31#define OPM_FLOW_PROBLEM_HPP
33#include <dune/common/version.hh>
34#include <dune/common/fvector.hh>
35#include <dune/common/fmatrix.hh>
37#include <opm/common/utility/TimeService.hpp>
39#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
40#include <opm/input/eclipse/Schedule/Schedule.hpp>
41#include <opm/input/eclipse/Units/Units.hpp>
43#include <opm/material/common/ConditionalStorage.hpp>
44#include <opm/material/common/Valgrind.hpp>
45#include <opm/material/densead/Evaluation.hpp>
46#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
47#include <opm/material/thermal/EclThermalLawManager.hpp>
53#include <opm/output/eclipse/EclipseIO.hpp>
62#include <opm/simulators/flow/FlowUtils.hpp>
65#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
66#include <opm/simulators/timestepping/SimulatorReport.hpp>
68#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
69#include <opm/simulators/utils/ParallelSerialization.hpp>
70#include <opm/simulators/utils/satfunc/RelpermDiagnostics.hpp>
72#include <opm/utility/CopyablePtr.hpp>
90template <
class TypeTag>
93 GetPropType<TypeTag, Properties::FluidSystem>>
110 enum { dim = GridView::dimension };
111 enum { dimWorld = GridView::dimensionworld };
115 enum { numPhases = FluidSystem::numPhases };
116 enum { numComponents = FluidSystem::numComponents };
134 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
135 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
136 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
140 enum { gasCompIdx = FluidSystem::gasCompIdx };
141 enum { oilCompIdx = FluidSystem::oilCompIdx };
142 enum { waterCompIdx = FluidSystem::waterCompIdx };
147 using Element =
typename GridView::template Codim<0>::Entity;
151 using MaterialLawParams =
typename EclMaterialLawManager::MaterialLawParams;
152 using SolidEnergyLawParams =
typename EclThermalLawManager::SolidEnergyLawParams;
153 using ThermalConductionLawParams =
typename EclThermalLawManager::ThermalConductionLawParams;
163 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
167 using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
182 ParentType::registerParameters();
191 const std::string&)>
addKey,
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 if (! Parameters::Get<Parameters::CheckSatfuncConsistency>()) {
233 simulator.vanguard().levelCartesianIndexMapper());
239 void prefetch(
const Element&
elem)
const
240 { this->pffDofData_.prefetch(
elem); }
253 template <
class Restarter>
260 wellModel_.deserialize(
res);
263 aquiferModel_.deserialize(
res);
272 template <
class Restarter>
275 wellModel_.serialize(
res);
277 aquiferModel_.serialize(
res);
280 int episodeIndex()
const
282 return std::max(this->simulator().episodeIndex(), 0);
292 auto& simulator = this->simulator();
294 auto& eclState = simulator.vanguard().eclState();
295 const auto& schedule = simulator.vanguard().schedule();
296 const auto& events = schedule[
episodeIdx].events();
298 if (
episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
306 const auto&
cc = simulator.vanguard().grid().comm();
307 eclState.apply_schedule_keywords(
miniDeck );
308 eclBroadcast(
cc, eclState.getTransMult() );
311 std::function<
unsigned int(
unsigned int)>
equilGridToGrid = [&simulator](
unsigned int i) {
312 return simulator.vanguard().gridEquilIdxToGridIdx(i);
316 using TransUpdateQuantities =
typename Vanguard::TransmissibilityType::TransUpdateQuantities;
317 transmissibilities_.update(
true, TransUpdateQuantities::All,
equilGridToGrid);
318 this->referencePorosity_[1] = this->referencePorosity_[0];
319 updateReferencePorosity_();
321 this->model().linearizer().updateDiscretizationParameters();
324 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
327 wellModel_.beginEpisode();
330 aquiferModel_.beginEpisode();
333 Scalar
dt = limitNextTimeStepSize_(simulator.episodeLength());
338 dt = std::min(
dt, this->initialTimeStepSize_);
339 simulator.setTimeStepSize(
dt);
349 const int timeStepSize = this->simulator().timeStepSize();
351 this->beginTimeStep_(enableExperiments,
353 this->simulator().timeStepIndex(),
354 this->simulator().startTime(),
355 this->simulator().time(),
357 this->simulator().endTime());
365 if (nonTrivialBoundaryConditions()) {
366 this->model().linearizer().updateBoundaryConditionData();
369 wellModel_.beginTimeStep();
370 aquiferModel_.beginTimeStep();
371 tracerModel_.beginTimeStep();
381 wellModel_.beginIteration();
382 aquiferModel_.beginIteration();
391 wellModel_.endIteration();
392 aquiferModel_.endIteration();
409 const int rank = this->simulator().gridView().comm().rank();
411 std::cout <<
"checking conservativeness of solution\n";
414 this->model().checkConservativeness(-1,
true);
416 std::cout <<
"solution is sufficiently conservative\n";
421 auto& simulator = this->simulator();
422 simulator.setTimeStepIndex(simulator.timeStepIndex()+1);
424 this->wellModel_.endTimeStep();
425 this->aquiferModel_.endTimeStep();
426 this->tracerModel_.endTimeStep();
429 this->model().linearizer().updateFlowsInfo();
431 if (this->enableDriftCompensation_) {
434 const auto& residual = this->model().linearizer().residual();
436 for (
unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
437 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(globalDofIdx);
454 this->wellModel_.endEpisode();
455 this->aquiferModel_.endEpisode();
457 const auto& schedule = this->simulator().vanguard().schedule();
460 if (
episodeIdx + 1 >=
static_cast<int>(schedule.size()) - 1) {
461 this->simulator().setFinished(
true);
466 this->simulator().startNextEpisode(schedule.stepLength(
episodeIdx + 1));
478 if (Parameters::Get<Parameters::EnableWriteAllSolutions>() ||
479 this->simulator().episodeWillBeOver()) {
480 ParentType::writeOutput(
verbose);
487 template <
class Context>
508 template <
class Context>
514 return pffDofData_.get(context.element(),
toDofLocalIdx).transmissibility;
528 template <
class Context>
534 return *pffDofData_.get(context.element(),
toDofLocalIdx).diffusivity;
566 template <
class Context>
570 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
571 return transmissibilities_.transmissibilityBoundary(elemIdx,
boundaryFaceIdx);
596 template <
class Context>
601 const auto&
face = context.stencil(
timeIdx).interiorFace(faceIdx);
603 return *pffDofData_.get(context.element(),
toDofLocalIdx).thermalHalfTransIn;
609 template <
class Context>
614 const auto&
face = context.stencil(
timeIdx).interiorFace(faceIdx);
616 return *pffDofData_.get(context.element(),
toDofLocalIdx).thermalHalfTransOut;
622 template <
class Context>
626 unsigned elemIdx = elemCtx.globalSpaceIndex(0, 0);
627 return transmissibilities_.thermalHalfTransBoundary(elemIdx,
boundaryFaceIdx);
634 {
return transmissibilities_; }
638 {
return tracerModel_; }
640 TracerModel& tracerModel()
641 {
return tracerModel_; }
651 template <
class Context>
664 template <
class Context>
679 return this->simulator().vanguard().cellCenterDepth(
globalSpaceIdx);
685 template <
class Context>
695 template <
class Context>
707 const auto&
rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
709 return asImp_().initialFluidState(
globalSpaceIdx).pressure(refPressurePhaseIdx_());
712 if (this->rockParams_.empty())
716 if (!this->rockTableIdx_.empty()) {
719 return this->rockParams_[
tableIdx].referencePressure;
726 template <
class Context>
736 return materialLawManager_->materialLawParams(globalDofIdx);
741 return materialLawManager_->materialLawParams(globalDofIdx,
facedir);
747 template <
class Context>
748 const SolidEnergyLawParams&
760 template <
class Context>
761 const ThermalConductionLawParams &
765 return thermalLawManager_->thermalConductionLawParams(
globalSpaceIdx);
775 {
return materialLawManager_; }
777 template <
class Flu
idState>
779 std::array<Evaluation,numPhases> &mobility,
780 DirectionalMobilityPtr &
dirMob,
781 FluidState &fluidState,
789 MaterialLaw::relativePermeabilities(mobility,
materialParams, fluidState);
790 Valgrind::CheckDefined(mobility);
792 if (materialLawManager_->hasDirectionalRelperms()
793 || materialLawManager_->hasDirectionalImbnum())
795 using Dir = FaceDir::DirEnum;
796 constexpr int ndim = 3;
797 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
798 Dir
facedirs[
ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
799 for (
int i = 0; i<
ndim; i++) {
811 {
return materialLawManager_; }
817 template <
class Context>
825 template <
class Context>
833 template <
class Context>
841 template <
class Context>
850 template <
class Context>
858 {
return this->simulator().vanguard().caseName(); }
863 template <
class Context>
869 return asImp_().initialFluidState(globalDofIdx).temperature(0);
873 Scalar
temperature(
unsigned globalDofIdx,
unsigned )
const
877 return asImp_().initialFluidState(globalDofIdx).temperature(0);
880 const SolidEnergyLawParams&
884 return this->thermalLawManager_->solidEnergyLawParams(
globalSpaceIdx);
886 const ThermalConductionLawParams &
890 return this->thermalLawManager_->thermalConductionLawParams(
globalSpaceIdx);
904 if (!this->vapparsActive(this->episodeIndex()))
907 return this->maxOilSaturation_[globalDofIdx];
921 if (!this->vapparsActive(this->episodeIndex()))
924 this->maxOilSaturation_[globalDofIdx] = value;
933 this->model().invalidateAndUpdateIntensiveQuantities(0);
939 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
940 this->model().invalidateAndUpdateIntensiveQuantities(1);
948 aquiferModel_.initialSolutionApplied();
953 this->model().invalidateAndUpdateIntensiveQuantities(0);
962 template <
class Context>
964 const Context& context,
968 const unsigned globalDofIdx = context.globalSpaceIndex(
spaceIdx,
timeIdx);
972 void source(RateVector& rate,
973 unsigned globalDofIdx,
980 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
985 rate[
eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
987 Valgrind::CheckDefined(rate[
eqIdx]);
992 addToSourceDense(rate, globalDofIdx,
timeIdx);
995 virtual void addToSourceDense(RateVector& rate,
996 unsigned globalDofIdx,
1005 {
return wellModel_; }
1008 {
return wellModel_; }
1010 const AquiferModel& aquiferModel()
const
1011 {
return aquiferModel_; }
1013 AquiferModel& mutableAquiferModel()
1014 {
return aquiferModel_; }
1016 bool nonTrivialBoundaryConditions()
const
1017 {
return nonTrivialBoundaryConditions_; }
1029 if (this->nextTimeStepSize_ > 0.0)
1030 return this->nextTimeStepSize_;
1032 const auto& simulator = this->simulator();
1037 return this->initialTimeStepSize_;
1042 const auto& newtonMethod = this->model().newtonMethod();
1043 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1051 template <
class LhsEval>
1055 if (this->rockCompPoroMult_.empty() &&
this->rockCompPoroMultWc_.empty())
1059 if (!this->rockTableIdx_.empty())
1062 const auto& fs = intQuants.fluidState();
1064 const auto&
rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
1065 if (!this->minRefPressure_.empty())
1071 if (!this->overburdenPressure_.empty())
1078 if (!this->rockCompPoroMult_.empty()) {
1083 assert(!this->rockCompPoroMultWc_.empty());
1095 template <
class LhsEval>
1098 const bool implicit = !this->explicitRockCompaction_;
1100 : this->simulator().problem().getRockCompTransMultVal(
elementIdx);
1107 template <
class LhsEval>
1112 const bool implicit = !this->explicitRockCompaction_;
1114 : this->simulator().problem().getRockCompTransMultVal(
elementIdx);
1123 if (!nonTrivialBoundaryConditions_) {
1124 return { BCType::NONE, RateVector(0.0) };
1126 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(
directionId);
1127 const auto& schedule = this->simulator().vanguard().schedule();
1129 return { BCType::NONE, RateVector(0.0) };
1131 if (schedule[this->episodeIndex()].
bcprop.size() == 0) {
1132 return { BCType::NONE, RateVector(0.0) };
1134 const auto&
bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[
globalSpaceIdx]];
1135 if (
bc.bctype!=BCType::RATE) {
1136 return {
bc.bctype, RateVector(0.0) };
1139 RateVector rate = 0.0;
1140 switch (
bc.component) {
1141 case BCComponent::OIL:
1142 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] =
bc.rate;
1144 case BCComponent::GAS:
1145 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] =
bc.rate;
1147 case BCComponent::WATER:
1148 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] =
bc.rate;
1150 case BCComponent::SOLVENT:
1151 this->handleSolventBC(
bc, rate);
1153 case BCComponent::POLYMER:
1154 this->handlePolymerBC(
bc, rate);
1156 case BCComponent::MICR:
1157 this->handleMicrBC(
bc, rate);
1159 case BCComponent::OXYG:
1160 this->handleOxygBC(
bc, rate);
1162 case BCComponent::UREA:
1163 this->handleUreaBC(
bc, rate);
1165 case BCComponent::NONE:
1166 throw std::logic_error(
"you need to specify the component when RATE type is set in BC");
1170 return {
bc.bctype, rate};
1174 template<
class Serializer>
1186 Implementation& asImp_()
1187 {
return *
static_cast<Implementation *
>(
this); }
1189 const Implementation& asImp_()
const
1190 {
return *
static_cast<const Implementation *
>(
this); }
1193 template<
class UpdateFunc>
1194 void updateProperty_(
const std::string&
failureMsg,
1198 const auto& model = this->simulator().model();
1199 const auto& primaryVars = model.solution(0);
1200 const auto& vanguard = this->simulator().vanguard();
1201 std::size_t numGridDof = primaryVars.size();
1202 OPM_BEGIN_PARALLEL_TRY_CATCH();
1204#pragma omp parallel for
1206 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1207 const auto&
iq = *model.cachedIntensiveQuantities(dofIdx, 0);
1210 OPM_END_PARALLEL_TRY_CATCH(
failureMsg, vanguard.grid().comm());
1213 bool updateMaxOilSaturation_()
1220 this->updateProperty_(
"FlowProblem::updateMaxOilSaturation_() failed:",
1231 bool updateMaxOilSaturation_(
unsigned compressedDofIdx,
const IntensiveQuantities&
iq)
1234 const auto& fs =
iq.fluidState();
1235 const Scalar So =
decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1236 auto&
mos = this->maxOilSaturation_;
1245 bool updateMaxWaterSaturation_()
1249 if (this->maxWaterSaturation_.empty())
1252 this->maxWaterSaturation_[1] = this->maxWaterSaturation_[0];
1253 this->updateProperty_(
"FlowProblem::updateMaxWaterSaturation_() failed:",
1262 bool updateMaxWaterSaturation_(
unsigned compressedDofIdx,
const IntensiveQuantities&
iq)
1265 const auto& fs =
iq.fluidState();
1266 const Scalar Sw =
decay<Scalar>(fs.saturation(waterPhaseIdx));
1267 auto&
mow = this->maxWaterSaturation_;
1276 bool updateMinPressure_()
1280 if (this->minRefPressure_.empty())
1283 this->updateProperty_(
"FlowProblem::updateMinPressure_() failed:",
1293 const auto& fs =
iq.fluidState();
1308 std::function<std::vector<double>(
const FieldPropsManager&,
const std::string&)>
1309 fieldPropDoubleOnLeafAssigner_()
1311 const auto&
lookup = this->lookUpData_;
1322 template<
typename IntType>
1323 std::function<std::vector<IntType>(
const FieldPropsManager&,
const std::string&,
bool)>
1324 fieldPropIntTypeOnLeafAssigner_()
1326 const auto&
lookup = this->lookUpData_;
1333 void readMaterialParameters_()
1336 const auto& simulator = this->simulator();
1337 const auto& vanguard = simulator.vanguard();
1338 const auto& eclState = vanguard.eclState();
1341 OPM_BEGIN_PARALLEL_TRY_CATCH();
1342 this->updatePvtnum_();
1343 this->updateSatnum_();
1346 this->updateMiscnum_();
1348 this->updatePlmixnum_();
1350 OPM_END_PARALLEL_TRY_CATCH(
"Invalid region numbers: ", vanguard.gridView().comm());
1353 updateReferencePorosity_();
1354 this->referencePorosity_[1] = this->referencePorosity_[0];
1359 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
1360 materialLawManager_->initFromState(eclState);
1361 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1363 this-> lookupIdxOnLevelZeroAssigner_());
1367 void readThermalParameters_()
1369 if constexpr (enableEnergy)
1371 const auto& simulator = this->simulator();
1372 const auto& vanguard = simulator.vanguard();
1373 const auto& eclState = vanguard.eclState();
1376 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
1377 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1378 this-> fieldPropDoubleOnLeafAssigner_(),
1383 void updateReferencePorosity_()
1385 const auto& simulator = this->simulator();
1386 const auto& vanguard = simulator.vanguard();
1387 const auto& eclState = vanguard.eclState();
1389 std::size_t numDof = this->model().numGridDof();
1391 this->referencePorosity_[0].resize(numDof);
1393 const auto&
fp = eclState.fieldProps();
1394 const std::vector<double>
porvData =
this -> fieldPropDoubleOnLeafAssigner_()(
fp,
"PORV");
1395 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1396 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(dofIdx);
1402 Scalar dofVolume = simulator.model().dofTotalVolume(
sfcdofIdx);
1408 virtual void readInitialCondition_()
1411 const auto& simulator = this->simulator();
1412 const auto& vanguard = simulator.vanguard();
1413 const auto& eclState = vanguard.eclState();
1415 if (eclState.getInitConfig().hasEquil())
1416 readEquilInitialCondition_();
1418 readExplicitInitialCondition_();
1421 std::size_t
numElems = this->model().numGridDof();
1422 for (std::size_t elemIdx = 0; elemIdx <
numElems; ++elemIdx) {
1423 const auto& fs = asImp_().initialFluidStates()[elemIdx];
1424 if (!this->maxWaterSaturation_.empty() && waterPhaseIdx > -1)
1425 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
1426 if (!this->maxOilSaturation_.empty() && oilPhaseIdx > -1)
1427 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
1428 if (!this->minRefPressure_.empty() && refPressurePhaseIdx_() > -1)
1429 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
1433 virtual void readEquilInitialCondition_() = 0;
1434 virtual void readExplicitInitialCondition_() = 0;
1437 bool updateHysteresis_()
1439 if (!materialLawManager_->enableHysteresis())
1444 this->updateProperty_(
"FlowProblem::updateHysteresis_() failed:",
1461 Scalar getRockCompTransMultVal(std::size_t dofIdx)
const
1463 if (this->rockCompTransMultVal_.empty())
1466 return this->rockCompTransMultVal_[dofIdx];
1476 Scalar transmissibility;
1480 void updatePffDofData_()
1484 const Stencil& stencil,
1488 const auto& elementMapper = this->model().elementMapper();
1495 if constexpr (enableEnergy) {
1499 if constexpr (enableDiffusion)
1501 if (enableDispersion)
1506 pffDofData_.update(
distFn);
1511 void readBoundaryConditions_()
1513 const auto& simulator = this->simulator();
1514 const auto& vanguard = simulator.vanguard();
1515 const auto&
bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
1517 nonTrivialBoundaryConditions_ =
true;
1519 std::size_t
numCartDof = vanguard.cartesianSize();
1520 unsigned numElems = vanguard.gridView().size(0);
1523 for (
unsigned elemIdx = 0; elemIdx <
numElems; ++elemIdx)
1528 &vanguard](
const auto&
bcface,
1534 std::array<int, 3> tmp = {i,j,
k};
1543 std::vector<int>& data = bcindex_(
bcface.dir);
1544 const int index =
bcface.index;
1546 [&data,index](
int elemIdx)
1547 { data[elemIdx] = index; });
1554 Scalar limitNextTimeStepSize_(Scalar
dtNext)
const
1556 if constexpr (enableExperiments) {
1557 const auto& simulator = this->simulator();
1558 const auto& schedule = simulator.vanguard().schedule();
1562 Scalar maxTimeStepSize = Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60;
1564 if (this->enableTuning_) {
1566 maxTimeStepSize =
tuning.TSMAXZ;
1572 simulator.episodeStartTime() + simulator.episodeLength()
1573 - (simulator.startTime() + simulator.time());
1583 if (simulator.episodeStarts()) {
1586 const auto& events = simulator.vanguard().schedule()[
reportStepIdx].events();
1588 events.hasEvent(ScheduleEvents::NEW_WELL)
1589 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
1590 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
1591 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
1593 dtNext = std::min(
dtNext, this->maxTimeStepAfterWellEvent_);
1600 int refPressurePhaseIdx_()
const {
1601 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1604 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1608 return waterPhaseIdx;
1612 void updateRockCompTransMultVal_()
1614 const auto& model = this->simulator().model();
1615 std::size_t numGridDof = this->model().numGridDof();
1616 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
1618 const auto&
iq = *model.cachedIntensiveQuantities(
elementIdx, 0);
1629 template <
class LhsEval>
1633 if (this->rockCompTransMult_.empty() &&
this->rockCompTransMultWc_.empty())
1637 if (!this->rockTableIdx_.empty())
1640 const auto& fs = intQuants.fluidState();
1642 const auto&
rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
1643 if (!this->minRefPressure_.empty())
1649 if (!this->overburdenPressure_.empty())
1656 if (!this->rockCompTransMult_.empty())
1660 assert(!this->rockCompTransMultWc_.empty());
1667 typename Vanguard::TransmissibilityType transmissibilities_;
1669 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
1670 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
1672 GlobalEqVector drift_;
1674 WellModel wellModel_;
1675 AquiferModel aquiferModel_;
1683 std::array<std::vector<T>,6> data;
1685 void resize(std::size_t size, T
defVal)
1687 for (
auto&
d : data)
1691 const std::vector<T>& operator()(FaceDir::DirEnum dir)
const
1693 if (dir == FaceDir::DirEnum::Unknown)
1694 throw std::runtime_error(
"Tried to access BC data for the 'Unknown' direction");
1696 int div =
static_cast<int>(dir);
1697 while ((
div /= 2) >= 1)
1703 std::vector<T>& operator()(FaceDir::DirEnum dir)
1705 return const_cast<std::vector<T>&
>(std::as_const(*
this)(dir));
1709 virtual void handleSolventBC(
const BCProp::BCFace&, RateVector&)
const = 0;
1711 virtual void handlePolymerBC(
const BCProp::BCFace&, RateVector&)
const = 0;
1713 virtual void handleMicrBC(
const BCProp::BCFace&, RateVector&)
const = 0;
1715 virtual void handleOxygBC(
const BCProp::BCFace&, RateVector&)
const = 0;
1717 virtual void handleUreaBC(
const BCProp::BCFace&, RateVector&)
const = 0;
1720 bool nonTrivialBoundaryConditions_ =
false;
1721 bool first_step_ =
true;
Helper class for grid instantiation of ECL file-format using problems.
This is a "dummy" gradient calculator which does not do anything.
Collects necessary output values and pass it to opm-common's ECL output.
Computes the initial condition based on the EQUIL keyword from ECL.
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
A class which handles tracers as specified in by ECL.
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowGenericProblem.hpp:61
Scalar maxPolymerAdsorption(unsigned elemIdx) const
Returns the max polymer adsorption value.
Definition FlowGenericProblem_impl.hpp:731
unsigned pvtRegionIndex(unsigned elemIdx) const
Returns the index the relevant PVT region given a cell index.
Definition FlowGenericProblem_impl.hpp:690
static std::string briefDescription()
Returns a human readable description of the problem for the help message.
Definition FlowGenericProblem_impl.hpp:126
Scalar porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
Direct indexed access to the porosity.
Definition FlowGenericProblem_impl.hpp:327
Scalar rockCompressibility(unsigned globalSpaceIdx) const
Direct access to rock compressibility.
Definition FlowGenericProblem_impl.hpp:312
unsigned miscnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant MISC region given a cell index.
Definition FlowGenericProblem_impl.hpp:710
unsigned satnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant saturation function region given a cell index.
Definition FlowGenericProblem_impl.hpp:700
unsigned plmixnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant PLMIXNUM (for polymer module) region given a cell index.
Definition FlowGenericProblem_impl.hpp:720
bool shouldWriteOutput() const
Always returns true.
Definition FlowGenericProblem.hpp:277
static std::string helpPreamble(int, const char **argv)
Returns the string that is printed before the list of command line parameters in the help message.
Definition FlowGenericProblem_impl.hpp:111
bool shouldWriteRestartFile() const
Returns true if an eWoms restart file should be written to disk.
Definition FlowGenericProblem.hpp:286
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:94
const WellModel & wellModel() const
Returns a reference to the ECL well manager used by the problem.
Definition FlowProblem.hpp:1004
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
Direct access to the transmissibility between two elements.
Definition FlowProblem.hpp:520
static int handlePositionalParameter(std::function< void(const std::string &, const std::string &)> addKey, std::set< std::string > &seenParams, std::string &errorMsg, int, const char **argv, int paramIdx, int)
Handles positional command line parameters.
Definition FlowProblem.hpp:190
virtual void writeOutput(bool verbose)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:473
const DimMatrix & intrinsicPermeability(unsigned globalElemIdx) const
This method returns the intrinsic permeability tensor given a global element index.
Definition FlowProblem.hpp:502
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:818
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:652
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition FlowProblem.hpp:378
unsigned satnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:826
Scalar thermalHalfTransmissibilityOut(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:610
Scalar transmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:567
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn, const unsigned globalSpaceIdxOut) const
Definition FlowProblem.hpp:587
Scalar rockCompressibility(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:686
Scalar maxOilSaturation(unsigned globalDofIdx) const
Returns an element's historic maximum oil phase saturation that was observed during the simulation.
Definition FlowProblem.hpp:902
std::string name() const
The problem name.
Definition FlowProblem.hpp:857
LhsEval rockCompTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1096
void endIteration()
Called by the simulator after each Newton-Raphson iteration.
Definition FlowProblem.hpp:388
FlowProblem(Simulator &simulator)
Definition FlowProblem.hpp:209
const Vanguard::TransmissibilityType & eclTransmissibilities() const
Return a reference to the object that handles the "raw" transmissibilities.
Definition FlowProblem.hpp:633
Scalar nextTimeStepSize() const
Propose the size of the next time step to the simulator.
Definition FlowProblem.hpp:1025
LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1630
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:762
Scalar transmissibility(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:509
virtual void beginEpisode()
Called by the simulator before an episode begins.
Definition FlowProblem.hpp:288
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
Direct indexed acces to the depth of an degree of freedom [m].
Definition FlowProblem.hpp:677
LhsEval wellTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Return the well transmissibility multiplier due to rock changes.
Definition FlowProblem.hpp:1108
std::shared_ptr< const EclMaterialLawManager > materialLawManager() const
Returns the ECL material law manager.
Definition FlowProblem.hpp:774
unsigned plmixnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:842
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:864
Scalar thermalHalfTransmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:623
unsigned miscnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:834
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:488
Scalar maxPolymerAdsorption(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the max polymer adsorption value.
Definition FlowProblem.hpp:851
void beginTimeStep()
Called by the simulator before each time integration.
Definition FlowProblem.hpp:345
Scalar thermalHalfTransmissibilityIn(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:597
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
Sets an element's maximum oil phase saturation observed during the simulation.
Definition FlowProblem.hpp:919
const SolidEnergyLawParams & solidEnergyLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the parameters for the energy storage law of the rock.
Definition FlowProblem.hpp:749
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the dispersivity for a face i.e.
Definition FlowProblem.hpp:547
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblem.hpp:180
virtual void endEpisode()
Called by the simulator after the end of an episode.
Definition FlowProblem.hpp:450
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition FlowProblem.hpp:963
virtual void endTimeStep()
Called by the simulator after each time integration.
Definition FlowProblem.hpp:398
virtual void initialSolutionApplied()
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition FlowProblem.hpp:930
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:727
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:554
Scalar diffusivity(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:529
Scalar rockReferencePressure(unsigned globalSpaceIdx) const
Definition FlowProblem.hpp:705
void deserialize(Restarter &res)
This method restores the complete state of the problem and its sub-objects from disk.
Definition FlowProblem.hpp:254
std::shared_ptr< EclMaterialLawManager > materialLawManager()
Definition FlowProblem.hpp:810
Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the depth of an degree of freedom [m].
Definition FlowProblem.hpp:665
LhsEval rockCompPoroMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the porosity multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1052
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:577
void serialize(Restarter &res)
This method writes the complete state of the problem and its subobjects to disk.
Definition FlowProblem.hpp:273
Scalar rockReferencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:696
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the transmissibility for a face i.e.
Definition FlowProblem.hpp:540
A random-access container which stores data attached to a grid's degrees of freedom in a prefetch fri...
Definition pffgridvector.hh:49
This class is intend to be a relperm diagnostics, to detect wrong input of relperm table and endpoint...
Definition RelpermDiagnostics.hpp:51
void diagnosis(const EclipseState &eclState, const LevelCartesianIndexMapper &levelCartesianIndexMapper)
This function is used to diagnosis relperm in eclipse data file.
Definition RelpermDiagnostics.cpp:829
A class which handles tracers as specified in by ECL.
Definition TracerModel.hpp:72
This file contains definitions related to directional mobilities.
The base class for the element-centered finite-volume discretization scheme.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition propertysystem.hh:226
Definition FlowProblem.hpp:1682
Definition FlowProblem.hpp:1471