FlowProblemBlackoil.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 Copyright 2023 INRIA
5 Copyright 2024 SINTEF Digital
6
7 This file is part of the Open Porous Media project (OPM).
8
9 OPM is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 2 of the License, or
12 (at your option) any later version.
13
14 OPM is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with OPM. If not, see <http://www.gnu.org/licenses/>.
21
22 Consult the COPYING file in the top-level source directory of this
23 module for the precise wording of the license and the list of
24 copyright holders.
25*/
31#ifndef OPM_FLOW_PROBLEM_BLACK_HPP
32#define OPM_FLOW_PROBLEM_BLACK_HPP
33
34#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
35#include <opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp>
36#include <opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp>
37#include <opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp>
38#include <opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp>
39#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp>
40#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp>
41
43
44#include <opm/output/eclipse/EclipseIO.hpp>
45
46#include <opm/input/eclipse/Units/Units.hpp>
47
55
57
58#if HAVE_DAMARIS
60#endif
61
62#include <algorithm>
63#include <cstddef>
64#include <functional>
65#include <limits>
66#include <memory>
67#include <stdexcept>
68#include <string>
69#include <string_view>
70#include <vector>
71
72namespace Opm {
73
80template <class TypeTag>
81class FlowProblemBlackoil : public FlowProblem<TypeTag>
82{
83 // TODO: the naming of the Types might be able to be adjusted
84public:
86
87private:
88 using typename FlowProblemType::Scalar;
89 using typename FlowProblemType::Simulator;
90 using typename FlowProblemType::GridView;
91 using typename FlowProblemType::FluidSystem;
92 using typename FlowProblemType::Vanguard;
94 using typename FlowProblemType::EqVector;
100
101 // TODO: potentially some cleaning up depending on the usage later here
117
121
125
127 using typename FlowProblemType::RateVector;
129 using typename FlowProblemType::Indices;
131 using typename FlowProblemType::ElementContext;
132
133 using typename FlowProblemType::MaterialLaw;
134 using typename FlowProblemType::DimMatrix;
135
136 enum { enableDissolvedGas = Indices::compositionSwitchIdx >= 0 };
137 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
138 enum { enableDisgasInWater = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
139
140 using SolventModule = BlackOilSolventModule<TypeTag>;
141 using PolymerModule = BlackOilPolymerModule<TypeTag>;
142 using FoamModule = BlackOilFoamModule<TypeTag>;
143 using BrineModule = BlackOilBrineModule<TypeTag>;
144 using ExtboModule = BlackOilExtboModule<TypeTag>;
145 using MICPModule = BlackOilMICPModule<TypeTag>;
146 using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
147 using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
148 using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
149 using ModuleParams = typename BlackOilLocalResidualTPFA<TypeTag>::ModuleParams;
150
151 using InitialFluidState = typename EquilInitializer<TypeTag>::ScalarFluidState;
152 using EclWriterType = EclWriter<TypeTag, OutputBlackOilModule<TypeTag> >;
153 using IndexTraits = typename FluidSystem::IndexTraitsType;
154#if HAVE_DAMARIS
155 using DamarisWriterType = DamarisWriter<TypeTag>;
156#endif
157
158
159public:
162
166 static void registerParameters()
167 {
169
171#if HAVE_DAMARIS
172 DamarisWriterType::registerParameters();
173#endif
175 }
176
180 explicit FlowProblemBlackoil(Simulator& simulator)
181 : FlowProblemType(simulator)
182 , thresholdPressures_(simulator)
183 , mixControls_(simulator.vanguard().schedule())
184 , actionHandler_(simulator.vanguard().eclState(),
185 simulator.vanguard().schedule(),
186 simulator.vanguard().actionState(),
187 simulator.vanguard().summaryState(),
188 this->wellModel_,
189 simulator.vanguard().grid().comm())
190 {
191 this->model().addOutputModule(std::make_unique<VtkTracerModule<TypeTag>>(simulator));
192
193 // Tell the black-oil extensions to initialize their internal data structures
194 const auto& vanguard = simulator.vanguard();
195
196 BlackOilBrineParams<Scalar> brineParams;
197 brineParams.template initFromState<enableBrine,
198 enableSaltPrecipitation>(vanguard.eclState());
199 BrineModule::setParams(std::move(brineParams));
200
201 DiffusionModule::initFromState(vanguard.eclState());
202 DispersionModule::initFromState(vanguard.eclState());
203
204 BlackOilExtboParams<Scalar> extboParams;
205 extboParams.template initFromState<enableExtbo>(vanguard.eclState());
206 ExtboModule::setParams(std::move(extboParams));
207
209 foamParams.template initFromState<enableFoam>(vanguard.eclState());
210 FoamModule::setParams(std::move(foamParams));
211
213 micpParams.template initFromState<enableMICP>(vanguard.eclState());
214 MICPModule::setParams(std::move(micpParams));
215
216 BlackOilPolymerParams<Scalar> polymerParams;
217 polymerParams.template initFromState<enablePolymer, enablePolymerMolarWeight>(vanguard.eclState());
218 PolymerModule::setParams(std::move(polymerParams));
219
220 BlackOilSolventParams<Scalar> solventParams;
221 solventParams.template initFromState<enableSolvent>(vanguard.eclState(), vanguard.schedule());
222 SolventModule::setParams(std::move(solventParams));
223
224 // create the ECL writer
225 eclWriter_ = std::make_unique<EclWriterType>(simulator);
226 enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
227
228#if HAVE_DAMARIS
229 // create Damaris writer
230 damarisWriter_ = std::make_unique<DamarisWriterType>(simulator);
231 enableDamarisOutput_ = Parameters::Get<Parameters::EnableDamarisOutput>();
232#endif
233 }
234
238 void beginEpisode() override
239 {
241
242 auto& simulator = this->simulator();
243
244 const int episodeIdx = simulator.episodeIndex();
245 const auto& schedule = simulator.vanguard().schedule();
246
247 // Evaluate UDQ assign statements to make sure the settings are
248 // available as UDA controls for the current report step.
249 this->actionHandler_
250 .evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
251
252 if (episodeIdx >= 0) {
253 const auto& oilVap = schedule[episodeIdx].oilvap();
254 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
255 FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2());
256 }
257 else {
258 FluidSystem::setVapPars(0.0, 0.0);
259 }
260 }
261
262 ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), schedule, episodeIdx,
263 this->moduleParams_.convectiveMixingModuleParam);
264 }
265
270 {
271 // TODO: there should be room to remove duplication for this
272 // function, but there is relatively complicated logic in the
273 // function calls here. Some refactoring is needed.
274 FlowProblemType::finishInit();
275
276 auto& simulator = this->simulator();
277
278 auto finishTransmissibilities = [updated = false, this]() mutable
279 {
280 if (updated) { return; }
281
282 this->transmissibilities_.finishInit([&vg = this->simulator().vanguard()](const unsigned int it) {
283 return vg.gridIdxToEquilGridIdx(it);
284 });
285
286 updated = true;
287 };
288
289 // calculating the TRANX, TRANY, TRANZ and NNC for output purpose
290 // for parallel running, it is based on global trans_
291 // for serial running, it is based on the transmissibilities_
292 // we try to avoid for the parallel running, has both global trans_ and transmissibilities_ allocated at the same time
293 if (enableEclOutput_) {
294 if (simulator.vanguard().grid().comm().size() > 1) {
295 if (simulator.vanguard().grid().comm().rank() == 0)
296 eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
297 } else {
298 finishTransmissibilities();
299 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
300 }
301
302 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
303 return simulator.vanguard().gridEquilIdxToGridIdx(i);
304 };
305
306 this->eclWriter_->extractOutputTransAndNNC(equilGridToGrid);
307 }
308 simulator.vanguard().releaseGlobalTransmissibilities();
309
310 const auto& eclState = simulator.vanguard().eclState();
311 const auto& schedule = simulator.vanguard().schedule();
312
313 // Set the start time of the simulation
314 simulator.setStartTime(schedule.getStartTime());
315 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
316
317 // We want the episode index to be the same as the report step index to make
318 // things simpler, so we have to set the episode index to -1 because it is
319 // incremented by endEpisode(). The size of the initial time step and
320 // length of the initial episode is set to zero for the same reason.
321 simulator.setEpisodeIndex(-1);
322 simulator.setEpisodeLength(0.0);
323
324 // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
325 // disables gravity, else the standard value of the gravity constant at sea level
326 // on earth is used
327 this->gravity_ = 0.0;
328 if (Parameters::Get<Parameters::EnableGravity>() &&
329 eclState.getInitConfig().hasGravity())
330 {
331 // unit::gravity is 9.80665 m^2/s--i.e., standard measure at Tellus equator.
332 this->gravity_[dim - 1] = unit::gravity;
333 }
334
335 if (this->enableTuning_) {
336 // if support for the TUNING keyword is enabled, we get the initial time
337 // steping parameters from it instead of from command line parameters
338 const auto& tuning = schedule[0].tuning();
339 this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
340 this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
341 }
342
343 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
344 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
345 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
346 }
347
348 this->readRockParameters_(simulator.vanguard().cellCenterDepths(),
349 [&simulator](const unsigned idx)
350 {
351 std::array<int,dim> coords;
352 simulator.vanguard().cartesianCoordinate(idx, coords);
353 std::transform(coords.begin(), coords.end(), coords.begin(),
354 [](const auto c) { return c + 1; });
355 return coords;
356 });
357
360
361 // write the static output files (EGRID, INIT)
362 if (enableEclOutput_) {
363 this->eclWriter_->writeInit();
364 }
365
366 finishTransmissibilities();
367
368 const auto& initconfig = eclState.getInitConfig();
369 this->tracerModel_.init(initconfig.restartRequested());
370 if (initconfig.restartRequested()) {
372 }
373 else {
374 this->readInitialCondition_();
375 }
376
377 this->tracerModel_.prepareTracerBatches();
378
379 this->updatePffDofData_();
380
381 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
382 const auto& vanguard = this->simulator().vanguard();
383 const auto& gridView = vanguard.gridView();
384 const int numElements = gridView.size(/*codim=*/0);
385 this->polymer_.maxAdsorption.resize(numElements, 0.0);
386 }
387
389
390 // compute and set eq weights based on initial b values
392
393 if (this->enableDriftCompensation_) {
394 this->drift_.resize(this->model().numGridDof());
395 this->drift_ = 0.0;
396 }
397
398 // after finishing the initialization and writing the initial solution, we move
399 // to the first "real" episode/report step
400 // for restart the episode index and start is already set
401 if (!initconfig.restartRequested() && !eclState.getIOConfig().initOnly()) {
402 simulator.startNextEpisode(schedule.seconds(1));
403 simulator.setEpisodeIndex(0);
404 simulator.setTimeStepIndex(0);
405 }
406
407 if (Parameters::Get<Parameters::CheckSatfuncConsistency>() &&
409 {
410 // User requested that saturation functions be checked for
411 // consistency and essential/critical requirements are not met.
412 // Abort simulation run.
413 //
414 // Note: We need synchronisation here lest ranks other than the
415 // I/O rank throw exceptions too early thereby risking an
416 // incomplete failure report being shown to the user.
417 this->simulator().vanguard().grid().comm().barrier();
418
419 throw std::domain_error {
420 "Saturation function end-points do not "
421 "meet requisite consistency conditions"
422 };
423 }
424
425 // TODO: move to the end for later refactoring of the function finishInit()
426 //
427 // deal with DRSDT
428 this->mixControls_.init(this->model().numGridDof(),
429 this->episodeIndex(),
430 eclState.runspec().tabdims().getNumPVTTables());
431
432 if (this->enableVtkOutput_() && eclState.getIOConfig().initOnly()) {
433 simulator.setTimeStepSize(0.0);
434 simulator.model().applyInitialSolution();
436 }
437
438 if (!eclState.getIOConfig().initOnly()) {
439 if (!this->enableTuning_ && eclState.getSimulationConfig().anyTUNING()) {
440 OpmLog::info("\nThe deck has TUNING in the SCHEDULE section, but "
441 "it is ignored due\nto the flag --enable-tuning=false. "
442 "Set this flag to true to activate it.\n"
443 "Manually tuning the simulator with the TUNING keyword may "
444 "increase run time.\nIt is recommended using the simulator's "
445 "default tuning (--enable-tuning=false).");
446 }
447 }
448 }
449
453 void endTimeStep() override
454 {
455 FlowProblemType::endTimeStep();
456 this->endStepApplyAction();
457 }
458
460 {
461 // After the solution is updated, the values in output module needs
462 // also updated.
463 this->eclWriter().mutableOutputModule().invalidateLocalData();
464
465 // For CpGrid with LGRs, ecl/vtk output is not supported yet.
466 const auto& grid = this->simulator().vanguard().gridView().grid();
467
468 using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
469 constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
470 if (!isCpGrid || (grid.maxLevel() == 0)) {
471 this->eclWriter_->evalSummaryState(!this->episodeWillBeOver());
472 }
473
474 {
475 OPM_TIMEBLOCK(applyActions);
476
477 const int episodeIdx = this->episodeIndex();
478 auto& simulator = this->simulator();
479
480 // Clear out any existing events as these have already been
481 // processed when we're running an action block
482 this->simulator().vanguard().schedule().clearEvents(episodeIdx);
483
484 // Re-ordering in case of Alugrid
485 this->actionHandler_
486 .applyActions(episodeIdx, simulator.time() + simulator.timeStepSize(),
487 [this](const bool global)
488 {
489 using TransUpdateQuantities = typename
490 Vanguard::TransmissibilityType::TransUpdateQuantities;
491
492 this->transmissibilities_
493 .update(global, TransUpdateQuantities::All,
494 [&vg = this->simulator().vanguard()]
495 (const unsigned int i)
496 {
497 return vg.gridIdxToEquilGridIdx(i);
498 });
499 });
500 }
501 }
502
506 void endEpisode() override
507 {
508 OPM_TIMEBLOCK(endEpisode);
509
510 // Rerun UDQ assignents following action processing on the final
511 // time step of this episode to make sure that any UDQ ASSIGN
512 // operations triggered in action blocks take effect. This is
513 // mainly to work around a shortcoming of the ScheduleState copy
514 // constructor which clears pending UDQ assignments under the
515 // assumption that all such assignments have been processed. If an
516 // action block happens to trigger on the final time step of an
517 // episode and that action block runs a UDQ assignment, then that
518 // assignment would be dropped and the rest of the simulator will
519 // never see its effect without this hack.
520 this->actionHandler_
521 .evalUDQAssignments(this->episodeIndex(), this->simulator().vanguard().udqState());
522
523 FlowProblemType::endEpisode();
524 }
525
526 void writeReports(const SimulatorTimer& timer)
527 {
528 if (this->enableEclOutput_) {
529 this->eclWriter_->writeReports(timer);
530 }
531 }
532
533
538 void writeOutput(const bool verbose) override
539 {
540 FlowProblemType::writeOutput(verbose);
541
542 const auto isSubStep = !this->episodeWillBeOver();
543
544 auto localCellData = data::Solution {};
545
546#if HAVE_DAMARIS
547 // N.B. the Damaris output has to be done before the ECL output as the ECL one
548 // does all kinds of std::move() relocation of data
549 if (this->enableDamarisOutput_ && (this->damarisWriter_ != nullptr)) {
550 this->damarisWriter_->writeOutput(localCellData, isSubStep);
551 }
552#endif
553
554 if (this->enableEclOutput_ && (this->eclWriter_ != nullptr)) {
555 this->eclWriter_->writeOutput(std::move(localCellData), isSubStep);
556 }
557 }
558
560 {
561 OPM_TIMEBLOCK(finalizeOutput);
562 // this will write all pending output to disk
563 // to avoid corruption of output files
564 eclWriter_.reset();
565 }
566
567
572 {
573 FlowProblemType::initialSolutionApplied();
574
575 // let the object for threshold pressures initialize itself. this is done only at
576 // this point, because determining the threshold pressures may require to access
577 // the initial solution.
578 this->thresholdPressures_.finishInit();
579
580 // For CpGrid with LGRs, ecl-output is not supported yet.
581 const auto& grid = this->simulator().vanguard().gridView().grid();
582
583 using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
584 constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
585 // Skip - for now - calculate the initial fip values for CpGrid with LGRs.
586 if (!isCpGrid || (grid.maxLevel() == 0)) {
587 if (this->simulator().episodeIndex() == 0) {
588 eclWriter_->writeInitialFIPReport();
589 }
590 }
591 }
592
594 unsigned globalDofIdx,
595 unsigned timeIdx) const override
596 {
597 this->aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
598
599 // Add source term from deck
600 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
601 std::array<int,3> ijk;
602 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
603
604 if (source.hasSource(ijk)) {
605 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
606 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
607 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
608 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
609
610 for (unsigned i = 0; i < phidx_map.size(); ++i) {
611 const auto phaseIdx = phidx_map[i];
612 const auto sourceComp = sc_map[i];
613 const auto compIdx = cidx_map[i];
614 if (!FluidSystem::phaseIsActive(phaseIdx)) {
615 continue;
616 }
617 Scalar mass_rate = source.rate(ijk, sourceComp) / this->model().dofTotalVolume(globalDofIdx);
618 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
619 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
620 }
621 rate[FluidSystem::canonicalToActiveCompIdx(compIdx)] += mass_rate;
622 }
623
624 if constexpr (enableSolvent) {
625 Scalar mass_rate = source.rate(ijk, SourceComponent::SOLVENT) / this->model().dofTotalVolume(globalDofIdx);
626 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
627 const auto& solventPvt = SolventModule::solventPvt();
628 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
629 }
630 rate[Indices::contiSolventEqIdx] += mass_rate;
631 }
632 if constexpr (enablePolymer) {
633 rate[Indices::polymerConcentrationIdx] += source.rate(ijk, SourceComponent::POLYMER) / this->model().dofTotalVolume(globalDofIdx);
634 }
635 if constexpr (enableMICP) {
636 rate[Indices::microbialConcentrationIdx] += source.rate(ijk, SourceComponent::MICR) / this->model().dofTotalVolume(globalDofIdx);
637 rate[Indices::oxygenConcentrationIdx] += source.rate(ijk, SourceComponent::OXYG) / this->model().dofTotalVolume(globalDofIdx);
638 rate[Indices::ureaConcentrationIdx] += source.rate(ijk, SourceComponent::UREA) / (this->model().dofTotalVolume(globalDofIdx));
639 }
640 if constexpr (enableEnergy) {
641 for (unsigned i = 0; i < phidx_map.size(); ++i) {
642 const auto phaseIdx = phidx_map[i];
643 if (!FluidSystem::phaseIsActive(phaseIdx)) {
644 continue;
645 }
646 const auto sourceComp = sc_map[i];
647 const auto source_hrate = source.hrate(ijk, sourceComp);
648 if (source_hrate) {
649 rate[Indices::contiEnergyEqIdx] += source_hrate.value() / this->model().dofTotalVolume(globalDofIdx);
650 } else {
651 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, /*timeIdx*/ 0);
652 auto fs = intQuants.fluidState();
653 // if temperature is not set, use cell temperature as default
654 const auto source_temp = source.temperature(ijk, sourceComp);
655 if (source_temp) {
656 Scalar temperature = source_temp.value();
657 fs.setTemperature(temperature);
658 }
659 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
660 Scalar mass_rate = source.rate(ijk, sourceComp)/ this->model().dofTotalVolume(globalDofIdx);
661 Scalar energy_rate = getValue(h)*mass_rate;
662 rate[Indices::contiEnergyEqIdx] += energy_rate;
663 }
664 }
665 }
666 }
667
668 // if requested, compensate systematic mass loss for cells which were "well
669 // behaved" in the last time step
670 if (this->enableDriftCompensation_) {
671 const auto& simulator = this->simulator();
672 const auto& model = this->model();
673
674 // we use a lower tolerance for the compensation too
675 // assure the added drift from the last step does not
676 // cause convergence issues on the current step
677 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
678 Scalar poro = this->porosity(globalDofIdx, timeIdx);
679 Scalar dt = simulator.timeStepSize();
680 EqVector dofDriftRate = this->drift_[globalDofIdx];
681 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
682
683 // restrict drift compensation to the CNV tolerance
684 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
685 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
686 if (cnv > maxCompensation) {
687 dofDriftRate[eqIdx] *= maxCompensation/cnv;
688 }
689 }
690
691 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
692 rate[eqIdx] -= dofDriftRate[eqIdx];
693 }
694 }
695
701 template <class LhsEval>
702 LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
703 {
704 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
705 if constexpr (enableSaltPrecipitation) {
706 const auto& fs = intQuants.fluidState();
707 unsigned tableIdx = this->simulator().problem().satnumRegionIndex(elementIdx);
708 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
709 porosityFactor = min(porosityFactor, 1.0);
710 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
711 return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
712 }
713 else if constexpr (enableMICP) {
714 return intQuants.permFactor().value();
715 }
716 else {
717 return 1.0;
718 }
719 }
720
721 // temporary solution to facilitate output of initial state from flow
722 const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
723 { return initialFluidStates_[globalDofIdx]; }
724
725 std::vector<InitialFluidState>& initialFluidStates()
726 { return initialFluidStates_; }
727
728 const std::vector<InitialFluidState>& initialFluidStates() const
729 { return initialFluidStates_; }
730
731 const EclipseIO& eclIO() const
732 { return eclWriter_->eclIO(); }
733
735 { return eclWriter_->setSubStepReport(report); }
736
738 { return eclWriter_->setSimulationReport(report); }
739
740 InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const
741 {
742 OPM_TIMEBLOCK_LOCAL(boundaryFluidState);
743 const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop;
744 if (bcprop.size() > 0) {
745 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
746
747 // index == 0: no boundary conditions for this
748 // global cell and direction
749 if (this->bcindex_(dir)[globalDofIdx] == 0)
750 return initialFluidStates_[globalDofIdx];
751
752 const auto& bc = bcprop[this->bcindex_(dir)[globalDofIdx]];
753 if (bc.bctype == BCType::DIRICHLET )
754 {
755 InitialFluidState fluidState;
756 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
757 fluidState.setPvtRegionIndex(pvtRegionIdx);
758
759 switch (bc.component) {
760 case BCComponent::OIL:
761 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
762 throw std::logic_error("oil is not active and you're trying to add oil BC");
763
764 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
765 break;
766 case BCComponent::GAS:
767 if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
768 throw std::logic_error("gas is not active and you're trying to add gas BC");
769
770 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
771 break;
772 case BCComponent::WATER:
773 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
774 throw std::logic_error("water is not active and you're trying to add water BC");
775
776 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
777 break;
778 case BCComponent::SOLVENT:
779 case BCComponent::POLYMER:
780 case BCComponent::MICR:
781 case BCComponent::OXYG:
782 case BCComponent::UREA:
784 throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
785 }
786 fluidState.setTotalSaturation(1.0);
787 double pressure = initialFluidStates_[globalDofIdx].pressure(this->refPressurePhaseIdx_());
788 const auto pressure_input = bc.pressure;
789 if (pressure_input) {
790 pressure = *pressure_input;
791 }
792
793 std::array<Scalar, numPhases> pc = {0};
794 const auto& matParams = this->materialLawParams(globalDofIdx);
795 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
796 Valgrind::CheckDefined(pressure);
797 Valgrind::CheckDefined(pc);
798 for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
799 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
800 if (Indices::oilEnabled)
801 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
802 else if (Indices::gasEnabled)
803 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
804 else if (Indices::waterEnabled)
805 //single (water) phase
806 fluidState.setPressure(phaseIdx, pressure);
807 }
808
809 double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature
810 const auto temperature_input = bc.temperature;
811 if(temperature_input)
812 temperature = *temperature_input;
813 fluidState.setTemperature(temperature);
814
815 if constexpr (enableDissolvedGas) {
816 if (FluidSystem::enableDissolvedGas()) {
817 fluidState.setRs(0.0);
818 fluidState.setRv(0.0);
819 }
820 }
821 if constexpr (enableDisgasInWater) {
822 if (FluidSystem::enableDissolvedGasInWater()) {
823 fluidState.setRsw(0.0);
824 }
825 }
826 if constexpr (enableVapwat) {
827 if (FluidSystem::enableVaporizedWater()) {
828 fluidState.setRvw(0.0);
829 }
830 }
831
832 for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
833 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
834
835 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
836 fluidState.setInvB(phaseIdx, b);
837
838 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
839 fluidState.setDensity(phaseIdx, rho);
840 if constexpr (enableEnergy) {
841 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
842 fluidState.setEnthalpy(phaseIdx, h);
843 }
844 }
845 fluidState.checkDefined();
846 return fluidState;
847 }
848 }
849 return initialFluidStates_[globalDofIdx];
850 }
851
852
854 { return *eclWriter_; }
855
857 { return *eclWriter_; }
858
863 Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
864 {
865 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
866 this->episodeIndex(),
867 this->pvtRegionIndex(globalDofIdx));
868 }
869
874 Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
875 {
876 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
877 this->episodeIndex(),
878 this->pvtRegionIndex(globalDofIdx));
879 }
880
890 {
891 int episodeIdx = this->episodeIndex();
892 return !this->mixControls_.drsdtActive(episodeIdx) &&
893 !this->mixControls_.drvdtActive(episodeIdx) &&
894 this->rockCompPoroMultWc_.empty() &&
895 this->rockCompPoroMult_.empty();
896 }
897
904 template <class Context>
905 void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
906 {
907 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
908
909 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
910 values.assignNaive(initialFluidStates_[globalDofIdx]);
911
912 SolventModule::assignPrimaryVars(values,
913 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
914 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
915
916 if constexpr (enablePolymer)
917 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
918
919 if constexpr (enablePolymerMolarWeight)
920 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
921
922 if constexpr (enableBrine) {
923 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
924 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
925 }
926 else {
927 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
928 }
929 }
930
931 if constexpr (enableMICP){
932 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
933 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
934 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
935 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
936 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
937 }
938
939 values.checkDefined();
940 }
941
942
943 Scalar drsdtcon(unsigned elemIdx, int episodeIdx) const
944 {
945 return this->mixControls_.drsdtcon(elemIdx, episodeIdx,
946 this->pvtRegionIndex(elemIdx));
947 }
948
949 bool drsdtconIsActive(unsigned elemIdx, int episodeIdx) const
950 {
951 return this->mixControls_.drsdtConvective(episodeIdx, this->pvtRegionIndex(elemIdx));
952 }
953
959 template <class Context>
960 void boundary(BoundaryRateVector& values,
961 const Context& context,
962 unsigned spaceIdx,
963 unsigned timeIdx) const
964 {
965 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
966 if (!context.intersection(spaceIdx).boundary())
967 return;
968
969 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
970 values.setNoFlow();
971 else {
972 // in the energy case we need to specify a non-trivial boundary condition
973 // because the geothermal gradient needs to be maintained. for this, we
974 // simply assume the initial temperature at the boundary and specify the
975 // thermal flow accordingly. in this context, "thermal flow" means energy
976 // flow due to a temerature gradient while assuming no-flow for mass
977 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
978 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
979 values.setThermalFlow(context, spaceIdx, timeIdx, this->initialFluidStates_[globalDofIdx] );
980 }
981
982 if (this->nonTrivialBoundaryConditions()) {
983 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
984 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
985 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
986 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
987 const auto [type, massrate] = this->boundaryCondition(globalDofIdx, indexInInside);
988 if (type == BCType::THERMAL)
989 values.setThermalFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
990 else if (type == BCType::FREE || type == BCType::DIRICHLET)
991 values.setFreeFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
992 else if (type == BCType::RATE)
993 values.setMassRate(massrate, pvtRegionIdx);
994 }
995 }
996
1001 void readSolutionFromOutputModule(const int restart_step, bool fip_init)
1002 {
1003 auto& simulator = this->simulator();
1004 const auto& eclState = simulator.vanguard().eclState();
1005
1006 std::size_t numElems = this->model().numGridDof();
1007 this->initialFluidStates_.resize(numElems);
1008 if constexpr (enableSolvent) {
1009 this->solventSaturation_.resize(numElems, 0.0);
1010 this->solventRsw_.resize(numElems, 0.0);
1011 }
1012
1013 if constexpr (enablePolymer)
1014 this->polymer_.concentration.resize(numElems, 0.0);
1015
1016 if constexpr (enablePolymerMolarWeight) {
1017 const std::string msg {"Support of the RESTART for polymer molecular weight "
1018 "is not implemented yet. The polymer weight value will be "
1019 "zero when RESTART begins"};
1020 OpmLog::warning("NO_POLYMW_RESTART", msg);
1021 this->polymer_.moleWeight.resize(numElems, 0.0);
1022 }
1023
1024 if constexpr (enableMICP) {
1025 this->micp_.resize(numElems);
1026 }
1027
1028 // Initialize mixing controls before trying to set any lastRx valuesx
1029 this->mixControls_.init(numElems, restart_step, eclState.runspec().tabdims().getNumPVTTables());
1030
1031 if constexpr (enableMICP) {
1032 this->micp_ = this->eclWriter_->outputModule().getMICP().getSolution();
1033 }
1034
1035 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1036 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1037 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
1038 this->eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
1039 this->eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
1040
1041 // Note: Function processRestartSaturations_() mutates the
1042 // 'ssol' argument--the value from the restart file--if solvent
1043 // is enabled. Then, store the updated solvent saturation into
1044 // 'solventSaturation_'. Otherwise, just pass a dummy value to
1045 // the function and discard the unchanged result. Do not index
1046 // into 'solventSaturation_' unless solvent is enabled.
1047 {
1048 auto ssol = enableSolvent
1049 ? this->eclWriter_->outputModule().getSolventSaturation(elemIdx)
1050 : Scalar(0);
1051
1052 this->processRestartSaturations_(elemFluidState, ssol);
1053
1054 if constexpr (enableSolvent) {
1055 this->solventSaturation_[elemIdx] = ssol;
1056 this->solventRsw_[elemIdx] = this->eclWriter_->outputModule().getSolventRsw(elemIdx);
1057 }
1058 }
1059
1060 // For CO2STORE and H2STORE we need to set the initial temperature for isothermal simulations
1061 bool isThermal = eclState.getSimulationConfig().isThermal();
1062 bool needTemperature = (eclState.runspec().co2Storage() || eclState.runspec().h2Storage());
1063 if (!isThermal && needTemperature) {
1064 const auto& fp = simulator.vanguard().eclState().fieldProps();
1065 elemFluidState.setTemperature(fp.get_double("TEMPI")[elemIdx]);
1066 }
1067
1068 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
1069
1070 if constexpr (enablePolymer)
1071 this->polymer_.concentration[elemIdx] = this->eclWriter_->outputModule().getPolymerConcentration(elemIdx);
1072 // if we need to restart for polymer molecular weight simulation, we need to add related here
1073 }
1074
1075 const int episodeIdx = this->episodeIndex();
1076 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
1077
1078 // assign the restart solution to the current solution. note that we still need
1079 // to compute real initial solution after this because the initial fluid states
1080 // need to be correct for stuff like boundary conditions.
1081 auto& sol = this->model().solution(/*timeIdx=*/0);
1082 const auto& gridView = this->gridView();
1083 ElementContext elemCtx(simulator);
1084 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1085 elemCtx.updatePrimaryStencil(elem);
1086 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1087 this->initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0);
1088 }
1089
1090 // make sure that the ghost and overlap entities exhibit the correct
1091 // solution. alternatively, this could be done in the loop above by also
1092 // considering non-interior elements. Since the initial() method might not work
1093 // 100% correctly for such elements, let's play safe and explicitly synchronize
1094 // using message passing.
1095 this->model().syncOverlap();
1096
1097 if (fip_init) {
1098 this->updateReferencePorosity_();
1099 this->mixControls_.init(this->model().numGridDof(),
1100 this->episodeIndex(),
1101 eclState.runspec().tabdims().getNumPVTTables());
1102 }
1103 }
1104
1108 Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
1109 { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
1110
1112 { return thresholdPressures_; }
1113
1115 { return thresholdPressures_; }
1116
1117 const ModuleParams& moduleParams() const
1118 {
1119 return moduleParams_;
1120 }
1121
1122 template<class Serializer>
1123 void serializeOp(Serializer& serializer)
1124 {
1125 serializer(static_cast<FlowProblemType&>(*this));
1126 serializer(mixControls_);
1127 serializer(*eclWriter_);
1128 }
1129
1130protected:
1131 void updateExplicitQuantities_(int episodeIdx, int timeStepSize, const bool first_step_after_restart) override
1132 {
1133 this->updateExplicitQuantities_(first_step_after_restart);
1134
1135 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
1136 updateMaxPolymerAdsorption_();
1137
1138 mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize);
1139 }
1140
1142 {
1143 // we need to update the max polymer adsoption data for all elements
1144 this->updateProperty_("FlowProblemBlackoil::updateMaxPolymerAdsorption_() failed:",
1145 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1146 {
1147 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
1148 });
1149 }
1150
1151 bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1152 {
1153 const Scalar pa = scalarValue(iq.polymerAdsorption());
1154 auto& mpa = this->polymer_.maxAdsorption;
1155 if (mpa[compressedDofIdx] < pa) {
1156 mpa[compressedDofIdx] = pa;
1157 return true;
1158 } else {
1159 return false;
1160 }
1161 }
1162
1164 {
1165 std::vector<Scalar> sumInvB(numPhases, 0.0);
1166 const auto& gridView = this->gridView();
1167 ElementContext elemCtx(this->simulator());
1168 for(const auto& elem: elements(gridView, Dune::Partitions::interior)) {
1169 elemCtx.updatePrimaryStencil(elem);
1170 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1171 const auto& dofFluidState = this->initialFluidStates_[elemIdx];
1172 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1173 if (!FluidSystem::phaseIsActive(phaseIdx))
1174 continue;
1175
1176 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
1177 }
1178 }
1179
1180 std::size_t numDof = this->model().numGridDof();
1181 const auto& comm = this->simulator().vanguard().grid().comm();
1182 comm.sum(sumInvB.data(),sumInvB.size());
1183 Scalar numTotalDof = comm.sum(numDof);
1184
1185 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1186 if (!FluidSystem::phaseIsActive(phaseIdx))
1187 continue;
1188
1189 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
1190 const unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
1191 const unsigned activeSolventCompIdx = FluidSystem::canonicalToActiveCompIdx(solventCompIdx);
1192 this->model().setEqWeight(activeSolventCompIdx, avgB);
1193 }
1194 }
1195
1196 // update the parameters needed for DRSDT and DRVDT
1198 {
1199 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1200 // update the "last Rs" values for all elements, including the ones in the ghost
1201 // and overlap regions
1202 int episodeIdx = this->episodeIndex();
1203 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1204 this->mixControls_.drsdtActive(episodeIdx),
1205 this->mixControls_.drvdtActive(episodeIdx)};
1206 if (!active[0] && !active[1] && !active[2]) {
1207 return false;
1208 }
1209
1210 this->updateProperty_("FlowProblemBlackoil::updateCompositionChangeLimits_()) failed:",
1211 [this,episodeIdx,active](unsigned compressedDofIdx,
1212 const IntensiveQuantities& iq)
1213 {
1214 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1215 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1216 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1217 this->mixControls_.update(compressedDofIdx,
1218 iq,
1219 episodeIdx,
1220 this->gravity_[dim - 1],
1221 perm[dim - 1][dim - 1],
1222 distZ,
1223 pvtRegionIdx);
1224 }
1225 );
1226
1227 return true;
1228 }
1229
1231 {
1232 // Throw an exception if the grid has LGRs. Refined grid are not supported for restart.
1233 if(this->simulator().vanguard().grid().maxLevel() > 0) {
1234 throw std::invalid_argument("Refined grids are not yet supported for restart ");
1235 }
1236
1237 // Set the start time of the simulation
1238 auto& simulator = this->simulator();
1239 const auto& schedule = simulator.vanguard().schedule();
1240 const auto& eclState = simulator.vanguard().eclState();
1241 const auto& initconfig = eclState.getInitConfig();
1242 const int restart_step = initconfig.getRestartStep();
1243 {
1244 simulator.setTime(schedule.seconds(restart_step));
1245
1246 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
1247 schedule.stepLength(restart_step));
1248 simulator.setEpisodeIndex(restart_step);
1249 }
1250 this->eclWriter_->beginRestart();
1251
1252 Scalar dt = std::min(this->eclWriter_->restartTimeStepSize(), simulator.episodeLength());
1253 simulator.setTimeStepSize(dt);
1254
1255 this->readSolutionFromOutputModule(restart_step, false);
1256
1257 this->eclWriter_->endRestart();
1258 }
1259
1261 {
1262 const auto& simulator = this->simulator();
1263
1264 // initial condition corresponds to hydrostatic conditions.
1265 EquilInitializer<TypeTag> equilInitializer(simulator, *(this->materialLawManager_));
1266
1267 std::size_t numElems = this->model().numGridDof();
1268 this->initialFluidStates_.resize(numElems);
1269 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1270 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1271 elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
1272 }
1273 }
1274
1276 {
1277 const auto& simulator = this->simulator();
1278 const auto& vanguard = simulator.vanguard();
1279 const auto& eclState = vanguard.eclState();
1280 const auto& fp = eclState.fieldProps();
1281 bool has_swat = fp.has_double("SWAT");
1282 bool has_sgas = fp.has_double("SGAS");
1283 bool has_rs = fp.has_double("RS");
1284 bool has_rsw = fp.has_double("RSW");
1285 bool has_rv = fp.has_double("RV");
1286 bool has_rvw = fp.has_double("RVW");
1287 bool has_pressure = fp.has_double("PRESSURE");
1288 bool has_salt = fp.has_double("SALT");
1289 bool has_saltp = fp.has_double("SALTP");
1290
1291 // make sure all required quantities are enables
1292 if (Indices::numPhases > 1) {
1293 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
1294 throw std::runtime_error("The ECL input file requires the presence of the SWAT keyword if "
1295 "the water phase is active");
1296 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
1297 throw std::runtime_error("The ECL input file requires the presence of the SGAS keyword if "
1298 "the gas phase is active");
1299 }
1300 if (!has_pressure)
1301 throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
1302 "keyword if the model is initialized explicitly");
1303 if (FluidSystem::enableDissolvedGas() && !has_rs)
1304 throw std::runtime_error("The ECL input file requires the RS keyword to be present if"
1305 " dissolved gas is enabled and the model is initialized explicitly");
1306 if (FluidSystem::enableDissolvedGasInWater() && !has_rsw)
1307 OpmLog::warning("The model is initialized explicitly and the RSW keyword is not present in the"
1308 " ECL input file. The RSW values are set equal to 0");
1309 if (FluidSystem::enableVaporizedOil() && !has_rv)
1310 throw std::runtime_error("The ECL input file requires the RV keyword to be present if"
1311 " vaporized oil is enabled and the model is initialized explicitly");
1312 if (FluidSystem::enableVaporizedWater() && !has_rvw)
1313 throw std::runtime_error("The ECL input file requires the RVW keyword to be present if"
1314 " vaporized water is enabled and the model is initialized explicitly");
1315 if (enableBrine && !has_salt)
1316 throw std::runtime_error("The ECL input file requires the SALT keyword to be present if"
1317 " brine is enabled and the model is initialized explicitly");
1318 if (enableSaltPrecipitation && !has_saltp)
1319 throw std::runtime_error("The ECL input file requires the SALTP keyword to be present if"
1320 " salt precipitation is enabled and the model is initialized explicitly");
1321
1322 std::size_t numDof = this->model().numGridDof();
1323
1324 initialFluidStates_.resize(numDof);
1325
1326 std::vector<double> waterSaturationData;
1327 std::vector<double> gasSaturationData;
1328 std::vector<double> pressureData;
1329 std::vector<double> rsData;
1330 std::vector<double> rswData;
1331 std::vector<double> rvData;
1332 std::vector<double> rvwData;
1333 std::vector<double> tempiData;
1334 std::vector<double> saltData;
1335 std::vector<double> saltpData;
1336
1337 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
1338 waterSaturationData = fp.get_double("SWAT");
1339 else
1340 waterSaturationData.resize(numDof);
1341
1342 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
1343 gasSaturationData = fp.get_double("SGAS");
1344 else
1345 gasSaturationData.resize(numDof);
1346
1347 pressureData = fp.get_double("PRESSURE");
1348 if (FluidSystem::enableDissolvedGas())
1349 rsData = fp.get_double("RS");
1350
1351 if (FluidSystem::enableDissolvedGasInWater() && has_rsw)
1352 rswData = fp.get_double("RSW");
1353
1354 if (FluidSystem::enableVaporizedOil())
1355 rvData = fp.get_double("RV");
1356
1357 if (FluidSystem::enableVaporizedWater())
1358 rvwData = fp.get_double("RVW");
1359
1360 // initial reservoir temperature
1361 tempiData = fp.get_double("TEMPI");
1362
1363 // initial salt concentration data
1364 if constexpr (enableBrine)
1365 saltData = fp.get_double("SALT");
1366
1367 // initial precipitated salt saturation data
1368 if constexpr (enableSaltPrecipitation)
1369 saltpData = fp.get_double("SALTP");
1370
1371 // calculate the initial fluid states
1372 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1373 auto& dofFluidState = initialFluidStates_[dofIdx];
1374
1375 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
1376
1378 // set temperature
1380 Scalar temperatureLoc = tempiData[dofIdx];
1381 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
1382 temperatureLoc = FluidSystem::surfaceTemperature;
1383 dofFluidState.setTemperature(temperatureLoc);
1384
1386 // set salt concentration
1388 if constexpr (enableBrine)
1389 dofFluidState.setSaltConcentration(saltData[dofIdx]);
1390
1392 // set precipitated salt saturation
1394 if constexpr (enableSaltPrecipitation)
1395 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
1396
1398 // set saturations
1400 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1401 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
1402 waterSaturationData[dofIdx]);
1403
1404 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
1405 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
1406 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1407 1.0
1408 - waterSaturationData[dofIdx]);
1409 }
1410 else
1411 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1412 gasSaturationData[dofIdx]);
1413 }
1414 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1415 const Scalar soil = 1.0 - waterSaturationData[dofIdx] - gasSaturationData[dofIdx];
1416 if (soil < smallSaturationTolerance_) {
1417 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
1418 }
1419 else {
1420 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx, soil);
1421 }
1422 }
1423
1425 // set phase pressures
1427 Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
1428
1429 // this assumes that capillary pressures only depend on the phase saturations
1430 // and possibly on temperature. (this is always the case for ECL problems.)
1431 std::array<Scalar, numPhases> pc = {0};
1432 const auto& matParams = this->materialLawParams(dofIdx);
1433 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
1434 Valgrind::CheckDefined(pressure);
1435 Valgrind::CheckDefined(pc);
1436 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1437 if (!FluidSystem::phaseIsActive(phaseIdx))
1438 continue;
1439
1440 if (Indices::oilEnabled)
1441 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1442 else if (Indices::gasEnabled)
1443 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1444 else if (Indices::waterEnabled)
1445 //single (water) phase
1446 dofFluidState.setPressure(phaseIdx, pressure);
1447 }
1448
1449 if constexpr (enableDissolvedGas) {
1450 if (FluidSystem::enableDissolvedGas())
1451 dofFluidState.setRs(rsData[dofIdx]);
1452 else if (Indices::gasEnabled && Indices::oilEnabled)
1453 dofFluidState.setRs(0.0);
1454 if (FluidSystem::enableVaporizedOil())
1455 dofFluidState.setRv(rvData[dofIdx]);
1456 else if (Indices::gasEnabled && Indices::oilEnabled)
1457 dofFluidState.setRv(0.0);
1458 }
1459
1460 if constexpr (enableDisgasInWater) {
1461 if (FluidSystem::enableDissolvedGasInWater() && has_rsw)
1462 dofFluidState.setRsw(rswData[dofIdx]);
1463 }
1464
1465 if constexpr (enableVapwat) {
1466 if (FluidSystem::enableVaporizedWater())
1467 dofFluidState.setRvw(rvwData[dofIdx]);
1468 }
1469
1471 // set invB_
1473 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1474 if (!FluidSystem::phaseIsActive(phaseIdx))
1475 continue;
1476
1477 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1478 dofFluidState.setInvB(phaseIdx, b);
1479
1480 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1481 dofFluidState.setDensity(phaseIdx, rho);
1482
1483 }
1484 }
1485 }
1486
1487
1488 void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
1489 {
1490 // each phase needs to be above certain value to be claimed to be existing
1491 // this is used to recover some RESTART running with the defaulted single-precision format
1492 Scalar sumSaturation = 0.0;
1493 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1494 if (FluidSystem::phaseIsActive(phaseIdx)) {
1495 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance_)
1496 elemFluidState.setSaturation(phaseIdx, 0.0);
1497
1498 sumSaturation += elemFluidState.saturation(phaseIdx);
1499 }
1500
1501 }
1502 if constexpr (enableSolvent) {
1503 if (solventSaturation < smallSaturationTolerance_)
1504 solventSaturation = 0.0;
1505
1506 sumSaturation += solventSaturation;
1507 }
1508
1509 assert(sumSaturation > 0.0);
1510
1511 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1512 if (FluidSystem::phaseIsActive(phaseIdx)) {
1513 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
1514 elemFluidState.setSaturation(phaseIdx, saturation);
1515 }
1516 }
1517 if constexpr (enableSolvent) {
1518 solventSaturation = solventSaturation / sumSaturation;
1519 }
1520 }
1521
1523 {
1524 FlowProblemType::readInitialCondition_();
1525
1526 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
1527 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
1528 enableSolvent,
1529 enablePolymer,
1530 enablePolymerMolarWeight,
1531 enableMICP);
1532
1533 }
1534
1535 void handleSolventBC(const BCProp::BCFace& bc, RateVector& rate) const override
1536 {
1537 if constexpr (!enableSolvent)
1538 throw std::logic_error("solvent is disabled and you're trying to add solvent to BC");
1539
1540 rate[Indices::solventSaturationIdx] = bc.rate;
1541 }
1542
1543 void handlePolymerBC(const BCProp::BCFace& bc, RateVector& rate) const override
1544 {
1545 if constexpr (!enablePolymer)
1546 throw std::logic_error("polymer is disabled and you're trying to add polymer to BC");
1547
1548 rate[Indices::polymerConcentrationIdx] = bc.rate;
1549 }
1550
1551 void handleMicrBC(const BCProp::BCFace& bc, RateVector& rate) const override
1552 {
1553 if constexpr (!enableMICP)
1554 throw std::logic_error("MICP is disabled and you're trying to add microbes to BC");
1555
1556 rate[Indices::microbialConcentrationIdx] = bc.rate;
1557 }
1558
1559 void handleOxygBC(const BCProp::BCFace& bc, RateVector& rate) const override
1560 {
1561 if constexpr (!enableMICP)
1562 throw std::logic_error("MICP is disabled and you're trying to add oxygen to BC");
1563
1564 rate[Indices::oxygenConcentrationIdx] = bc.rate;
1565 }
1566
1567 void handleUreaBC(const BCProp::BCFace& bc, RateVector& rate) const override
1568 {
1569 if constexpr (!enableMICP)
1570 throw std::logic_error("MICP is disabled and you're trying to add urea to BC");
1571
1572 rate[Indices::ureaConcentrationIdx] = bc.rate;
1573 // since the urea concentration can be much larger than 1, then we apply a scaling factor
1574 rate[Indices::ureaConcentrationIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
1575 }
1576
1577 void updateExplicitQuantities_(const bool first_step_after_restart)
1578 {
1579 OPM_TIMEBLOCK(updateExplicitQuantities);
1580 const bool invalidateFromMaxWaterSat = this->updateMaxWaterSaturation_();
1581 const bool invalidateFromMinPressure = this->updateMinPressure_();
1582
1583 // update hysteresis and max oil saturation used in vappars
1584 const bool invalidateFromHyst = this->updateHysteresis_();
1585 const bool invalidateFromMaxOilSat = this->updateMaxOilSaturation_();
1586
1587 // deal with DRSDT and DRVDT
1588 const bool invalidateDRDT = !first_step_after_restart && this->updateCompositionChangeLimits_();
1589
1590 // the derivatives may have changed
1591 const bool invalidateIntensiveQuantities
1592 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat || invalidateDRDT;
1593 if (invalidateIntensiveQuantities) {
1594 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1595 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
1596 }
1597
1598 this->updateRockCompTransMultVal_();
1599 }
1600
1602 {
1603 if (const auto nph = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
1604 + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
1605 + FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1606 nph < 2)
1607 {
1608 // Single phase runs don't need saturation functions and there's
1609 // nothing to do here. Return 'true' to tell caller that the
1610 // consistency requirements are Met.
1611 return true;
1612 }
1613
1614 const auto numSamplePoints = static_cast<std::size_t>
1615 (Parameters::Get<Parameters::NumSatfuncConsistencySamplePoints>());
1616
1617 auto sfuncConsistencyChecks =
1619 numSamplePoints, this->simulator().vanguard().eclState(),
1620 [&cmap = this->simulator().vanguard().cartesianIndexMapper()](const int elemIdx)
1621 { return cmap.cartesianIndex(elemIdx); }
1622 };
1623
1624 const auto ioRank = 0;
1625 const auto isIoRank = this->simulator().vanguard()
1626 .grid().comm().rank() == ioRank;
1627
1628 // Note: Run saturation function consistency checks on main grid
1629 // only (i.e., levelGridView(0)). These checks are not supported
1630 // for LGRs at this time.
1631 sfuncConsistencyChecks.collectFailuresTo(ioRank)
1632 .run(this->simulator().vanguard().grid().levelGridView(0),
1633 [&vg = this->simulator().vanguard(),
1634 &emap = this->simulator().model().elementMapper()]
1635 (const auto& elem)
1636 { return vg.gridIdxToEquilGridIdx(emap.index(elem)); });
1637
1638 using ViolationLevel = typename Satfunc::PhaseChecks::
1639 SatfuncConsistencyCheckManager<Scalar>::ViolationLevel;
1640
1641 auto reportFailures = [&sfuncConsistencyChecks]
1642 (const ViolationLevel level)
1643 {
1644 sfuncConsistencyChecks.reportFailures
1645 (level, [](std::string_view record)
1646 { OpmLog::info(std::string { record }); });
1647 };
1648
1649 if (sfuncConsistencyChecks.anyFailedStandardChecks()) {
1650 if (isIoRank) {
1651 OpmLog::warning("Saturation Function "
1652 "End-point Consistency Problems");
1653
1654 reportFailures(ViolationLevel::Standard);
1655 }
1656 }
1657
1658 if (sfuncConsistencyChecks.anyFailedCriticalChecks()) {
1659 if (isIoRank) {
1660 OpmLog::error("Saturation Function "
1661 "End-point Consistency Failures");
1662
1663 reportFailures(ViolationLevel::Critical);
1664 }
1665
1666 // There are "critical" check failures. Report that consistency
1667 // requirements are not Met.
1668 return false;
1669 }
1670
1671 // If we get here then there are no critical failures. Report
1672 // Met = true, i.e., that the consistency requirements ARE met.
1673 return true;
1674 }
1675
1677
1678 std::vector<InitialFluidState> initialFluidStates_;
1679
1681 std::unique_ptr<EclWriterType> eclWriter_;
1682
1683 const Scalar smallSaturationTolerance_ = 1.e-6;
1684#if HAVE_DAMARIS
1685 bool enableDamarisOutput_ = false ;
1686 std::unique_ptr<DamarisWriterType> damarisWriter_;
1687#endif
1689
1691
1692 ModuleParams moduleParams_;
1693
1694private:
1705 bool episodeWillBeOver() const override
1706 {
1707 const auto currTime = this->simulator().time()
1708 + this->simulator().timeStepSize();
1709
1710 const auto nextReportStep =
1711 this->simulator().vanguard().schedule()
1712 .seconds(this->simulator().episodeIndex() + 1);
1713
1714 const auto isSubStep = (nextReportStep - currTime)
1715 > (2 * std::numeric_limits<float>::epsilon()) * nextReportStep;
1716
1717 return !isSubStep;
1718 }
1719};
1720
1721} // namespace Opm
1722
1723#endif // OPM_FLOW_PROBLEM_BLACK_HPP
Class handling Action support in simulator.
Definition: ActionHandler.hpp:52
static void setParams(BlackOilBrineParams< Scalar > &&params)
Set parameters.
Definition: blackoilbrinemodules.hh:88
static void setParams(BlackOilExtboParams< Scalar > &&params)
Set parameters.
Definition: blackoilextbomodules.hh:93
static void setParams(BlackOilFoamParams< Scalar > &&params)
Set parameters.
Definition: blackoilfoammodules.hh:90
static void setParams(BlackOilMICPParams< Scalar > &&params)
Set parameters.
Definition: blackoilmicpmodules.hh:93
static void setParams(BlackOilPolymerParams< Scalar > &&params)
Set parameters.
Definition: blackoilpolymermodules.hh:97
static void setParams(BlackOilSolventParams< Scalar > &&params)
Set parameters.
Definition: blackoilsolventmodules.hh:101
Collects necessary output values and pass it to opm-common's ECL output.
Definition: EclWriter.hpp:114
static void registerParameters()
Definition: EclWriter.hpp:137
Computes the initial condition based on the EQUIL keyword from ECL.
Definition: EquilInitializer.hpp:58
BlackOilFluidState< Scalar, FluidSystem, enableTemperature, enableEnergy, enableDissolution, enableVapwat, enableBrine, enableSaltPrecipitation, enableDisgasInWater, Indices::numPhases > ScalarFluidState
Definition: EquilInitializer.hpp:98
const ScalarFluidState & initialFluidState(unsigned elemIdx) const
Return the initial thermodynamic state which should be used as the initial condition.
Definition: EquilInitializer.hpp:198
void readRockParameters_(const std::vector< Scalar > &cellCenterDepths, std::function< std::array< int, 3 >(const unsigned)> ijkIndex)
Definition: FlowGenericProblem_impl.hpp:136
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowProblemBlackoil.hpp:82
void updateExplicitQuantities_(int episodeIdx, int timeStepSize, const bool first_step_after_restart) override
Definition: FlowProblemBlackoil.hpp:1131
bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities &iq)
Definition: FlowProblemBlackoil.hpp:1151
void handlePolymerBC(const BCProp::BCFace &bc, RateVector &rate) const override
Definition: FlowProblemBlackoil.hpp:1543
void writeOutput(const bool verbose) override
Write the requested quantities of the current solution into the output files.
Definition: FlowProblemBlackoil.hpp:538
void readInitialCondition_() override
Definition: FlowProblemBlackoil.hpp:1522
void handleOxygBC(const BCProp::BCFace &bc, RateVector &rate) const override
Definition: FlowProblemBlackoil.hpp:1559
void readEquilInitialCondition_() override
Definition: FlowProblemBlackoil.hpp:1260
void handleSolventBC(const BCProp::BCFace &bc, RateVector &rate) const override
Definition: FlowProblemBlackoil.hpp:1535
Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
Returns the maximum value of the gas dissolution factor at the current time for a given degree of fre...
Definition: FlowProblemBlackoil.hpp:863
const std::vector< InitialFluidState > & initialFluidStates() const
Definition: FlowProblemBlackoil.hpp:728
void processRestartSaturations_(InitialFluidState &elemFluidState, Scalar &solventSaturation)
Definition: FlowProblemBlackoil.hpp:1488
std::vector< InitialFluidState > & initialFluidStates()
Definition: FlowProblemBlackoil.hpp:725
FlowProblemBlackoil(Simulator &simulator)
Definition: FlowProblemBlackoil.hpp:180
bool enableEclOutput_
Definition: FlowProblemBlackoil.hpp:1680
Scalar drsdtcon(unsigned elemIdx, int episodeIdx) const
Definition: FlowProblemBlackoil.hpp:943
void endStepApplyAction()
Definition: FlowProblemBlackoil.hpp:459
bool drsdtconIsActive(unsigned elemIdx, int episodeIdx) const
Definition: FlowProblemBlackoil.hpp:949
Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
Returns the maximum value of the oil vaporization factor at the current time for a given degree of fr...
Definition: FlowProblemBlackoil.hpp:874
std::vector< InitialFluidState > initialFluidStates_
Definition: FlowProblemBlackoil.hpp:1678
void endTimeStep() override
Called by the simulator after each time integration.
Definition: FlowProblemBlackoil.hpp:453
void updateMaxPolymerAdsorption_()
Definition: FlowProblemBlackoil.hpp:1141
const InitialFluidState & initialFluidState(unsigned globalDofIdx) const
Definition: FlowProblemBlackoil.hpp:722
void endEpisode() override
Called by the simulator after the end of an episode.
Definition: FlowProblemBlackoil.hpp:506
void setSubStepReport(const SimulatorReportSingle &report)
Definition: FlowProblemBlackoil.hpp:734
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: FlowProblemBlackoil.hpp:905
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: FlowProblemBlackoil.hpp:269
void handleMicrBC(const BCProp::BCFace &bc, RateVector &rate) const override
Definition: FlowProblemBlackoil.hpp:1551
const FlowThresholdPressure< TypeTag > & thresholdPressure() const
Definition: FlowProblemBlackoil.hpp:1111
void finalizeOutput()
Definition: FlowProblemBlackoil.hpp:559
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: FlowProblemBlackoil.hpp:960
InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const
Definition: FlowProblemBlackoil.hpp:740
std::unique_ptr< EclWriterType > eclWriter_
Definition: FlowProblemBlackoil.hpp:1681
void initialSolutionApplied() override
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition: FlowProblemBlackoil.hpp:571
const ModuleParams & moduleParams() const
Definition: FlowProblemBlackoil.hpp:1117
void readEclRestartSolution_()
Definition: FlowProblemBlackoil.hpp:1230
FlowThresholdPressure< TypeTag > & thresholdPressure()
Definition: FlowProblemBlackoil.hpp:1114
FlowThresholdPressure< TypeTag > thresholdPressures_
Definition: FlowProblemBlackoil.hpp:1676
void readExplicitInitialCondition_() override
Definition: FlowProblemBlackoil.hpp:1275
void beginEpisode() override
Called by the simulator before an episode begins.
Definition: FlowProblemBlackoil.hpp:238
bool recycleFirstIterationStorage() const
Return if the storage term of the first iteration is identical to the storage term for the solution o...
Definition: FlowProblemBlackoil.hpp:889
LhsEval permFactTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to porosity reduction.
Definition: FlowProblemBlackoil.hpp:702
void serializeOp(Serializer &serializer)
Definition: FlowProblemBlackoil.hpp:1123
MixingRateControls< FluidSystem > mixControls_
Definition: FlowProblemBlackoil.hpp:1688
void writeReports(const SimulatorTimer &timer)
Definition: FlowProblemBlackoil.hpp:526
ModuleParams moduleParams_
Definition: FlowProblemBlackoil.hpp:1692
const EclWriterType & eclWriter() const
Definition: FlowProblemBlackoil.hpp:853
void setSimulationReport(const SimulatorReport &report)
Definition: FlowProblemBlackoil.hpp:737
void addToSourceDense(RateVector &rate, unsigned globalDofIdx, unsigned timeIdx) const override
Definition: FlowProblemBlackoil.hpp:593
void computeAndSetEqWeights_()
Definition: FlowProblemBlackoil.hpp:1163
void updateExplicitQuantities_(const bool first_step_after_restart)
Definition: FlowProblemBlackoil.hpp:1577
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition: FlowProblemBlackoil.hpp:166
ActionHandler< Scalar, IndexTraits > actionHandler_
Definition: FlowProblemBlackoil.hpp:1690
void readSolutionFromOutputModule(const int restart_step, bool fip_init)
Read simulator solution state from the outputmodule (used with restart)
Definition: FlowProblemBlackoil.hpp:1001
const EclipseIO & eclIO() const
Definition: FlowProblemBlackoil.hpp:731
Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
Definition: FlowProblemBlackoil.hpp:1108
void handleUreaBC(const BCProp::BCFace &bc, RateVector &rate) const override
Definition: FlowProblemBlackoil.hpp:1567
EclWriterType & eclWriter()
Definition: FlowProblemBlackoil.hpp:856
bool satfuncConsistencyRequirementsMet() const
Definition: FlowProblemBlackoil.hpp:1601
bool updateCompositionChangeLimits_()
Definition: FlowProblemBlackoil.hpp:1197
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowProblem.hpp:94
@ dim
Definition: FlowProblem.hpp:110
virtual void writeOutput(bool verbose)
Write the requested quantities of the current solution into the output files.
Definition: FlowProblem.hpp:479
@ gasCompIdx
Definition: FlowProblem.hpp:140
@ enableMICP
Definition: FlowProblem.hpp:126
@ enableExperiments
Definition: FlowProblem.hpp:123
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition: FlowProblem.hpp:826
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:659
GetPropType< TypeTag, Properties::Vanguard > Vanguard
Definition: FlowProblem.hpp:107
@ enableEnergy
Definition: FlowProblem.hpp:122
@ oilPhaseIdx
Definition: FlowProblem.hpp:135
@ numPhases
Definition: FlowProblem.hpp:115
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: FlowProblem.hpp:101
GetPropType< TypeTag, Properties::EqVector > EqVector
Definition: FlowProblem.hpp:106
@ enableConvectiveMixing
Definition: FlowProblem.hpp:118
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: FlowProblem.hpp:148
GlobalEqVector drift_
Definition: FlowProblem.hpp:1680
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: FlowProblem.hpp:145
@ enablePolymerMolarWeight
Definition: FlowProblem.hpp:128
Dune::FieldMatrix< Scalar, dimWorld, dimWorld > DimMatrix
Definition: FlowProblem.hpp:163
int episodeIndex() const
Definition: FlowProblem.hpp:286
GetPropType< TypeTag, Properties::Indices > Indices
Definition: FlowProblem.hpp:157
@ enablePolymer
Definition: FlowProblem.hpp:127
@ enableSaltPrecipitation
Definition: FlowProblem.hpp:129
@ gasPhaseIdx
Definition: FlowProblem.hpp:134
@ enableBrine
Definition: FlowProblem.hpp:119
GetPropType< TypeTag, Properties::GlobalEqVector > GlobalEqVector
Definition: FlowProblem.hpp:105
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: FlowProblem.hpp:146
@ enableFoam
Definition: FlowProblem.hpp:125
@ enableExtbo
Definition: FlowProblem.hpp:124
@ dimWorld
Definition: FlowProblem.hpp:111
TracerModel tracerModel_
Definition: FlowProblem.hpp:1686
@ waterPhaseIdx
Definition: FlowProblem.hpp:136
WellModel wellModel_
Definition: FlowProblem.hpp:1682
virtual void beginEpisode()
Called by the simulator before an episode begins.
Definition: FlowProblem.hpp:294
void readThermalParameters_()
Definition: FlowProblem.hpp:1375
@ enableDispersion
Definition: FlowProblem.hpp:121
@ numEq
Definition: FlowProblem.hpp:114
@ enableTemperature
Definition: FlowProblem.hpp:131
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: FlowProblem.hpp:158
@ oilCompIdx
Definition: FlowProblem.hpp:141
GetPropType< TypeTag, Properties::GridView > GridView
Definition: FlowProblem.hpp:102
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition: FlowProblem.hpp:180
void updatePffDofData_()
Definition: FlowProblem.hpp:1488
GetPropType< TypeTag, Properties::PrimaryVariables > PrimaryVariables
Definition: FlowProblem.hpp:144
@ enableDiffusion
Definition: FlowProblem.hpp:120
@ enableThermalFluxBoundaries
Definition: FlowProblem.hpp:132
void readBoundaryConditions_()
Definition: FlowProblem.hpp:1519
Vanguard::TransmissibilityType transmissibilities_
Definition: FlowProblem.hpp:1675
@ enableSolvent
Definition: FlowProblem.hpp:130
@ waterCompIdx
Definition: FlowProblem.hpp:142
@ numComponents
Definition: FlowProblem.hpp:116
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: FlowProblem.hpp:104
GetPropType< TypeTag, Properties::MaterialLaw > MaterialLaw
Definition: FlowProblem.hpp:154
void readMaterialParameters_()
Definition: FlowProblem.hpp:1341
This class calculates the threshold pressure for grid faces according to the Eclipse Reference Manual...
Definition: FlowThresholdPressure.hpp:59
Class handling mixing rate controls for a FlowProblemBlackoil.
Definition: MixingRateControls.hpp:46
Definition: SatfuncConsistencyCheckManager.hpp:58
SatfuncConsistencyCheckManager & collectFailuresTo(const int root)
Definition: SatfuncConsistencyCheckManager.hpp:99
void run(const GridView &gv, GetCellIndex &&getCellIndex)
Definition: SatfuncConsistencyCheckManager.hpp:128
Definition: SimulatorTimer.hpp:39
VTK output module for the tracer model's parameters.
Definition: VtkTracerModule.hpp:58
static void registerParameters()
Register all run-time parameters for the tracer VTK output module.
Definition: VtkTracerModule.hpp:84
@ NONE
Definition: DeferredLogger.hpp:46
static constexpr int dim
Definition: structuredgridvanguard.hh:68
Definition: blackoilboundaryratevector.hh:39
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:233
Struct holding the parameters for the BlackoilBrineModule class.
Definition: blackoilbrineparams.hpp:44
Struct holding the parameters for the BlackoilExtboModule class.
Definition: blackoilextboparams.hpp:49
Struct holding the parameters for the BlackoilFoamModule class.
Definition: blackoilfoamparams.hpp:46
Struct holding the parameters for the BlackOilMICPModule class.
Definition: blackoilmicpparams.hpp:44
Struct holding the parameters for the BlackOilPolymerModule class.
Definition: blackoilpolymerparams.hpp:45
Struct holding the parameters for the BlackOilSolventModule class.
Definition: blackoilsolventparams.hpp:49
Definition: SimulatorReport.hpp:122
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:34