64 using typename FlowProblemType::Scalar;
65 using typename FlowProblemType::Simulator;
66 using typename FlowProblemType::GridView;
67 using typename FlowProblemType::FluidSystem;
68 using typename FlowProblemType::Vanguard;
71 using FlowProblemType::dim;
72 using FlowProblemType::dimWorld;
74 using FlowProblemType::numPhases;
75 using FlowProblemType::numComponents;
77 using FlowProblemType::gasPhaseIdx;
78 using FlowProblemType::oilPhaseIdx;
79 using FlowProblemType::waterPhaseIdx;
81 using typename FlowProblemType::Indices;
82 using typename FlowProblemType::PrimaryVariables;
84 using typename FlowProblemType::Evaluation;
85 using typename FlowProblemType::MaterialLaw;
86 using typename FlowProblemType::RateVector;
102 EclWriterType::registerParameters();
105 Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1
e-7);
108 Opm::CompositionalConfig::EOSType getEosType()
const
110 auto& simulator = this->simulator();
111 const auto& eclState = simulator.vanguard().eclState();
112 return eclState.compositionalConfig().eosType(0);
120 , thresholdPressures_(simulator)
122 eclWriter_ = std::make_unique<EclWriterType>(simulator);
123 enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
134 FlowProblemType::finishInit();
136 auto& simulator = this->simulator();
142 this->transmissibilities_.finishInit(
143 [&
vg = this->simulator().vanguard()](
const unsigned int it) {
return vg.gridIdxToEquilGridIdx(it); });
150 if (enableEclOutput_) {
151 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
152 std::function<
unsigned int(
unsigned int)>
equilGridToGrid = [&simulator](
unsigned int i) {
153 return simulator.vanguard().gridEquilIdxToGridIdx(i);
158 const auto& eclState = simulator.vanguard().eclState();
159 const auto& schedule = simulator.vanguard().schedule();
162 simulator.setStartTime(schedule.getStartTime());
163 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
169 simulator.setEpisodeIndex(-1);
170 simulator.setEpisodeLength(0.0);
175 this->gravity_ = 0.0;
176 if (Parameters::Get<Parameters::EnableGravity>())
177 this->gravity_[dim - 1] = 9.80665;
178 if (!eclState.getInitConfig().hasGravity())
179 this->gravity_[dim - 1] = 0.0;
181 if (this->enableTuning_) {
184 const auto&
tuning = schedule[0].tuning();
185 this->initialTimeStepSize_ =
tuning.TSINIT.has_value() ?
tuning.TSINIT.value() : -1.0;
186 this->maxTimeStepAfterWellEvent_ =
tuning.TMAXWC;
189 this->initFluidSystem_();
191 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
192 && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
193 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
196 this->readRockParameters_(simulator.vanguard().cellCenterDepths(), [&simulator](
const unsigned idx) {
197 std::array<int, dim> coords;
198 simulator.vanguard().cartesianCoordinate(idx, coords);
199 for (auto& c : coords) {
204 FlowProblemType::readMaterialParameters_();
205 FlowProblemType::readThermalParameters_();
208 if (enableEclOutput_) {
209 eclWriter_->writeInit();
212 const auto&
initconfig = eclState.getInitConfig();
214 readEclRestartSolution_();
216 this->readInitialCondition_();
218 FlowProblemType::updatePffDofData_();
221 const auto& vanguard = this->simulator().vanguard();
222 const auto& gridView = vanguard.gridView();
224 this->polymer_.maxAdsorption.resize(
numElements, 0.0);
238 if (this->enableVtkOutput_() && eclState.getIOConfig().initOnly()) {
239 simulator.setTimeStepSize(0.0);
247 simulator.startNextEpisode(schedule.seconds(1));
248 simulator.setEpisodeIndex(0);
249 simulator.setTimeStepIndex(0);
258 FlowProblemType::endTimeStep();
260 const bool isSubStep = !this->simulator().episodeWillBeOver();
263 this->eclWriter_->mutableOutputModule().invalidateLocalData();
266 const auto& grid = this->simulator().vanguard().gridView().grid();
268 using GridType = std::remove_cv_t<std::remove_reference_t<
decltype(grid)>>;
269 constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
270 if (!
isCpGrid || (grid.maxLevel() == 0)) {
271 this->eclWriter_->evalSummaryState(
isSubStep);
277 if (enableEclOutput_){
278 eclWriter_->writeReports(timer);
288 FlowProblemType::writeOutput(
verbose);
290 const bool isSubStep = !this->simulator().episodeWillBeOver();
293 if (enableEclOutput_) {
294 if (Parameters::Get<Parameters::EnableWriteAllSolutions>() || !
isSubStep) {
305 template <
class Context>
307 const Context& context,
312 if (!context.intersection(
spaceIdx).boundary())
317 if (this->nonTrivialBoundaryConditions()) {
318 throw std::logic_error(
"boundary condition is not supported by compostional modeling yet");
328 template <
class Context>
331 const unsigned globalDofIdx = context.globalSpaceIndex(
spaceIdx,
timeIdx);
332 const auto&
initial_fs = initialFluidStates_[globalDofIdx];
334 for (
unsigned p = 0;
p < numPhases; ++
p) {
346 if (!zmf_initialization_) {
347 for (
unsigned p = 0;
p < numPhases; ++
p) {
354 const auto&
eos_type = getEosType();
356 paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
357 paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
358 fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs,
paramCache, FluidSystem::oilPhaseIdx));
359 fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs,
paramCache, FluidSystem::gasPhaseIdx));
362 Dune::FieldVector<Scalar, numComponents>
z(0.0);
365 if (Indices::waterEnabled &&
phaseIdx ==
static_cast<unsigned int>(waterPhaseIdx)){
368 const auto saturation = fs.saturation(
phaseIdx);
371 tmp = max(tmp, 1
e-8);
393 const Scalar&
Ltmp = -1.0;
396 values.assignNaive(fs);
399 void addToSourceDense(RateVector&,
unsigned,
unsigned)
const override
404 const InitialFluidState& initialFluidState(
unsigned globalDofIdx)
const
405 {
return initialFluidStates_[globalDofIdx]; }
407 std::vector<InitialFluidState>& initialFluidStates()
408 {
return initialFluidStates_; }
410 const std::vector<InitialFluidState>& initialFluidStates()
const
411 {
return initialFluidStates_; }
415 assert( !thresholdPressures_.enableThresholdPressure() &&
416 " Threshold Pressures are not supported by compostional simulation ");
417 return thresholdPressures_;
420 const EclWriterType& eclWriter()
const
421 {
return *eclWriter_; }
423 EclWriterType& eclWriter()
424 {
return *eclWriter_; }
427 template<
class Serializer>
430 serializer(
static_cast<FlowProblemType&
>(*
this));
435 void updateExplicitQuantities_(
int ,
int ,
bool )
override
440 void readEquilInitialCondition_()
override
442 throw std::logic_error(
"Equilibration is not supported by compositional modeling yet");
445 void readEclRestartSolution_()
447 throw std::logic_error(
"Restarting is not supported by compositional modeling yet");
450 void readExplicitInitialCondition_()
override
452 readExplicitInitialConditionCompositional_();
455 void readExplicitInitialConditionCompositional_()
457 const auto& simulator = this->simulator();
458 const auto& vanguard = simulator.vanguard();
459 const auto& eclState = vanguard.eclState();
460 const auto&
fp = eclState.fieldProps();
463 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
464 "keyword if the model is initialized explicitly");
466 const bool has_xmf =
fp.has_double(
"XMF");
467 const bool has_ymf =
fp.has_double(
"YMF");
468 const bool has_zmf =
fp.has_double(
"ZMF");
470 throw std::runtime_error(
"The ECL input file requires the presence of ZMF or XMF and YMF "
471 "keyword if the model is initialized explicitly");
475 throw std::runtime_error(
"The ECL input file can not handle explicit initialization "
476 "with both ZMF and XMF or YMF");
480 throw std::runtime_error(
"The ECL input file needs XMF and YMF combined to do the explicit "
481 "initializtion when using XMF or YMF");
489 std::size_t numDof = this->model().numGridDof();
491 initialFluidStates_.resize(numDof);
499 const bool water_active = FluidSystem::phaseIsActive(waterPhaseIdx);
500 const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
501 const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
521 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
550 const std::array<Scalar, numPhases>
pc = {0};
552 if (!FluidSystem::phaseIsActive(
phaseIdx))
555 if (Indices::oilEnabled)
557 else if (Indices::gasEnabled)
559 else if (Indices::waterEnabled)
565 const auto&
xmfData =
fp.get_double(
"XMF");
566 const auto&
ymfData =
fp.get_double(
"YMF");
578 zmf_initialization_ =
true;
579 const auto&
zmfData =
fp.get_double(
"ZMF");
586 const auto ymf = (
dofFluidState.saturation(FluidSystem::gasPhaseIdx) > 0.) ?
zmf : Scalar{0};
590 const auto xmf = (
dofFluidState.saturation(FluidSystem::oilPhaseIdx) > 0.) ?
zmf : Scalar{0};
600 void handleSolventBC(
const BCProp::BCFace& , RateVector& )
const override
602 throw std::logic_error(
"solvent is disabled for compositional modeling and you're trying to add solvent to BC");
605 void handlePolymerBC(
const BCProp::BCFace& , RateVector& )
const override
607 throw std::logic_error(
"polymer is disabled for compositional modeling and you're trying to add polymer to BC");
610 void handleMicrBC(
const BCProp::BCFace& , RateVector& )
const override
612 throw std::logic_error(
"MICP is disabled for compositional modeling and you're trying to add microbes to BC");
615 void handleOxygBC(
const BCProp::BCFace& , RateVector& )
const override
617 throw std::logic_error(
"MICP is disabled for compositional modeling and you're trying to add oxygen to BC");
620 void handleUreaBC(
const BCProp::BCFace& , RateVector& )
const override
622 throw std::logic_error(
"MICP is disabled for compositional modeling and you're trying to add urea to BC");
627 std::vector<InitialFluidState> initialFluidStates_;
629 bool zmf_initialization_ {
false};
631 bool enableEclOutput_{
false};
632 std::unique_ptr<EclWriterType> eclWriter_;