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 // Deal with "clogging" for the MICP model
491 if constexpr (enableMICP) {
492 auto& model = this->model();
493 const auto& residual = model.linearizer().residual();
494
495 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); ++globalDofIdx) {
496 auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx];
497 MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
498 }
499 }
500 }
501
505 void endEpisode() override
506 {
507 OPM_TIMEBLOCK(endEpisode);
508
509 // Rerun UDQ assignents following action processing on the final
510 // time step of this episode to make sure that any UDQ ASSIGN
511 // operations triggered in action blocks take effect. This is
512 // mainly to work around a shortcoming of the ScheduleState copy
513 // constructor which clears pending UDQ assignments under the
514 // assumption that all such assignments have been processed. If an
515 // action block happens to trigger on the final time step of an
516 // episode and that action block runs a UDQ assignment, then that
517 // assignment would be dropped and the rest of the simulator will
518 // never see its effect without this hack.
519 this->actionHandler_
520 .evalUDQAssignments(this->episodeIndex(), this->simulator().vanguard().udqState());
521
522 FlowProblemType::endEpisode();
523 }
524
525 void writeReports(const SimulatorTimer& timer)
526 {
527 if (this->enableEclOutput_) {
528 this->eclWriter_->writeReports(timer);
529 }
530 }
531
532
537 void writeOutput(bool verbose) override
538 {
539 FlowProblemType::writeOutput(verbose);
540
541 bool isSubStep = !this->simulator().episodeWillBeOver();
542
543 data::Solution localCellData = {};
544#if HAVE_DAMARIS
545 // N.B. the Damaris output has to be done before the ECL output as the ECL one
546 // does all kinds of std::move() relocation of data
547 if (enableDamarisOutput_) {
548 damarisWriter_->writeOutput(localCellData, isSubStep) ;
549 }
550#endif
551 if (enableEclOutput_){
552 eclWriter_->writeOutput(std::move(localCellData), isSubStep);
553 }
554 }
555
556 void finalizeOutput()
557 {
558 OPM_TIMEBLOCK(finalizeOutput);
559 // this will write all pending output to disk
560 // to avoid corruption of output files
561 eclWriter_.reset();
562 }
563
564
569 {
570 FlowProblemType::initialSolutionApplied();
571
572 // let the object for threshold pressures initialize itself. this is done only at
573 // this point, because determining the threshold pressures may require to access
574 // the initial solution.
575 this->thresholdPressures_.finishInit();
576
577 // For CpGrid with LGRs, ecl-output is not supported yet.
578 const auto& grid = this->simulator().vanguard().gridView().grid();
579
580 using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
581 constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
582 // Skip - for now - calculate the initial fip values for CpGrid with LGRs.
583 if (!isCpGrid || (grid.maxLevel() == 0)) {
584 if (this->simulator().episodeIndex() == 0) {
585 eclWriter_->writeInitialFIPReport();
586 }
587 }
588 }
589
590 void addToSourceDense(RateVector& rate,
591 unsigned globalDofIdx,
592 unsigned timeIdx) const override
593 {
594 this->aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
595
596 // Add source term from deck
597 const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source();
598 std::array<int,3> ijk;
599 this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
600
601 if (source.hasSource(ijk)) {
602 const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
603 static std::array<SourceComponent, 3> sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS};
604 static std::array<int, 3> phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
605 static std::array<int, 3> cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx};
606
607 for (unsigned i = 0; i < phidx_map.size(); ++i) {
608 const auto phaseIdx = phidx_map[i];
609 const auto sourceComp = sc_map[i];
610 const auto compIdx = cidx_map[i];
611 if (!FluidSystem::phaseIsActive(phaseIdx)) {
612 continue;
613 }
614 Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
615 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
616 mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
617 }
618 rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate;
619 }
620
621 if constexpr (enableSolvent) {
622 Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
623 if constexpr (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
624 const auto& solventPvt = SolventModule::solventPvt();
625 mass_rate /= solventPvt.referenceDensity(pvtRegionIdx);
626 }
627 rate[Indices::contiSolventEqIdx] += mass_rate;
628 }
629 if constexpr (enablePolymer) {
630 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
631 }
632 if constexpr (enableEnergy) {
633 for (unsigned i = 0; i < phidx_map.size(); ++i) {
634 const auto phaseIdx = phidx_map[i];
635 if (!FluidSystem::phaseIsActive(phaseIdx)) {
636 continue;
637 }
638 const auto sourceComp = sc_map[i];
639 if (source.hasHrate({ijk, sourceComp})) {
640 rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx);
641 } else {
642 const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, /*timeIdx*/ 0);
643 auto fs = intQuants.fluidState();
644 // if temperature is not set, use cell temperature as default
645 if (source.hasTemperature({ijk, sourceComp})) {
646 Scalar temperature = source.temperature({ijk, sourceComp});
647 fs.setTemperature(temperature);
648 }
649 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx);
650 Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx);
651 Scalar energy_rate = getValue(h)*mass_rate;
652 rate[Indices::contiEnergyEqIdx] += energy_rate;
653 }
654 }
655 }
656 }
657
658 // if requested, compensate systematic mass loss for cells which were "well
659 // behaved" in the last time step
660 if (this->enableDriftCompensation_) {
661 const auto& simulator = this->simulator();
662 const auto& model = this->model();
663
664 // we use a lower tolerance for the compensation too
665 // assure the added drift from the last step does not
666 // cause convergence issues on the current step
667 Scalar maxCompensation = model.newtonMethod().tolerance()/10;
668 Scalar poro = this->porosity(globalDofIdx, timeIdx);
669 Scalar dt = simulator.timeStepSize();
670 EqVector dofDriftRate = this->drift_[globalDofIdx];
671 dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
672
673 // restrict drift compensation to the CNV tolerance
674 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
675 Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
676 if (cnv > maxCompensation) {
677 dofDriftRate[eqIdx] *= maxCompensation/cnv;
678 }
679 }
680
681 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
682 rate[eqIdx] -= dofDriftRate[eqIdx];
683 }
684 }
685
691 template <class LhsEval>
692 LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
693 {
694 OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
695 if (!enableSaltPrecipitation)
696 return 1.0;
697
698 const auto& fs = intQuants.fluidState();
699 unsigned tableIdx = this->simulator().problem().satnumRegionIndex(elementIdx);
700
701 LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
702 porosityFactor = min(porosityFactor, 1.0);
703 const auto& permfactTable = BrineModule::permfactTable(tableIdx);
704 return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
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::NONE:
767 throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
768 }
769 fluidState.setTotalSaturation(1.0);
770 double pressure = initialFluidStates_[globalDofIdx].pressure(this->refPressurePhaseIdx_());
771 const auto pressure_input = bc.pressure;
772 if (pressure_input) {
773 pressure = *pressure_input;
774 }
775
776 std::array<Scalar, numPhases> pc = {0};
777 const auto& matParams = this->materialLawParams(globalDofIdx);
778 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
779 Valgrind::CheckDefined(pressure);
780 Valgrind::CheckDefined(pc);
781 for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
782 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
783
784 fluidState.setPc(phaseIdx, pc[phaseIdx]);
785 if (Indices::oilEnabled)
786 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
787 else if (Indices::gasEnabled)
788 fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
789 else if (Indices::waterEnabled)
790 //single (water) phase
791 fluidState.setPressure(phaseIdx, pressure);
792 }
793
794 double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature
795 const auto temperature_input = bc.temperature;
796 if(temperature_input)
797 temperature = *temperature_input;
798 fluidState.setTemperature(temperature);
799
800 if (FluidSystem::enableDissolvedGas()) {
801 fluidState.setRs(0.0);
802 fluidState.setRv(0.0);
803 }
804 if (FluidSystem::enableDissolvedGasInWater()) {
805 fluidState.setRsw(0.0);
806 }
807 if (FluidSystem::enableVaporizedWater())
808 fluidState.setRvw(0.0);
809
810 for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) {
811 const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx);
812
813 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
814 fluidState.setInvB(phaseIdx, b);
815
816 const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
817 fluidState.setDensity(phaseIdx, rho);
818 if (enableEnergy) {
819 const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
820 fluidState.setEnthalpy(phaseIdx, h);
821 }
822 }
823 fluidState.checkDefined();
824 return fluidState;
825 }
826 }
827 return initialFluidStates_[globalDofIdx];
828 }
829
830
831 const EclWriterType& eclWriter() const
832 { return *eclWriter_; }
833
834 EclWriterType& eclWriter()
835 { return *eclWriter_; }
836
841 Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const
842 {
843 return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx,
844 this->episodeIndex(),
845 this->pvtRegionIndex(globalDofIdx));
846 }
847
852 Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const
853 {
854 return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx,
855 this->episodeIndex(),
856 this->pvtRegionIndex(globalDofIdx));
857 }
858
868 {
869 int episodeIdx = this->episodeIndex();
870 return !this->mixControls_.drsdtActive(episodeIdx) &&
871 !this->mixControls_.drvdtActive(episodeIdx) &&
872 this->rockCompPoroMultWc_.empty() &&
873 this->rockCompPoroMult_.empty();
874 }
875
882 template <class Context>
883 void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
884 {
885 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
886
887 values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx));
888 values.assignNaive(initialFluidStates_[globalDofIdx]);
889
890 SolventModule::assignPrimaryVars(values,
891 enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0,
892 enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0);
893
894 if constexpr (enablePolymer)
895 values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx];
896
897 if constexpr (enablePolymerMolarWeight)
898 values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx];
899
900 if constexpr (enableBrine) {
901 if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
902 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
903 }
904 else {
905 values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
906 }
907 }
908
909 if constexpr (enableMICP){
910 values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx];
911 values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx];
912 values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx];
913 values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx];
914 values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx];
915 }
916
917 values.checkDefined();
918 }
919
920
921 Scalar drsdtcon(unsigned elemIdx, int episodeIdx) const
922 {
923 return this->mixControls_.drsdtcon(elemIdx, episodeIdx,
924 this->pvtRegionIndex(elemIdx));
925 }
926
932 template <class Context>
933 void boundary(BoundaryRateVector& values,
934 const Context& context,
935 unsigned spaceIdx,
936 unsigned timeIdx) const
937 {
938 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
939 if (!context.intersection(spaceIdx).boundary())
940 return;
941
942 if constexpr (!enableEnergy || !enableThermalFluxBoundaries)
943 values.setNoFlow();
944 else {
945 // in the energy case we need to specify a non-trivial boundary condition
946 // because the geothermal gradient needs to be maintained. for this, we
947 // simply assume the initial temperature at the boundary and specify the
948 // thermal flow accordingly. in this context, "thermal flow" means energy
949 // flow due to a temerature gradient while assuming no-flow for mass
950 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
951 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
952 values.setThermalFlow(context, spaceIdx, timeIdx, this->initialFluidStates_[globalDofIdx] );
953 }
954
955 if (this->nonTrivialBoundaryConditions()) {
956 unsigned indexInInside = context.intersection(spaceIdx).indexInInside();
957 unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx);
958 unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx);
959 unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx);
960 const auto [type, massrate] = this->boundaryCondition(globalDofIdx, indexInInside);
961 if (type == BCType::THERMAL)
962 values.setThermalFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
963 else if (type == BCType::FREE || type == BCType::DIRICHLET)
964 values.setFreeFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside));
965 else if (type == BCType::RATE)
966 values.setMassRate(massrate, pvtRegionIdx);
967 }
968 }
969
974 void readSolutionFromOutputModule(const int restart_step, bool fip_init)
975 {
976 auto& simulator = this->simulator();
977 const auto& eclState = simulator.vanguard().eclState();
978
979 std::size_t numElems = this->model().numGridDof();
980 this->initialFluidStates_.resize(numElems);
981 if constexpr (enableSolvent) {
982 this->solventSaturation_.resize(numElems, 0.0);
983 this->solventRsw_.resize(numElems, 0.0);
984 }
985
986 if constexpr (enablePolymer)
987 this->polymer_.concentration.resize(numElems, 0.0);
988
989 if constexpr (enablePolymerMolarWeight) {
990 const std::string msg {"Support of the RESTART for polymer molecular weight "
991 "is not implemented yet. The polymer weight value will be "
992 "zero when RESTART begins"};
993 OpmLog::warning("NO_POLYMW_RESTART", msg);
994 this->polymer_.moleWeight.resize(numElems, 0.0);
995 }
996
997 if constexpr (enableMICP) {
998 this->micp_.resize(numElems);
999 }
1000
1001 // Initialize mixing controls before trying to set any lastRx valuesx
1002 this->mixControls_.init(numElems, restart_step, eclState.runspec().tabdims().getNumPVTTables());
1003
1004 if constexpr (enableMICP) {
1005 this->micp_ = this->eclWriter_->outputModule().getMICP().getSolution();
1006 }
1007
1008 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1009 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1010 elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
1011 this->eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx);
1012 this->eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx);
1013
1014 // Note: Function processRestartSaturations_() mutates the
1015 // 'ssol' argument--the value from the restart file--if solvent
1016 // is enabled. Then, store the updated solvent saturation into
1017 // 'solventSaturation_'. Otherwise, just pass a dummy value to
1018 // the function and discard the unchanged result. Do not index
1019 // into 'solventSaturation_' unless solvent is enabled.
1020 {
1021 auto ssol = enableSolvent
1022 ? this->eclWriter_->outputModule().getSolventSaturation(elemIdx)
1023 : Scalar(0);
1024
1025 this->processRestartSaturations_(elemFluidState, ssol);
1026
1027 if constexpr (enableSolvent) {
1028 this->solventSaturation_[elemIdx] = ssol;
1029 this->solventRsw_[elemIdx] = this->eclWriter_->outputModule().getSolventRsw(elemIdx);
1030 }
1031 }
1032
1033 // For CO2STORE and H2STORE we need to set the initial temperature for isothermal simulations
1034 bool isThermal = eclState.getSimulationConfig().isThermal();
1035 bool needTemperature = (eclState.runspec().co2Storage() || eclState.runspec().h2Storage());
1036 if (!isThermal && needTemperature) {
1037 const auto& fp = simulator.vanguard().eclState().fieldProps();
1038 elemFluidState.setTemperature(fp.get_double("TEMPI")[elemIdx]);
1039 }
1040
1041 this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv());
1042
1043 if constexpr (enablePolymer)
1044 this->polymer_.concentration[elemIdx] = this->eclWriter_->outputModule().getPolymerConcentration(elemIdx);
1045 // if we need to restart for polymer molecular weight simulation, we need to add related here
1046 }
1047
1048 const int episodeIdx = this->episodeIndex();
1049 this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize());
1050
1051 // assign the restart solution to the current solution. note that we still need
1052 // to compute real initial solution after this because the initial fluid states
1053 // need to be correct for stuff like boundary conditions.
1054 auto& sol = this->model().solution(/*timeIdx=*/0);
1055 const auto& gridView = this->gridView();
1056 ElementContext elemCtx(simulator);
1057 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1058 elemCtx.updatePrimaryStencil(elem);
1059 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1060 this->initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0);
1061 }
1062
1063 // make sure that the ghost and overlap entities exhibit the correct
1064 // solution. alternatively, this could be done in the loop above by also
1065 // considering non-interior elements. Since the initial() method might not work
1066 // 100% correctly for such elements, let's play safe and explicitly synchronize
1067 // using message passing.
1068 this->model().syncOverlap();
1069
1070 if (fip_init) {
1071 this->updateReferencePorosity_();
1072 this->mixControls_.init(this->model().numGridDof(),
1073 this->episodeIndex(),
1074 eclState.runspec().tabdims().getNumPVTTables());
1075 }
1076 }
1077
1081 Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
1082 { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
1083
1084 const FlowThresholdPressure<TypeTag>& thresholdPressure() const
1085 { return thresholdPressures_; }
1086
1087 FlowThresholdPressure<TypeTag>& thresholdPressure()
1088 { return thresholdPressures_; }
1089
1090 const ModuleParams& moduleParams() const
1091 {
1092 return moduleParams_;
1093 }
1094
1095 template<class Serializer>
1096 void serializeOp(Serializer& serializer)
1097 {
1098 serializer(static_cast<FlowProblemType&>(*this));
1099 serializer(mixControls_);
1100 serializer(*eclWriter_);
1101 }
1102
1103protected:
1104 void updateExplicitQuantities_(int episodeIdx, int timeStepSize, const bool first_step_after_restart) override
1105 {
1106 this->updateExplicitQuantities_(first_step_after_restart);
1107
1108 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
1109 updateMaxPolymerAdsorption_();
1110
1111 mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize);
1112 }
1113
1114 void updateMaxPolymerAdsorption_()
1115 {
1116 // we need to update the max polymer adsoption data for all elements
1117 this->updateProperty_("FlowProblemBlackoil::updateMaxPolymerAdsorption_() failed:",
1118 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1119 {
1120 this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
1121 });
1122 }
1123
1124 bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1125 {
1126 const Scalar pa = scalarValue(iq.polymerAdsorption());
1127 auto& mpa = this->polymer_.maxAdsorption;
1128 if (mpa[compressedDofIdx] < pa) {
1129 mpa[compressedDofIdx] = pa;
1130 return true;
1131 } else {
1132 return false;
1133 }
1134 }
1135
1136 void computeAndSetEqWeights_()
1137 {
1138 std::vector<Scalar> sumInvB(numPhases, 0.0);
1139 const auto& gridView = this->gridView();
1140 ElementContext elemCtx(this->simulator());
1141 for(const auto& elem: elements(gridView, Dune::Partitions::interior)) {
1142 elemCtx.updatePrimaryStencil(elem);
1143 int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
1144 const auto& dofFluidState = this->initialFluidStates_[elemIdx];
1145 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1146 if (!FluidSystem::phaseIsActive(phaseIdx))
1147 continue;
1148
1149 sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
1150 }
1151 }
1152
1153 std::size_t numDof = this->model().numGridDof();
1154 const auto& comm = this->simulator().vanguard().grid().comm();
1155 comm.sum(sumInvB.data(),sumInvB.size());
1156 Scalar numTotalDof = comm.sum(numDof);
1157
1158 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1159 if (!FluidSystem::phaseIsActive(phaseIdx))
1160 continue;
1161
1162 Scalar avgB = numTotalDof / sumInvB[phaseIdx];
1163 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
1164 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
1165 this->model().setEqWeight(activeSolventCompIdx, avgB);
1166 }
1167 }
1168
1169 // update the parameters needed for DRSDT and DRVDT
1170 bool updateCompositionChangeLimits_()
1171 {
1172 OPM_TIMEBLOCK(updateCompositionChangeLimits);
1173 // update the "last Rs" values for all elements, including the ones in the ghost
1174 // and overlap regions
1175 int episodeIdx = this->episodeIndex();
1176 std::array<bool,3> active{this->mixControls_.drsdtConvective(episodeIdx),
1177 this->mixControls_.drsdtActive(episodeIdx),
1178 this->mixControls_.drvdtActive(episodeIdx)};
1179 if (!active[0] && !active[1] && !active[2]) {
1180 return false;
1181 }
1182
1183 this->updateProperty_("FlowProblemBlackoil::updateCompositionChangeLimits_()) failed:",
1184 [this,episodeIdx,active](unsigned compressedDofIdx,
1185 const IntensiveQuantities& iq)
1186 {
1187 const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx);
1188 const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0;
1189 const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
1190 this->mixControls_.update(compressedDofIdx,
1191 iq,
1192 episodeIdx,
1193 this->gravity_[dim - 1],
1194 perm[dim - 1][dim - 1],
1195 distZ,
1196 pvtRegionIdx);
1197 }
1198 );
1199
1200 return true;
1201 }
1202
1203 void readEclRestartSolution_()
1204 {
1205 // Throw an exception if the grid has LGRs. Refined grid are not supported for restart.
1206 if(this->simulator().vanguard().grid().maxLevel() > 0) {
1207 throw std::invalid_argument("Refined grids are not yet supported for restart ");
1208 }
1209
1210 // Set the start time of the simulation
1211 auto& simulator = this->simulator();
1212 const auto& schedule = simulator.vanguard().schedule();
1213 const auto& eclState = simulator.vanguard().eclState();
1214 const auto& initconfig = eclState.getInitConfig();
1215 const int restart_step = initconfig.getRestartStep();
1216 {
1217 simulator.setTime(schedule.seconds(restart_step));
1218
1219 simulator.startNextEpisode(simulator.startTime() + simulator.time(),
1220 schedule.stepLength(restart_step));
1221 simulator.setEpisodeIndex(restart_step);
1222 }
1223 this->eclWriter_->beginRestart();
1224
1225 Scalar dt = std::min(this->eclWriter_->restartTimeStepSize(), simulator.episodeLength());
1226 simulator.setTimeStepSize(dt);
1227
1228 this->readSolutionFromOutputModule(restart_step, false);
1229
1230 this->eclWriter_->endRestart();
1231 }
1232
1233 void readEquilInitialCondition_() override
1234 {
1235 const auto& simulator = this->simulator();
1236
1237 // initial condition corresponds to hydrostatic conditions.
1238 EquilInitializer<TypeTag> equilInitializer(simulator, *(this->materialLawManager_));
1239
1240 std::size_t numElems = this->model().numGridDof();
1241 this->initialFluidStates_.resize(numElems);
1242 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1243 auto& elemFluidState = this->initialFluidStates_[elemIdx];
1244 elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
1245 }
1246 }
1247
1248 void readExplicitInitialCondition_() override
1249 {
1250 const auto& simulator = this->simulator();
1251 const auto& vanguard = simulator.vanguard();
1252 const auto& eclState = vanguard.eclState();
1253 const auto& fp = eclState.fieldProps();
1254 bool has_swat = fp.has_double("SWAT");
1255 bool has_sgas = fp.has_double("SGAS");
1256 bool has_rs = fp.has_double("RS");
1257 bool has_rsw = fp.has_double("RSW");
1258 bool has_rv = fp.has_double("RV");
1259 bool has_rvw = fp.has_double("RVW");
1260 bool has_pressure = fp.has_double("PRESSURE");
1261 bool has_salt = fp.has_double("SALT");
1262 bool has_saltp = fp.has_double("SALTP");
1263
1264 // make sure all required quantities are enables
1265 if (Indices::numPhases > 1) {
1266 if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat)
1267 throw std::runtime_error("The ECL input file requires the presence of the SWAT keyword if "
1268 "the water phase is active");
1269 if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx))
1270 throw std::runtime_error("The ECL input file requires the presence of the SGAS keyword if "
1271 "the gas phase is active");
1272 }
1273 if (!has_pressure)
1274 throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
1275 "keyword if the model is initialized explicitly");
1276 if (FluidSystem::enableDissolvedGas() && !has_rs)
1277 throw std::runtime_error("The ECL input file requires the RS keyword to be present if"
1278 " dissolved gas is enabled and the model is initialized explicitly");
1279 if (FluidSystem::enableDissolvedGasInWater() && !has_rsw)
1280 OpmLog::warning("The model is initialized explicitly and the RSW keyword is not present in the"
1281 " ECL input file. The RSW values are set equal to 0");
1282 if (FluidSystem::enableVaporizedOil() && !has_rv)
1283 throw std::runtime_error("The ECL input file requires the RV keyword to be present if"
1284 " vaporized oil is enabled and the model is initialized explicitly");
1285 if (FluidSystem::enableVaporizedWater() && !has_rvw)
1286 throw std::runtime_error("The ECL input file requires the RVW keyword to be present if"
1287 " vaporized water is enabled and the model is initialized explicitly");
1288 if (enableBrine && !has_salt)
1289 throw std::runtime_error("The ECL input file requires the SALT keyword to be present if"
1290 " brine is enabled and the model is initialized explicitly");
1291 if (enableSaltPrecipitation && !has_saltp)
1292 throw std::runtime_error("The ECL input file requires the SALTP keyword to be present if"
1293 " salt precipitation is enabled and the model is initialized explicitly");
1294
1295 std::size_t numDof = this->model().numGridDof();
1296
1297 initialFluidStates_.resize(numDof);
1298
1299 std::vector<double> waterSaturationData;
1300 std::vector<double> gasSaturationData;
1301 std::vector<double> pressureData;
1302 std::vector<double> rsData;
1303 std::vector<double> rswData;
1304 std::vector<double> rvData;
1305 std::vector<double> rvwData;
1306 std::vector<double> tempiData;
1307 std::vector<double> saltData;
1308 std::vector<double> saltpData;
1309
1310 if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1)
1311 waterSaturationData = fp.get_double("SWAT");
1312 else
1313 waterSaturationData.resize(numDof);
1314
1315 if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx))
1316 gasSaturationData = fp.get_double("SGAS");
1317 else
1318 gasSaturationData.resize(numDof);
1319
1320 pressureData = fp.get_double("PRESSURE");
1321 if (FluidSystem::enableDissolvedGas())
1322 rsData = fp.get_double("RS");
1323
1324 if (FluidSystem::enableDissolvedGasInWater() && has_rsw)
1325 rswData = fp.get_double("RSW");
1326
1327 if (FluidSystem::enableVaporizedOil())
1328 rvData = fp.get_double("RV");
1329
1330 if (FluidSystem::enableVaporizedWater())
1331 rvwData = fp.get_double("RVW");
1332
1333 // initial reservoir temperature
1334 tempiData = fp.get_double("TEMPI");
1335
1336 // initial salt concentration data
1337 if constexpr (enableBrine)
1338 saltData = fp.get_double("SALT");
1339
1340 // initial precipitated salt saturation data
1341 if constexpr (enableSaltPrecipitation)
1342 saltpData = fp.get_double("SALTP");
1343
1344 // calculate the initial fluid states
1345 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1346 auto& dofFluidState = initialFluidStates_[dofIdx];
1347
1348 dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
1349
1351 // set temperature
1353 Scalar temperatureLoc = tempiData[dofIdx];
1354 if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0)
1355 temperatureLoc = FluidSystem::surfaceTemperature;
1356 dofFluidState.setTemperature(temperatureLoc);
1357
1359 // set salt concentration
1361 if constexpr (enableBrine)
1362 dofFluidState.setSaltConcentration(saltData[dofIdx]);
1363
1365 // set precipitated salt saturation
1367 if constexpr (enableSaltPrecipitation)
1368 dofFluidState.setSaltSaturation(saltpData[dofIdx]);
1369
1371 // set saturations
1373 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))
1374 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
1375 waterSaturationData[dofIdx]);
1376
1377 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
1378 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
1379 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1380 1.0
1381 - waterSaturationData[dofIdx]);
1382 }
1383 else
1384 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
1385 gasSaturationData[dofIdx]);
1386 }
1387 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
1388 const Scalar soil = 1.0 - waterSaturationData[dofIdx] - gasSaturationData[dofIdx];
1389 if (soil < smallSaturationTolerance_) {
1390 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
1391 }
1392 else {
1393 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx, soil);
1394 }
1395 }
1396
1398 // set phase pressures
1400 Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
1401
1402 // this assumes that capillary pressures only depend on the phase saturations
1403 // and possibly on temperature. (this is always the case for ECL problems.)
1404 std::array<Scalar, numPhases> pc = {0};
1405 const auto& matParams = this->materialLawParams(dofIdx);
1406 MaterialLaw::capillaryPressures(pc, matParams, dofFluidState);
1407 Valgrind::CheckDefined(pressure);
1408 Valgrind::CheckDefined(pc);
1409 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1410 if (!FluidSystem::phaseIsActive(phaseIdx))
1411 continue;
1412
1413 if (Indices::oilEnabled)
1414 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1415 else if (Indices::gasEnabled)
1416 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1417 else if (Indices::waterEnabled)
1418 //single (water) phase
1419 dofFluidState.setPressure(phaseIdx, pressure);
1420 }
1421
1422 if (FluidSystem::enableDissolvedGas())
1423 dofFluidState.setRs(rsData[dofIdx]);
1424 else if (Indices::gasEnabled && Indices::oilEnabled)
1425 dofFluidState.setRs(0.0);
1426
1427 if (FluidSystem::enableDissolvedGasInWater() && has_rsw)
1428 dofFluidState.setRsw(rswData[dofIdx]);
1429
1430 if (FluidSystem::enableVaporizedOil())
1431 dofFluidState.setRv(rvData[dofIdx]);
1432 else if (Indices::gasEnabled && Indices::oilEnabled)
1433 dofFluidState.setRv(0.0);
1434
1435 if (FluidSystem::enableVaporizedWater())
1436 dofFluidState.setRvw(rvwData[dofIdx]);
1437
1439 // set invB_
1441 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1442 if (!FluidSystem::phaseIsActive(phaseIdx))
1443 continue;
1444
1445 const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1446 dofFluidState.setInvB(phaseIdx, b);
1447
1448 const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx));
1449 dofFluidState.setDensity(phaseIdx, rho);
1450
1451 }
1452 }
1453 }
1454
1455
1456 void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation)
1457 {
1458 // each phase needs to be above certain value to be claimed to be existing
1459 // this is used to recover some RESTART running with the defaulted single-precision format
1460 Scalar sumSaturation = 0.0;
1461 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1462 if (FluidSystem::phaseIsActive(phaseIdx)) {
1463 if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance_)
1464 elemFluidState.setSaturation(phaseIdx, 0.0);
1465
1466 sumSaturation += elemFluidState.saturation(phaseIdx);
1467 }
1468
1469 }
1470 if constexpr (enableSolvent) {
1471 if (solventSaturation < smallSaturationTolerance_)
1472 solventSaturation = 0.0;
1473
1474 sumSaturation += solventSaturation;
1475 }
1476
1477 assert(sumSaturation > 0.0);
1478
1479 for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1480 if (FluidSystem::phaseIsActive(phaseIdx)) {
1481 const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation;
1482 elemFluidState.setSaturation(phaseIdx, saturation);
1483 }
1484 }
1485 if constexpr (enableSolvent) {
1486 solventSaturation = solventSaturation / sumSaturation;
1487 }
1488 }
1489
1490 void readInitialCondition_() override
1491 {
1492 FlowProblemType::readInitialCondition_();
1493
1494 if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP)
1495 this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(),
1496 enableSolvent,
1497 enablePolymer,
1498 enablePolymerMolarWeight,
1499 enableMICP);
1500
1501 }
1502
1503 void handleSolventBC(const BCProp::BCFace& bc, RateVector& rate) const override
1504 {
1505 if constexpr (!enableSolvent)
1506 throw std::logic_error("solvent is disabled and you're trying to add solvent to BC");
1507
1508 rate[Indices::solventSaturationIdx] = bc.rate;
1509 }
1510
1511 void handlePolymerBC(const BCProp::BCFace& bc, RateVector& rate) const override
1512 {
1513 if constexpr (!enablePolymer)
1514 throw std::logic_error("polymer is disabled and you're trying to add polymer to BC");
1515
1516 rate[Indices::polymerConcentrationIdx] = bc.rate;
1517 }
1518
1519 void updateExplicitQuantities_(const bool first_step_after_restart)
1520 {
1521 OPM_TIMEBLOCK(updateExplicitQuantities);
1522 const bool invalidateFromMaxWaterSat = this->updateMaxWaterSaturation_();
1523 const bool invalidateFromMinPressure = this->updateMinPressure_();
1524
1525 // update hysteresis and max oil saturation used in vappars
1526 const bool invalidateFromHyst = this->updateHysteresis_();
1527 const bool invalidateFromMaxOilSat = this->updateMaxOilSaturation_();
1528
1529 // deal with DRSDT and DRVDT
1530 const bool invalidateDRDT = !first_step_after_restart && this->updateCompositionChangeLimits_();
1531
1532 // the derivatives may have changed
1533 const bool invalidateIntensiveQuantities
1534 = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat || invalidateDRDT;
1535 if (invalidateIntensiveQuantities) {
1536 OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
1537 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
1538 }
1539
1540 this->updateRockCompTransMultVal_();
1541 }
1542
1543 bool satfuncConsistencyRequirementsMet() const
1544 {
1545 if (const auto nph = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
1546 + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
1547 + FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1548 nph < 2)
1549 {
1550 // Single phase runs don't need saturation functions and there's
1551 // nothing to do here. Return 'true' to tell caller that the
1552 // consistency requirements are Met.
1553 return true;
1554 }
1555
1556 const auto numSamplePoints = static_cast<std::size_t>
1557 (Parameters::Get<Parameters::NumSatfuncConsistencySamplePoints>());
1558
1559 auto sfuncConsistencyChecks =
1560 Satfunc::PhaseChecks::SatfuncConsistencyCheckManager<Scalar> {
1561 numSamplePoints, this->simulator().vanguard().eclState(),
1562 [&cmap = this->simulator().vanguard().cartesianIndexMapper()](const int elemIdx)
1563 { return cmap.cartesianIndex(elemIdx); }
1564 };
1565
1566 const auto ioRank = 0;
1567 const auto isIoRank = this->simulator().vanguard()
1568 .grid().comm().rank() == ioRank;
1569
1570 sfuncConsistencyChecks.collectFailuresTo(ioRank)
1571 .run(this->simulator().vanguard().grid().leafGridView(),
1572 [&vg = this->simulator().vanguard(),
1573 &emap = this->simulator().model().elementMapper()]
1574 (const auto& elem)
1575 { return vg.gridIdxToEquilGridIdx(emap.index(elem)); });
1576
1577 using ViolationLevel = typename Satfunc::PhaseChecks::
1578 SatfuncConsistencyCheckManager<Scalar>::ViolationLevel;
1579
1580 auto reportFailures = [&sfuncConsistencyChecks]
1581 (const ViolationLevel level)
1582 {
1583 sfuncConsistencyChecks.reportFailures
1584 (level, [](std::string_view record)
1585 { OpmLog::info(std::string { record }); });
1586 };
1587
1588 if (sfuncConsistencyChecks.anyFailedStandardChecks()) {
1589 if (isIoRank) {
1590 OpmLog::warning("Saturation Function "
1591 "End-point Consistency Problems");
1592
1593 reportFailures(ViolationLevel::Standard);
1594 }
1595 }
1596
1597 if (sfuncConsistencyChecks.anyFailedCriticalChecks()) {
1598 if (isIoRank) {
1599 OpmLog::error("Saturation Function "
1600 "End-point Consistency Failures");
1601
1602 reportFailures(ViolationLevel::Critical);
1603 }
1604
1605 // There are "critical" check failures. Report that consistency
1606 // requirements are not Met.
1607 return false;
1608 }
1609
1610 // If we get here then there are no critical failures. Report
1611 // Met = true, i.e., that the consistency requirements ARE met.
1612 return true;
1613 }
1614
1615 FlowThresholdPressure<TypeTag> thresholdPressures_;
1616
1617 std::vector<InitialFluidState> initialFluidStates_;
1618
1619 bool enableEclOutput_;
1620 std::unique_ptr<EclWriterType> eclWriter_;
1621
1622 const Scalar smallSaturationTolerance_ = 1.e-6;
1623#if HAVE_DAMARIS
1624 bool enableDamarisOutput_ = false ;
1625 std::unique_ptr<DamarisWriterType> damarisWriter_;
1626#endif
1627 MixingRateControls<FluidSystem> mixControls_;
1628
1629 ActionHandler<Scalar> actionHandler_;
1630
1631 ModuleParams moduleParams_;
1632};
1633
1634} // namespace Opm
1635
1636#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:48
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:83
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:841
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:852
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:505
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition FlowProblemBlackoil.hpp:883
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:537
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition FlowProblemBlackoil.hpp:933
void initialSolutionApplied() override
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition FlowProblemBlackoil.hpp:568
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:867
LhsEval permFactTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to porosity reduction.
Definition FlowProblemBlackoil.hpp:692
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:974
Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
Definition FlowProblemBlackoil.hpp:1081
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:42
Struct holding the parameters for the BlackOilPolymerModule class.
Definition blackoilpolymerparams.hpp:45
Struct holding the parameters for the BlackOilSolventModule class.
Definition blackoilsolventparams.hpp:49