My Project
Loading...
Searching...
No Matches
FlowProblem.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
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 2 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20
21 Consult the COPYING file in the top-level source directory of this
22 module for the precise wording of the license and the list of
23 copyright holders.
24*/
30#ifndef OPM_FLOW_PROBLEM_HPP
31#define OPM_FLOW_PROBLEM_HPP
32
33#include <dune/common/version.hh>
34#include <dune/common/fvector.hh>
35#include <dune/common/fmatrix.hh>
36
37#include <opm/common/utility/TimeService.hpp>
38
39#include <opm/core/props/satfunc/RelpermDiagnostics.hpp>
40
41#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
42#include <opm/input/eclipse/Parser/ParserKeywords/E.hpp>
43#include <opm/input/eclipse/Schedule/Schedule.hpp>
44
45#include <opm/material/common/ConditionalStorage.hpp>
46#include <opm/material/common/Valgrind.hpp>
47#include <opm/material/densead/Evaluation.hpp>
48#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
49#include <opm/material/fluidstates/CompositionalFluidState.hpp>
50#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
51#include <opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp>
52#include <opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp>
53#include <opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp>
54#include <opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp>
55#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp>
56#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp>
57#include <opm/material/thermal/EclThermalLawManager.hpp>
58
59#include <opm/models/common/directionalmobility.hh>
60#include <opm/models/utils/pffgridvector.hh>
61#include <opm/models/blackoil/blackoilmodel.hh>
62#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
63
64#include <opm/output/eclipse/EclipseIO.hpp>
65
66#include <opm/simulators/flow/ActionHandler.hpp>
81#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
82#include <opm/simulators/timestepping/SimulatorReport.hpp>
83#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
84#include <opm/simulators/utils/ParallelSerialization.hpp>
85
86#include <opm/utility/CopyablePtr.hpp>
87
88#include <opm/common/OpmLog/OpmLog.hpp>
89
90#if HAVE_DAMARIS
92#endif
93
94#include <algorithm>
95#include <functional>
96#include <set>
97#include <string>
98#include <vector>
99
100namespace Opm {
101
108template <class TypeTag>
109class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
110 , public FlowGenericProblem<GetPropType<TypeTag, Properties::GridView>,
111 GetPropType<TypeTag, Properties::FluidSystem>>
112{
114 GetPropType<TypeTag, Properties::FluidSystem>>;
115 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
116 using Implementation = GetPropType<TypeTag, Properties::Problem>;
117
118 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
119 using GridView = GetPropType<TypeTag, Properties::GridView>;
120 using Stencil = GetPropType<TypeTag, Properties::Stencil>;
121 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
122 using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
123 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
124 using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
125
126 // Grid and world dimension
127 enum { dim = GridView::dimension };
128 enum { dimWorld = GridView::dimensionworld };
129
130 // copy some indices for convenience
131 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
132 enum { numPhases = FluidSystem::numPhases };
133 enum { numComponents = FluidSystem::numComponents };
134 enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
135 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
136 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
137 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
138 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
139 enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
140 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
141 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
142 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
143 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
144 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
145 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
146 enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
147 enum { enableApiTracking = getPropValue<TypeTag, Properties::EnableApiTracking>() };
148 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
149 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
150 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
151 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
152 enum { gasCompIdx = FluidSystem::gasCompIdx };
153 enum { oilCompIdx = FluidSystem::oilCompIdx };
154 enum { waterCompIdx = FluidSystem::waterCompIdx };
155
156 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
157 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
158 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
159 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
160 using Element = typename GridView::template Codim<0>::Entity;
161 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
162 using EclMaterialLawManager = typename GetProp<TypeTag, Properties::MaterialLaw>::EclMaterialLawManager;
163 using EclThermalLawManager = typename GetProp<TypeTag, Properties::SolidEnergyLaw>::EclThermalLawManager;
164 using MaterialLawParams = typename EclMaterialLawManager::MaterialLawParams;
165 using SolidEnergyLawParams = typename EclThermalLawManager::SolidEnergyLawParams;
166 using ThermalConductionLawParams = typename EclThermalLawManager::ThermalConductionLawParams;
167 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
168 using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
169 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
170 using Indices = GetPropType<TypeTag, Properties::Indices>;
171 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
172 using WellModel = GetPropType<TypeTag, Properties::WellModel>;
173 using AquiferModel = GetPropType<TypeTag, Properties::AquiferModel>;
174
175 using SolventModule = BlackOilSolventModule<TypeTag>;
176 using PolymerModule = BlackOilPolymerModule<TypeTag>;
177 using FoamModule = BlackOilFoamModule<TypeTag>;
178 using BrineModule = BlackOilBrineModule<TypeTag>;
179 using ExtboModule = BlackOilExtboModule<TypeTag>;
180 using MICPModule = BlackOilMICPModule<TypeTag>;
181 using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
182 using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
183
184 using InitialFluidState = typename EquilInitializer<TypeTag>::ScalarFluidState;
185
186 using Toolbox = MathToolbox<Evaluation>;
187 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
188
189 using EclWriterType = EclWriter<TypeTag>;
190#if HAVE_DAMARIS
191 using DamarisWriterType = DamarisWriter<TypeTag>;
192#endif
193
195 using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
196
197public:
204 using BaseType::porosity;
205
209 static void registerParameters()
210 {
211 ParentType::registerParameters();
212 EclWriterType::registerParameters();
213#if HAVE_DAMARIS
214 DamarisWriterType::registerParameters();
215#endif
216
218
219 Parameters::registerParam<TypeTag, Properties::EnableWriteAllSolutions>
220 ("Write all solutions to disk instead of only the ones for the "
221 "report steps");
222 Parameters::registerParam<TypeTag, Properties::EnableEclOutput>
223 ("Write binary output which is compatible with the commercial "
224 "Eclipse simulator");
225#if HAVE_DAMARIS
226 Parameters::registerParam<TypeTag, Properties::EnableDamarisOutput>
227 ("Write a specific variable using Damaris in a separate core");
228#endif
229 Parameters::registerParam<TypeTag, Properties::EclOutputDoublePrecision>
230 ("Tell the output writer to use double precision. Useful for 'perfect' restarts");
231 Parameters::registerParam<TypeTag, Properties::RestartWritingInterval>
232 ("The frequencies of which time steps are serialized to disk");
233 Parameters::registerParam<TypeTag, Properties::EnableDriftCompensation>
234 ("Enable partial compensation of systematic mass losses via "
235 "the source term of the next time step");
236 Parameters::registerParam<TypeTag, Properties::OutputMode>
237 ("Specify which messages are going to be printed. "
238 "Valid values are: none, log, all (default)");
239 Parameters::registerParam<TypeTag, Properties::NumPressurePointsEquil>
240 ("Number of pressure points (in each direction) in tables used for equilibration");
241 Parameters::hideParam<TypeTag, Properties::NumPressurePointsEquil>(); // Users will typically not need to modify this parameter..
242 Parameters::registerParam<TypeTag, Properties::ExplicitRockCompaction>
243 ("Use pressure from end of the last time step when evaluating rock compaction");
244 Parameters::hideParam<TypeTag, Properties::ExplicitRockCompaction>(); // Users will typically not need to modify this parameter..
245 }
246
247
251 static int handlePositionalParameter(std::set<std::string>& seenParams,
252 std::string& errorMsg,
253 int,
254 const char** argv,
255 int paramIdx,
256 int)
257 {
258 using ParamsMeta = GetProp<TypeTag, Properties::ParameterMetaData>;
259 Dune::ParameterTree& tree = ParamsMeta::tree();
260 return eclPositionalParameter(tree,
261 seenParams,
262 errorMsg,
263 argv,
264 paramIdx);
265 }
266
270 FlowProblem(Simulator& simulator)
271 : ParentType(simulator)
272 , BaseType(simulator.vanguard().eclState(),
273 simulator.vanguard().schedule(),
274 simulator.vanguard().gridView())
275 , transmissibilities_(simulator.vanguard().eclState(),
276 simulator.vanguard().gridView(),
277 simulator.vanguard().cartesianIndexMapper(),
278 simulator.vanguard().grid(),
279 simulator.vanguard().cellCentroids(),
280 enableEnergy,
281 enableDiffusion,
282 enableDispersion)
283 , thresholdPressures_(simulator)
284 , wellModel_(simulator)
285 , aquiferModel_(simulator)
286 , pffDofData_(simulator.gridView(), this->elementMapper())
287 , tracerModel_(simulator)
288 , actionHandler_(simulator.vanguard().eclState(),
289 simulator.vanguard().schedule(),
290 simulator.vanguard().actionState(),
291 simulator.vanguard().summaryState(),
292 wellModel_,
293 simulator.vanguard().grid().comm())
294 {
295 this->model().addOutputModule(new VtkTracerModule<TypeTag>(simulator));
296 // Tell the black-oil extensions to initialize their internal data structures
297 const auto& vanguard = simulator.vanguard();
298 SolventModule::initFromState(vanguard.eclState(), vanguard.schedule());
299 PolymerModule::initFromState(vanguard.eclState());
300 FoamModule::initFromState(vanguard.eclState());
301 BrineModule::initFromState(vanguard.eclState());
302 ExtboModule::initFromState(vanguard.eclState());
303 MICPModule::initFromState(vanguard.eclState());
304 DispersionModule::initFromState(vanguard.eclState());
305 DiffusionModule::initFromState(vanguard.eclState());
306
307 // create the ECL writer
308 eclWriter_ = std::make_unique<EclWriterType>(simulator);
309#if HAVE_DAMARIS
310 // create Damaris writer
311 damarisWriter_ = std::make_unique<DamarisWriterType>(simulator);
312 enableDamarisOutput_ = Parameters::get<TypeTag, Properties::EnableDamarisOutput>();
313#endif
314 enableDriftCompensation_ = Parameters::get<TypeTag, Properties::EnableDriftCompensation>();
315
316 enableEclOutput_ = Parameters::get<TypeTag, Properties::EnableEclOutput>();
317
318 this->enableTuning_ = Parameters::get<TypeTag, Properties::EnableTuning>();
319 this->initialTimeStepSize_ = Parameters::get<TypeTag, Properties::InitialTimeStepSize>();
320 this->maxTimeStepAfterWellEvent_ = Parameters::get<TypeTag, Properties::TimeStepAfterEventInDays>() * 24 * 60 * 60;
321
322 // The value N for this parameter is defined in the following order of presedence:
323 // 1. Command line value (--num-pressure-points-equil=N)
324 // 2. EQLDIMS item 2
325 // Default value is defined in opm-common/src/opm/input/eclipse/share/keywords/000_Eclipse100/E/EQLDIMS
326 if (Parameters::isSet<TypeTag, Properties::NumPressurePointsEquil>())
327 {
328 this->numPressurePointsEquil_ = Parameters::get<TypeTag, Properties::NumPressurePointsEquil>();
329 } else {
330 this->numPressurePointsEquil_ = simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP();
331 }
332
333 RelpermDiagnostics relpermDiagnostics;
334 relpermDiagnostics.diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
335 }
336
341 {
342 ParentType::finishInit();
343
344 auto& simulator = this->simulator();
345 const auto& eclState = simulator.vanguard().eclState();
346 const auto& schedule = simulator.vanguard().schedule();
347
348 // Set the start time of the simulation
349 simulator.setStartTime(schedule.getStartTime());
350 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
351
352 // We want the episode index to be the same as the report step index to make
353 // things simpler, so we have to set the episode index to -1 because it is
354 // incremented by endEpisode(). The size of the initial time step and
355 // length of the initial episode is set to zero for the same reason.
356 simulator.setEpisodeIndex(-1);
357 simulator.setEpisodeLength(0.0);
358
359 // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
360 // disables gravity, else the standard value of the gravity constant at sea level
361 // on earth is used
362 this->gravity_ = 0.0;
363 if (Parameters::get<TypeTag, Properties::EnableGravity>())
364 this->gravity_[dim - 1] = 9.80665;
365 if (!eclState.getInitConfig().hasGravity())
366 this->gravity_[dim - 1] = 0.0;
367
368 if (this->enableTuning_) {
369 // if support for the TUNING keyword is enabled, we get the initial time
370 // steping parameters from it instead of from command line parameters
371 const auto& tuning = schedule[0].tuning();
372 this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
373 this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
374 }
375
376 this->initFluidSystem_();
377
378 // deal with DRSDT
379 this->mixControls_.init(this->model().numGridDof(),
380 this->episodeIndex(),
381 eclState.runspec().tabdims().getNumPVTTables());
382
383 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
384 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
385 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
386 }
387
388 this->readRockParameters_(simulator.vanguard().cellCenterDepths(),
389 [&simulator](const unsigned idx)
390 {
391 std::array<int,dim> coords;
392 simulator.vanguard().cartesianCoordinate(idx, coords);
393 for (auto& c : coords) {
394 ++c;
395 }
396 return coords;
397 });
398 readMaterialParameters_();
399 readThermalParameters_();
400
401 // Re-ordering in case of ALUGrid
402 std::function<unsigned int(unsigned int)> gridToEquilGrid = [&simulator](unsigned int i) {
403 return simulator.vanguard().gridIdxToEquilGridIdx(i);
404 };
405 transmissibilities_.finishInit(gridToEquilGrid);
406
407 const auto& initconfig = eclState.getInitConfig();
408 tracerModel_.init(initconfig.restartRequested());
409 if (initconfig.restartRequested())
410 readEclRestartSolution_();
411 else
412 readInitialCondition_();
413
414 tracerModel_.prepareTracerBatches();
415
416 updatePffDofData_();
417
418 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
419 const auto& vanguard = this->simulator().vanguard();
420 const auto& gridView = vanguard.gridView();
421 int numElements = gridView.size(/*codim=*/0);
422 this->polymer_.maxAdsorption.resize(numElements, 0.0);
423 }
424
425 readBoundaryConditions_();
426
427 // compute and set eq weights based on initial b values
428 computeAndSetEqWeights_();
429
430 if (enableDriftCompensation_) {
431 drift_.resize(this->model().numGridDof());
432 drift_ = 0.0;
433 }
434
435 // write the static output files (EGRID, INIT, SMSPEC, etc.)
436 if (enableEclOutput_) {
437 if (simulator.vanguard().grid().comm().size() > 1) {
438 if (simulator.vanguard().grid().comm().rank() == 0)
439 eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
440 } else
441 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
442
443 // Re-ordering in case of ALUGrid
444 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
445 return simulator.vanguard().gridEquilIdxToGridIdx(i);
446 };
447 eclWriter_->writeInit(equilGridToGrid);
448 }
449
450 simulator.vanguard().releaseGlobalTransmissibilities();
451
452 // after finishing the initialization and writing the initial solution, we move
453 // to the first "real" episode/report step
454 // for restart the episode index and start is already set
455 if (!initconfig.restartRequested()) {
456 simulator.startNextEpisode(schedule.seconds(1));
457 simulator.setEpisodeIndex(0);
458 }
459 }
460
461 void prefetch(const Element& elem) const
462 { pffDofData_.prefetch(elem); }
463
475 template <class Restarter>
476 void deserialize(Restarter& res)
477 {
478 // reload the current episode/report step from the deck
479 beginEpisode();
480
481 // deserialize the wells
482 wellModel_.deserialize(res);
483
484 // deserialize the aquifer
485 aquiferModel_.deserialize(res);
486 }
487
494 template <class Restarter>
495 void serialize(Restarter& res)
496 {
497 wellModel_.serialize(res);
498
499 aquiferModel_.serialize(res);
500 }
501
502 int episodeIndex() const
503 {
504 return std::max(this->simulator().episodeIndex(), 0);
505 }
506
511 {
512 OPM_TIMEBLOCK(beginEpisode);
513 // Proceed to the next report step
514 auto& simulator = this->simulator();
515 int episodeIdx = simulator.episodeIndex();
516 auto& eclState = simulator.vanguard().eclState();
517 const auto& schedule = simulator.vanguard().schedule();
518 const auto& events = schedule[episodeIdx].events();
519
520 if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
521 // bring the contents of the keywords to the current state of the SCHEDULE
522 // section.
523 //
524 // TODO (?): make grid topology changes possible (depending on what exactly
525 // has changed, the grid may need be re-created which has some serious
526 // implications on e.g., the solution of the simulation.)
527 const auto& miniDeck = schedule[episodeIdx].geo_keywords();
528 const auto& cc = simulator.vanguard().grid().comm();
529 eclState.apply_schedule_keywords( miniDeck );
530 eclBroadcast(cc, eclState.getTransMult() );
531
532 // Re-ordering in case of ALUGrid
533 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
534 return simulator.vanguard().gridEquilIdxToGridIdx(i);
535 };
536
537 // re-compute all quantities which may possibly be affected.
538 transmissibilities_.update(true, equilGridToGrid);
539 this->referencePorosity_[1] = this->referencePorosity_[0];
540 updateReferencePorosity_();
541 updatePffDofData_();
542 this->model().linearizer().updateDiscretizationParameters();
543 }
544
545 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
546
547 // set up the wells for the next episode.
548 wellModel_.beginEpisode();
549
550 // set up the aquifers for the next episode.
551 aquiferModel_.beginEpisode();
552
553 // set the size of the initial time step of the episode
554 Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
555 // negative value of initialTimeStepSize_ indicates no active limit from TSINIT or NEXTSTEP
556 if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
557 // allow the size of the initial time step to be set via an external parameter
558 // if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
559 dt = std::min(dt, this->initialTimeStepSize_);
560 simulator.setTimeStepSize(dt);
561
562 // Evaluate UDQ assign statements to make sure the settings are
563 // available as UDA controls for the current report step.
564 actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
565
566 if (episodeIdx >= 0) {
567 const auto& oilVap = schedule[episodeIdx].oilvap();
568 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
569 FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2());
570 } else {
571 FluidSystem::setVapPars(0.0, 0.0);
572 }
573 }
574 }
575
580 {
581 OPM_TIMEBLOCK(beginTimeStep);
582 int episodeIdx = this->episodeIndex();
583
584 this->beginTimeStep_(enableExperiments,
585 episodeIdx,
586 this->simulator().timeStepIndex(),
587 this->simulator().startTime(),
588 this->simulator().time(),
589 this->simulator().timeStepSize(),
590 this->simulator().endTime());
591
592 // update maximum water saturation and minimum pressure
593 // used when ROCKCOMP is activated
594 asImp_().updateExplicitQuantities_();
595
596 if (nonTrivialBoundaryConditions()) {
597 this->model().linearizer().updateBoundaryConditionData();
598 }
599
600 wellModel_.beginTimeStep();
601 aquiferModel_.beginTimeStep();
602 tracerModel_.beginTimeStep();
603
604 }
605
610 {
611 OPM_TIMEBLOCK(beginIteration);
612 wellModel_.beginIteration();
613 aquiferModel_.beginIteration();
614 }
615
620 {
621 OPM_TIMEBLOCK(endIteration);
622 wellModel_.endIteration();
623 aquiferModel_.endIteration();
624 }
625
630 {
631 OPM_TIMEBLOCK(endTimeStep);
632#ifndef NDEBUG
633 if constexpr (getPropValue<TypeTag, Properties::EnableDebuggingChecks>()) {
634 // in debug mode, we don't care about performance, so we check if the model does
635 // the right thing (i.e., the mass change inside the whole reservoir must be
636 // equivalent to the fluxes over the grid's boundaries plus the source rates
637 // specified by the problem)
638 int rank = this->simulator().gridView().comm().rank();
639 if (rank == 0)
640 std::cout << "checking conservativeness of solution\n";
641 this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
642 if (rank == 0)
643 std::cout << "solution is sufficiently conservative\n";
644 }
645#endif // NDEBUG
646
647 auto& simulator = this->simulator();
648 wellModel_.endTimeStep();
649 aquiferModel_.endTimeStep();
650 tracerModel_.endTimeStep();
651
652
653 // Compute flux for output
654 this->model().linearizer().updateFlowsInfo();
655
656 // deal with DRSDT and DRVDT
657 asImp_().updateCompositionChangeLimits_();
658 {
659 OPM_TIMEBLOCK(driftCompansation);
660 if (enableDriftCompensation_) {
661 const auto& residual = this->model().linearizer().residual();
662 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
663 drift_[globalDofIdx] = residual[globalDofIdx];
664 drift_[globalDofIdx] *= simulator.timeStepSize();
665 if constexpr (getPropValue<TypeTag, Properties::UseVolumetricResidual>())
666 drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
667 }
668 }
669 }
670 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
671 !this->simulator().episodeWillBeOver();
672 // For CpGrid with LGRs, ecl/vtk output is not supported yet.
673 const auto& grid = this->simulator().vanguard().gridView().grid();
674 using GridType = std::remove_cv_t< typename std::remove_reference<decltype(grid)>::type>;
675 bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
676 if ( !isCpGrid || (this->simulator().vanguard().gridView().grid().maxLevel()==0)) {
677 eclWriter_->evalSummaryState(isSubStep);
678 }
679
680 int episodeIdx = this->episodeIndex();
681
682 // Re-ordering in case of Alugrid
683 std::function<unsigned int(unsigned int)> gridToEquilGrid = [&simulator](unsigned int i) {
684 return simulator.vanguard().gridIdxToEquilGridIdx(i);
685 };
686
687 std::function<void(bool)> transUp =
688 [this,gridToEquilGrid](bool global) {
689 this->transmissibilities_.update(global,gridToEquilGrid);
690 };
691 {
692 OPM_TIMEBLOCK(applyActions);
693 actionHandler_.applyActions(episodeIdx,
694 simulator.time() + simulator.timeStepSize(),
695 transUp);
696 }
697 // deal with "clogging" for the MICP model
698 if constexpr (enableMICP){
699 auto& model = this->model();
700 const auto& residual = this->model().linearizer().residual();
701 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
702 auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx];
703 MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
704 }
705 }
706 }
707
712 {
713 OPM_TIMEBLOCK(endEpisode);
714 auto& simulator = this->simulator();
715 auto& schedule = simulator.vanguard().schedule();
716
717 wellModel_.endEpisode();
718 aquiferModel_.endEpisode();
719
720 int episodeIdx = this->episodeIndex();
721 // check if we're finished ...
722 if (episodeIdx + 1 >= static_cast<int>(schedule.size() - 1)) {
723 simulator.setFinished(true);
724 return;
725 }
726
727 // .. if we're not yet done, start the next episode (report step)
728 simulator.startNextEpisode(schedule.stepLength(episodeIdx + 1));
729 }
730
735 void writeOutput(const SimulatorTimer& timer, bool verbose = true)
736 {
737 OPM_TIMEBLOCK(problemWriteOutput);
738 // use the generic code to prepare the output fields and to
739 // write the desired VTK files.
740 if (Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() ||
741 this->simulator().episodeWillBeOver()) {
742 ParentType::writeOutput(verbose);
743 }
744
745 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
746 !this->simulator().episodeWillBeOver();
747
748 data::Solution localCellData = {};
749#if HAVE_DAMARIS
750 // N.B. the Damaris output has to be done before the ECL output as the ECL one
751 // does all kinds of std::move() relocation of data
752 if (enableDamarisOutput_) {
753 damarisWriter_->writeOutput(localCellData, isSubStep) ;
754 }
755#endif
756 if (enableEclOutput_){
757 eclWriter_->writeOutput(std::move(localCellData), timer, isSubStep);
758 }
759 }
760
761 void finalizeOutput() {
762 OPM_TIMEBLOCK(finalizeOutput);
763 // this will write all pending output to disk
764 // to avoid corruption of output files
765 eclWriter_.reset();
766 }
767
768
772 template <class Context>
773 const DimMatrix& intrinsicPermeability(const Context& context,
774 unsigned spaceIdx,
775 unsigned timeIdx) const
776 {
777 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
778 return transmissibilities_.permeability(globalSpaceIdx);
779 }
780
787 const DimMatrix& intrinsicPermeability(unsigned globalElemIdx) const
788 { return transmissibilities_.permeability(globalElemIdx); }
789
793 template <class Context>
794 Scalar transmissibility(const Context& context,
795 [[maybe_unused]] unsigned fromDofLocalIdx,
796 unsigned toDofLocalIdx) const
797 {
798 assert(fromDofLocalIdx == 0);
799 return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
800 }
801
805 Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
806 {
807 return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
808 }
809
813 template <class Context>
814 Scalar diffusivity(const Context& context,
815 [[maybe_unused]] unsigned fromDofLocalIdx,
816 unsigned toDofLocalIdx) const
817 {
818 assert(fromDofLocalIdx == 0);
819 return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
820 }
821
825 Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
826 return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
827 }
828
832 Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
833 return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
834 }
835
839 Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx,
840 const unsigned boundaryFaceIdx) const
841 {
842 return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
843 }
844
845
846
847
851 template <class Context>
852 Scalar transmissibilityBoundary(const Context& elemCtx,
853 unsigned boundaryFaceIdx) const
854 {
855 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
856 return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
857 }
858
862 Scalar transmissibilityBoundary(const unsigned globalSpaceIdx,
863 const unsigned boundaryFaceIdx) const
864 {
865 return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
866 }
867
868
872 Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn,
873 const unsigned globalSpaceIdxOut) const
874 {
875 return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
876 }
877
881 template <class Context>
882 Scalar thermalHalfTransmissibilityIn(const Context& context,
883 unsigned faceIdx,
884 unsigned timeIdx) const
885 {
886 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
887 unsigned toDofLocalIdx = face.exteriorIndex();
888 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
889 }
890
894 template <class Context>
895 Scalar thermalHalfTransmissibilityOut(const Context& context,
896 unsigned faceIdx,
897 unsigned timeIdx) const
898 {
899 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
900 unsigned toDofLocalIdx = face.exteriorIndex();
901 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
902 }
903
907 template <class Context>
908 Scalar thermalHalfTransmissibilityBoundary(const Context& elemCtx,
909 unsigned boundaryFaceIdx) const
910 {
911 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
912 return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
913 }
914
918 const typename Vanguard::TransmissibilityType& eclTransmissibilities() const
919 { return transmissibilities_; }
920
924 Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
925 { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
926
927 const FlowThresholdPressure<TypeTag>& thresholdPressure() const
928 { return thresholdPressures_; }
929
930 FlowThresholdPressure<TypeTag>& thresholdPressure()
931 { return thresholdPressures_; }
932
933 const TracerModel& tracerModel() const
934 { return tracerModel_; }
935
936 TracerModel& tracerModel()
937 { return tracerModel_; }
938
947 template <class Context>
948 Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
949 {
950 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
951 return this->porosity(globalSpaceIdx, timeIdx);
952 }
953
960 template <class Context>
961 Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
962 {
963 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
964 return this->dofCenterDepth(globalSpaceIdx);
965 }
966
973 Scalar dofCenterDepth(unsigned globalSpaceIdx) const
974 {
975 return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
976 }
977
981 template <class Context>
982 Scalar rockCompressibility(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
983 {
984 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
985 return this->rockCompressibility(globalSpaceIdx);
986 }
987
991 template <class Context>
992 Scalar rockReferencePressure(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
993 {
994 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
995 return this->rockReferencePressure(globalSpaceIdx);
996 }
997
1001 template <class Context>
1002 const MaterialLawParams& materialLawParams(const Context& context,
1003 unsigned spaceIdx, unsigned timeIdx) const
1004 {
1005 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1006 return this->materialLawParams(globalSpaceIdx);
1007 }
1008
1009 const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const
1010 {
1011 return materialLawManager_->materialLawParams(globalDofIdx);
1012 }
1013
1014 const MaterialLawParams& materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const
1015 {
1016 return materialLawManager_->materialLawParams(globalDofIdx, facedir);
1017 }
1018
1022 template <class Context>
1023 const SolidEnergyLawParams&
1024 solidEnergyLawParams(const Context& context,
1025 unsigned spaceIdx,
1026 unsigned timeIdx) const
1027 {
1028 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1029 return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1030 }
1031
1035 template <class Context>
1036 const ThermalConductionLawParams &
1037 thermalConductionLawParams(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1038 {
1039 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1040 return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1041 }
1042
1049 std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
1050 { return materialLawManager_; }
1051
1052 template <class FluidState>
1053 void updateRelperms(
1054 std::array<Evaluation,numPhases> &mobility,
1055 DirectionalMobilityPtr &dirMob,
1056 FluidState &fluidState,
1057 unsigned globalSpaceIdx) const
1058 {
1059 OPM_TIMEBLOCK_LOCAL(updateRelperms);
1060 {
1061 // calculate relative permeabilities. note that we store the result into the
1062 // mobility_ class attribute. the division by the phase viscosity happens later.
1063 const auto& materialParams = materialLawParams(globalSpaceIdx);
1064 MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
1065 Valgrind::CheckDefined(mobility);
1066 }
1067 if (materialLawManager_->hasDirectionalRelperms()
1068 || materialLawManager_->hasDirectionalImbnum())
1069 {
1070 using Dir = FaceDir::DirEnum;
1071 constexpr int ndim = 3;
1072 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
1073 Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
1074 for (int i = 0; i<ndim; i++) {
1075 const auto& materialParams = materialLawParams(globalSpaceIdx, facedirs[i]);
1076 auto& mob_array = dirMob->getArray(i);
1077 MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
1078 }
1079 }
1080 }
1081
1085 std::shared_ptr<EclMaterialLawManager> materialLawManager()
1086 { return materialLawManager_; }
1087
1088 using BaseType::pvtRegionIndex;
1092 template <class Context>
1093 unsigned pvtRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1094 { return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1095
1096 using BaseType::satnumRegionIndex;
1100 template <class Context>
1101 unsigned satnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1102 { return this->satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1103
1104 using BaseType::miscnumRegionIndex;
1108 template <class Context>
1109 unsigned miscnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1110 { return this->miscnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1111
1112 using BaseType::plmixnumRegionIndex;
1116 template <class Context>
1117 unsigned plmixnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1118 { return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1119
1120 using BaseType::maxPolymerAdsorption;
1124 template <class Context>
1125 Scalar maxPolymerAdsorption(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1126 { return this->maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1127
1131 std::string name() const
1132 { return this->simulator().vanguard().caseName(); }
1133
1137 template <class Context>
1138 Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1139 {
1140 // use the initial temperature of the DOF if temperature is not a primary
1141 // variable
1142 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1143 return initialFluidStates_[globalDofIdx].temperature(/*phaseIdx=*/0);
1144 }
1145
1146
1147 Scalar temperature(unsigned globalDofIdx, unsigned /*timeIdx*/) const
1148 {
1149 // use the initial temperature of the DOF if temperature is not a primary
1150 // variable
1151 return initialFluidStates_[globalDofIdx].temperature(/*phaseIdx=*/0);
1152 }
1153
1154 const SolidEnergyLawParams&
1155 solidEnergyLawParams(unsigned globalSpaceIdx,
1156 unsigned /*timeIdx*/) const
1157 {
1158 return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1159 }
1160 const ThermalConductionLawParams &
1161 thermalConductionLawParams(unsigned globalSpaceIdx,
1162 unsigned /*timeIdx*/)const
1163 {
1164 return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1165 }
1166
1172 template <class Context>
1173 void boundary(BoundaryRateVector& values,
1174 const Context& context,
1175 unsigned spaceIdx,
1176 unsigned timeIdx) const
1177 {
1178 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
1179 if (!context.intersection(spaceIdx).boundary())
1180 return;
1181
1182 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
1183 values.setNoFlow();
1184 else {
1185 // in the energy case we need to specify a non-trivial boundary condition
1186 // because the geothermal gradient needs to be maintained. for this, we
1187 // simply assume the initial temperature at the boundary and specify the
1188 // thermal flow accordingly. in this context, "thermal flow" means energy
1189 // flow due to a temerature gradient while assuming no-flow for mass
1190 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1191 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1192 values.setThermalFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
1193 }
1194
1195 if (nonTrivialBoundaryConditions()) {
1196 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
1197 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1198 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1199 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
1200 const auto [type, massrate] = boundaryCondition(globalDofIdx, indexInInside);
1201 if (type == BCType::THERMAL)
1202 values.setThermalFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1203 else if (type == BCType::FREE || type == BCType::DIRICHLET)
1204 values.setFreeFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1205 else if (type == BCType::RATE)
1206 values.setMassRate(massrate, pvtRegionIdx);
1207 }
1208 }
1209
1219 Scalar maxOilSaturation(unsigned globalDofIdx) const
1220 {
1221 if (!this->vapparsActive(this->episodeIndex()))
1222 return 0.0;
1223
1224 return this->maxOilSaturation_[globalDofIdx];
1225 }
1226
1236 void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
1237 {
1238 if (!this->vapparsActive(this->episodeIndex()))
1239 return;
1240
1241 this->maxOilSaturation_[globalDofIdx] = value;
1242 }
1243
1248 Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
1249 {
1250 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
1251 this->episodeIndex(),
1252 this->pvtRegionIndex(globalDofIdx));
1253 }
1254
1259 Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
1260 {
1261 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
1262 this->episodeIndex(),
1263 this->pvtRegionIndex(globalDofIdx));
1264 }
1265
1275 {
1276 int episodeIdx = this->episodeIndex();
1277 return !this->mixControls_.drsdtActive(episodeIdx) &&
1278 !this->mixControls_.drvdtActive(episodeIdx) &&
1279 this->rockCompPoroMultWc_.empty() &&
1280 this->rockCompPoroMult_.empty();
1281 }
1282
1289 template <class Context>
1290 void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1291 {
1292 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1293
1294 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
1295 values.assignNaive(initialFluidStates_[globalDofIdx]);
1296
1297 SolventModule::assignPrimaryVars(values,
1298 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
1299 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
1300
1301 if constexpr (enablePolymer)
1302 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
1303
1304 if constexpr (enablePolymerMolarWeight)
1305 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
1306
1307 if constexpr (enableBrine) {
1308 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
1309 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
1310 }
1311 else {
1312 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
1313 }
1314 }
1315
1316 if constexpr (enableMICP){
1317 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
1318 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
1319 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
1320 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
1321 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
1322 }
1323
1324 values.checkDefined();
1325 }
1326
1331 {
1332 // Calculate all intensive quantities.
1333 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/0);
1334
1335 // We also need the intensive quantities for timeIdx == 1
1336 // corresponding to the start of the current timestep, if we
1337 // do not use the storage cache, or if we cannot recycle the
1338 // first iteration storage.
1339 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
1340 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/1);
1341 }
1342
1343 // initialize the wells. Note that this needs to be done after initializing the
1344 // intrinsic permeabilities and the after applying the initial solution because
1345 // the well model uses these...
1346 wellModel_.init();
1347
1348 // let the object for threshold pressures initialize itself. this is done only at
1349 // this point, because determining the threshold pressures may require to access
1350 // the initial solution.
1351 thresholdPressures_.finishInit();
1352
1353 updateCompositionChangeLimits_();
1354
1355 aquiferModel_.initialSolutionApplied();
1356
1357 if (this->simulator().episodeIndex() == 0) {
1358 eclWriter_->writeInitialFIPReport();
1359 }
1360 }
1361
1367 template <class Context>
1368 void source(RateVector& rate,
1369 const Context& context,
1370 unsigned spaceIdx,
1371 unsigned timeIdx) const
1372 {
1373 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1374 source(rate, globalDofIdx, timeIdx);
1375 }
1376
1377 void source(RateVector& rate,
1378 unsigned globalDofIdx,
1379 unsigned timeIdx) const
1380 {
1381 OPM_TIMEBLOCK_LOCAL(eclProblemSource);
1382 rate = 0.0;
1383
1384 // Add well contribution to source here.
1385 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
1386
1387 // convert the source term from the total mass rate of the
1388 // cell to the one per unit of volume as used by the model.
1389 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1390 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
1391
1392 Valgrind::CheckDefined(rate[eqIdx]);
1393 assert(isfinite(rate[eqIdx]));
1394 }
1395
1396 // Add non-well sources.
1397 addToSourceDense(rate, globalDofIdx, timeIdx);
1398 }
1399
1400 void addToSourceDense(RateVector& rate,
1401 unsigned globalDofIdx,
1402 unsigned timeIdx) const
1403 {
1404 aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
1405
1406 // Add source term from deck
1407 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
1408 std::array<int,3> ijk;
1409 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
1410
1411 if (source.hasSource(ijk)) {
1412 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1413 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
1414 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
1415 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
1416
1417 for (unsigned i = 0; i < phidx_map.size(); ++i) {
1418 const auto phaseIdx = phidx_map[i];
1419 const auto sourceComp = sc_map[i];
1420 const auto compIdx = cidx_map[i];
1421 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1422 continue;
1423 }
1424 Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1425 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1426 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
1427 }
1428 rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate;
1429 }
1430
1431 if constexpr (enableSolvent) {
1432 Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
1433 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1434 const auto& solventPvt = SolventModule::solventPvt();
1435 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
1436 }
1437 rate[Indices::contiSolventEqIdx] += mass_rate;
1438 }
1439 if constexpr (enablePolymer) {
1440 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
1441 }
1442 if constexpr (enableEnergy) {
1443 for (unsigned i = 0; i < phidx_map.size(); ++i) {
1444 const auto phaseIdx = phidx_map[i];
1445 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1446 continue;
1447 }
1448 const auto sourceComp = sc_map[i];
1449 if (source.hasHrate({ijk, sourceComp})) {
1450 rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1451 } else {
1452 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, /*timeIdx*/ 0);
1453 auto fs = intQuants.fluidState();
1454 // if temperature is not set, use cell temperature as default
1455 if (source.hasTemperature({ijk, sourceComp})) {
1456 Scalar temperature = source.temperature({ijk, sourceComp});
1457 fs.setTemperature(temperature);
1458 }
1459 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
1460 Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx);
1461 Scalar energy_rate = getValue(h)*mass_rate;
1462 rate[Indices::contiEnergyEqIdx] += energy_rate;
1463 }
1464 }
1465 }
1466 }
1467
1468 // if requested, compensate systematic mass loss for cells which were "well
1469 // behaved" in the last time step
1470 // Note that we don't allow for drift compensation if there are no active wells.
1471 const bool compensateDrift = wellModel_.wellsActive();
1472 if (enableDriftCompensation_ && compensateDrift) {
1473 const auto& simulator = this->simulator();
1474 const auto& model = this->model();
1475
1476 // we use a lower tolerance for the compensation too
1477 // assure the added drift from the last step does not
1478 // cause convergence issues on the current step
1479 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
1480 Scalar poro = this->porosity(globalDofIdx, timeIdx);
1481 Scalar dt = simulator.timeStepSize();
1482 EqVector dofDriftRate = drift_[globalDofIdx];
1483 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
1484
1485 // restrict drift compensation to the CNV tolerance
1486 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1487 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
1488 if (cnv > maxCompensation) {
1489 dofDriftRate[eqIdx] *= maxCompensation/cnv;
1490 }
1491 }
1492
1493 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
1494 rate[eqIdx] -= dofDriftRate[eqIdx];
1495 }
1496 }
1497
1503 const WellModel& wellModel() const
1504 { return wellModel_; }
1505
1506 WellModel& wellModel()
1507 { return wellModel_; }
1508
1509 const AquiferModel& aquiferModel() const
1510 { return aquiferModel_; }
1511
1512 AquiferModel& mutableAquiferModel()
1513 { return aquiferModel_; }
1514
1515 // temporary solution to facilitate output of initial state from flow
1516 const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
1517 { return initialFluidStates_[globalDofIdx]; }
1518
1519 const EclipseIO& eclIO() const
1520 { return eclWriter_->eclIO(); }
1521
1522 void setSubStepReport(const SimulatorReportSingle& report)
1523 { return eclWriter_->setSubStepReport(report); }
1524
1525 void setSimulationReport(const SimulatorReport& report)
1526 { return eclWriter_->setSimulationReport(report); }
1527
1528 bool nonTrivialBoundaryConditions() const
1529 { return nonTrivialBoundaryConditions_; }
1530
1531 const InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const
1532 {
1533 OPM_TIMEBLOCK_LOCAL(boundaryFluidState);
1534 const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop;
1535 if (bcprop.size() > 0) {
1536 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1537
1538 // index == 0: no boundary conditions for this
1539 // global cell and direction
1540 if (bcindex_(dir)[globalDofIdx] == 0)
1541 return initialFluidStates_[globalDofIdx];
1542
1543 const auto& bc = bcprop[bcindex_(dir)[globalDofIdx]];
1544 if (bc.bctype == BCType::DIRICHLET )
1545 {
1546 InitialFluidState fluidState;
1547 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1548 fluidState.setPvtRegionIndex(pvtRegionIdx);
1549
1550 switch (bc.component) {
1551 case BCComponent::OIL:
1552 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
1553 throw std::logic_error("oil is not active and you're trying to add oil BC");
1554
1555 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
1556 break;
1557 case BCComponent::GAS:
1558 if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1559 throw std::logic_error("gas is not active and you're trying to add gas BC");
1560
1561 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
1562 break;
1563 case BCComponent::WATER:
1564 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1565 throw std::logic_error("water is not active and you're trying to add water BC");
1566
1567 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
1568 break;
1569 case BCComponent::SOLVENT:
1570 case BCComponent::POLYMER:
1571 case BCComponent::NONE:
1572 throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
1573 break;
1574 }
1575 double pressure = initialFluidStates_[globalDofIdx].pressure(refPressurePhaseIdx_());
1576 const auto pressure_input = bc.pressure;
1577 if (pressure_input) {
1578 pressure = *pressure_input;
1579 }
1580
1581 std::array<Scalar, numPhases> pc = {0};
1582 const auto& matParams = materialLawParams(globalDofIdx);
1583 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
1584 Valgrind::CheckDefined(pressure);
1585 Valgrind::CheckDefined(pc);
1586 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1587 if (!FluidSystem::phaseIsActive(phaseIdx))
1588 continue;
1589
1590 if (Indices::oilEnabled)
1591 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1592 else if (Indices::gasEnabled)
1593 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1594 else if (Indices::waterEnabled)
1595 //single (water) phase
1596 fluidState.setPressure(phaseIdx, pressure);
1597 }
1598
1599 double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature
1600 const auto temperature_input = bc.temperature;
1601 if(temperature_input)
1602 temperature = *temperature_input;
1603 fluidState.setTemperature(temperature);
1604
1605 if (FluidSystem::enableDissolvedGas()) {
1606 fluidState.setRs(0.0);
1607 fluidState.setRv(0.0);
1608 }
1609 if (FluidSystem::enableDissolvedGasInWater()) {
1610 fluidState.setRsw(0.0);
1611 }
1612 if (FluidSystem::enableVaporizedWater())
1613 fluidState.setRvw(0.0);
1614
1615 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1616 if (!FluidSystem::phaseIsActive(phaseIdx))
1617 continue;
1618
1619 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
1620 fluidState.setInvB(phaseIdx, b);
1621
1622 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
1623 fluidState.setDensity(phaseIdx, rho);
1624 if (enableEnergy) {
1625 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
1626 fluidState.setEnthalpy(phaseIdx, h);
1627 }
1628 }
1629 fluidState.checkDefined();
1630 return fluidState;
1631 }
1632 }
1633 return initialFluidStates_[globalDofIdx];
1634 }
1635
1642 Scalar nextTimeStepSize() const
1643 {
1644 OPM_TIMEBLOCK(nexTimeStepSize);
1645 // allow external code to do the timestepping
1646 if (this->nextTimeStepSize_ > 0.0)
1647 return this->nextTimeStepSize_;
1648
1649 const auto& simulator = this->simulator();
1650 int episodeIdx = simulator.episodeIndex();
1651
1652 // for the initial episode, we use a fixed time step size
1653 if (episodeIdx < 0)
1654 return this->initialTimeStepSize_;
1655
1656 // ask the newton method for a suggestion. This suggestion will be based on how
1657 // well the previous time step converged. After that, apply the runtime time
1658 // stepping constraints.
1659 const auto& newtonMethod = this->model().newtonMethod();
1660 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1661 }
1662
1668 template <class LhsEval>
1669 LhsEval rockCompPoroMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1670 {
1671 OPM_TIMEBLOCK_LOCAL(rockCompPoroMultiplier);
1672 if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
1673 return 1.0;
1674
1675 unsigned tableIdx = 0;
1676 if (!this->rockTableIdx_.empty())
1677 tableIdx = this->rockTableIdx_[elementIdx];
1678
1679 const auto& fs = intQuants.fluidState();
1680 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1681 if (!this->minRefPressure_.empty())
1682 // The pore space change is irreversible
1683 effectivePressure =
1684 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1685 this->minRefPressure_[elementIdx]);
1686
1687 if (!this->overburdenPressure_.empty())
1688 effectivePressure -= this->overburdenPressure_[elementIdx];
1689
1690
1691 if (!this->rockCompPoroMult_.empty()) {
1692 return this->rockCompPoroMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
1693 }
1694
1695 // water compaction
1696 assert(!this->rockCompPoroMultWc_.empty());
1697 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1698 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
1699
1700 return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
1701 }
1702
1708 template <class LhsEval>
1709 LhsEval rockCompTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1710 {
1711 bool implicit = !Parameters::get<TypeTag, Properties::ExplicitRockCompaction>();
1712 return implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<LhsEval>(intQuants, elementIdx)
1713 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1714 }
1715
1721 template <class LhsEval>
1722 LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants) const
1723 {
1724 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
1725 if (!enableSaltPrecipitation)
1726 return 1.0;
1727
1728 const auto& fs = intQuants.fluidState();
1729 unsigned tableIdx = fs.pvtRegionIndex();
1730 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
1731 porosityFactor = min(porosityFactor, 1.0);
1732 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
1733 return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
1734 }
1735
1739 template <class LhsEval>
1740 LhsEval wellTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1741 {
1742 OPM_TIMEBLOCK_LOCAL(wellTransMultiplier);
1743
1744 bool implicit = !Parameters::get<TypeTag, Properties::ExplicitRockCompaction>();
1745 double trans_mult = implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<double>(intQuants, elementIdx)
1746 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1747 trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants);
1748
1749 return trans_mult;
1750 }
1751
1752 std::pair<BCType, RateVector> boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) const
1753 {
1754 OPM_TIMEBLOCK_LOCAL(boundaryCondition);
1755 if (!nonTrivialBoundaryConditions_) {
1756 return { BCType::NONE, RateVector(0.0) };
1757 }
1758 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1759 const auto& schedule = this->simulator().vanguard().schedule();
1760 if (bcindex_(dir)[globalSpaceIdx] == 0) {
1761 return { BCType::NONE, RateVector(0.0) };
1762 }
1763 if (schedule[this->episodeIndex()].bcprop.size() == 0) {
1764 return { BCType::NONE, RateVector(0.0) };
1765 }
1766 const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
1767 if (bc.bctype!=BCType::RATE) {
1768 return { bc.bctype, RateVector(0.0) };
1769 }
1770
1771 RateVector rate = 0.0;
1772 switch (bc.component) {
1773 case BCComponent::OIL:
1774 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = bc.rate;
1775 break;
1776 case BCComponent::GAS:
1777 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = bc.rate;
1778 break;
1779 case BCComponent::WATER:
1780 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate;
1781 break;
1782 case BCComponent::SOLVENT:
1783 if constexpr (!enableSolvent)
1784 throw std::logic_error("solvent is disabled and you're trying to add solvent to BC");
1785
1786 rate[Indices::solventSaturationIdx] = bc.rate;
1787 break;
1788 case BCComponent::POLYMER:
1789 if constexpr (!enablePolymer)
1790 throw std::logic_error("polymer is disabled and you're trying to add polymer to BC");
1791
1792 rate[Indices::polymerConcentrationIdx] = bc.rate;
1793 break;
1794 case BCComponent::NONE:
1795 throw std::logic_error("you need to specify the component when RATE type is set in BC");
1796 break;
1797 }
1798 //TODO add support for enthalpy rate
1799 return {bc.bctype, rate};
1800 }
1801
1802 const std::unique_ptr<EclWriterType>& eclWriter() const
1803 {
1804 return eclWriter_;
1805 }
1806
1807 void setConvData(const std::vector<std::vector<int>>& data)
1808 {
1809 eclWriter_->mutableOutputModule().setCnvData(data);
1810 }
1811
1812 template<class Serializer>
1813 void serializeOp(Serializer& serializer)
1814 {
1815 serializer(static_cast<BaseType&>(*this));
1816 serializer(drift_);
1817 serializer(wellModel_);
1818 serializer(aquiferModel_);
1819 serializer(tracerModel_);
1820 serializer(*materialLawManager_);
1821 serializer(*eclWriter_);
1822 }
1823private:
1824 Implementation& asImp_()
1825 { return *static_cast<Implementation *>(this); }
1826protected:
1827 void updateExplicitQuantities_()
1828 {
1829 OPM_TIMEBLOCK(updateExplicitQuantities);
1830 const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
1831 const bool invalidateFromMinPressure = updateMinPressure_();
1832
1833 // update hysteresis and max oil saturation used in vappars
1834 const bool invalidateFromHyst = updateHysteresis_();
1835 const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
1836
1837 // the derivatives may have change
1838 bool invalidateIntensiveQuantities
1839 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
1840 if (invalidateIntensiveQuantities) {
1841 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1842 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
1843 }
1844
1845 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
1846 updateMaxPolymerAdsorption_();
1847
1848 updateRockCompTransMultVal_();
1849 }
1850
1851 template<class UpdateFunc>
1852 void updateProperty_(const std::string& failureMsg,
1853 UpdateFunc func)
1854 {
1855 OPM_TIMEBLOCK(updateProperty);
1856 const auto& model = this->simulator().model();
1857 const auto& primaryVars = model.solution(/*timeIdx*/0);
1858 const auto& vanguard = this->simulator().vanguard();
1859 std::size_t numGridDof = primaryVars.size();
1860 OPM_BEGIN_PARALLEL_TRY_CATCH();
1861#ifdef _OPENMP
1862#pragma omp parallel for
1863#endif
1864 for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1865 const auto& iq = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/ 0);
1866 func(dofIdx, iq);
1867 }
1868 OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
1869 }
1870
1871 // update the parameters needed for DRSDT and DRVDT
1872 void updateCompositionChangeLimits_()
1873 {
1874 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1875 // update the "last Rs" values for all elements, including the ones in the ghost
1876 // and overlap regions
1877 int episodeIdx = this->episodeIndex();
1878 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1879 this->mixControls_.drsdtActive(episodeIdx),
1880 this->mixControls_.drvdtActive(episodeIdx)};
1881 if (!active[0] && !active[1] && !active[2]) {
1882 return;
1883 }
1884
1885 this->updateProperty_("FlowProblem::updateCompositionChangeLimits_()) failed:",
1886 [this,episodeIdx,active](unsigned compressedDofIdx,
1887 const IntensiveQuantities& iq)
1888 {
1889 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1890 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1891 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1892 this->mixControls_.update(compressedDofIdx,
1893 iq,
1894 episodeIdx,
1895 this->gravity_[dim - 1],
1896 perm[dim - 1][dim - 1],
1897 distZ,
1898 pvtRegionIdx,
1899 active);
1900 }
1901 );
1902 }
1903
1904 bool updateMaxOilSaturation_()
1905 {
1906 OPM_TIMEBLOCK(updateMaxOilSaturation);
1907 int episodeIdx = this->episodeIndex();
1908
1909 // we use VAPPARS
1910 if (this->vapparsActive(episodeIdx)) {
1911 this->updateProperty_("FlowProblem::updateMaxOilSaturation_() failed:",
1912 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1913 {
1914 this->updateMaxOilSaturation_(compressedDofIdx,iq);
1915 });
1916 return true;
1917 }
1918
1919 return false;
1920 }
1921
1922 bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1923 {
1924 OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
1925 const auto& fs = iq.fluidState();
1926 const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1927 auto& mos = this->maxOilSaturation_;
1928 if(mos[compressedDofIdx] < So){
1929 mos[compressedDofIdx] = So;
1930 return true;
1931 }else{
1932 return false;
1933 }
1934 }
1935
1936 bool updateMaxWaterSaturation_()
1937 {
1938 OPM_TIMEBLOCK(updateMaxWaterSaturation);
1939 // water compaction is activated in ROCKCOMP
1940 if (this->maxWaterSaturation_.empty())
1941 return false;
1942
1943 this->maxWaterSaturation_[/*timeIdx=*/1] = this->maxWaterSaturation_[/*timeIdx=*/0];
1944 this->updateProperty_("FlowProblem::updateMaxWaterSaturation_() failed:",
1945 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1946 {
1947 this->updateMaxWaterSaturation_(compressedDofIdx,iq);
1948 });
1949 return true;
1950 }
1951
1952
1953 bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1954 {
1955 OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation);
1956 const auto& fs = iq.fluidState();
1957 const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
1958 auto& mow = this->maxWaterSaturation_;
1959 if(mow[compressedDofIdx]< Sw){
1960 mow[compressedDofIdx] = Sw;
1961 return true;
1962 }else{
1963 return false;
1964 }
1965 }
1966
1967 bool updateMinPressure_()
1968 {
1969 OPM_TIMEBLOCK(updateMinPressure);
1970 // IRREVERS option is used in ROCKCOMP
1971 if (this->minRefPressure_.empty())
1972 return false;
1973
1974 this->updateProperty_("FlowProblem::updateMinPressure_() failed:",
1975 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1976 {
1977 this->updateMinPressure_(compressedDofIdx,iq);
1978 });
1979 return true;
1980 }
1981
1982 bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
1983 OPM_TIMEBLOCK_LOCAL(updateMinPressure);
1984 const auto& fs = iq.fluidState();
1985 const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
1986 auto& min_pressures = this->minRefPressure_;
1987 if(min_pressures[compressedDofIdx]> min_pressure){
1988 min_pressures[compressedDofIdx] = min_pressure;
1989 return true;
1990 }else{
1991 return false;
1992 }
1993 }
1994
1995 // \brief Function to assign field properties of type double, on the leaf grid view.
1996 //
1997 // For CpGrid with local grid refinement, the field property of a cell on the leaf
1998 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
1999 std::function<std::vector<double>(const FieldPropsManager&, const std::string&)>
2000 fieldPropDoubleOnLeafAssigner_()
2001 {
2002 const auto& lookup = this->lookUpData_;
2003 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString)
2004 {
2005 return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
2006 };
2007 }
2008
2009 // \brief Function to assign field properties of type int, unsigned int, ..., on the leaf grid view.
2010 //
2011 // For CpGrid with local grid refinement, the field property of a cell on the leaf
2012 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
2013 template<typename IntType>
2014 std::function<std::vector<IntType>(const FieldPropsManager&, const std::string&, bool)>
2015 fieldPropIntTypeOnLeafAssigner_()
2016 {
2017 const auto& lookup = this->lookUpData_;
2018 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString, bool needsTranslation)
2019 {
2020 return lookup.template assignFieldPropsIntOnLeaf<IntType>(fieldPropManager, propString, needsTranslation);
2021 };
2022 }
2023
2024 void readMaterialParameters_()
2025 {
2026 OPM_TIMEBLOCK(readMaterialParameters);
2027 const auto& simulator = this->simulator();
2028 const auto& vanguard = simulator.vanguard();
2029 const auto& eclState = vanguard.eclState();
2030
2031 // the PVT and saturation region numbers
2032 OPM_BEGIN_PARALLEL_TRY_CATCH();
2033 this->updatePvtnum_();
2034 this->updateSatnum_();
2035
2036 // the MISC region numbers (solvent model)
2037 this->updateMiscnum_();
2038 // the PLMIX region numbers (polymer model)
2039 this->updatePlmixnum_();
2040
2041 // directional relative permeabilities
2042 this->updateKrnum_();
2043 OPM_END_PARALLEL_TRY_CATCH("Invalid region numbers: ", vanguard.gridView().comm());
2045 // porosity
2046 updateReferencePorosity_();
2047 this->referencePorosity_[1] = this->referencePorosity_[0];
2049
2051 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
2052 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
2053 materialLawManager_->initFromState(eclState);
2054 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2055 this-> template fieldPropIntTypeOnLeafAssigner_<int>(),
2056 this-> lookupIdxOnLevelZeroAssigner_());
2058 }
2059
2060 void readThermalParameters_()
2061 {
2062 if constexpr (enableEnergy)
2063 {
2064 const auto& simulator = this->simulator();
2065 const auto& vanguard = simulator.vanguard();
2066 const auto& eclState = vanguard.eclState();
2067
2068 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
2069 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
2070 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2071 this-> fieldPropDoubleOnLeafAssigner_(),
2072 this-> template fieldPropIntTypeOnLeafAssigner_<unsigned int>());
2073 }
2074 }
2075
2076 void updateReferencePorosity_()
2077 {
2078 const auto& simulator = this->simulator();
2079 const auto& vanguard = simulator.vanguard();
2080 const auto& eclState = vanguard.eclState();
2081
2082 std::size_t numDof = this->model().numGridDof();
2083
2084 this->referencePorosity_[/*timeIdx=*/0].resize(numDof);
2085
2086 const auto& fp = eclState.fieldProps();
2087 const std::vector<double> porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV");
2088 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2089 Scalar poreVolume = porvData[dofIdx];
2090
2091 // we define the porosity as the accumulated pore volume divided by the
2092 // geometric volume of the element. Note that -- in pathetic cases -- it can
2093 // be larger than 1.0!
2094 Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx);
2095 assert(dofVolume > 0.0);
2096 this->referencePorosity_[/*timeIdx=*/0][dofIdx] = poreVolume/dofVolume;
2097 }
2098 }
2099
2100 void readInitialCondition_()
2101 {
2102 const auto& simulator = this->simulator();
2103 const auto& vanguard = simulator.vanguard();
2104 const auto& eclState = vanguard.eclState();
2105
2106 if (eclState.getInitConfig().hasEquil())
2107 readEquilInitialCondition_();
2108 else
2109 readExplicitInitialCondition_();
2110
2111 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
2112 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
2113 enableSolvent,
2114 enablePolymer,
2115 enablePolymerMolarWeight,
2116 enableMICP);
2117
2118 //initialize min/max values
2119 std::size_t numElems = this->model().numGridDof();
2120 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2121 const auto& fs = initialFluidStates_[elemIdx];
2122 if (!this->maxWaterSaturation_.empty())
2123 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
2124 if (!this->maxOilSaturation_.empty())
2125 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
2126 if (!this->minRefPressure_.empty())
2127 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
2128 }
2129
2130
2131 }
2132
2133 void readEquilInitialCondition_()
2134 {
2135 const auto& simulator = this->simulator();
2136
2137 // initial condition corresponds to hydrostatic conditions.
2138 EquilInitializer<TypeTag> equilInitializer(simulator, *materialLawManager_);
2139
2140 std::size_t numElems = this->model().numGridDof();
2141 initialFluidStates_.resize(numElems);
2142 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2143 auto& elemFluidState = initialFluidStates_[elemIdx];
2144 elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
2145 }
2146 }
2147
2148 void readEclRestartSolution_()
2149 {
2150 // Throw an exception if the grid has LGRs. Refined grid are not supported for restart.
2151 if(this->simulator().vanguard().grid().maxLevel() > 0) {
2152 throw std::invalid_argument("Refined grids are not yet supported for restart ");
2153 }
2154
2155 // Set the start time of the simulation
2156 auto& simulator = this->simulator();
2157 const auto& schedule = simulator.vanguard().schedule();
2158 const auto& eclState = simulator.vanguard().eclState();
2159 const auto& initconfig = eclState.getInitConfig();
2160 {
2161 int restart_step = initconfig.getRestartStep();
2162
2163 simulator.setTime(schedule.seconds(restart_step));
2164
2165 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
2166 schedule.stepLength(restart_step));
2167 simulator.setEpisodeIndex(restart_step);
2168 }
2169 eclWriter_->beginRestart();
2170
2171 Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength());
2172 simulator.setTimeStepSize(dt);
2173
2174 std::size_t numElems = this->model().numGridDof();
2175 initialFluidStates_.resize(numElems);
2176 if constexpr (enableSolvent) {
2177 this->solventSaturation_.resize(numElems, 0.0);
2178 this->solventRsw_.resize(numElems, 0.0);
2179 }
2180
2181 if constexpr (enablePolymer)
2182 this->polymer_.concentration.resize(numElems, 0.0);
2183
2184 if constexpr (enablePolymerMolarWeight) {
2185 const std::string msg {"Support of the RESTART for polymer molecular weight "
2186 "is not implemented yet. The polymer weight value will be "
2187 "zero when RESTART begins"};
2188 OpmLog::warning("NO_POLYMW_RESTART", msg);
2189 this->polymer_.moleWeight.resize(numElems, 0.0);
2190 }
2191
2192 if constexpr (enableMICP) {
2193 this->micp_.resize(numElems);
2194 }
2195
2196 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2197 auto& elemFluidState = initialFluidStates_[elemIdx];
2198 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
2199 eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
2200 eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
2201
2202 // Note: Function processRestartSaturations_() mutates the
2203 // 'ssol' argument--the value from the restart file--if solvent
2204 // is enabled. Then, store the updated solvent saturation into
2205 // 'solventSaturation_'. Otherwise, just pass a dummy value to
2206 // the function and discard the unchanged result. Do not index
2207 // into 'solventSaturation_' unless solvent is enabled.
2208 {
2209 auto ssol = enableSolvent
2210 ? eclWriter_->outputModule().getSolventSaturation(elemIdx)
2211 : Scalar(0);
2212
2213 processRestartSaturations_(elemFluidState, ssol);
2214
2215 if constexpr (enableSolvent) {
2216 this->solventSaturation_[elemIdx] = ssol;
2217 this->solventRsw_[elemIdx] = eclWriter_->outputModule().getSolventRsw(elemIdx);
2218 }
2219 }
2220
2221 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
2222
2223 if constexpr (enablePolymer)
2224 this->polymer_.concentration[elemIdx] = eclWriter_->outputModule().getPolymerConcentration(elemIdx);
2225 if constexpr (enableMICP){
2226 this->micp_.microbialConcentration[elemIdx] = eclWriter_->outputModule().getMicrobialConcentration(elemIdx);
2227 this->micp_.oxygenConcentration[elemIdx] = eclWriter_->outputModule().getOxygenConcentration(elemIdx);
2228 this->micp_.ureaConcentration[elemIdx] = eclWriter_->outputModule().getUreaConcentration(elemIdx);
2229 this->micp_.biofilmConcentration[elemIdx] = eclWriter_->outputModule().getBiofilmConcentration(elemIdx);
2230 this->micp_.calciteConcentration[elemIdx] = eclWriter_->outputModule().getCalciteConcentration(elemIdx);
2231 }
2232 // if we need to restart for polymer molecular weight simulation, we need to add related here
2233 }
2234
2235 const int episodeIdx = this->episodeIndex();
2236 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
2237
2238 // assign the restart solution to the current solution. note that we still need
2239 // to compute real initial solution after this because the initial fluid states
2240 // need to be correct for stuff like boundary conditions.
2241 auto& sol = this->model().solution(/*timeIdx=*/0);
2242 const auto& gridView = this->gridView();
2243 ElementContext elemCtx(simulator);
2244 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2245 elemCtx.updatePrimaryStencil(elem);
2246 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
2247 initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0);
2248 }
2249
2250 // make sure that the ghost and overlap entities exhibit the correct
2251 // solution. alternatively, this could be done in the loop above by also
2252 // considering non-interior elements. Since the initial() method might not work
2253 // 100% correctly for such elements, let's play safe and explicitly synchronize
2254 // using message passing.
2255 this->model().syncOverlap();
2256
2257 eclWriter_->endRestart();
2258 }
2259
2260 void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
2261 {
2262 // each phase needs to be above certain value to be claimed to be existing
2263 // this is used to recover some RESTART running with the defaulted single-precision format
2264 const Scalar smallSaturationTolerance = 1.e-6;
2265 Scalar sumSaturation = 0.0;
2266 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2267 if (FluidSystem::phaseIsActive(phaseIdx)) {
2268 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
2269 elemFluidState.setSaturation(phaseIdx, 0.0);
2270
2271 sumSaturation += elemFluidState.saturation(phaseIdx);
2272 }
2273
2274 }
2275 if constexpr (enableSolvent) {
2276 if (solventSaturation < smallSaturationTolerance)
2277 solventSaturation = 0.0;
2278
2279 sumSaturation += solventSaturation;
2280 }
2281
2282 assert(sumSaturation > 0.0);
2283
2284 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2285 if (FluidSystem::phaseIsActive(phaseIdx)) {
2286 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
2287 elemFluidState.setSaturation(phaseIdx, saturation);
2288 }
2289 }
2290 if constexpr (enableSolvent) {
2291 solventSaturation = solventSaturation / sumSaturation;
2292 }
2293 }
2294
2295 void readExplicitInitialCondition_()
2296 {
2297 const auto& simulator = this->simulator();
2298 const auto& vanguard = simulator.vanguard();
2299 const auto& eclState = vanguard.eclState();
2300 const auto& fp = eclState.fieldProps();
2301 bool has_swat = fp.has_double("SWAT");
2302 bool has_sgas = fp.has_double("SGAS");
2303 bool has_rs = fp.has_double("RS");
2304 bool has_rv = fp.has_double("RV");
2305 bool has_rvw = fp.has_double("RVW");
2306 bool has_pressure = fp.has_double("PRESSURE");
2307 bool has_salt = fp.has_double("SALT");
2308 bool has_saltp = fp.has_double("SALTP");
2309
2310 // make sure all required quantities are enables
2311 if (Indices::numPhases > 1) {
2312 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
2313 throw std::runtime_error("The ECL input file requires the presence of the SWAT keyword if "
2314 "the water phase is active");
2315 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
2316 throw std::runtime_error("The ECL input file requires the presence of the SGAS keyword if "
2317 "the gas phase is active");
2318 }
2319 if (!has_pressure)
2320 throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
2321 "keyword if the model is initialized explicitly");
2322 if (FluidSystem::enableDissolvedGas() && !has_rs)
2323 throw std::runtime_error("The ECL input file requires the RS keyword to be present if"
2324 " dissolved gas is enabled");
2325 if (FluidSystem::enableVaporizedOil() && !has_rv)
2326 throw std::runtime_error("The ECL input file requires the RV keyword to be present if"
2327 " vaporized oil is enabled");
2328 if (FluidSystem::enableVaporizedWater() && !has_rvw)
2329 throw std::runtime_error("The ECL input file requires the RVW keyword to be present if"
2330 " vaporized water is enabled");
2331 if (enableBrine && !has_salt)
2332 throw std::runtime_error("The ECL input file requires the SALT keyword to be present if"
2333 " brine is enabled and the model is initialized explicitly");
2334 if (enableSaltPrecipitation && !has_saltp)
2335 throw std::runtime_error("The ECL input file requires the SALTP keyword to be present if"
2336 " salt precipitation is enabled and the model is initialized explicitly");
2337
2338 std::size_t numDof = this->model().numGridDof();
2339
2340 initialFluidStates_.resize(numDof);
2341
2342 std::vector<double> waterSaturationData;
2343 std::vector<double> gasSaturationData;
2344 std::vector<double> pressureData;
2345 std::vector<double> rsData;
2346 std::vector<double> rvData;
2347 std::vector<double> rvwData;
2348 std::vector<double> tempiData;
2349 std::vector<double> saltData;
2350 std::vector<double> saltpData;
2351
2352 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
2353 waterSaturationData = fp.get_double("SWAT");
2354 else
2355 waterSaturationData.resize(numDof);
2356
2357 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
2358 gasSaturationData = fp.get_double("SGAS");
2359 else
2360 gasSaturationData.resize(numDof);
2361
2362 pressureData = fp.get_double("PRESSURE");
2363 if (FluidSystem::enableDissolvedGas())
2364 rsData = fp.get_double("RS");
2365
2366 if (FluidSystem::enableVaporizedOil())
2367 rvData = fp.get_double("RV");
2368
2369 if (FluidSystem::enableVaporizedWater())
2370 rvwData = fp.get_double("RVW");
2371
2372 // initial reservoir temperature
2373 tempiData = fp.get_double("TEMPI");
2374
2375 // initial salt concentration data
2376 if constexpr (enableBrine)
2377 saltData = fp.get_double("SALT");
2378
2379 // initial precipitated salt saturation data
2380 if constexpr (enableSaltPrecipitation)
2381 saltpData = fp.get_double("SALTP");
2382
2383 // calculate the initial fluid states
2384 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2385 auto& dofFluidState = initialFluidStates_[dofIdx];
2386
2387 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
2388
2390 // set temperature
2392 Scalar temperatureLoc = tempiData[dofIdx];
2393 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
2394 temperatureLoc = FluidSystem::surfaceTemperature;
2395 dofFluidState.setTemperature(temperatureLoc);
2396
2398 // set salt concentration
2400 if constexpr (enableBrine)
2401 dofFluidState.setSaltConcentration(saltData[dofIdx]);
2402
2404 // set precipitated salt saturation
2406 if constexpr (enableSaltPrecipitation)
2407 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
2408
2410 // set saturations
2412 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
2413 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
2414 waterSaturationData[dofIdx]);
2415
2416 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
2417 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
2418 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2419 1.0
2420 - waterSaturationData[dofIdx]);
2421 }
2422 else
2423 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2424 gasSaturationData[dofIdx]);
2425 }
2426 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
2427 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
2428 1.0
2429 - waterSaturationData[dofIdx]
2430 - gasSaturationData[dofIdx]);
2431
2433 // set phase pressures
2435 Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
2436
2437 // this assumes that capillary pressures only depend on the phase saturations
2438 // and possibly on temperature. (this is always the case for ECL problems.)
2439 std::array<Scalar, numPhases> pc = {0};
2440 const auto& matParams = materialLawParams(dofIdx);
2441 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
2442 Valgrind::CheckDefined(pressure);
2443 Valgrind::CheckDefined(pc);
2444 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2445 if (!FluidSystem::phaseIsActive(phaseIdx))
2446 continue;
2447
2448 if (Indices::oilEnabled)
2449 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
2450 else if (Indices::gasEnabled)
2451 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
2452 else if (Indices::waterEnabled)
2453 //single (water) phase
2454 dofFluidState.setPressure(phaseIdx, pressure);
2455 }
2456
2457 if (FluidSystem::enableDissolvedGas())
2458 dofFluidState.setRs(rsData[dofIdx]);
2459 else if (Indices::gasEnabled && Indices::oilEnabled)
2460 dofFluidState.setRs(0.0);
2461
2462 if (FluidSystem::enableVaporizedOil())
2463 dofFluidState.setRv(rvData[dofIdx]);
2464 else if (Indices::gasEnabled && Indices::oilEnabled)
2465 dofFluidState.setRv(0.0);
2466
2467 if (FluidSystem::enableVaporizedWater())
2468 dofFluidState.setRvw(rvwData[dofIdx]);
2469
2471 // set invB_
2473 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2474 if (!FluidSystem::phaseIsActive(phaseIdx))
2475 continue;
2476
2477 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2478 dofFluidState.setInvB(phaseIdx, b);
2479
2480 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2481 dofFluidState.setDensity(phaseIdx, rho);
2482
2483 }
2484 }
2485 }
2486
2487 // update the hysteresis parameters of the material laws for the whole grid
2488 bool updateHysteresis_()
2489 {
2490 if (!materialLawManager_->enableHysteresis())
2491 return false;
2492
2493 // we need to update the hysteresis data for _all_ elements (i.e., not just the
2494 // interior ones) to avoid desynchronization of the processes in the parallel case!
2495 this->updateProperty_("FlowProblem::updateHysteresis_() failed:",
2496 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
2497 {
2498 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2499 });
2500 return true;
2501 }
2502
2503
2504 bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
2505 {
2506 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
2507 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2508 //TODO change materials to give a bool
2509 return true;
2510 }
2511
2512 void updateMaxPolymerAdsorption_()
2513 {
2514 // we need to update the max polymer adsoption data for all elements
2515 this->updateProperty_("FlowProblem::updateMaxPolymerAdsorption_() failed:",
2516 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
2517 {
2518 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
2519 });
2520 }
2521
2522 bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
2523 {
2524 const Scalar pa = scalarValue(iq.polymerAdsorption());
2525 auto& mpa = this->polymer_.maxAdsorption;
2526 if (mpa[compressedDofIdx] < pa) {
2527 mpa[compressedDofIdx] = pa;
2528 return true;
2529 } else {
2530 return false;
2531 }
2532 }
2533
2534 Scalar getRockCompTransMultVal(std::size_t dofIdx) const
2535 {
2536 if (this->rockCompTransMultVal_.empty())
2537 return 1.0;
2538
2539 return this->rockCompTransMultVal_[dofIdx];
2540 }
2541
2542
2543private:
2544 struct PffDofData_
2545 {
2546 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
2547 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
2548 ConditionalStorage<enableDiffusion, Scalar> diffusivity;
2549 ConditionalStorage<enableDispersion, Scalar> dispersivity;
2550 Scalar transmissibility;
2551 };
2552
2553 // update the prefetch friendly data object
2554 void updatePffDofData_()
2555 {
2556 const auto& distFn =
2557 [this](PffDofData_& dofData,
2558 const Stencil& stencil,
2559 unsigned localDofIdx)
2560 -> void
2561 {
2562 const auto& elementMapper = this->model().elementMapper();
2563
2564 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
2565 if (localDofIdx != 0) {
2566 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(/*dofIdx=*/0));
2567 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
2568
2569 if constexpr (enableEnergy) {
2570 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
2571 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
2572 }
2573 if constexpr (enableDiffusion)
2574 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
2575 if (enableDispersion)
2576 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
2577 }
2578 };
2579
2580 pffDofData_.update(distFn);
2581 }
2582
2583 void readBoundaryConditions_()
2584 {
2585 const auto& simulator = this->simulator();
2586 const auto& vanguard = simulator.vanguard();
2587 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
2588 if (bcconfig.size() > 0) {
2589 nonTrivialBoundaryConditions_ = true;
2590
2591 std::size_t numCartDof = vanguard.cartesianSize();
2592 unsigned numElems = vanguard.gridView().size(/*codim=*/0);
2593 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
2594
2595 for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
2596 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
2597
2598 bcindex_.resize(numElems, 0);
2599 auto loopAndApply = [&cartesianToCompressedElemIdx,
2600 &vanguard](const auto& bcface,
2601 auto apply)
2602 {
2603 for (int i = bcface.i1; i <= bcface.i2; ++i) {
2604 for (int j = bcface.j1; j <= bcface.j2; ++j) {
2605 for (int k = bcface.k1; k <= bcface.k2; ++k) {
2606 std::array<int, 3> tmp = {i,j,k};
2607 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
2608 if (elemIdx >= 0)
2609 apply(elemIdx);
2610 }
2611 }
2612 }
2613 };
2614 for (const auto& bcface : bcconfig) {
2615 std::vector<int>& data = bcindex_(bcface.dir);
2616 const int index = bcface.index;
2617 loopAndApply(bcface,
2618 [&data,index](int elemIdx)
2619 { data[elemIdx] = index; });
2620 }
2621 }
2622 }
2623
2624 // this method applies the runtime constraints specified via the deck and/or command
2625 // line parameters for the size of the next time step.
2626 Scalar limitNextTimeStepSize_(Scalar dtNext) const
2627 {
2628 if constexpr (enableExperiments) {
2629 const auto& simulator = this->simulator();
2630 const auto& schedule = simulator.vanguard().schedule();
2631 int episodeIdx = simulator.episodeIndex();
2632
2633 // first thing in the morning, limit the time step size to the maximum size
2634 Scalar maxTimeStepSize = Parameters::get<TypeTag, Properties::SolverMaxTimeStepInDays>() * 24 * 60 * 60;
2635 int reportStepIdx = std::max(episodeIdx, 0);
2636 if (this->enableTuning_) {
2637 const auto& tuning = schedule[reportStepIdx].tuning();
2638 maxTimeStepSize = tuning.TSMAXZ;
2639 }
2640
2641 dtNext = std::min(dtNext, maxTimeStepSize);
2642
2643 Scalar remainingEpisodeTime =
2644 simulator.episodeStartTime() + simulator.episodeLength()
2645 - (simulator.startTime() + simulator.time());
2646 assert(remainingEpisodeTime >= 0.0);
2647
2648 // if we would have a small amount of time left over in the current episode, make
2649 // two equal time steps instead of a big and a small one
2650 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
2651 // note: limiting to the maximum time step size here is probably not strictly
2652 // necessary, but it should not hurt and is more fool-proof
2653 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
2654
2655 if (simulator.episodeStarts()) {
2656 // if a well event occurred, respect the limit for the maximum time step after
2657 // that, too
2658 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
2659 bool wellEventOccured =
2660 events.hasEvent(ScheduleEvents::NEW_WELL)
2661 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
2662 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
2663 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
2664 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
2665 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
2666 }
2667 }
2668
2669 return dtNext;
2670 }
2671
2672 void computeAndSetEqWeights_()
2673 {
2674 std::vector<Scalar> sumInvB(numPhases, 0.0);
2675 const auto& gridView = this->gridView();
2676 ElementContext elemCtx(this->simulator());
2677 for(const auto& elem: elements(gridView, Dune::Partitions::interior)) {
2678 elemCtx.updatePrimaryStencil(elem);
2679 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
2680 const auto& dofFluidState = initialFluidStates_[elemIdx];
2681 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2682 if (!FluidSystem::phaseIsActive(phaseIdx))
2683 continue;
2684
2685 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
2686 }
2687 }
2688
2689 std::size_t numDof = this->model().numGridDof();
2690 const auto& comm = this->simulator().vanguard().grid().comm();
2691 comm.sum(sumInvB.data(),sumInvB.size());
2692 Scalar numTotalDof = comm.sum(numDof);
2693
2694 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2695 if (!FluidSystem::phaseIsActive(phaseIdx))
2696 continue;
2697
2698 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
2699 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
2700 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
2701 this->model().setEqWeight(activeSolventCompIdx, avgB);
2702 }
2703 }
2704
2705 int refPressurePhaseIdx_() const {
2706 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2707 return oilPhaseIdx;
2708 }
2709 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2710 return gasPhaseIdx;
2711 }
2712 else {
2713 return waterPhaseIdx;
2714 }
2715 }
2716
2717 void updateRockCompTransMultVal_()
2718 {
2719 const auto& model = this->simulator().model();
2720 std::size_t numGridDof = this->model().numGridDof();
2721 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
2722 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
2723 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, /*timeIdx=*/ 0);
2724 Scalar trans_mult = computeRockCompTransMultiplier_<Scalar>(iq, elementIdx);
2725 this->rockCompTransMultVal_[elementIdx] = trans_mult;
2726 }
2727 }
2728
2734 template <class LhsEval>
2735 LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities& intQuants, unsigned elementIdx) const
2736 {
2737 OPM_TIMEBLOCK_LOCAL(computeRockCompTransMultiplier);
2738 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
2739 return 1.0;
2740
2741 unsigned tableIdx = 0;
2742 if (!this->rockTableIdx_.empty())
2743 tableIdx = this->rockTableIdx_[elementIdx];
2744
2745 const auto& fs = intQuants.fluidState();
2746 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
2747
2748 if (!this->minRefPressure_.empty())
2749 // The pore space change is irreversible
2750 effectivePressure =
2751 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
2752 this->minRefPressure_[elementIdx]);
2753
2754 if (!this->overburdenPressure_.empty())
2755 effectivePressure -= this->overburdenPressure_[elementIdx];
2756
2757 if (!this->rockCompTransMult_.empty())
2758 return this->rockCompTransMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
2759
2760 // water compaction
2761 assert(!this->rockCompTransMultWc_.empty());
2762 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
2763 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
2764
2765 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
2766 }
2767
2768 typename Vanguard::TransmissibilityType transmissibilities_;
2769
2770 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
2771 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
2772
2773 FlowThresholdPressure<TypeTag> thresholdPressures_;
2774
2775 std::vector<InitialFluidState> initialFluidStates_;
2776
2777 bool enableDriftCompensation_;
2778 GlobalEqVector drift_;
2779
2780 WellModel wellModel_;
2781 AquiferModel aquiferModel_;
2782
2783 bool enableEclOutput_;
2784 std::unique_ptr<EclWriterType> eclWriter_;
2785
2786#if HAVE_DAMARIS
2787 bool enableDamarisOutput_ = false ;
2788 std::unique_ptr<DamarisWriterType> damarisWriter_;
2789#endif
2790
2791 PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
2792 TracerModel tracerModel_;
2793
2794 ActionHandler actionHandler_;
2795
2796 template<class T>
2797 struct BCData
2798 {
2799 std::array<std::vector<T>,6> data;
2800
2801 void resize(std::size_t size, T defVal)
2802 {
2803 for (auto& d : data)
2804 d.resize(size, defVal);
2805 }
2806
2807 const std::vector<T>& operator()(FaceDir::DirEnum dir) const
2808 {
2809 if (dir == FaceDir::DirEnum::Unknown)
2810 throw std::runtime_error("Tried to access BC data for the 'Unknown' direction");
2811 int idx = 0;
2812 int div = static_cast<int>(dir);
2813 while ((div /= 2) >= 1)
2814 ++idx;
2815 assert(idx >= 0 && idx <= 5);
2816 return data[idx];
2817 }
2818
2819 std::vector<T>& operator()(FaceDir::DirEnum dir)
2820 {
2821 return const_cast<std::vector<T>&>(std::as_const(*this)(dir));
2822 }
2823 };
2824
2825 BCData<int> bcindex_;
2826 bool nonTrivialBoundaryConditions_ = false;
2827};
2828
2829} // namespace Opm
2830
2831#endif // OPM_FLOW_PROBLEM_HPP
The base class which specifies the API of aquifer models.
Helper class for grid instantiation of ECL file-format using problems.
Collects necessary output values and pass them to Damaris server processes.
This is a "dummy" gradient calculator which does not do anything.
Collects necessary output values and pass it to opm-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...
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
This class calculates the threshold pressure for grid faces according to the Eclipse Reference Manual...
This file contains the flux module which is used for flow problems.
Output module for the results black oil model writing in ECL binary format.
A class which handles tracers as specified in by ECL.
VTK output module for the tracer model's parameters.
Collects necessary output values and pass them to Damaris server processes.
Definition DamarisWriter.hpp:88
Collects necessary output values and pass it to opm-output.
Definition EclWriter.hpp:100
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowGenericProblem.hpp:70
Scalar rockReferencePressure(unsigned globalSpaceIdx) const
Direct access to rock reference pressure.
Definition FlowGenericProblem_impl.hpp:324
Scalar porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
Direct indexed access to the porosity.
Definition FlowGenericProblem_impl.hpp:339
Scalar rockCompressibility(unsigned globalSpaceIdx) const
Direct access to rock compressibility.
Definition FlowGenericProblem_impl.hpp:309
static std::string helpPreamble(int, const char **argv)
Definition FlowGenericProblem_impl.hpp:118
bool shouldWriteRestartFile() const
Returns true if an eWoms restart file should be written to disk.
Definition FlowGenericProblem.hpp:310
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:112
const WellModel & wellModel() const
Returns a reference to the ECL well manager used by the problem.
Definition FlowProblem.hpp:1503
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
Direct access to the transmissibility between two elements.
Definition FlowProblem.hpp:805
Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
Definition FlowProblem.hpp:924
const DimMatrix & intrinsicPermeability(unsigned globalElemIdx) const
This method returns the intrinsic permeability tensor given a global element index.
Definition FlowProblem.hpp:787
void endEpisode()
Called by the simulator after the end of an episode.
Definition FlowProblem.hpp:711
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1093
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 FlowProblem.hpp:1259
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:948
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition FlowProblem.hpp:609
unsigned satnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1101
Scalar thermalHalfTransmissibilityOut(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:895
Scalar transmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:852
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn, const unsigned globalSpaceIdxOut) const
Definition FlowProblem.hpp:872
void endTimeStep()
Called by the simulator after each time integration.
Definition FlowProblem.hpp:629
Scalar rockCompressibility(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:982
bool recycleFirstIterationStorage() const
Return if the storage term of the first iteration is identical to the storage term for the solution o...
Definition FlowProblem.hpp:1274
Scalar maxOilSaturation(unsigned globalDofIdx) const
Returns an element's historic maximum oil phase saturation that was observed during the simulation.
Definition FlowProblem.hpp:1219
std::string name() const
Definition FlowProblem.hpp:1131
LhsEval rockCompTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1709
void endIteration()
Called by the simulator after each Newton-Raphson iteration.
Definition FlowProblem.hpp:619
FlowProblem(Simulator &simulator)
Definition FlowProblem.hpp:270
const Vanguard::TransmissibilityType & eclTransmissibilities() const
Return a reference to the object that handles the "raw" transmissibilities.
Definition FlowProblem.hpp:918
Scalar nextTimeStepSize() const
Propose the size of the next time step to the simulator.
Definition FlowProblem.hpp:1642
void writeOutput(const SimulatorTimer &timer, bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:735
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1037
Scalar transmissibility(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:794
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
Direct indexed acces to the depth of an degree of freedom [m].
Definition FlowProblem.hpp:973
LhsEval wellTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Return the well transmissibility multiplier due to rock changues.
Definition FlowProblem.hpp:1740
std::shared_ptr< const EclMaterialLawManager > materialLawManager() const
Returns the ECL material law manager.
Definition FlowProblem.hpp:1049
unsigned plmixnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1117
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1138
Scalar thermalHalfTransmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:908
unsigned miscnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1109
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 FlowProblem.hpp:1248
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:773
Scalar maxPolymerAdsorption(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the max polymer adsorption value.
Definition FlowProblem.hpp:1125
void beginTimeStep()
Called by the simulator before each time integration.
Definition FlowProblem.hpp:579
Scalar thermalHalfTransmissibilityIn(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:882
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
Sets an element's maximum oil phase saturation observed during the simulation.
Definition FlowProblem.hpp:1236
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:1024
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the dispersivity for a face i.e.
Definition FlowProblem.hpp:832
void finishInit()
Definition FlowProblem.hpp:340
void initialSolutionApplied()
Definition FlowProblem.hpp:1330
static void registerParameters()
Definition FlowProblem.hpp:209
static int handlePositionalParameter(std::set< std::string > &seenParams, std::string &errorMsg, int, const char **argv, int paramIdx, int)
Definition FlowProblem.hpp:251
void beginEpisode()
Called by the simulator before an episode begins.
Definition FlowProblem.hpp:510
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1368
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1002
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:839
Scalar diffusivity(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:814
void deserialize(Restarter &res)
This method restores the complete state of the problem and its sub-objects from disk.
Definition FlowProblem.hpp:476
std::shared_ptr< EclMaterialLawManager > materialLawManager()
Definition FlowProblem.hpp:1085
Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the depth of an degree of freedom [m].
Definition FlowProblem.hpp:961
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1290
LhsEval permFactTransMultiplier(const IntensiveQuantities &intQuants) const
Calculate the transmissibility multiplier due to porosity reduction.
Definition FlowProblem.hpp:1722
LhsEval rockCompPoroMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the porosity multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1669
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:862
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1173
void serialize(Restarter &res)
This method writes the complete state of the problem and its subobjects to disk.
Definition FlowProblem.hpp:495
Scalar rockReferencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:992
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the transmissibility for a face i.e.
Definition FlowProblem.hpp:825
This class calculates the threshold pressure for grid faces according to the Eclipse Reference Manual...
Definition FlowThresholdPressure.hpp:59
This class is intend to be a relperm diagnostics, to detect wrong input of relperm table and endpoint...
Definition RelpermDiagnostics.hpp:50
void diagnosis(const EclipseState &eclState, const CartesianIndexMapper &cartesianIndexMapper)
This function is used to diagnosis relperm in eclipse data file.
Definition RelpermDiagnostics.cpp:824
Definition SimulatorTimer.hpp:39
A class which handles tracers as specified in by ECL.
Definition TracerModel.hpp:67
VTK output module for the tracer model's parameters.
Definition VtkTracerModule.hpp:70
static void registerParameters()
Register all run-time parameters for the tracer VTK output module.
Definition VtkTracerModule.hpp:96
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27