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 explicitRockCompaction_ = Parameters::get<TypeTag, Properties::ExplicitRockCompaction>();
334
335
336 RelpermDiagnostics relpermDiagnostics;
337 relpermDiagnostics.diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
338 }
339
344 {
345 ParentType::finishInit();
346
347 auto& simulator = this->simulator();
348 const auto& eclState = simulator.vanguard().eclState();
349 const auto& schedule = simulator.vanguard().schedule();
350
351 // Set the start time of the simulation
352 simulator.setStartTime(schedule.getStartTime());
353 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
354
355 // We want the episode index to be the same as the report step index to make
356 // things simpler, so we have to set the episode index to -1 because it is
357 // incremented by endEpisode(). The size of the initial time step and
358 // length of the initial episode is set to zero for the same reason.
359 simulator.setEpisodeIndex(-1);
360 simulator.setEpisodeLength(0.0);
361
362 // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
363 // disables gravity, else the standard value of the gravity constant at sea level
364 // on earth is used
365 this->gravity_ = 0.0;
366 if (Parameters::get<TypeTag, Properties::EnableGravity>())
367 this->gravity_[dim - 1] = 9.80665;
368 if (!eclState.getInitConfig().hasGravity())
369 this->gravity_[dim - 1] = 0.0;
370
371 if (this->enableTuning_) {
372 // if support for the TUNING keyword is enabled, we get the initial time
373 // steping parameters from it instead of from command line parameters
374 const auto& tuning = schedule[0].tuning();
375 this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
376 this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
377 }
378
379 this->initFluidSystem_();
380
381 // deal with DRSDT
382 this->mixControls_.init(this->model().numGridDof(),
383 this->episodeIndex(),
384 eclState.runspec().tabdims().getNumPVTTables());
385
386 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
387 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
388 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
389 }
390
391 this->readRockParameters_(simulator.vanguard().cellCenterDepths(),
392 [&simulator](const unsigned idx)
393 {
394 std::array<int,dim> coords;
395 simulator.vanguard().cartesianCoordinate(idx, coords);
396 for (auto& c : coords) {
397 ++c;
398 }
399 return coords;
400 });
401 readMaterialParameters_();
402 readThermalParameters_();
403
404 // Re-ordering in case of ALUGrid
405 std::function<unsigned int(unsigned int)> gridToEquilGrid = [&simulator](unsigned int i) {
406 return simulator.vanguard().gridIdxToEquilGridIdx(i);
407 };
408 transmissibilities_.finishInit(gridToEquilGrid);
409
410 const auto& initconfig = eclState.getInitConfig();
411 tracerModel_.init(initconfig.restartRequested());
412 if (initconfig.restartRequested())
413 readEclRestartSolution_();
414 else
415 readInitialCondition_();
416
417 tracerModel_.prepareTracerBatches();
418
419 updatePffDofData_();
420
421 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
422 const auto& vanguard = this->simulator().vanguard();
423 const auto& gridView = vanguard.gridView();
424 int numElements = gridView.size(/*codim=*/0);
425 this->polymer_.maxAdsorption.resize(numElements, 0.0);
426 }
427
428 readBoundaryConditions_();
429
430 // compute and set eq weights based on initial b values
431 computeAndSetEqWeights_();
432
433 if (enableDriftCompensation_) {
434 drift_.resize(this->model().numGridDof());
435 drift_ = 0.0;
436 }
437
438 // write the static output files (EGRID, INIT, SMSPEC, etc.)
439 if (enableEclOutput_) {
440 if (simulator.vanguard().grid().comm().size() > 1) {
441 if (simulator.vanguard().grid().comm().rank() == 0)
442 eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
443 } else
444 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
445
446 // Re-ordering in case of ALUGrid
447 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
448 return simulator.vanguard().gridEquilIdxToGridIdx(i);
449 };
450 eclWriter_->writeInit(equilGridToGrid);
451 }
452
453 simulator.vanguard().releaseGlobalTransmissibilities();
454
455 // after finishing the initialization and writing the initial solution, we move
456 // to the first "real" episode/report step
457 // for restart the episode index and start is already set
458 if (!initconfig.restartRequested()) {
459 simulator.startNextEpisode(schedule.seconds(1));
460 simulator.setEpisodeIndex(0);
461 }
462 }
463
464 void prefetch(const Element& elem) const
465 { pffDofData_.prefetch(elem); }
466
478 template <class Restarter>
479 void deserialize(Restarter& res)
480 {
481 // reload the current episode/report step from the deck
482 beginEpisode();
483
484 // deserialize the wells
485 wellModel_.deserialize(res);
486
487 // deserialize the aquifer
488 aquiferModel_.deserialize(res);
489 }
490
497 template <class Restarter>
498 void serialize(Restarter& res)
499 {
500 wellModel_.serialize(res);
501
502 aquiferModel_.serialize(res);
503 }
504
505 int episodeIndex() const
506 {
507 return std::max(this->simulator().episodeIndex(), 0);
508 }
509
514 {
515 OPM_TIMEBLOCK(beginEpisode);
516 // Proceed to the next report step
517 auto& simulator = this->simulator();
518 int episodeIdx = simulator.episodeIndex();
519 auto& eclState = simulator.vanguard().eclState();
520 const auto& schedule = simulator.vanguard().schedule();
521 const auto& events = schedule[episodeIdx].events();
522
523 if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
524 // bring the contents of the keywords to the current state of the SCHEDULE
525 // section.
526 //
527 // TODO (?): make grid topology changes possible (depending on what exactly
528 // has changed, the grid may need be re-created which has some serious
529 // implications on e.g., the solution of the simulation.)
530 const auto& miniDeck = schedule[episodeIdx].geo_keywords();
531 const auto& cc = simulator.vanguard().grid().comm();
532 eclState.apply_schedule_keywords( miniDeck );
533 eclBroadcast(cc, eclState.getTransMult() );
534
535 // Re-ordering in case of ALUGrid
536 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
537 return simulator.vanguard().gridEquilIdxToGridIdx(i);
538 };
539
540 // re-compute all quantities which may possibly be affected.
541 transmissibilities_.update(true, equilGridToGrid);
542 this->referencePorosity_[1] = this->referencePorosity_[0];
543 updateReferencePorosity_();
544 updatePffDofData_();
545 this->model().linearizer().updateDiscretizationParameters();
546 }
547
548 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
549
550 // set up the wells for the next episode.
551 wellModel_.beginEpisode();
552
553 // set up the aquifers for the next episode.
554 aquiferModel_.beginEpisode();
555
556 // set the size of the initial time step of the episode
557 Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
558 // negative value of initialTimeStepSize_ indicates no active limit from TSINIT or NEXTSTEP
559 if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
560 // allow the size of the initial time step to be set via an external parameter
561 // if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
562 dt = std::min(dt, this->initialTimeStepSize_);
563 simulator.setTimeStepSize(dt);
564
565 // Evaluate UDQ assign statements to make sure the settings are
566 // available as UDA controls for the current report step.
567 actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
568
569 if (episodeIdx >= 0) {
570 const auto& oilVap = schedule[episodeIdx].oilvap();
571 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
572 FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2());
573 } else {
574 FluidSystem::setVapPars(0.0, 0.0);
575 }
576 }
577 }
578
583 {
584 OPM_TIMEBLOCK(beginTimeStep);
585 int episodeIdx = this->episodeIndex();
586
587 this->beginTimeStep_(enableExperiments,
588 episodeIdx,
589 this->simulator().timeStepIndex(),
590 this->simulator().startTime(),
591 this->simulator().time(),
592 this->simulator().timeStepSize(),
593 this->simulator().endTime());
594
595 // update maximum water saturation and minimum pressure
596 // used when ROCKCOMP is activated
597 asImp_().updateExplicitQuantities_();
598
599 if (nonTrivialBoundaryConditions()) {
600 this->model().linearizer().updateBoundaryConditionData();
601 }
602
603 wellModel_.beginTimeStep();
604 aquiferModel_.beginTimeStep();
605 tracerModel_.beginTimeStep();
606
607 }
608
613 {
614 OPM_TIMEBLOCK(beginIteration);
615 wellModel_.beginIteration();
616 aquiferModel_.beginIteration();
617 }
618
623 {
624 OPM_TIMEBLOCK(endIteration);
625 wellModel_.endIteration();
626 aquiferModel_.endIteration();
627 }
628
633 {
634 OPM_TIMEBLOCK(endTimeStep);
635#ifndef NDEBUG
636 if constexpr (getPropValue<TypeTag, Properties::EnableDebuggingChecks>()) {
637 // in debug mode, we don't care about performance, so we check if the model does
638 // the right thing (i.e., the mass change inside the whole reservoir must be
639 // equivalent to the fluxes over the grid's boundaries plus the source rates
640 // specified by the problem)
641 int rank = this->simulator().gridView().comm().rank();
642 if (rank == 0)
643 std::cout << "checking conservativeness of solution\n";
644 this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
645 if (rank == 0)
646 std::cout << "solution is sufficiently conservative\n";
647 }
648#endif // NDEBUG
649
650 auto& simulator = this->simulator();
651 wellModel_.endTimeStep();
652 aquiferModel_.endTimeStep();
653 tracerModel_.endTimeStep();
654
655
656 // Compute flux for output
657 this->model().linearizer().updateFlowsInfo();
658
659 // deal with DRSDT and DRVDT
660 asImp_().updateCompositionChangeLimits_();
661 {
662 OPM_TIMEBLOCK(driftCompansation);
663 if (enableDriftCompensation_) {
664 const auto& residual = this->model().linearizer().residual();
665 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
666 drift_[globalDofIdx] = residual[globalDofIdx];
667 drift_[globalDofIdx] *= simulator.timeStepSize();
668 if constexpr (getPropValue<TypeTag, Properties::UseVolumetricResidual>())
669 drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
670 }
671 }
672 }
673 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
674 !this->simulator().episodeWillBeOver();
675 // For CpGrid with LGRs, ecl/vtk output is not supported yet.
676 const auto& grid = this->simulator().vanguard().gridView().grid();
677 using GridType = std::remove_cv_t< typename std::remove_reference<decltype(grid)>::type>;
678 bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
679 if ( !isCpGrid || (this->simulator().vanguard().gridView().grid().maxLevel()==0)) {
680 eclWriter_->evalSummaryState(isSubStep);
681 }
682
683 int episodeIdx = this->episodeIndex();
684
685 // Re-ordering in case of Alugrid
686 std::function<unsigned int(unsigned int)> gridToEquilGrid = [&simulator](unsigned int i) {
687 return simulator.vanguard().gridIdxToEquilGridIdx(i);
688 };
689
690 std::function<void(bool)> transUp =
691 [this,gridToEquilGrid](bool global) {
692 this->transmissibilities_.update(global,gridToEquilGrid);
693 };
694 {
695 OPM_TIMEBLOCK(applyActions);
696 actionHandler_.applyActions(episodeIdx,
697 simulator.time() + simulator.timeStepSize(),
698 transUp);
699 }
700 // deal with "clogging" for the MICP model
701 if constexpr (enableMICP){
702 auto& model = this->model();
703 const auto& residual = this->model().linearizer().residual();
704 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
705 auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx];
706 MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
707 }
708 }
709 }
710
715 {
716 OPM_TIMEBLOCK(endEpisode);
717 auto& simulator = this->simulator();
718 auto& schedule = simulator.vanguard().schedule();
719
720 wellModel_.endEpisode();
721 aquiferModel_.endEpisode();
722
723 int episodeIdx = this->episodeIndex();
724 // check if we're finished ...
725 if (episodeIdx + 1 >= static_cast<int>(schedule.size() - 1)) {
726 simulator.setFinished(true);
727 return;
728 }
729
730 // .. if we're not yet done, start the next episode (report step)
731 simulator.startNextEpisode(schedule.stepLength(episodeIdx + 1));
732 }
733
738 void writeOutput(const SimulatorTimer& timer, bool verbose = true)
739 {
740 OPM_TIMEBLOCK(problemWriteOutput);
741 // use the generic code to prepare the output fields and to
742 // write the desired VTK files.
743 if (Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() ||
744 this->simulator().episodeWillBeOver()) {
745 ParentType::writeOutput(verbose);
746 }
747
748 bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
749 !this->simulator().episodeWillBeOver();
750
751 data::Solution localCellData = {};
752#if HAVE_DAMARIS
753 // N.B. the Damaris output has to be done before the ECL output as the ECL one
754 // does all kinds of std::move() relocation of data
755 if (enableDamarisOutput_) {
756 damarisWriter_->writeOutput(localCellData, isSubStep) ;
757 }
758#endif
759 if (enableEclOutput_){
760 eclWriter_->writeOutput(std::move(localCellData), timer, isSubStep);
761 }
762 }
763
764 void finalizeOutput() {
765 OPM_TIMEBLOCK(finalizeOutput);
766 // this will write all pending output to disk
767 // to avoid corruption of output files
768 eclWriter_.reset();
769 }
770
771
775 template <class Context>
776 const DimMatrix& intrinsicPermeability(const Context& context,
777 unsigned spaceIdx,
778 unsigned timeIdx) const
779 {
780 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
781 return transmissibilities_.permeability(globalSpaceIdx);
782 }
783
790 const DimMatrix& intrinsicPermeability(unsigned globalElemIdx) const
791 { return transmissibilities_.permeability(globalElemIdx); }
792
796 template <class Context>
797 Scalar transmissibility(const Context& context,
798 [[maybe_unused]] unsigned fromDofLocalIdx,
799 unsigned toDofLocalIdx) const
800 {
801 assert(fromDofLocalIdx == 0);
802 return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
803 }
804
808 Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
809 {
810 return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
811 }
812
816 template <class Context>
817 Scalar diffusivity(const Context& context,
818 [[maybe_unused]] unsigned fromDofLocalIdx,
819 unsigned toDofLocalIdx) const
820 {
821 assert(fromDofLocalIdx == 0);
822 return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
823 }
824
828 Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
829 return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
830 }
831
835 Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
836 return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
837 }
838
842 Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx,
843 const unsigned boundaryFaceIdx) const
844 {
845 return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
846 }
847
848
849
850
854 template <class Context>
855 Scalar transmissibilityBoundary(const Context& elemCtx,
856 unsigned boundaryFaceIdx) const
857 {
858 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
859 return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
860 }
861
865 Scalar transmissibilityBoundary(const unsigned globalSpaceIdx,
866 const unsigned boundaryFaceIdx) const
867 {
868 return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
869 }
870
871
875 Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn,
876 const unsigned globalSpaceIdxOut) const
877 {
878 return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
879 }
880
884 template <class Context>
885 Scalar thermalHalfTransmissibilityIn(const Context& context,
886 unsigned faceIdx,
887 unsigned timeIdx) const
888 {
889 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
890 unsigned toDofLocalIdx = face.exteriorIndex();
891 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
892 }
893
897 template <class Context>
898 Scalar thermalHalfTransmissibilityOut(const Context& context,
899 unsigned faceIdx,
900 unsigned timeIdx) const
901 {
902 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
903 unsigned toDofLocalIdx = face.exteriorIndex();
904 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
905 }
906
910 template <class Context>
911 Scalar thermalHalfTransmissibilityBoundary(const Context& elemCtx,
912 unsigned boundaryFaceIdx) const
913 {
914 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
915 return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
916 }
917
921 const typename Vanguard::TransmissibilityType& eclTransmissibilities() const
922 { return transmissibilities_; }
923
927 Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
928 { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
929
930 const FlowThresholdPressure<TypeTag>& thresholdPressure() const
931 { return thresholdPressures_; }
932
933 FlowThresholdPressure<TypeTag>& thresholdPressure()
934 { return thresholdPressures_; }
935
936 const TracerModel& tracerModel() const
937 { return tracerModel_; }
938
939 TracerModel& tracerModel()
940 { return tracerModel_; }
941
950 template <class Context>
951 Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
952 {
953 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
954 return this->porosity(globalSpaceIdx, timeIdx);
955 }
956
963 template <class Context>
964 Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
965 {
966 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
967 return this->dofCenterDepth(globalSpaceIdx);
968 }
969
976 Scalar dofCenterDepth(unsigned globalSpaceIdx) const
977 {
978 return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
979 }
980
984 template <class Context>
985 Scalar rockCompressibility(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
986 {
987 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
988 return this->rockCompressibility(globalSpaceIdx);
989 }
990
994 template <class Context>
995 Scalar rockReferencePressure(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
996 {
997 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
998 return this->rockReferencePressure(globalSpaceIdx);
999 }
1000
1004 template <class Context>
1005 const MaterialLawParams& materialLawParams(const Context& context,
1006 unsigned spaceIdx, unsigned timeIdx) const
1007 {
1008 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1009 return this->materialLawParams(globalSpaceIdx);
1010 }
1011
1012 const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const
1013 {
1014 return materialLawManager_->materialLawParams(globalDofIdx);
1015 }
1016
1017 const MaterialLawParams& materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const
1018 {
1019 return materialLawManager_->materialLawParams(globalDofIdx, facedir);
1020 }
1021
1025 template <class Context>
1026 const SolidEnergyLawParams&
1027 solidEnergyLawParams(const Context& context,
1028 unsigned spaceIdx,
1029 unsigned timeIdx) const
1030 {
1031 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1032 return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1033 }
1034
1038 template <class Context>
1039 const ThermalConductionLawParams &
1040 thermalConductionLawParams(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1041 {
1042 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1043 return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1044 }
1045
1052 std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
1053 { return materialLawManager_; }
1054
1055 template <class FluidState>
1056 void updateRelperms(
1057 std::array<Evaluation,numPhases> &mobility,
1058 DirectionalMobilityPtr &dirMob,
1059 FluidState &fluidState,
1060 unsigned globalSpaceIdx) const
1061 {
1062 OPM_TIMEBLOCK_LOCAL(updateRelperms);
1063 {
1064 // calculate relative permeabilities. note that we store the result into the
1065 // mobility_ class attribute. the division by the phase viscosity happens later.
1066 const auto& materialParams = materialLawParams(globalSpaceIdx);
1067 MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
1068 Valgrind::CheckDefined(mobility);
1069 }
1070 if (materialLawManager_->hasDirectionalRelperms()
1071 || materialLawManager_->hasDirectionalImbnum())
1072 {
1073 using Dir = FaceDir::DirEnum;
1074 constexpr int ndim = 3;
1075 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
1076 Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
1077 for (int i = 0; i<ndim; i++) {
1078 const auto& materialParams = materialLawParams(globalSpaceIdx, facedirs[i]);
1079 auto& mob_array = dirMob->getArray(i);
1080 MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
1081 }
1082 }
1083 }
1084
1088 std::shared_ptr<EclMaterialLawManager> materialLawManager()
1089 { return materialLawManager_; }
1090
1091 using BaseType::pvtRegionIndex;
1095 template <class Context>
1096 unsigned pvtRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1097 { return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1098
1099 using BaseType::satnumRegionIndex;
1103 template <class Context>
1104 unsigned satnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1105 { return this->satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1106
1107 using BaseType::miscnumRegionIndex;
1111 template <class Context>
1112 unsigned miscnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1113 { return this->miscnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1114
1115 using BaseType::plmixnumRegionIndex;
1119 template <class Context>
1120 unsigned plmixnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1121 { return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1122
1123 using BaseType::maxPolymerAdsorption;
1127 template <class Context>
1128 Scalar maxPolymerAdsorption(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1129 { return this->maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); }
1130
1134 std::string name() const
1135 { return this->simulator().vanguard().caseName(); }
1136
1140 template <class Context>
1141 Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1142 {
1143 // use the initial temperature of the DOF if temperature is not a primary
1144 // variable
1145 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1146 return initialFluidStates_[globalDofIdx].temperature(/*phaseIdx=*/0);
1147 }
1148
1149
1150 Scalar temperature(unsigned globalDofIdx, unsigned /*timeIdx*/) const
1151 {
1152 // use the initial temperature of the DOF if temperature is not a primary
1153 // variable
1154 return initialFluidStates_[globalDofIdx].temperature(/*phaseIdx=*/0);
1155 }
1156
1157 const SolidEnergyLawParams&
1158 solidEnergyLawParams(unsigned globalSpaceIdx,
1159 unsigned /*timeIdx*/) const
1160 {
1161 return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
1162 }
1163 const ThermalConductionLawParams &
1164 thermalConductionLawParams(unsigned globalSpaceIdx,
1165 unsigned /*timeIdx*/)const
1166 {
1167 return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
1168 }
1169
1175 template <class Context>
1176 void boundary(BoundaryRateVector& values,
1177 const Context& context,
1178 unsigned spaceIdx,
1179 unsigned timeIdx) const
1180 {
1181 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
1182 if (!context.intersection(spaceIdx).boundary())
1183 return;
1184
1185 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
1186 values.setNoFlow();
1187 else {
1188 // in the energy case we need to specify a non-trivial boundary condition
1189 // because the geothermal gradient needs to be maintained. for this, we
1190 // simply assume the initial temperature at the boundary and specify the
1191 // thermal flow accordingly. in this context, "thermal flow" means energy
1192 // flow due to a temerature gradient while assuming no-flow for mass
1193 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1194 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1195 values.setThermalFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]);
1196 }
1197
1198 if (nonTrivialBoundaryConditions()) {
1199 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
1200 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
1201 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
1202 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
1203 const auto [type, massrate] = boundaryCondition(globalDofIdx, indexInInside);
1204 if (type == BCType::THERMAL)
1205 values.setThermalFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1206 else if (type == BCType::FREE || type == BCType::DIRICHLET)
1207 values.setFreeFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside));
1208 else if (type == BCType::RATE)
1209 values.setMassRate(massrate, pvtRegionIdx);
1210 }
1211 }
1212
1222 Scalar maxOilSaturation(unsigned globalDofIdx) const
1223 {
1224 if (!this->vapparsActive(this->episodeIndex()))
1225 return 0.0;
1226
1227 return this->maxOilSaturation_[globalDofIdx];
1228 }
1229
1239 void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
1240 {
1241 if (!this->vapparsActive(this->episodeIndex()))
1242 return;
1243
1244 this->maxOilSaturation_[globalDofIdx] = value;
1245 }
1246
1251 Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
1252 {
1253 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
1254 this->episodeIndex(),
1255 this->pvtRegionIndex(globalDofIdx));
1256 }
1257
1262 Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
1263 {
1264 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
1265 this->episodeIndex(),
1266 this->pvtRegionIndex(globalDofIdx));
1267 }
1268
1278 {
1279 int episodeIdx = this->episodeIndex();
1280 return !this->mixControls_.drsdtActive(episodeIdx) &&
1281 !this->mixControls_.drvdtActive(episodeIdx) &&
1282 this->rockCompPoroMultWc_.empty() &&
1283 this->rockCompPoroMult_.empty();
1284 }
1285
1292 template <class Context>
1293 void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
1294 {
1295 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1296
1297 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
1298 values.assignNaive(initialFluidStates_[globalDofIdx]);
1299
1300 SolventModule::assignPrimaryVars(values,
1301 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
1302 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
1303
1304 if constexpr (enablePolymer)
1305 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
1306
1307 if constexpr (enablePolymerMolarWeight)
1308 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
1309
1310 if constexpr (enableBrine) {
1311 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
1312 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
1313 }
1314 else {
1315 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
1316 }
1317 }
1318
1319 if constexpr (enableMICP){
1320 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
1321 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
1322 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
1323 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
1324 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
1325 }
1326
1327 values.checkDefined();
1328 }
1329
1334 {
1335 // Calculate all intensive quantities.
1336 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/0);
1337
1338 // We also need the intensive quantities for timeIdx == 1
1339 // corresponding to the start of the current timestep, if we
1340 // do not use the storage cache, or if we cannot recycle the
1341 // first iteration storage.
1342 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
1343 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/1);
1344 }
1345
1346 // initialize the wells. Note that this needs to be done after initializing the
1347 // intrinsic permeabilities and the after applying the initial solution because
1348 // the well model uses these...
1349 wellModel_.init();
1350
1351 // let the object for threshold pressures initialize itself. this is done only at
1352 // this point, because determining the threshold pressures may require to access
1353 // the initial solution.
1354 thresholdPressures_.finishInit();
1355
1356 updateCompositionChangeLimits_();
1357
1358 aquiferModel_.initialSolutionApplied();
1359
1360 if (this->simulator().episodeIndex() == 0) {
1361 eclWriter_->writeInitialFIPReport();
1362 }
1363 }
1364
1370 template <class Context>
1371 void source(RateVector& rate,
1372 const Context& context,
1373 unsigned spaceIdx,
1374 unsigned timeIdx) const
1375 {
1376 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
1377 source(rate, globalDofIdx, timeIdx);
1378 }
1379
1380 void source(RateVector& rate,
1381 unsigned globalDofIdx,
1382 unsigned timeIdx) const
1383 {
1384 OPM_TIMEBLOCK_LOCAL(eclProblemSource);
1385 rate = 0.0;
1386
1387 // Add well contribution to source here.
1388 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
1389
1390 // convert the source term from the total mass rate of the
1391 // cell to the one per unit of volume as used by the model.
1392 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1393 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
1394
1395 Valgrind::CheckDefined(rate[eqIdx]);
1396 assert(isfinite(rate[eqIdx]));
1397 }
1398
1399 // Add non-well sources.
1400 addToSourceDense(rate, globalDofIdx, timeIdx);
1401 }
1402
1403 void addToSourceDense(RateVector& rate,
1404 unsigned globalDofIdx,
1405 unsigned timeIdx) const
1406 {
1407 aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
1408
1409 // Add source term from deck
1410 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
1411 std::array<int,3> ijk;
1412 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
1413
1414 if (source.hasSource(ijk)) {
1415 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1416 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
1417 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
1418 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
1419
1420 for (unsigned i = 0; i < phidx_map.size(); ++i) {
1421 const auto phaseIdx = phidx_map[i];
1422 const auto sourceComp = sc_map[i];
1423 const auto compIdx = cidx_map[i];
1424 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1425 continue;
1426 }
1427 Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1428 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1429 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
1430 }
1431 rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate;
1432 }
1433
1434 if constexpr (enableSolvent) {
1435 Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
1436 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
1437 const auto& solventPvt = SolventModule::solventPvt();
1438 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
1439 }
1440 rate[Indices::contiSolventEqIdx] += mass_rate;
1441 }
1442 if constexpr (enablePolymer) {
1443 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
1444 }
1445 if constexpr (enableEnergy) {
1446 for (unsigned i = 0; i < phidx_map.size(); ++i) {
1447 const auto phaseIdx = phidx_map[i];
1448 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1449 continue;
1450 }
1451 const auto sourceComp = sc_map[i];
1452 if (source.hasHrate({ijk, sourceComp})) {
1453 rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
1454 } else {
1455 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, /*timeIdx*/ 0);
1456 auto fs = intQuants.fluidState();
1457 // if temperature is not set, use cell temperature as default
1458 if (source.hasTemperature({ijk, sourceComp})) {
1459 Scalar temperature = source.temperature({ijk, sourceComp});
1460 fs.setTemperature(temperature);
1461 }
1462 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
1463 Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx);
1464 Scalar energy_rate = getValue(h)*mass_rate;
1465 rate[Indices::contiEnergyEqIdx] += energy_rate;
1466 }
1467 }
1468 }
1469 }
1470
1471 // if requested, compensate systematic mass loss for cells which were "well
1472 // behaved" in the last time step
1473 // Note that we don't allow for drift compensation if there are no active wells.
1474 const bool compensateDrift = wellModel_.wellsActive();
1475 if (enableDriftCompensation_ && compensateDrift) {
1476 const auto& simulator = this->simulator();
1477 const auto& model = this->model();
1478
1479 // we use a lower tolerance for the compensation too
1480 // assure the added drift from the last step does not
1481 // cause convergence issues on the current step
1482 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
1483 Scalar poro = this->porosity(globalDofIdx, timeIdx);
1484 Scalar dt = simulator.timeStepSize();
1485 EqVector dofDriftRate = drift_[globalDofIdx];
1486 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
1487
1488 // restrict drift compensation to the CNV tolerance
1489 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
1490 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
1491 if (cnv > maxCompensation) {
1492 dofDriftRate[eqIdx] *= maxCompensation/cnv;
1493 }
1494 }
1495
1496 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
1497 rate[eqIdx] -= dofDriftRate[eqIdx];
1498 }
1499 }
1500
1506 const WellModel& wellModel() const
1507 { return wellModel_; }
1508
1509 WellModel& wellModel()
1510 { return wellModel_; }
1511
1512 const AquiferModel& aquiferModel() const
1513 { return aquiferModel_; }
1514
1515 AquiferModel& mutableAquiferModel()
1516 { return aquiferModel_; }
1517
1518 // temporary solution to facilitate output of initial state from flow
1519 const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
1520 { return initialFluidStates_[globalDofIdx]; }
1521
1522 const EclipseIO& eclIO() const
1523 { return eclWriter_->eclIO(); }
1524
1525 void setSubStepReport(const SimulatorReportSingle& report)
1526 { return eclWriter_->setSubStepReport(report); }
1527
1528 void setSimulationReport(const SimulatorReport& report)
1529 { return eclWriter_->setSimulationReport(report); }
1530
1531 bool nonTrivialBoundaryConditions() const
1532 { return nonTrivialBoundaryConditions_; }
1533
1534 const InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const
1535 {
1536 OPM_TIMEBLOCK_LOCAL(boundaryFluidState);
1537 const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop;
1538 if (bcprop.size() > 0) {
1539 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1540
1541 // index == 0: no boundary conditions for this
1542 // global cell and direction
1543 if (bcindex_(dir)[globalDofIdx] == 0)
1544 return initialFluidStates_[globalDofIdx];
1545
1546 const auto& bc = bcprop[bcindex_(dir)[globalDofIdx]];
1547 if (bc.bctype == BCType::DIRICHLET )
1548 {
1549 InitialFluidState fluidState;
1550 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
1551 fluidState.setPvtRegionIndex(pvtRegionIdx);
1552
1553 switch (bc.component) {
1554 case BCComponent::OIL:
1555 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
1556 throw std::logic_error("oil is not active and you're trying to add oil BC");
1557
1558 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
1559 break;
1560 case BCComponent::GAS:
1561 if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1562 throw std::logic_error("gas is not active and you're trying to add gas BC");
1563
1564 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
1565 break;
1566 case BCComponent::WATER:
1567 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1568 throw std::logic_error("water is not active and you're trying to add water BC");
1569
1570 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
1571 break;
1572 case BCComponent::SOLVENT:
1573 case BCComponent::POLYMER:
1574 case BCComponent::NONE:
1575 throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
1576 break;
1577 }
1578 double pressure = initialFluidStates_[globalDofIdx].pressure(refPressurePhaseIdx_());
1579 const auto pressure_input = bc.pressure;
1580 if (pressure_input) {
1581 pressure = *pressure_input;
1582 }
1583
1584 std::array<Scalar, numPhases> pc = {0};
1585 const auto& matParams = materialLawParams(globalDofIdx);
1586 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
1587 Valgrind::CheckDefined(pressure);
1588 Valgrind::CheckDefined(pc);
1589 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1590 if (!FluidSystem::phaseIsActive(phaseIdx))
1591 continue;
1592
1593 if (Indices::oilEnabled)
1594 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1595 else if (Indices::gasEnabled)
1596 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1597 else if (Indices::waterEnabled)
1598 //single (water) phase
1599 fluidState.setPressure(phaseIdx, pressure);
1600 }
1601
1602 double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature
1603 const auto temperature_input = bc.temperature;
1604 if(temperature_input)
1605 temperature = *temperature_input;
1606 fluidState.setTemperature(temperature);
1607
1608 if (FluidSystem::enableDissolvedGas()) {
1609 fluidState.setRs(0.0);
1610 fluidState.setRv(0.0);
1611 }
1612 if (FluidSystem::enableDissolvedGasInWater()) {
1613 fluidState.setRsw(0.0);
1614 }
1615 if (FluidSystem::enableVaporizedWater())
1616 fluidState.setRvw(0.0);
1617
1618 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1619 if (!FluidSystem::phaseIsActive(phaseIdx))
1620 continue;
1621
1622 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
1623 fluidState.setInvB(phaseIdx, b);
1624
1625 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
1626 fluidState.setDensity(phaseIdx, rho);
1627 if (enableEnergy) {
1628 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
1629 fluidState.setEnthalpy(phaseIdx, h);
1630 }
1631 }
1632 fluidState.checkDefined();
1633 return fluidState;
1634 }
1635 }
1636 return initialFluidStates_[globalDofIdx];
1637 }
1638
1645 Scalar nextTimeStepSize() const
1646 {
1647 OPM_TIMEBLOCK(nexTimeStepSize);
1648 // allow external code to do the timestepping
1649 if (this->nextTimeStepSize_ > 0.0)
1650 return this->nextTimeStepSize_;
1651
1652 const auto& simulator = this->simulator();
1653 int episodeIdx = simulator.episodeIndex();
1654
1655 // for the initial episode, we use a fixed time step size
1656 if (episodeIdx < 0)
1657 return this->initialTimeStepSize_;
1658
1659 // ask the newton method for a suggestion. This suggestion will be based on how
1660 // well the previous time step converged. After that, apply the runtime time
1661 // stepping constraints.
1662 const auto& newtonMethod = this->model().newtonMethod();
1663 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1664 }
1665
1671 template <class LhsEval>
1672 LhsEval rockCompPoroMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1673 {
1674 OPM_TIMEBLOCK_LOCAL(rockCompPoroMultiplier);
1675 if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
1676 return 1.0;
1677
1678 unsigned tableIdx = 0;
1679 if (!this->rockTableIdx_.empty())
1680 tableIdx = this->rockTableIdx_[elementIdx];
1681
1682 const auto& fs = intQuants.fluidState();
1683 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1684 if (!this->minRefPressure_.empty())
1685 // The pore space change is irreversible
1686 effectivePressure =
1687 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1688 this->minRefPressure_[elementIdx]);
1689
1690 if (!this->overburdenPressure_.empty())
1691 effectivePressure -= this->overburdenPressure_[elementIdx];
1692
1693
1694 if (!this->rockCompPoroMult_.empty()) {
1695 return this->rockCompPoroMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
1696 }
1697
1698 // water compaction
1699 assert(!this->rockCompPoroMultWc_.empty());
1700 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1701 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
1702
1703 return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
1704 }
1705
1711 template <class LhsEval>
1712 LhsEval rockCompTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1713 {
1714 const bool implicit = !this->explicitRockCompaction_;
1715 return implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<LhsEval>(intQuants, elementIdx)
1716 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1717 }
1718
1724 template <class LhsEval>
1725 LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants) const
1726 {
1727 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
1728 if (!enableSaltPrecipitation)
1729 return 1.0;
1730
1731 const auto& fs = intQuants.fluidState();
1732 unsigned tableIdx = fs.pvtRegionIndex();
1733 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
1734 porosityFactor = min(porosityFactor, 1.0);
1735 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
1736 return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
1737 }
1738
1742 template <class LhsEval>
1743 LhsEval wellTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1744 {
1745 OPM_TIMEBLOCK_LOCAL(wellTransMultiplier);
1746
1747 const bool implicit = !this->explicitRockCompaction_;
1748 double trans_mult = implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<double>(intQuants, elementIdx)
1749 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1750 trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants);
1751
1752 return trans_mult;
1753 }
1754
1755 std::pair<BCType, RateVector> boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) const
1756 {
1757 OPM_TIMEBLOCK_LOCAL(boundaryCondition);
1758 if (!nonTrivialBoundaryConditions_) {
1759 return { BCType::NONE, RateVector(0.0) };
1760 }
1761 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1762 const auto& schedule = this->simulator().vanguard().schedule();
1763 if (bcindex_(dir)[globalSpaceIdx] == 0) {
1764 return { BCType::NONE, RateVector(0.0) };
1765 }
1766 if (schedule[this->episodeIndex()].bcprop.size() == 0) {
1767 return { BCType::NONE, RateVector(0.0) };
1768 }
1769 const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
1770 if (bc.bctype!=BCType::RATE) {
1771 return { bc.bctype, RateVector(0.0) };
1772 }
1773
1774 RateVector rate = 0.0;
1775 switch (bc.component) {
1776 case BCComponent::OIL:
1777 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = bc.rate;
1778 break;
1779 case BCComponent::GAS:
1780 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = bc.rate;
1781 break;
1782 case BCComponent::WATER:
1783 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate;
1784 break;
1785 case BCComponent::SOLVENT:
1786 if constexpr (!enableSolvent)
1787 throw std::logic_error("solvent is disabled and you're trying to add solvent to BC");
1788
1789 rate[Indices::solventSaturationIdx] = bc.rate;
1790 break;
1791 case BCComponent::POLYMER:
1792 if constexpr (!enablePolymer)
1793 throw std::logic_error("polymer is disabled and you're trying to add polymer to BC");
1794
1795 rate[Indices::polymerConcentrationIdx] = bc.rate;
1796 break;
1797 case BCComponent::NONE:
1798 throw std::logic_error("you need to specify the component when RATE type is set in BC");
1799 break;
1800 }
1801 //TODO add support for enthalpy rate
1802 return {bc.bctype, rate};
1803 }
1804
1805 const std::unique_ptr<EclWriterType>& eclWriter() const
1806 {
1807 return eclWriter_;
1808 }
1809
1810 void setConvData(const std::vector<std::vector<int>>& data)
1811 {
1812 eclWriter_->mutableOutputModule().setCnvData(data);
1813 }
1814
1815 template<class Serializer>
1816 void serializeOp(Serializer& serializer)
1817 {
1818 serializer(static_cast<BaseType&>(*this));
1819 serializer(drift_);
1820 serializer(wellModel_);
1821 serializer(aquiferModel_);
1822 serializer(tracerModel_);
1823 serializer(*materialLawManager_);
1824 serializer(*eclWriter_);
1825 }
1826private:
1827 Implementation& asImp_()
1828 { return *static_cast<Implementation *>(this); }
1829protected:
1830 void updateExplicitQuantities_()
1831 {
1832 OPM_TIMEBLOCK(updateExplicitQuantities);
1833 const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
1834 const bool invalidateFromMinPressure = updateMinPressure_();
1835
1836 // update hysteresis and max oil saturation used in vappars
1837 const bool invalidateFromHyst = updateHysteresis_();
1838 const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
1839
1840 // the derivatives may have change
1841 bool invalidateIntensiveQuantities
1842 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
1843 if (invalidateIntensiveQuantities) {
1844 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1845 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
1846 }
1847
1848 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
1849 updateMaxPolymerAdsorption_();
1850
1851 updateRockCompTransMultVal_();
1852 }
1853
1854 template<class UpdateFunc>
1855 void updateProperty_(const std::string& failureMsg,
1856 UpdateFunc func)
1857 {
1858 OPM_TIMEBLOCK(updateProperty);
1859 const auto& model = this->simulator().model();
1860 const auto& primaryVars = model.solution(/*timeIdx*/0);
1861 const auto& vanguard = this->simulator().vanguard();
1862 std::size_t numGridDof = primaryVars.size();
1863 OPM_BEGIN_PARALLEL_TRY_CATCH();
1864#ifdef _OPENMP
1865#pragma omp parallel for
1866#endif
1867 for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1868 const auto& iq = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/ 0);
1869 func(dofIdx, iq);
1870 }
1871 OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
1872 }
1873
1874 // update the parameters needed for DRSDT and DRVDT
1875 void updateCompositionChangeLimits_()
1876 {
1877 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1878 // update the "last Rs" values for all elements, including the ones in the ghost
1879 // and overlap regions
1880 int episodeIdx = this->episodeIndex();
1881 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1882 this->mixControls_.drsdtActive(episodeIdx),
1883 this->mixControls_.drvdtActive(episodeIdx)};
1884 if (!active[0] && !active[1] && !active[2]) {
1885 return;
1886 }
1887
1888 this->updateProperty_("FlowProblem::updateCompositionChangeLimits_()) failed:",
1889 [this,episodeIdx,active](unsigned compressedDofIdx,
1890 const IntensiveQuantities& iq)
1891 {
1892 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1893 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1894 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1895 this->mixControls_.update(compressedDofIdx,
1896 iq,
1897 episodeIdx,
1898 this->gravity_[dim - 1],
1899 perm[dim - 1][dim - 1],
1900 distZ,
1901 pvtRegionIdx,
1902 active);
1903 }
1904 );
1905 }
1906
1907 bool updateMaxOilSaturation_()
1908 {
1909 OPM_TIMEBLOCK(updateMaxOilSaturation);
1910 int episodeIdx = this->episodeIndex();
1911
1912 // we use VAPPARS
1913 if (this->vapparsActive(episodeIdx)) {
1914 this->updateProperty_("FlowProblem::updateMaxOilSaturation_() failed:",
1915 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1916 {
1917 this->updateMaxOilSaturation_(compressedDofIdx,iq);
1918 });
1919 return true;
1920 }
1921
1922 return false;
1923 }
1924
1925 bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1926 {
1927 OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
1928 const auto& fs = iq.fluidState();
1929 const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1930 auto& mos = this->maxOilSaturation_;
1931 if(mos[compressedDofIdx] < So){
1932 mos[compressedDofIdx] = So;
1933 return true;
1934 }else{
1935 return false;
1936 }
1937 }
1938
1939 bool updateMaxWaterSaturation_()
1940 {
1941 OPM_TIMEBLOCK(updateMaxWaterSaturation);
1942 // water compaction is activated in ROCKCOMP
1943 if (this->maxWaterSaturation_.empty())
1944 return false;
1945
1946 this->maxWaterSaturation_[/*timeIdx=*/1] = this->maxWaterSaturation_[/*timeIdx=*/0];
1947 this->updateProperty_("FlowProblem::updateMaxWaterSaturation_() failed:",
1948 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1949 {
1950 this->updateMaxWaterSaturation_(compressedDofIdx,iq);
1951 });
1952 return true;
1953 }
1954
1955
1956 bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1957 {
1958 OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation);
1959 const auto& fs = iq.fluidState();
1960 const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
1961 auto& mow = this->maxWaterSaturation_;
1962 if(mow[compressedDofIdx]< Sw){
1963 mow[compressedDofIdx] = Sw;
1964 return true;
1965 }else{
1966 return false;
1967 }
1968 }
1969
1970 bool updateMinPressure_()
1971 {
1972 OPM_TIMEBLOCK(updateMinPressure);
1973 // IRREVERS option is used in ROCKCOMP
1974 if (this->minRefPressure_.empty())
1975 return false;
1976
1977 this->updateProperty_("FlowProblem::updateMinPressure_() failed:",
1978 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1979 {
1980 this->updateMinPressure_(compressedDofIdx,iq);
1981 });
1982 return true;
1983 }
1984
1985 bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
1986 OPM_TIMEBLOCK_LOCAL(updateMinPressure);
1987 const auto& fs = iq.fluidState();
1988 const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
1989 auto& min_pressures = this->minRefPressure_;
1990 if(min_pressures[compressedDofIdx]> min_pressure){
1991 min_pressures[compressedDofIdx] = min_pressure;
1992 return true;
1993 }else{
1994 return false;
1995 }
1996 }
1997
1998 // \brief Function to assign field properties of type double, on the leaf grid view.
1999 //
2000 // For CpGrid with local grid refinement, the field property of a cell on the leaf
2001 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
2002 std::function<std::vector<double>(const FieldPropsManager&, const std::string&)>
2003 fieldPropDoubleOnLeafAssigner_()
2004 {
2005 const auto& lookup = this->lookUpData_;
2006 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString)
2007 {
2008 return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
2009 };
2010 }
2011
2012 // \brief Function to assign field properties of type int, unsigned int, ..., on the leaf grid view.
2013 //
2014 // For CpGrid with local grid refinement, the field property of a cell on the leaf
2015 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
2016 template<typename IntType>
2017 std::function<std::vector<IntType>(const FieldPropsManager&, const std::string&, bool)>
2018 fieldPropIntTypeOnLeafAssigner_()
2019 {
2020 const auto& lookup = this->lookUpData_;
2021 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString, bool needsTranslation)
2022 {
2023 return lookup.template assignFieldPropsIntOnLeaf<IntType>(fieldPropManager, propString, needsTranslation);
2024 };
2025 }
2026
2027 void readMaterialParameters_()
2028 {
2029 OPM_TIMEBLOCK(readMaterialParameters);
2030 const auto& simulator = this->simulator();
2031 const auto& vanguard = simulator.vanguard();
2032 const auto& eclState = vanguard.eclState();
2033
2034 // the PVT and saturation region numbers
2035 OPM_BEGIN_PARALLEL_TRY_CATCH();
2036 this->updatePvtnum_();
2037 this->updateSatnum_();
2038
2039 // the MISC region numbers (solvent model)
2040 this->updateMiscnum_();
2041 // the PLMIX region numbers (polymer model)
2042 this->updatePlmixnum_();
2043
2044 // directional relative permeabilities
2045 this->updateKrnum_();
2046 OPM_END_PARALLEL_TRY_CATCH("Invalid region numbers: ", vanguard.gridView().comm());
2048 // porosity
2049 updateReferencePorosity_();
2050 this->referencePorosity_[1] = this->referencePorosity_[0];
2052
2054 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
2055 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
2056 materialLawManager_->initFromState(eclState);
2057 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2058 this-> template fieldPropIntTypeOnLeafAssigner_<int>(),
2059 this-> lookupIdxOnLevelZeroAssigner_());
2061 }
2062
2063 void readThermalParameters_()
2064 {
2065 if constexpr (enableEnergy)
2066 {
2067 const auto& simulator = this->simulator();
2068 const auto& vanguard = simulator.vanguard();
2069 const auto& eclState = vanguard.eclState();
2070
2071 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
2072 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
2073 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
2074 this-> fieldPropDoubleOnLeafAssigner_(),
2075 this-> template fieldPropIntTypeOnLeafAssigner_<unsigned int>());
2076 }
2077 }
2078
2079 void updateReferencePorosity_()
2080 {
2081 const auto& simulator = this->simulator();
2082 const auto& vanguard = simulator.vanguard();
2083 const auto& eclState = vanguard.eclState();
2084
2085 std::size_t numDof = this->model().numGridDof();
2086
2087 this->referencePorosity_[/*timeIdx=*/0].resize(numDof);
2088
2089 const auto& fp = eclState.fieldProps();
2090 const std::vector<double> porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV");
2091 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2092 Scalar poreVolume = porvData[dofIdx];
2093
2094 // we define the porosity as the accumulated pore volume divided by the
2095 // geometric volume of the element. Note that -- in pathetic cases -- it can
2096 // be larger than 1.0!
2097 Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx);
2098 assert(dofVolume > 0.0);
2099 this->referencePorosity_[/*timeIdx=*/0][dofIdx] = poreVolume/dofVolume;
2100 }
2101 }
2102
2103 void readInitialCondition_()
2104 {
2105 const auto& simulator = this->simulator();
2106 const auto& vanguard = simulator.vanguard();
2107 const auto& eclState = vanguard.eclState();
2108
2109 if (eclState.getInitConfig().hasEquil())
2110 readEquilInitialCondition_();
2111 else
2112 readExplicitInitialCondition_();
2113
2114 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
2115 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
2116 enableSolvent,
2117 enablePolymer,
2118 enablePolymerMolarWeight,
2119 enableMICP);
2120
2121 //initialize min/max values
2122 std::size_t numElems = this->model().numGridDof();
2123 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2124 const auto& fs = initialFluidStates_[elemIdx];
2125 if (!this->maxWaterSaturation_.empty())
2126 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
2127 if (!this->maxOilSaturation_.empty())
2128 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
2129 if (!this->minRefPressure_.empty())
2130 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
2131 }
2132
2133
2134 }
2135
2136 void readEquilInitialCondition_()
2137 {
2138 const auto& simulator = this->simulator();
2139
2140 // initial condition corresponds to hydrostatic conditions.
2141 EquilInitializer<TypeTag> equilInitializer(simulator, *materialLawManager_);
2142
2143 std::size_t numElems = this->model().numGridDof();
2144 initialFluidStates_.resize(numElems);
2145 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2146 auto& elemFluidState = initialFluidStates_[elemIdx];
2147 elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
2148 }
2149 }
2150
2151 void readEclRestartSolution_()
2152 {
2153 // Throw an exception if the grid has LGRs. Refined grid are not supported for restart.
2154 if(this->simulator().vanguard().grid().maxLevel() > 0) {
2155 throw std::invalid_argument("Refined grids are not yet supported for restart ");
2156 }
2157
2158 // Set the start time of the simulation
2159 auto& simulator = this->simulator();
2160 const auto& schedule = simulator.vanguard().schedule();
2161 const auto& eclState = simulator.vanguard().eclState();
2162 const auto& initconfig = eclState.getInitConfig();
2163 {
2164 int restart_step = initconfig.getRestartStep();
2165
2166 simulator.setTime(schedule.seconds(restart_step));
2167
2168 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
2169 schedule.stepLength(restart_step));
2170 simulator.setEpisodeIndex(restart_step);
2171 }
2172 eclWriter_->beginRestart();
2173
2174 Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength());
2175 simulator.setTimeStepSize(dt);
2176
2177 std::size_t numElems = this->model().numGridDof();
2178 initialFluidStates_.resize(numElems);
2179 if constexpr (enableSolvent) {
2180 this->solventSaturation_.resize(numElems, 0.0);
2181 this->solventRsw_.resize(numElems, 0.0);
2182 }
2183
2184 if constexpr (enablePolymer)
2185 this->polymer_.concentration.resize(numElems, 0.0);
2186
2187 if constexpr (enablePolymerMolarWeight) {
2188 const std::string msg {"Support of the RESTART for polymer molecular weight "
2189 "is not implemented yet. The polymer weight value will be "
2190 "zero when RESTART begins"};
2191 OpmLog::warning("NO_POLYMW_RESTART", msg);
2192 this->polymer_.moleWeight.resize(numElems, 0.0);
2193 }
2194
2195 if constexpr (enableMICP) {
2196 this->micp_.resize(numElems);
2197 }
2198
2199 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
2200 auto& elemFluidState = initialFluidStates_[elemIdx];
2201 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
2202 eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
2203 eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
2204
2205 // Note: Function processRestartSaturations_() mutates the
2206 // 'ssol' argument--the value from the restart file--if solvent
2207 // is enabled. Then, store the updated solvent saturation into
2208 // 'solventSaturation_'. Otherwise, just pass a dummy value to
2209 // the function and discard the unchanged result. Do not index
2210 // into 'solventSaturation_' unless solvent is enabled.
2211 {
2212 auto ssol = enableSolvent
2213 ? eclWriter_->outputModule().getSolventSaturation(elemIdx)
2214 : Scalar(0);
2215
2216 processRestartSaturations_(elemFluidState, ssol);
2217
2218 if constexpr (enableSolvent) {
2219 this->solventSaturation_[elemIdx] = ssol;
2220 this->solventRsw_[elemIdx] = eclWriter_->outputModule().getSolventRsw(elemIdx);
2221 }
2222 }
2223
2224 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
2225
2226 if constexpr (enablePolymer)
2227 this->polymer_.concentration[elemIdx] = eclWriter_->outputModule().getPolymerConcentration(elemIdx);
2228 if constexpr (enableMICP){
2229 this->micp_.microbialConcentration[elemIdx] = eclWriter_->outputModule().getMicrobialConcentration(elemIdx);
2230 this->micp_.oxygenConcentration[elemIdx] = eclWriter_->outputModule().getOxygenConcentration(elemIdx);
2231 this->micp_.ureaConcentration[elemIdx] = eclWriter_->outputModule().getUreaConcentration(elemIdx);
2232 this->micp_.biofilmConcentration[elemIdx] = eclWriter_->outputModule().getBiofilmConcentration(elemIdx);
2233 this->micp_.calciteConcentration[elemIdx] = eclWriter_->outputModule().getCalciteConcentration(elemIdx);
2234 }
2235 // if we need to restart for polymer molecular weight simulation, we need to add related here
2236 }
2237
2238 const int episodeIdx = this->episodeIndex();
2239 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
2240
2241 // assign the restart solution to the current solution. note that we still need
2242 // to compute real initial solution after this because the initial fluid states
2243 // need to be correct for stuff like boundary conditions.
2244 auto& sol = this->model().solution(/*timeIdx=*/0);
2245 const auto& gridView = this->gridView();
2246 ElementContext elemCtx(simulator);
2247 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2248 elemCtx.updatePrimaryStencil(elem);
2249 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
2250 initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0);
2251 }
2252
2253 // make sure that the ghost and overlap entities exhibit the correct
2254 // solution. alternatively, this could be done in the loop above by also
2255 // considering non-interior elements. Since the initial() method might not work
2256 // 100% correctly for such elements, let's play safe and explicitly synchronize
2257 // using message passing.
2258 this->model().syncOverlap();
2259
2260 eclWriter_->endRestart();
2261 }
2262
2263 void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
2264 {
2265 // each phase needs to be above certain value to be claimed to be existing
2266 // this is used to recover some RESTART running with the defaulted single-precision format
2267 const Scalar smallSaturationTolerance = 1.e-6;
2268 Scalar sumSaturation = 0.0;
2269 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2270 if (FluidSystem::phaseIsActive(phaseIdx)) {
2271 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
2272 elemFluidState.setSaturation(phaseIdx, 0.0);
2273
2274 sumSaturation += elemFluidState.saturation(phaseIdx);
2275 }
2276
2277 }
2278 if constexpr (enableSolvent) {
2279 if (solventSaturation < smallSaturationTolerance)
2280 solventSaturation = 0.0;
2281
2282 sumSaturation += solventSaturation;
2283 }
2284
2285 assert(sumSaturation > 0.0);
2286
2287 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2288 if (FluidSystem::phaseIsActive(phaseIdx)) {
2289 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
2290 elemFluidState.setSaturation(phaseIdx, saturation);
2291 }
2292 }
2293 if constexpr (enableSolvent) {
2294 solventSaturation = solventSaturation / sumSaturation;
2295 }
2296 }
2297
2298 void readExplicitInitialCondition_()
2299 {
2300 const auto& simulator = this->simulator();
2301 const auto& vanguard = simulator.vanguard();
2302 const auto& eclState = vanguard.eclState();
2303 const auto& fp = eclState.fieldProps();
2304 bool has_swat = fp.has_double("SWAT");
2305 bool has_sgas = fp.has_double("SGAS");
2306 bool has_rs = fp.has_double("RS");
2307 bool has_rv = fp.has_double("RV");
2308 bool has_rvw = fp.has_double("RVW");
2309 bool has_pressure = fp.has_double("PRESSURE");
2310 bool has_salt = fp.has_double("SALT");
2311 bool has_saltp = fp.has_double("SALTP");
2312
2313 // make sure all required quantities are enables
2314 if (Indices::numPhases > 1) {
2315 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
2316 throw std::runtime_error("The ECL input file requires the presence of the SWAT keyword if "
2317 "the water phase is active");
2318 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
2319 throw std::runtime_error("The ECL input file requires the presence of the SGAS keyword if "
2320 "the gas phase is active");
2321 }
2322 if (!has_pressure)
2323 throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
2324 "keyword if the model is initialized explicitly");
2325 if (FluidSystem::enableDissolvedGas() && !has_rs)
2326 throw std::runtime_error("The ECL input file requires the RS keyword to be present if"
2327 " dissolved gas is enabled");
2328 if (FluidSystem::enableVaporizedOil() && !has_rv)
2329 throw std::runtime_error("The ECL input file requires the RV keyword to be present if"
2330 " vaporized oil is enabled");
2331 if (FluidSystem::enableVaporizedWater() && !has_rvw)
2332 throw std::runtime_error("The ECL input file requires the RVW keyword to be present if"
2333 " vaporized water is enabled");
2334 if (enableBrine && !has_salt)
2335 throw std::runtime_error("The ECL input file requires the SALT keyword to be present if"
2336 " brine is enabled and the model is initialized explicitly");
2337 if (enableSaltPrecipitation && !has_saltp)
2338 throw std::runtime_error("The ECL input file requires the SALTP keyword to be present if"
2339 " salt precipitation is enabled and the model is initialized explicitly");
2340
2341 std::size_t numDof = this->model().numGridDof();
2342
2343 initialFluidStates_.resize(numDof);
2344
2345 std::vector<double> waterSaturationData;
2346 std::vector<double> gasSaturationData;
2347 std::vector<double> pressureData;
2348 std::vector<double> rsData;
2349 std::vector<double> rvData;
2350 std::vector<double> rvwData;
2351 std::vector<double> tempiData;
2352 std::vector<double> saltData;
2353 std::vector<double> saltpData;
2354
2355 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
2356 waterSaturationData = fp.get_double("SWAT");
2357 else
2358 waterSaturationData.resize(numDof);
2359
2360 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
2361 gasSaturationData = fp.get_double("SGAS");
2362 else
2363 gasSaturationData.resize(numDof);
2364
2365 pressureData = fp.get_double("PRESSURE");
2366 if (FluidSystem::enableDissolvedGas())
2367 rsData = fp.get_double("RS");
2368
2369 if (FluidSystem::enableVaporizedOil())
2370 rvData = fp.get_double("RV");
2371
2372 if (FluidSystem::enableVaporizedWater())
2373 rvwData = fp.get_double("RVW");
2374
2375 // initial reservoir temperature
2376 tempiData = fp.get_double("TEMPI");
2377
2378 // initial salt concentration data
2379 if constexpr (enableBrine)
2380 saltData = fp.get_double("SALT");
2381
2382 // initial precipitated salt saturation data
2383 if constexpr (enableSaltPrecipitation)
2384 saltpData = fp.get_double("SALTP");
2385
2386 // calculate the initial fluid states
2387 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
2388 auto& dofFluidState = initialFluidStates_[dofIdx];
2389
2390 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
2391
2393 // set temperature
2395 Scalar temperatureLoc = tempiData[dofIdx];
2396 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
2397 temperatureLoc = FluidSystem::surfaceTemperature;
2398 dofFluidState.setTemperature(temperatureLoc);
2399
2401 // set salt concentration
2403 if constexpr (enableBrine)
2404 dofFluidState.setSaltConcentration(saltData[dofIdx]);
2405
2407 // set precipitated salt saturation
2409 if constexpr (enableSaltPrecipitation)
2410 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
2411
2413 // set saturations
2415 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
2416 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
2417 waterSaturationData[dofIdx]);
2418
2419 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
2420 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
2421 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2422 1.0
2423 - waterSaturationData[dofIdx]);
2424 }
2425 else
2426 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
2427 gasSaturationData[dofIdx]);
2428 }
2429 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
2430 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
2431 1.0
2432 - waterSaturationData[dofIdx]
2433 - gasSaturationData[dofIdx]);
2434
2436 // set phase pressures
2438 Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
2439
2440 // this assumes that capillary pressures only depend on the phase saturations
2441 // and possibly on temperature. (this is always the case for ECL problems.)
2442 std::array<Scalar, numPhases> pc = {0};
2443 const auto& matParams = materialLawParams(dofIdx);
2444 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
2445 Valgrind::CheckDefined(pressure);
2446 Valgrind::CheckDefined(pc);
2447 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2448 if (!FluidSystem::phaseIsActive(phaseIdx))
2449 continue;
2450
2451 if (Indices::oilEnabled)
2452 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
2453 else if (Indices::gasEnabled)
2454 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
2455 else if (Indices::waterEnabled)
2456 //single (water) phase
2457 dofFluidState.setPressure(phaseIdx, pressure);
2458 }
2459
2460 if (FluidSystem::enableDissolvedGas())
2461 dofFluidState.setRs(rsData[dofIdx]);
2462 else if (Indices::gasEnabled && Indices::oilEnabled)
2463 dofFluidState.setRs(0.0);
2464
2465 if (FluidSystem::enableVaporizedOil())
2466 dofFluidState.setRv(rvData[dofIdx]);
2467 else if (Indices::gasEnabled && Indices::oilEnabled)
2468 dofFluidState.setRv(0.0);
2469
2470 if (FluidSystem::enableVaporizedWater())
2471 dofFluidState.setRvw(rvwData[dofIdx]);
2472
2474 // set invB_
2476 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2477 if (!FluidSystem::phaseIsActive(phaseIdx))
2478 continue;
2479
2480 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2481 dofFluidState.setInvB(phaseIdx, b);
2482
2483 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
2484 dofFluidState.setDensity(phaseIdx, rho);
2485
2486 }
2487 }
2488 }
2489
2490 // update the hysteresis parameters of the material laws for the whole grid
2491 bool updateHysteresis_()
2492 {
2493 if (!materialLawManager_->enableHysteresis())
2494 return false;
2495
2496 // we need to update the hysteresis data for _all_ elements (i.e., not just the
2497 // interior ones) to avoid desynchronization of the processes in the parallel case!
2498 this->updateProperty_("FlowProblem::updateHysteresis_() failed:",
2499 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
2500 {
2501 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2502 });
2503 return true;
2504 }
2505
2506
2507 bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
2508 {
2509 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
2510 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
2511 //TODO change materials to give a bool
2512 return true;
2513 }
2514
2515 void updateMaxPolymerAdsorption_()
2516 {
2517 // we need to update the max polymer adsoption data for all elements
2518 this->updateProperty_("FlowProblem::updateMaxPolymerAdsorption_() failed:",
2519 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
2520 {
2521 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
2522 });
2523 }
2524
2525 bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
2526 {
2527 const Scalar pa = scalarValue(iq.polymerAdsorption());
2528 auto& mpa = this->polymer_.maxAdsorption;
2529 if (mpa[compressedDofIdx] < pa) {
2530 mpa[compressedDofIdx] = pa;
2531 return true;
2532 } else {
2533 return false;
2534 }
2535 }
2536
2537 Scalar getRockCompTransMultVal(std::size_t dofIdx) const
2538 {
2539 if (this->rockCompTransMultVal_.empty())
2540 return 1.0;
2541
2542 return this->rockCompTransMultVal_[dofIdx];
2543 }
2544
2545
2546private:
2547 struct PffDofData_
2548 {
2549 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
2550 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
2551 ConditionalStorage<enableDiffusion, Scalar> diffusivity;
2552 ConditionalStorage<enableDispersion, Scalar> dispersivity;
2553 Scalar transmissibility;
2554 };
2555
2556 // update the prefetch friendly data object
2557 void updatePffDofData_()
2558 {
2559 const auto& distFn =
2560 [this](PffDofData_& dofData,
2561 const Stencil& stencil,
2562 unsigned localDofIdx)
2563 -> void
2564 {
2565 const auto& elementMapper = this->model().elementMapper();
2566
2567 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
2568 if (localDofIdx != 0) {
2569 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(/*dofIdx=*/0));
2570 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
2571
2572 if constexpr (enableEnergy) {
2573 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
2574 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
2575 }
2576 if constexpr (enableDiffusion)
2577 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
2578 if (enableDispersion)
2579 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
2580 }
2581 };
2582
2583 pffDofData_.update(distFn);
2584 }
2585
2586 void readBoundaryConditions_()
2587 {
2588 const auto& simulator = this->simulator();
2589 const auto& vanguard = simulator.vanguard();
2590 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
2591 if (bcconfig.size() > 0) {
2592 nonTrivialBoundaryConditions_ = true;
2593
2594 std::size_t numCartDof = vanguard.cartesianSize();
2595 unsigned numElems = vanguard.gridView().size(/*codim=*/0);
2596 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
2597
2598 for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
2599 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
2600
2601 bcindex_.resize(numElems, 0);
2602 auto loopAndApply = [&cartesianToCompressedElemIdx,
2603 &vanguard](const auto& bcface,
2604 auto apply)
2605 {
2606 for (int i = bcface.i1; i <= bcface.i2; ++i) {
2607 for (int j = bcface.j1; j <= bcface.j2; ++j) {
2608 for (int k = bcface.k1; k <= bcface.k2; ++k) {
2609 std::array<int, 3> tmp = {i,j,k};
2610 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
2611 if (elemIdx >= 0)
2612 apply(elemIdx);
2613 }
2614 }
2615 }
2616 };
2617 for (const auto& bcface : bcconfig) {
2618 std::vector<int>& data = bcindex_(bcface.dir);
2619 const int index = bcface.index;
2620 loopAndApply(bcface,
2621 [&data,index](int elemIdx)
2622 { data[elemIdx] = index; });
2623 }
2624 }
2625 }
2626
2627 // this method applies the runtime constraints specified via the deck and/or command
2628 // line parameters for the size of the next time step.
2629 Scalar limitNextTimeStepSize_(Scalar dtNext) const
2630 {
2631 if constexpr (enableExperiments) {
2632 const auto& simulator = this->simulator();
2633 const auto& schedule = simulator.vanguard().schedule();
2634 int episodeIdx = simulator.episodeIndex();
2635
2636 // first thing in the morning, limit the time step size to the maximum size
2637 Scalar maxTimeStepSize = Parameters::get<TypeTag, Properties::SolverMaxTimeStepInDays>() * 24 * 60 * 60;
2638 int reportStepIdx = std::max(episodeIdx, 0);
2639 if (this->enableTuning_) {
2640 const auto& tuning = schedule[reportStepIdx].tuning();
2641 maxTimeStepSize = tuning.TSMAXZ;
2642 }
2643
2644 dtNext = std::min(dtNext, maxTimeStepSize);
2645
2646 Scalar remainingEpisodeTime =
2647 simulator.episodeStartTime() + simulator.episodeLength()
2648 - (simulator.startTime() + simulator.time());
2649 assert(remainingEpisodeTime >= 0.0);
2650
2651 // if we would have a small amount of time left over in the current episode, make
2652 // two equal time steps instead of a big and a small one
2653 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
2654 // note: limiting to the maximum time step size here is probably not strictly
2655 // necessary, but it should not hurt and is more fool-proof
2656 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
2657
2658 if (simulator.episodeStarts()) {
2659 // if a well event occurred, respect the limit for the maximum time step after
2660 // that, too
2661 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
2662 bool wellEventOccured =
2663 events.hasEvent(ScheduleEvents::NEW_WELL)
2664 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
2665 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
2666 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
2667 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
2668 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
2669 }
2670 }
2671
2672 return dtNext;
2673 }
2674
2675 void computeAndSetEqWeights_()
2676 {
2677 std::vector<Scalar> sumInvB(numPhases, 0.0);
2678 const auto& gridView = this->gridView();
2679 ElementContext elemCtx(this->simulator());
2680 for(const auto& elem: elements(gridView, Dune::Partitions::interior)) {
2681 elemCtx.updatePrimaryStencil(elem);
2682 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
2683 const auto& dofFluidState = initialFluidStates_[elemIdx];
2684 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2685 if (!FluidSystem::phaseIsActive(phaseIdx))
2686 continue;
2687
2688 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
2689 }
2690 }
2691
2692 std::size_t numDof = this->model().numGridDof();
2693 const auto& comm = this->simulator().vanguard().grid().comm();
2694 comm.sum(sumInvB.data(),sumInvB.size());
2695 Scalar numTotalDof = comm.sum(numDof);
2696
2697 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
2698 if (!FluidSystem::phaseIsActive(phaseIdx))
2699 continue;
2700
2701 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
2702 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
2703 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
2704 this->model().setEqWeight(activeSolventCompIdx, avgB);
2705 }
2706 }
2707
2708 int refPressurePhaseIdx_() const {
2709 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
2710 return oilPhaseIdx;
2711 }
2712 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
2713 return gasPhaseIdx;
2714 }
2715 else {
2716 return waterPhaseIdx;
2717 }
2718 }
2719
2720 void updateRockCompTransMultVal_()
2721 {
2722 const auto& model = this->simulator().model();
2723 std::size_t numGridDof = this->model().numGridDof();
2724 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
2725 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
2726 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, /*timeIdx=*/ 0);
2727 Scalar trans_mult = computeRockCompTransMultiplier_<Scalar>(iq, elementIdx);
2728 this->rockCompTransMultVal_[elementIdx] = trans_mult;
2729 }
2730 }
2731
2737 template <class LhsEval>
2738 LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities& intQuants, unsigned elementIdx) const
2739 {
2740 OPM_TIMEBLOCK_LOCAL(computeRockCompTransMultiplier);
2741 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
2742 return 1.0;
2743
2744 unsigned tableIdx = 0;
2745 if (!this->rockTableIdx_.empty())
2746 tableIdx = this->rockTableIdx_[elementIdx];
2747
2748 const auto& fs = intQuants.fluidState();
2749 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
2750
2751 if (!this->minRefPressure_.empty())
2752 // The pore space change is irreversible
2753 effectivePressure =
2754 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
2755 this->minRefPressure_[elementIdx]);
2756
2757 if (!this->overburdenPressure_.empty())
2758 effectivePressure -= this->overburdenPressure_[elementIdx];
2759
2760 if (!this->rockCompTransMult_.empty())
2761 return this->rockCompTransMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
2762
2763 // water compaction
2764 assert(!this->rockCompTransMultWc_.empty());
2765 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
2766 LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
2767
2768 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
2769 }
2770
2771 typename Vanguard::TransmissibilityType transmissibilities_;
2772
2773 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
2774 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
2775
2776 FlowThresholdPressure<TypeTag> thresholdPressures_;
2777
2778 std::vector<InitialFluidState> initialFluidStates_;
2779
2780 bool enableDriftCompensation_;
2781 GlobalEqVector drift_;
2782
2783 WellModel wellModel_;
2784 AquiferModel aquiferModel_;
2785
2786 bool enableEclOutput_;
2787 std::unique_ptr<EclWriterType> eclWriter_;
2788
2789#if HAVE_DAMARIS
2790 bool enableDamarisOutput_ = false ;
2791 std::unique_ptr<DamarisWriterType> damarisWriter_;
2792#endif
2793
2794 PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
2795 TracerModel tracerModel_;
2796
2797 ActionHandler actionHandler_;
2798
2799 template<class T>
2800 struct BCData
2801 {
2802 std::array<std::vector<T>,6> data;
2803
2804 void resize(std::size_t size, T defVal)
2805 {
2806 for (auto& d : data)
2807 d.resize(size, defVal);
2808 }
2809
2810 const std::vector<T>& operator()(FaceDir::DirEnum dir) const
2811 {
2812 if (dir == FaceDir::DirEnum::Unknown)
2813 throw std::runtime_error("Tried to access BC data for the 'Unknown' direction");
2814 int idx = 0;
2815 int div = static_cast<int>(dir);
2816 while ((div /= 2) >= 1)
2817 ++idx;
2818 assert(idx >= 0 && idx <= 5);
2819 return data[idx];
2820 }
2821
2822 std::vector<T>& operator()(FaceDir::DirEnum dir)
2823 {
2824 return const_cast<std::vector<T>&>(std::as_const(*this)(dir));
2825 }
2826 };
2827
2828 BCData<int> bcindex_;
2829 bool nonTrivialBoundaryConditions_ = false;
2830 bool explicitRockCompaction_ = false;
2831};
2832
2833} // namespace Opm
2834
2835#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:1506
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
Direct access to the transmissibility between two elements.
Definition FlowProblem.hpp:808
Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
Definition FlowProblem.hpp:927
const DimMatrix & intrinsicPermeability(unsigned globalElemIdx) const
This method returns the intrinsic permeability tensor given a global element index.
Definition FlowProblem.hpp:790
void endEpisode()
Called by the simulator after the end of an episode.
Definition FlowProblem.hpp:714
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1096
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:1262
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:951
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition FlowProblem.hpp:612
unsigned satnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1104
Scalar thermalHalfTransmissibilityOut(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:898
Scalar transmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:855
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn, const unsigned globalSpaceIdxOut) const
Definition FlowProblem.hpp:875
void endTimeStep()
Called by the simulator after each time integration.
Definition FlowProblem.hpp:632
Scalar rockCompressibility(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:985
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:1277
Scalar maxOilSaturation(unsigned globalDofIdx) const
Returns an element's historic maximum oil phase saturation that was observed during the simulation.
Definition FlowProblem.hpp:1222
std::string name() const
Definition FlowProblem.hpp:1134
LhsEval rockCompTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1712
void endIteration()
Called by the simulator after each Newton-Raphson iteration.
Definition FlowProblem.hpp:622
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:921
Scalar nextTimeStepSize() const
Propose the size of the next time step to the simulator.
Definition FlowProblem.hpp:1645
void writeOutput(const SimulatorTimer &timer, bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:738
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1040
Scalar transmissibility(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:797
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
Direct indexed acces to the depth of an degree of freedom [m].
Definition FlowProblem.hpp:976
LhsEval wellTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Return the well transmissibility multiplier due to rock changues.
Definition FlowProblem.hpp:1743
std::shared_ptr< const EclMaterialLawManager > materialLawManager() const
Returns the ECL material law manager.
Definition FlowProblem.hpp:1052
unsigned plmixnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1120
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1141
Scalar thermalHalfTransmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:911
unsigned miscnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:1112
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:1251
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:776
Scalar maxPolymerAdsorption(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the max polymer adsorption value.
Definition FlowProblem.hpp:1128
void beginTimeStep()
Called by the simulator before each time integration.
Definition FlowProblem.hpp:582
Scalar thermalHalfTransmissibilityIn(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:885
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
Sets an element's maximum oil phase saturation observed during the simulation.
Definition FlowProblem.hpp:1239
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:1027
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the dispersivity for a face i.e.
Definition FlowProblem.hpp:835
void finishInit()
Definition FlowProblem.hpp:343
void initialSolutionApplied()
Definition FlowProblem.hpp:1333
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:513
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1371
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1005
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:842
Scalar diffusivity(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:817
void deserialize(Restarter &res)
This method restores the complete state of the problem and its sub-objects from disk.
Definition FlowProblem.hpp:479
std::shared_ptr< EclMaterialLawManager > materialLawManager()
Definition FlowProblem.hpp:1088
Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the depth of an degree of freedom [m].
Definition FlowProblem.hpp:964
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1293
LhsEval permFactTransMultiplier(const IntensiveQuantities &intQuants) const
Calculate the transmissibility multiplier due to porosity reduction.
Definition FlowProblem.hpp:1725
LhsEval rockCompPoroMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the porosity multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1672
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:865
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:1176
void serialize(Restarter &res)
This method writes the complete state of the problem and its subobjects to disk.
Definition FlowProblem.hpp:498
Scalar rockReferencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:995
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the transmissibility for a face i.e.
Definition FlowProblem.hpp:828
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