My Project
Loading...
Searching...
No Matches
FlowProblemBlackoil.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 Copyright 2023 INRIA
5 Copyright 2024 SINTEF Digital
6
7 This file is part of the Open Porous Media project (OPM).
8
9 OPM is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 2 of the License, or
12 (at your option) any later version.
13
14 OPM is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with OPM. If not, see <http://www.gnu.org/licenses/>.
21
22 Consult the COPYING file in the top-level source directory of this
23 module for the precise wording of the license and the list of
24 copyright holders.
25*/
31#ifndef OPM_FLOW_PROBLEM_BLACK_HPP
32#define OPM_FLOW_PROBLEM_BLACK_HPP
33
34#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
35#include <opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp>
36#include <opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp>
37#include <opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp>
38#include <opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp>
39#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp>
40#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp>
41
43
44#include <opm/output/eclipse/EclipseIO.hpp>
45
50
52
53#include <opm/simulators/flow/ActionHandler.hpp>
54#if HAVE_DAMARIS
56#endif
57
58#include <algorithm>
59#include <functional>
60#include <set>
61#include <string>
62#include <vector>
63
64namespace Opm {
65
72template <class TypeTag>
73class FlowProblemBlackoil : public FlowProblem<TypeTag>
74{
75 // TODO: the naming of the Types might be able to be adjusted
76 using FlowProblemType = FlowProblem<TypeTag>;
77
78 using typename FlowProblemType::Scalar;
79 using typename FlowProblemType::Simulator;
80 using typename FlowProblemType::GridView;
81 using typename FlowProblemType::FluidSystem;
82 using typename FlowProblemType::Vanguard;
83 using typename FlowProblemType::GlobalEqVector;
84 using typename FlowProblemType::EqVector;
85 using FlowProblemType::dim;
86 using FlowProblemType::dimWorld;
87 using FlowProblemType::numEq;
88 using FlowProblemType::numPhases;
89 using FlowProblemType::numComponents;
90
91 // TODO: potentially some cleaning up depending on the usage later here
92 using FlowProblemType::enableConvectiveMixing;
93 using FlowProblemType::enableBrine;
94 using FlowProblemType::enableDiffusion;
95 using FlowProblemType::enableDispersion;
96 using FlowProblemType::enableEnergy;
97 using FlowProblemType::enableExperiments;
98 using FlowProblemType::enableExtbo;
99 using FlowProblemType::enableFoam;
100 using FlowProblemType::enableMICP;
101 using FlowProblemType::enablePolymer;
102 using FlowProblemType::enablePolymerMolarWeight;
103 using FlowProblemType::enableSaltPrecipitation;
104 using FlowProblemType::enableSolvent;
105 using FlowProblemType::enableTemperature;
106 using FlowProblemType::enableThermalFluxBoundaries;
107
108 using FlowProblemType::gasPhaseIdx;
109 using FlowProblemType::oilPhaseIdx;
110 using FlowProblemType::waterPhaseIdx;
111
112 using FlowProblemType::waterCompIdx;
113 using FlowProblemType::oilCompIdx;
114 using FlowProblemType::gasCompIdx;
115
116 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
117 using typename FlowProblemType::RateVector;
118 using typename FlowProblemType::PrimaryVariables;
119 using typename FlowProblemType::Indices;
120 using typename FlowProblemType::IntensiveQuantities;
121 using typename FlowProblemType::ElementContext;
122
123 using typename FlowProblemType::MaterialLaw;
124 using typename FlowProblemType::DimMatrix;
125
126 using SolventModule = BlackOilSolventModule<TypeTag>;
127 using PolymerModule = BlackOilPolymerModule<TypeTag>;
128 using FoamModule = BlackOilFoamModule<TypeTag>;
129 using BrineModule = BlackOilBrineModule<TypeTag>;
130 using ExtboModule = BlackOilExtboModule<TypeTag>;
131 using MICPModule = BlackOilMICPModule<TypeTag>;
132 using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
133 using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
134 using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
135 using ModuleParams = typename BlackOilLocalResidualTPFA<TypeTag>::ModuleParams;
136
137 using InitialFluidState = typename EquilInitializer<TypeTag>::ScalarFluidState;
138 using EclWriterType = EclWriter<TypeTag>;
139#if HAVE_DAMARIS
140 using DamarisWriterType = DamarisWriter<TypeTag>;
141#endif
142
143
144public:
147
151 static void registerParameters()
152 {
154
155 EclWriterType::registerParameters();
156#if HAVE_DAMARIS
157 DamarisWriterType::registerParameters();
158#endif
160 }
161
165 explicit FlowProblemBlackoil(Simulator& simulator)
166 : FlowProblemType(simulator)
167 , thresholdPressures_(simulator)
168 , mixControls_(simulator.vanguard().schedule())
169 , actionHandler_(simulator.vanguard().eclState(),
170 simulator.vanguard().schedule(),
171 simulator.vanguard().actionState(),
172 simulator.vanguard().summaryState(),
173 this->wellModel_,
174 simulator.vanguard().grid().comm())
175 {
176 this->model().addOutputModule(new VtkTracerModule<TypeTag>(simulator));
177
178 // Tell the black-oil extensions to initialize their internal data structures
179 const auto& vanguard = simulator.vanguard();
180
181 BlackOilBrineParams<Scalar> brineParams;
182 brineParams.template initFromState<enableBrine,
183 enableSaltPrecipitation>(vanguard.eclState());
184 BrineModule::setParams(std::move(brineParams));
185
186 DiffusionModule::initFromState(vanguard.eclState());
187 DispersionModule::initFromState(vanguard.eclState());
188
189 BlackOilExtboParams<Scalar> extboParams;
190 extboParams.template initFromState<enableExtbo>(vanguard.eclState());
191 ExtboModule::setParams(std::move(extboParams));
192
194 foamParams.template initFromState<enableFoam>(vanguard.eclState());
195 FoamModule::setParams(std::move(foamParams));
196
198 micpParams.template initFromState<enableMICP>(vanguard.eclState());
199 MICPModule::setParams(std::move(micpParams));
200
201 BlackOilPolymerParams<Scalar> polymerParams;
202 polymerParams.template initFromState<enablePolymer, enablePolymerMolarWeight>(vanguard.eclState());
203 PolymerModule::setParams(std::move(polymerParams));
204
205 BlackOilSolventParams<Scalar> solventParams;
206 solventParams.template initFromState<enableSolvent>(vanguard.eclState(), vanguard.schedule());
207 SolventModule::setParams(std::move(solventParams));
208
209 // create the ECL writer
210 eclWriter_ = std::make_unique<EclWriterType>(simulator);
211 enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
212#if HAVE_DAMARIS
213 // create Damaris writer
214 damarisWriter_ = std::make_unique<DamarisWriterType>(simulator);
215 enableDamarisOutput_ = Parameters::Get<Parameters::EnableDamarisOutput>();
216#endif
217 }
218
222 void beginEpisode() override
223 {
225
226 auto& simulator = this->simulator();
227 int episodeIdx = simulator.episodeIndex();
228 const auto& schedule = simulator.vanguard().schedule();
229
230 // Evaluate UDQ assign statements to make sure the settings are
231 // available as UDA controls for the current report step.
232 actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState());
233
234 if (episodeIdx >= 0) {
235 const auto& oilVap = schedule[episodeIdx].oilvap();
236 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
237 FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2());
238 } else {
239 FluidSystem::setVapPars(0.0, 0.0);
240 }
241 }
242
243 ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), simulator.vanguard().schedule(), episodeIdx, moduleParams_.convectiveMixingModuleParam);
244 }
245
250 {
251 // TODO: there should be room to remove duplication for this function,
252 // but there is relatively complicated logic in the function calls in this function
253 // some refactoring is needed for this function
254 FlowProblemType::finishInit();
255 auto& simulator = this->simulator();
256
257 auto finishTransmissibilities = [updated = false, this]() mutable
258 {
259 if (updated) { return; }
260 this->transmissibilities_.finishInit([&vg = this->simulator().vanguard()](const unsigned int it) {
261 return vg.gridIdxToEquilGridIdx(it);
262 });
263 updated = true;
264 };
265
266 // calculating the TRANX, TRANY, TRANZ and NNC for output purpose
267 // for parallel running, it is based on global trans_
268 // for serial running, it is based on the transmissibilities_
269 // we try to avoid for the parallel running, has both global trans_ and transmissibilities_ allocated at the same time
270 if (enableEclOutput_) {
271 if (simulator.vanguard().grid().comm().size() > 1) {
272 if (simulator.vanguard().grid().comm().rank() == 0)
273 eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
274 } else {
275 finishTransmissibilities();
276 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
277 }
278
279 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
280 return simulator.vanguard().gridEquilIdxToGridIdx(i);
281 };
282 eclWriter_->extractOutputTransAndNNC(equilGridToGrid);
283 }
284 simulator.vanguard().releaseGlobalTransmissibilities();
285
286 const auto& eclState = simulator.vanguard().eclState();
287 const auto& schedule = simulator.vanguard().schedule();
288
289 // Set the start time of the simulation
290 simulator.setStartTime(schedule.getStartTime());
291 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
292
293 // We want the episode index to be the same as the report step index to make
294 // things simpler, so we have to set the episode index to -1 because it is
295 // incremented by endEpisode(). The size of the initial time step and
296 // length of the initial episode is set to zero for the same reason.
297 simulator.setEpisodeIndex(-1);
298 simulator.setEpisodeLength(0.0);
299
300 // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
301 // disables gravity, else the standard value of the gravity constant at sea level
302 // on earth is used
303 this->gravity_ = 0.0;
304 if (Parameters::Get<Parameters::EnableGravity>())
305 this->gravity_[dim - 1] = 9.80665;
306 if (!eclState.getInitConfig().hasGravity())
307 this->gravity_[dim - 1] = 0.0;
308
309 if (this->enableTuning_) {
310 // if support for the TUNING keyword is enabled, we get the initial time
311 // steping parameters from it instead of from command line parameters
312 const auto& tuning = schedule[0].tuning();
313 this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
314 this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
315 }
316
317 this->initFluidSystem_();
318
319
320
321 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
322 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
323 this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
324 }
325
326 this->readRockParameters_(simulator.vanguard().cellCenterDepths(),
327 [&simulator](const unsigned idx)
328 {
329 std::array<int,dim> coords;
330 simulator.vanguard().cartesianCoordinate(idx, coords);
331 for (auto& c : coords) {
332 ++c;
333 }
334 return coords;
335 });
336 this->readMaterialParameters_();
337 this->readThermalParameters_();
338
339 // write the static output files (EGRID, INIT)
340 if (enableEclOutput_) {
341 eclWriter_->writeInit();
342 }
343
344 finishTransmissibilities();
345
346 const auto& initconfig = eclState.getInitConfig();
347 this->tracerModel_.init(initconfig.restartRequested());
348 if (initconfig.restartRequested())
349 readEclRestartSolution_();
350 else
351 readInitialCondition_();
352
353 this->tracerModel_.prepareTracerBatches();
354
355 this->updatePffDofData_();
356
358 const auto& vanguard = this->simulator().vanguard();
359 const auto& gridView = vanguard.gridView();
360 int numElements = gridView.size(/*codim=*/0);
361 this->polymer_.maxAdsorption.resize(numElements, 0.0);
362 }
363
364 this->readBoundaryConditions_();
365
366 // compute and set eq weights based on initial b values
367 computeAndSetEqWeights_();
368
369 if (this->enableDriftCompensation_) {
370 this->drift_.resize(this->model().numGridDof());
371 this->drift_ = 0.0;
372 }
373
374 if (this->enableVtkOutput_ && eclState.getIOConfig().initOnly()) {
375 simulator.setTimeStepSize(0.0);
377 }
378
379 // after finishing the initialization and writing the initial solution, we move
380 // to the first "real" episode/report step
381 // for restart the episode index and start is already set
382 if (!initconfig.restartRequested()) {
383 simulator.startNextEpisode(schedule.seconds(1));
384 simulator.setEpisodeIndex(0);
385 simulator.setTimeStepIndex(0);
386 }
387
388 // TODO: move to the end for later refactoring of the function finishInit()
389 // deal with DRSDT
390 this->mixControls_.init(this->model().numGridDof(),
391 this->episodeIndex(),
392 eclState.runspec().tabdims().getNumPVTTables());
393 }
394
398 void endTimeStep () override
399 {
400 FlowProblemType::endTimeStep();
401
402 // after the solution is updated, the values in output module needs also updated
403 this->eclWriter()->mutableOutputModule().invalidateLocalData();
404
405 const bool isSubStep = !this->simulator().episodeWillBeOver();
406
407 // For CpGrid with LGRs, ecl/vtk output is not supported yet.
408 const auto& grid = this->simulator().vanguard().gridView().grid();
409
410 using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
411 constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
412 if (!isCpGrid || (grid.maxLevel() == 0)) {
413 this->eclWriter_->evalSummaryState(isSubStep);
414 }
415
416 {
417 OPM_TIMEBLOCK(applyActions);
418
419 const int episodeIdx = this->episodeIndex();
420 auto& simulator = this->simulator();
421
422 // Re-ordering in case of Alugrid
423 this->actionHandler_
424 .applyActions(episodeIdx, simulator.time() + simulator.timeStepSize(),
425 [this](const bool global)
426 {
427 using TransUpdateQuantities = typename Vanguard::TransmissibilityType::TransUpdateQuantities;
428 this->transmissibilities_
429 .update(global, TransUpdateQuantities::All, [&vg = this->simulator().vanguard()]
430 (const unsigned int i)
431 {
432 return vg.gridIdxToEquilGridIdx(i);
433 });
434 });
435 }
436
437 // Deal with "clogging" for the MICP model
438 if constexpr (enableMICP) {
439 auto& model = this->model();
440 const auto& residual = model.linearizer().residual();
441
442 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); ++globalDofIdx) {
443 auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx];
444 MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
445 }
446 }
447
448 }
452 void endEpisode() override
453 {
454 OPM_TIMEBLOCK(endEpisode);
455 const int episodeIdx = this->episodeIndex();
456 // Rerun UDQ assignents following action processing on the final
457 // time step of this episode to make sure that any UDQ ASSIGN
458 // operations triggered in action blocks take effect. This is
459 // mainly to work around a shortcoming of the ScheduleState copy
460 // constructor which clears pending UDQ assignments under the
461 // assumption that all such assignments have been processed. If an
462 // action block happens to trigger on the final time step of an
463 // episode and that action block runs a UDQ assignment, then that
464 // assignment would be dropped and the rest of the simulator will
465 // never see its effect without this hack.
466 this->actionHandler_
467 .evalUDQAssignments(episodeIdx, this->simulator().vanguard().udqState());
468
469 FlowProblemType::endEpisode();
470 }
471
472 void writeReports(const SimulatorTimer& timer) {
473 if (enableEclOutput_){
474 eclWriter_->writeReports(timer);
475 }
476 }
477
478
483 void writeOutput(bool verbose = true)
484 {
485 FlowProblemType::writeOutput(verbose);
486
487 bool isSubStep = !this->simulator().episodeWillBeOver();
488
489 data::Solution localCellData = {};
490#if HAVE_DAMARIS
491 // N.B. the Damaris output has to be done before the ECL output as the ECL one
492 // does all kinds of std::move() relocation of data
493 if (enableDamarisOutput_) {
494 damarisWriter_->writeOutput(localCellData, isSubStep) ;
495 }
496#endif
497 if (enableEclOutput_){
498 eclWriter_->writeOutput(std::move(localCellData), isSubStep);
499 }
500 }
501
502 void finalizeOutput()
503 {
504 OPM_TIMEBLOCK(finalizeOutput);
505 // this will write all pending output to disk
506 // to avoid corruption of output files
507 eclWriter_.reset();
508 }
509
510
515 {
516 FlowProblemType::initialSolutionApplied();
517
518 // let the object for threshold pressures initialize itself. this is done only at
519 // this point, because determining the threshold pressures may require to access
520 // the initial solution.
521 this->thresholdPressures_.finishInit();
522
523 if (this->simulator().episodeIndex() == 0) {
524 eclWriter_->writeInitialFIPReport();
525 }
526 }
527
528 void addToSourceDense(RateVector& rate,
529 unsigned globalDofIdx,
530 unsigned timeIdx) const override
531 {
532 this->aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
533
534 // Add source term from deck
535 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
536 std::array<int,3> ijk;
537 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
538
539 if (source.hasSource(ijk)) {
540 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
541 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
542 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
543 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
544
545 for (unsigned i = 0; i < phidx_map.size(); ++i) {
546 const auto phaseIdx = phidx_map[i];
547 const auto sourceComp = sc_map[i];
548 const auto compIdx = cidx_map[i];
549 if (!FluidSystem::phaseIsActive(phaseIdx)) {
550 continue;
551 }
552 Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
553 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
554 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
555 }
556 rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate;
557 }
558
559 if constexpr (enableSolvent) {
560 Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
561 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
562 const auto& solventPvt = SolventModule::solventPvt();
563 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
564 }
565 rate[Indices::contiSolventEqIdx] += mass_rate;
566 }
567 if constexpr (enablePolymer) {
568 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
569 }
570 if constexpr (enableEnergy) {
571 for (unsigned i = 0; i < phidx_map.size(); ++i) {
572 const auto phaseIdx = phidx_map[i];
573 if (!FluidSystem::phaseIsActive(phaseIdx)) {
574 continue;
575 }
576 const auto sourceComp = sc_map[i];
577 if (source.hasHrate({ijk, sourceComp})) {
578 rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
579 } else {
580 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, /*timeIdx*/ 0);
581 auto fs = intQuants.fluidState();
582 // if temperature is not set, use cell temperature as default
583 if (source.hasTemperature({ijk, sourceComp})) {
584 Scalar temperature = source.temperature({ijk, sourceComp});
585 fs.setTemperature(temperature);
586 }
587 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
588 Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx);
589 Scalar energy_rate = getValue(h)*mass_rate;
590 rate[Indices::contiEnergyEqIdx] += energy_rate;
591 }
592 }
593 }
594 }
595
596 // if requested, compensate systematic mass loss for cells which were "well
597 // behaved" in the last time step
598 if (this->enableDriftCompensation_) {
599 const auto& simulator = this->simulator();
600 const auto& model = this->model();
601
602 // we use a lower tolerance for the compensation too
603 // assure the added drift from the last step does not
604 // cause convergence issues on the current step
605 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
606 Scalar poro = this->porosity(globalDofIdx, timeIdx);
607 Scalar dt = simulator.timeStepSize();
608 EqVector dofDriftRate = this->drift_[globalDofIdx];
609 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
610
611 // restrict drift compensation to the CNV tolerance
612 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
613 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
614 if (cnv > maxCompensation) {
615 dofDriftRate[eqIdx] *= maxCompensation/cnv;
616 }
617 }
618
619 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
620 rate[eqIdx] -= dofDriftRate[eqIdx];
621 }
622 }
623
629 template <class LhsEval>
630 LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants) const
631 {
632 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
633 if (!enableSaltPrecipitation)
634 return 1.0;
635
636 const auto& fs = intQuants.fluidState();
637 unsigned tableIdx = fs.pvtRegionIndex();
638 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
639 porosityFactor = min(porosityFactor, 1.0);
640 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
641 return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
642 }
643
644 // temporary solution to facilitate output of initial state from flow
645 const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
646 { return initialFluidStates_[globalDofIdx]; }
647
648 std::vector<InitialFluidState>& initialFluidStates()
649 { return initialFluidStates_; }
650
651 const std::vector<InitialFluidState>& initialFluidStates() const
652 { return initialFluidStates_; }
653
654 const EclipseIO& eclIO() const
655 { return eclWriter_->eclIO(); }
656
657 void setSubStepReport(const SimulatorReportSingle& report)
658 { return eclWriter_->setSubStepReport(report); }
659
660 void setSimulationReport(const SimulatorReport& report)
661 { return eclWriter_->setSimulationReport(report); }
662
663 InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const
664 {
665 OPM_TIMEBLOCK_LOCAL(boundaryFluidState);
666 const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop;
667 if (bcprop.size() > 0) {
668 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
669
670 // index == 0: no boundary conditions for this
671 // global cell and direction
672 if (this->bcindex_(dir)[globalDofIdx] == 0)
673 return initialFluidStates_[globalDofIdx];
674
675 const auto& bc = bcprop[this->bcindex_(dir)[globalDofIdx]];
676 if (bc.bctype == BCType::DIRICHLET )
677 {
678 InitialFluidState fluidState;
679 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
680 fluidState.setPvtRegionIndex(pvtRegionIdx);
681
682 switch (bc.component) {
683 case BCComponent::OIL:
684 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
685 throw std::logic_error("oil is not active and you're trying to add oil BC");
686
687 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
688 break;
689 case BCComponent::GAS:
690 if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
691 throw std::logic_error("gas is not active and you're trying to add gas BC");
692
693 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0);
694 break;
695 case BCComponent::WATER:
696 if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
697 throw std::logic_error("water is not active and you're trying to add water BC");
698
699 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0);
700 break;
701 case BCComponent::SOLVENT:
702 case BCComponent::POLYMER:
703 case BCComponent::NONE:
704 throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
705 }
706 fluidState.setTotalSaturation(1.0);
707 double pressure = initialFluidStates_[globalDofIdx].pressure(this->refPressurePhaseIdx_());
708 const auto pressure_input = bc.pressure;
709 if (pressure_input) {
710 pressure = *pressure_input;
711 }
712
713 std::array<Scalar, numPhases> pc = {0};
714 const auto& matParams = this->materialLawParams(globalDofIdx);
715 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
716 Valgrind::CheckDefined(pressure);
717 Valgrind::CheckDefined(pc);
718 for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
719 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
720
721 fluidState.setPc(phaseIdx, pc[phaseIdx]);
722 if (Indices::oilEnabled)
723 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
724 else if (Indices::gasEnabled)
725 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
726 else if (Indices::waterEnabled)
727 //single (water) phase
728 fluidState.setPressure(phaseIdx, pressure);
729 }
730
731 double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature
732 const auto temperature_input = bc.temperature;
733 if(temperature_input)
734 temperature = *temperature_input;
735 fluidState.setTemperature(temperature);
736
737 if (FluidSystem::enableDissolvedGas()) {
738 fluidState.setRs(0.0);
739 fluidState.setRv(0.0);
740 }
741 if (FluidSystem::enableDissolvedGasInWater()) {
742 fluidState.setRsw(0.0);
743 }
744 if (FluidSystem::enableVaporizedWater())
745 fluidState.setRvw(0.0);
746
747 for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
748 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
749
750 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
751 fluidState.setInvB(phaseIdx, b);
752
753 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
754 fluidState.setDensity(phaseIdx, rho);
755 if (enableEnergy) {
756 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
757 fluidState.setEnthalpy(phaseIdx, h);
758 }
759 }
760 fluidState.checkDefined();
761 return fluidState;
762 }
763 }
764 return initialFluidStates_[globalDofIdx];
765 }
766
767
768 const std::unique_ptr<EclWriterType>& eclWriter() const
769 {
770 return eclWriter_;
771 }
772
773 void setConvData(const std::vector<std::vector<int>>& data)
774 {
775 eclWriter_->mutableOutputModule().setCnvData(data);
776 }
777
782 Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
783 {
784 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
785 this->episodeIndex(),
786 this->pvtRegionIndex(globalDofIdx));
787 }
788
793 Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
794 {
795 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
796 this->episodeIndex(),
797 this->pvtRegionIndex(globalDofIdx));
798 }
799
809 {
810 int episodeIdx = this->episodeIndex();
811 return !this->mixControls_.drsdtActive(episodeIdx) &&
812 !this->mixControls_.drvdtActive(episodeIdx) &&
813 this->rockCompPoroMultWc_.empty() &&
814 this->rockCompPoroMult_.empty();
815 }
816
823 template <class Context>
824 void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
825 {
826 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
827
828 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
829 values.assignNaive(initialFluidStates_[globalDofIdx]);
830
831 SolventModule::assignPrimaryVars(values,
832 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
833 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
834
835 if constexpr (enablePolymer)
836 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
837
838 if constexpr (enablePolymerMolarWeight)
839 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
840
841 if constexpr (enableBrine) {
842 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
843 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
844 }
845 else {
846 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
847 }
848 }
849
850 if constexpr (enableMICP){
851 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
852 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
853 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
854 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
855 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
856 }
857
858 values.checkDefined();
859 }
860
861
862 Scalar drsdtcon(unsigned elemIdx, int episodeIdx) const
863 {
864 return this->mixControls_.drsdtcon(elemIdx, episodeIdx,
865 this->pvtRegionIndex(elemIdx));
866 }
867
873 template <class Context>
874 void boundary(BoundaryRateVector& values,
875 const Context& context,
876 unsigned spaceIdx,
877 unsigned timeIdx) const
878 {
879 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
880 if (!context.intersection(spaceIdx).boundary())
881 return;
882
883 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
884 values.setNoFlow();
885 else {
886 // in the energy case we need to specify a non-trivial boundary condition
887 // because the geothermal gradient needs to be maintained. for this, we
888 // simply assume the initial temperature at the boundary and specify the
889 // thermal flow accordingly. in this context, "thermal flow" means energy
890 // flow due to a temerature gradient while assuming no-flow for mass
891 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
892 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
893 values.setThermalFlow(context, spaceIdx, timeIdx, this->initialFluidStates_[globalDofIdx] );
894 }
895
896 if (this->nonTrivialBoundaryConditions()) {
897 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
898 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
899 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
900 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
901 const auto [type, massrate] = this->boundaryCondition(globalDofIdx, indexInInside);
902 if (type == BCType::THERMAL)
903 values.setThermalFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
904 else if (type == BCType::FREE || type == BCType::DIRICHLET)
905 values.setFreeFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
906 else if (type == BCType::RATE)
907 values.setMassRate(massrate, pvtRegionIdx);
908 }
909 }
910
911
915 Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
916 { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
917
918 const FlowThresholdPressure<TypeTag>& thresholdPressure() const
919 { return thresholdPressures_; }
920
921 FlowThresholdPressure<TypeTag>& thresholdPressure()
922 { return thresholdPressures_; }
923
924 const ModuleParams& moduleParams() const
925 {
926 return moduleParams_;
927 }
928
929 template<class Serializer>
930 void serializeOp(Serializer& serializer)
931 {
932 serializer(static_cast<FlowProblemType&>(*this));
933 serializer(mixControls_);
934 serializer(*eclWriter_);
935 }
936
937protected:
938 void updateExplicitQuantities_(int episodeIdx, int timeStepSize, const bool first_step_after_restart) override
939 {
940 this->updateExplicitQuantities_(first_step_after_restart);
941
942 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
943 updateMaxPolymerAdsorption_();
944
945 mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize);
946 }
947
948 void updateMaxPolymerAdsorption_()
949 {
950 // we need to update the max polymer adsoption data for all elements
951 this->updateProperty_("FlowProblemBlackoil::updateMaxPolymerAdsorption_() failed:",
952 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
953 {
954 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
955 });
956 }
957
958 bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
959 {
960 const Scalar pa = scalarValue(iq.polymerAdsorption());
961 auto& mpa = this->polymer_.maxAdsorption;
962 if (mpa[compressedDofIdx] < pa) {
963 mpa[compressedDofIdx] = pa;
964 return true;
965 } else {
966 return false;
967 }
968 }
969
970 void computeAndSetEqWeights_()
971 {
972 std::vector<Scalar> sumInvB(numPhases, 0.0);
973 const auto& gridView = this->gridView();
974 ElementContext elemCtx(this->simulator());
975 for(const auto& elem: elements(gridView, Dune::Partitions::interior)) {
976 elemCtx.updatePrimaryStencil(elem);
977 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
978 const auto& dofFluidState = this->initialFluidStates_[elemIdx];
979 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
980 if (!FluidSystem::phaseIsActive(phaseIdx))
981 continue;
982
983 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
984 }
985 }
986
987 std::size_t numDof = this->model().numGridDof();
988 const auto& comm = this->simulator().vanguard().grid().comm();
989 comm.sum(sumInvB.data(),sumInvB.size());
990 Scalar numTotalDof = comm.sum(numDof);
991
992 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
993 if (!FluidSystem::phaseIsActive(phaseIdx))
994 continue;
995
996 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
997 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
998 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
999 this->model().setEqWeight(activeSolventCompIdx, avgB);
1000 }
1001 }
1002
1003 // update the parameters needed for DRSDT and DRVDT
1004 bool updateCompositionChangeLimits_()
1005 {
1006 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1007 // update the "last Rs" values for all elements, including the ones in the ghost
1008 // and overlap regions
1009 int episodeIdx = this->episodeIndex();
1010 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1011 this->mixControls_.drsdtActive(episodeIdx),
1012 this->mixControls_.drvdtActive(episodeIdx)};
1013 if (!active[0] && !active[1] && !active[2]) {
1014 return false;
1015 }
1016
1017 this->updateProperty_("FlowProblemBlackoil::updateCompositionChangeLimits_()) failed:",
1018 [this,episodeIdx,active](unsigned compressedDofIdx,
1019 const IntensiveQuantities& iq)
1020 {
1021 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1022 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1023 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1024 this->mixControls_.update(compressedDofIdx,
1025 iq,
1026 episodeIdx,
1027 this->gravity_[dim - 1],
1028 perm[dim - 1][dim - 1],
1029 distZ,
1030 pvtRegionIdx);
1031 }
1032 );
1033
1034 return true;
1035 }
1036
1037
1038 void readEclRestartSolution_()
1039 {
1040 // Throw an exception if the grid has LGRs. Refined grid are not supported for restart.
1041 if(this->simulator().vanguard().grid().maxLevel() > 0) {
1042 throw std::invalid_argument("Refined grids are not yet supported for restart ");
1043 }
1044
1045 // Set the start time of the simulation
1046 auto& simulator = this->simulator();
1047 const auto& schedule = simulator.vanguard().schedule();
1048 const auto& eclState = simulator.vanguard().eclState();
1049 const auto& initconfig = eclState.getInitConfig();
1050 const int restart_step = initconfig.getRestartStep();
1051 {
1052 simulator.setTime(schedule.seconds(restart_step));
1053
1054 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
1055 schedule.stepLength(restart_step));
1056 simulator.setEpisodeIndex(restart_step);
1057 }
1058 this->eclWriter_->beginRestart();
1059
1060 Scalar dt = std::min(this->eclWriter_->restartTimeStepSize(), simulator.episodeLength());
1061 simulator.setTimeStepSize(dt);
1062
1063 std::size_t numElems = this->model().numGridDof();
1064 this->initialFluidStates_.resize(numElems);
1065 if constexpr (enableSolvent) {
1066 this->solventSaturation_.resize(numElems, 0.0);
1067 this->solventRsw_.resize(numElems, 0.0);
1068 }
1069
1070 if constexpr (enablePolymer)
1071 this->polymer_.concentration.resize(numElems, 0.0);
1072
1073 if constexpr (enablePolymerMolarWeight) {
1074 const std::string msg {"Support of the RESTART for polymer molecular weight "
1075 "is not implemented yet. The polymer weight value will be "
1076 "zero when RESTART begins"};
1077 OpmLog::warning("NO_POLYMW_RESTART", msg);
1078 this->polymer_.moleWeight.resize(numElems, 0.0);
1079 }
1080
1081 if constexpr (enableMICP) {
1082 this->micp_.resize(numElems);
1083 }
1084
1085 // Initialize mixing controls before trying to set any lastRx valuesx
1086 this->mixControls_.init(numElems, restart_step, eclState.runspec().tabdims().getNumPVTTables());
1087
1088 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1089 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1090 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
1091 this->eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
1092 this->eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
1093
1094 // Note: Function processRestartSaturations_() mutates the
1095 // 'ssol' argument--the value from the restart file--if solvent
1096 // is enabled. Then, store the updated solvent saturation into
1097 // 'solventSaturation_'. Otherwise, just pass a dummy value to
1098 // the function and discard the unchanged result. Do not index
1099 // into 'solventSaturation_' unless solvent is enabled.
1100 {
1101 auto ssol = enableSolvent
1102 ? this->eclWriter_->outputModule().getSolventSaturation(elemIdx)
1103 : Scalar(0);
1104
1105 this->processRestartSaturations_(elemFluidState, ssol);
1106
1107 if constexpr (enableSolvent) {
1108 this->solventSaturation_[elemIdx] = ssol;
1109 this->solventRsw_[elemIdx] = this->eclWriter_->outputModule().getSolventRsw(elemIdx);
1110 }
1111 }
1112
1113 // For CO2STORE and H2STORE we need to set the initial temperature for isothermal simulations
1114 bool isThermal = eclState.getSimulationConfig().isThermal();
1115 bool needTemperature = (eclState.runspec().co2Storage() || eclState.runspec().h2Storage());
1116 if (!isThermal && needTemperature) {
1117 const auto& fp = simulator.vanguard().eclState().fieldProps();
1118 elemFluidState.setTemperature(fp.get_double("TEMPI")[elemIdx]);
1119 }
1120
1121 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
1122
1123 if constexpr (enablePolymer)
1124 this->polymer_.concentration[elemIdx] = this->eclWriter_->outputModule().getPolymerConcentration(elemIdx);
1125 if constexpr (enableMICP){
1126 this->micp_.microbialConcentration[elemIdx] = this->eclWriter_->outputModule().getMicrobialConcentration(elemIdx);
1127 this->micp_.oxygenConcentration[elemIdx] = this->eclWriter_->outputModule().getOxygenConcentration(elemIdx);
1128 this->micp_.ureaConcentration[elemIdx] = this->eclWriter_->outputModule().getUreaConcentration(elemIdx);
1129 this->micp_.biofilmConcentration[elemIdx] = this->eclWriter_->outputModule().getBiofilmConcentration(elemIdx);
1130 this->micp_.calciteConcentration[elemIdx] = this->eclWriter_->outputModule().getCalciteConcentration(elemIdx);
1131 }
1132 // if we need to restart for polymer molecular weight simulation, we need to add related here
1133 }
1134
1135 const int episodeIdx = this->episodeIndex();
1136 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
1137
1138 // assign the restart solution to the current solution. note that we still need
1139 // to compute real initial solution after this because the initial fluid states
1140 // need to be correct for stuff like boundary conditions.
1141 auto& sol = this->model().solution(/*timeIdx=*/0);
1142 const auto& gridView = this->gridView();
1143 ElementContext elemCtx(simulator);
1144 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1145 elemCtx.updatePrimaryStencil(elem);
1146 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1147 this->initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0);
1148 }
1149
1150 // make sure that the ghost and overlap entities exhibit the correct
1151 // solution. alternatively, this could be done in the loop above by also
1152 // considering non-interior elements. Since the initial() method might not work
1153 // 100% correctly for such elements, let's play safe and explicitly synchronize
1154 // using message passing.
1155 this->model().syncOverlap();
1156
1157 this->eclWriter_->endRestart();
1158 }
1159
1160 void readEquilInitialCondition_() override
1161 {
1162 const auto& simulator = this->simulator();
1163
1164 // initial condition corresponds to hydrostatic conditions.
1165 EquilInitializer<TypeTag> equilInitializer(simulator, *(this->materialLawManager_));
1166
1167 std::size_t numElems = this->model().numGridDof();
1168 this->initialFluidStates_.resize(numElems);
1169 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1170 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1171 elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
1172 }
1173 }
1174
1175 void readExplicitInitialCondition_() override
1176 {
1177 const auto& simulator = this->simulator();
1178 const auto& vanguard = simulator.vanguard();
1179 const auto& eclState = vanguard.eclState();
1180 const auto& fp = eclState.fieldProps();
1181 bool has_swat = fp.has_double("SWAT");
1182 bool has_sgas = fp.has_double("SGAS");
1183 bool has_rs = fp.has_double("RS");
1184 bool has_rv = fp.has_double("RV");
1185 bool has_rvw = fp.has_double("RVW");
1186 bool has_pressure = fp.has_double("PRESSURE");
1187 bool has_salt = fp.has_double("SALT");
1188 bool has_saltp = fp.has_double("SALTP");
1189
1190 // make sure all required quantities are enables
1191 if (Indices::numPhases > 1) {
1192 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
1193 throw std::runtime_error("The ECL input file requires the presence of the SWAT keyword if "
1194 "the water phase is active");
1195 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
1196 throw std::runtime_error("The ECL input file requires the presence of the SGAS keyword if "
1197 "the gas phase is active");
1198 }
1199 if (!has_pressure)
1200 throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
1201 "keyword if the model is initialized explicitly");
1202 if (FluidSystem::enableDissolvedGas() && !has_rs)
1203 throw std::runtime_error("The ECL input file requires the RS keyword to be present if"
1204 " dissolved gas is enabled");
1205 if (FluidSystem::enableVaporizedOil() && !has_rv)
1206 throw std::runtime_error("The ECL input file requires the RV keyword to be present if"
1207 " vaporized oil is enabled");
1208 if (FluidSystem::enableVaporizedWater() && !has_rvw)
1209 throw std::runtime_error("The ECL input file requires the RVW keyword to be present if"
1210 " vaporized water is enabled");
1211 if (enableBrine && !has_salt)
1212 throw std::runtime_error("The ECL input file requires the SALT keyword to be present if"
1213 " brine is enabled and the model is initialized explicitly");
1214 if (enableSaltPrecipitation && !has_saltp)
1215 throw std::runtime_error("The ECL input file requires the SALTP keyword to be present if"
1216 " salt precipitation is enabled and the model is initialized explicitly");
1217
1218 std::size_t numDof = this->model().numGridDof();
1219
1220 initialFluidStates_.resize(numDof);
1221
1222 std::vector<double> waterSaturationData;
1223 std::vector<double> gasSaturationData;
1224 std::vector<double> pressureData;
1225 std::vector<double> rsData;
1226 std::vector<double> rvData;
1227 std::vector<double> rvwData;
1228 std::vector<double> tempiData;
1229 std::vector<double> saltData;
1230 std::vector<double> saltpData;
1231
1232 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
1233 waterSaturationData = fp.get_double("SWAT");
1234 else
1235 waterSaturationData.resize(numDof);
1236
1237 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
1238 gasSaturationData = fp.get_double("SGAS");
1239 else
1240 gasSaturationData.resize(numDof);
1241
1242 pressureData = fp.get_double("PRESSURE");
1243 if (FluidSystem::enableDissolvedGas())
1244 rsData = fp.get_double("RS");
1245
1246 if (FluidSystem::enableVaporizedOil())
1247 rvData = fp.get_double("RV");
1248
1249 if (FluidSystem::enableVaporizedWater())
1250 rvwData = fp.get_double("RVW");
1251
1252 // initial reservoir temperature
1253 tempiData = fp.get_double("TEMPI");
1254
1255 // initial salt concentration data
1256 if constexpr (enableBrine)
1257 saltData = fp.get_double("SALT");
1258
1259 // initial precipitated salt saturation data
1260 if constexpr (enableSaltPrecipitation)
1261 saltpData = fp.get_double("SALTP");
1262
1263 // calculate the initial fluid states
1264 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1265 auto& dofFluidState = initialFluidStates_[dofIdx];
1266
1267 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
1268
1270 // set temperature
1272 Scalar temperatureLoc = tempiData[dofIdx];
1273 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
1274 temperatureLoc = FluidSystem::surfaceTemperature;
1275 dofFluidState.setTemperature(temperatureLoc);
1276
1278 // set salt concentration
1280 if constexpr (enableBrine)
1281 dofFluidState.setSaltConcentration(saltData[dofIdx]);
1282
1284 // set precipitated salt saturation
1286 if constexpr (enableSaltPrecipitation)
1287 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
1288
1290 // set saturations
1292 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1293 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
1294 waterSaturationData[dofIdx]);
1295
1296 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
1297 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
1298 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1299 1.0
1300 - waterSaturationData[dofIdx]);
1301 }
1302 else
1303 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1304 gasSaturationData[dofIdx]);
1305 }
1306 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx))
1307 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
1308 1.0
1309 - waterSaturationData[dofIdx]
1310 - gasSaturationData[dofIdx]);
1311
1313 // set phase pressures
1315 Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
1316
1317 // this assumes that capillary pressures only depend on the phase saturations
1318 // and possibly on temperature. (this is always the case for ECL problems.)
1319 std::array<Scalar, numPhases> pc = {0};
1320 const auto& matParams = this->materialLawParams(dofIdx);
1321 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
1322 Valgrind::CheckDefined(pressure);
1323 Valgrind::CheckDefined(pc);
1324 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1325 if (!FluidSystem::phaseIsActive(phaseIdx))
1326 continue;
1327
1328 if (Indices::oilEnabled)
1329 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1330 else if (Indices::gasEnabled)
1331 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1332 else if (Indices::waterEnabled)
1333 //single (water) phase
1334 dofFluidState.setPressure(phaseIdx, pressure);
1335 }
1336
1337 if (FluidSystem::enableDissolvedGas())
1338 dofFluidState.setRs(rsData[dofIdx]);
1339 else if (Indices::gasEnabled && Indices::oilEnabled)
1340 dofFluidState.setRs(0.0);
1341
1342 if (FluidSystem::enableVaporizedOil())
1343 dofFluidState.setRv(rvData[dofIdx]);
1344 else if (Indices::gasEnabled && Indices::oilEnabled)
1345 dofFluidState.setRv(0.0);
1346
1347 if (FluidSystem::enableVaporizedWater())
1348 dofFluidState.setRvw(rvwData[dofIdx]);
1349
1351 // set invB_
1353 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1354 if (!FluidSystem::phaseIsActive(phaseIdx))
1355 continue;
1356
1357 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1358 dofFluidState.setInvB(phaseIdx, b);
1359
1360 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1361 dofFluidState.setDensity(phaseIdx, rho);
1362
1363 }
1364 }
1365 }
1366
1367
1368 void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
1369 {
1370 // each phase needs to be above certain value to be claimed to be existing
1371 // this is used to recover some RESTART running with the defaulted single-precision format
1372 const Scalar smallSaturationTolerance = 1.e-6;
1373 Scalar sumSaturation = 0.0;
1374 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1375 if (FluidSystem::phaseIsActive(phaseIdx)) {
1376 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance)
1377 elemFluidState.setSaturation(phaseIdx, 0.0);
1378
1379 sumSaturation += elemFluidState.saturation(phaseIdx);
1380 }
1381
1382 }
1383 if constexpr (enableSolvent) {
1384 if (solventSaturation < smallSaturationTolerance)
1385 solventSaturation = 0.0;
1386
1387 sumSaturation += solventSaturation;
1388 }
1389
1390 assert(sumSaturation > 0.0);
1391
1392 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1393 if (FluidSystem::phaseIsActive(phaseIdx)) {
1394 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
1395 elemFluidState.setSaturation(phaseIdx, saturation);
1396 }
1397 }
1398 if constexpr (enableSolvent) {
1399 solventSaturation = solventSaturation / sumSaturation;
1400 }
1401 }
1402
1403 void readInitialCondition_() override
1404 {
1405 FlowProblemType::readInitialCondition_();
1406
1407 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
1408 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
1409 enableSolvent,
1410 enablePolymer,
1411 enablePolymerMolarWeight,
1412 enableMICP);
1413
1414 }
1415
1416 void handleSolventBC(const BCProp::BCFace& bc, RateVector& rate) const override
1417 {
1418 if constexpr (!enableSolvent)
1419 throw std::logic_error("solvent is disabled and you're trying to add solvent to BC");
1420
1421 rate[Indices::solventSaturationIdx] = bc.rate;
1422 }
1423
1424 void handlePolymerBC(const BCProp::BCFace& bc, RateVector& rate) const override
1425 {
1426 if constexpr (!enablePolymer)
1427 throw std::logic_error("polymer is disabled and you're trying to add polymer to BC");
1428
1429 rate[Indices::polymerConcentrationIdx] = bc.rate;
1430 }
1431
1432 void updateExplicitQuantities_(const bool first_step_after_restart)
1433 {
1434 OPM_TIMEBLOCK(updateExplicitQuantities);
1435 const bool invalidateFromMaxWaterSat = this->updateMaxWaterSaturation_();
1436 const bool invalidateFromMinPressure = this->updateMinPressure_();
1437
1438 // update hysteresis and max oil saturation used in vappars
1439 const bool invalidateFromHyst = this->updateHysteresis_();
1440 const bool invalidateFromMaxOilSat = this->updateMaxOilSaturation_();
1441
1442 // deal with DRSDT and DRVDT
1443 const bool invalidateDRDT = !first_step_after_restart && this->updateCompositionChangeLimits_();
1444
1445 // the derivatives may have changed
1446 const bool invalidateIntensiveQuantities
1447 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat || invalidateDRDT;
1448 if (invalidateIntensiveQuantities) {
1449 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1450 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
1451 }
1452
1453 this->updateRockCompTransMultVal_();
1454 }
1455
1456 FlowThresholdPressure<TypeTag> thresholdPressures_;
1457
1458 std::vector<InitialFluidState> initialFluidStates_;
1459
1460 bool enableEclOutput_;
1461 std::unique_ptr<EclWriterType> eclWriter_;
1462#if HAVE_DAMARIS
1463 bool enableDamarisOutput_ = false ;
1464 std::unique_ptr<DamarisWriterType> damarisWriter_;
1465#endif
1466 MixingRateControls<FluidSystem> mixControls_;
1467
1468 ActionHandler<Scalar> actionHandler_;
1469
1470 ModuleParams moduleParams_;
1471};
1472
1473} // namespace Opm
1474
1475#endif // OPM_FLOW_PROBLEM_BLACK_HPP
Collects necessary output values and pass them to Damaris server processes.
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 problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
VTK output module for the tracer model's parameters.
Calculates the local residual of the black oil model.
static void setParams(BlackOilBrineParams< Scalar > &&params)
Set parameters.
Definition blackoilbrinemodules.hh:81
static void setParams(BlackOilExtboParams< Scalar > &&params)
Set parameters.
Definition blackoilextbomodules.hh:90
static void setParams(BlackOilFoamParams< Scalar > &&params)
Set parameters.
Definition blackoilfoammodules.hh:92
static void setParams(BlackOilMICPParams< Scalar > &&params)
Set parameters.
Definition blackoilmicpmodules.hh:83
static void setParams(BlackOilPolymerParams< Scalar > &&params)
Set parameters.
Definition blackoilpolymermodules.hh:88
static void setParams(BlackOilSolventParams< Scalar > &&params)
Set parameters.
Definition blackoilsolventmodules.hh:90
Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
Returns the maximum value of the gas dissolution factor at the current time for a given degree of fre...
Definition FlowProblemBlackoil.hpp:782
FlowProblemBlackoil(Simulator &simulator)
Definition FlowProblemBlackoil.hpp:165
Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
Returns the maximum value of the oil vaporization factor at the current time for a given degree of fr...
Definition FlowProblemBlackoil.hpp:793
void endTimeStep() override
Called by the simulator after each time integration.
Definition FlowProblemBlackoil.hpp:398
void endEpisode() override
Called by the simulator after the end of an episode.
Definition FlowProblemBlackoil.hpp:452
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition FlowProblemBlackoil.hpp:824
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition FlowProblemBlackoil.hpp:249
LhsEval permFactTransMultiplier(const IntensiveQuantities &intQuants) const
Calculate the transmissibility multiplier due to porosity reduction.
Definition FlowProblemBlackoil.hpp:630
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition FlowProblemBlackoil.hpp:874
void initialSolutionApplied() override
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition FlowProblemBlackoil.hpp:514
void beginEpisode() override
Called by the simulator before an episode begins.
Definition FlowProblemBlackoil.hpp:222
bool recycleFirstIterationStorage() const
Return if the storage term of the first iteration is identical to the storage term for the solution o...
Definition FlowProblemBlackoil.hpp:808
void writeOutput(bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition FlowProblemBlackoil.hpp:483
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblemBlackoil.hpp:151
Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
Definition FlowProblemBlackoil.hpp:915
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:92
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:815
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:670
FlowProblem(Simulator &simulator)
Definition FlowProblem.hpp:209
virtual void beginEpisode()
Called by the simulator before an episode begins.
Definition FlowProblem.hpp:306
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblem.hpp:179
void writeOutput(bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:491
This class calculates the threshold pressure for grid faces according to the Eclipse Reference Manual...
Definition FlowThresholdPressure.hpp:59
Definition SimulatorTimer.hpp:39
VTK output module for the tracer model's parameters.
Definition VtkTracerModule.hpp:57
static void registerParameters()
Register all run-time parameters for the tracer VTK output module.
Definition VtkTracerModule.hpp:83
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Struct holding the parameters for the BlackoilBrineModule class.
Definition blackoilbrineparams.hpp:44
Struct holding the parameters for the BlackoilExtboModule class.
Definition blackoilextboparams.hpp:49
Struct holding the parameters for the BlackoilFoamModule class.
Definition blackoilfoamparams.hpp:46
Struct holding the parameters for the BlackOilMICPModule class.
Definition blackoilmicpparams.hpp:42
Struct holding the parameters for the BlackOilPolymerModule class.
Definition blackoilpolymerparams.hpp:45
Struct holding the parameters for the BlackOilSolventModule class.
Definition blackoilsolventparams.hpp:49