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