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