My Project
Loading...
Searching...
No Matches
FlowProblem.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 Copyright 2023 INRIA
5
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 2 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20
21 Consult the COPYING file in the top-level source directory of this
22 module for the precise wording of the license and the list of
23 copyright holders.
24*/
30#ifndef OPM_FLOW_PROBLEM_HPP
31#define OPM_FLOW_PROBLEM_HPP
32
33#include <dune/common/version.hh>
34#include <dune/common/fvector.hh>
35#include <dune/common/fmatrix.hh>
36
37#include <opm/common/utility/TimeService.hpp>
38
39#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
40#include <opm/input/eclipse/Parser/ParserKeywords/E.hpp>
41#include <opm/input/eclipse/Schedule/Schedule.hpp>
42
43#include <opm/material/common/ConditionalStorage.hpp>
44#include <opm/material/common/Valgrind.hpp>
45#include <opm/material/densead/Evaluation.hpp>
46#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
47#include <opm/material/thermal/EclThermalLawManager.hpp>
48
52
53#include <opm/output/eclipse/EclipseIO.hpp>
54
60// TODO: maybe we can name it FlowProblemProperties.hpp
62#include <opm/simulators/flow/FlowUtils.hpp>
65#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
66#include <opm/simulators/timestepping/SimulatorReport.hpp>
67
68#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
69#include <opm/simulators/utils/ParallelSerialization.hpp>
70#include <opm/simulators/utils/satfunc/RelpermDiagnostics.hpp>
71
72#include <opm/utility/CopyablePtr.hpp>
73
74#include <algorithm>
75#include <functional>
76#include <set>
77#include <string>
78#include <vector>
79
80namespace Opm {
81
88template <class TypeTag>
89class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
90 , public FlowGenericProblem<GetPropType<TypeTag, Properties::GridView>,
91 GetPropType<TypeTag, Properties::FluidSystem>>
92{
93protected:
97 using Implementation = GetPropType<TypeTag, Properties::Problem>;
98
106
107 // Grid and world dimension
108 enum { dim = GridView::dimension };
109 enum { dimWorld = GridView::dimensionworld };
110
111 // copy some indices for convenience
113 enum { numPhases = FluidSystem::numPhases };
114 enum { numComponents = FluidSystem::numComponents };
115
116 enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
117 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
118 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
119 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
120 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
121 enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
122 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
123 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
124 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
125 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
126 enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
127 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
128 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
129 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
130 enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
131
132 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
133 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
134 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
135
136 // TODO: later, gasCompIdx, oilCompIdx and waterCompIdx should go to the FlowProblemBlackoil in the future
137 // we do not want them in the compositional setting
138 enum { gasCompIdx = FluidSystem::gasCompIdx };
139 enum { oilCompIdx = FluidSystem::oilCompIdx };
140 enum { waterCompIdx = FluidSystem::waterCompIdx };
141
145 using Element = typename GridView::template Codim<0>::Entity;
147 using EclMaterialLawManager = typename GetProp<TypeTag, Properties::MaterialLaw>::EclMaterialLawManager;
148 using EclThermalLawManager = typename GetProp<TypeTag, Properties::SolidEnergyLaw>::EclThermalLawManager;
149 using MaterialLawParams = typename EclMaterialLawManager::MaterialLawParams;
150 using SolidEnergyLawParams = typename EclThermalLawManager::SolidEnergyLawParams;
151 using ThermalConductionLawParams = typename EclThermalLawManager::ThermalConductionLawParams;
159
160 using Toolbox = MathToolbox<Evaluation>;
161 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
162
163
165 using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
166
167public:
174 using BaseType::porosity;
175
179 static void registerParameters()
180 {
181 ParentType::registerParameters();
182
184 }
185
186
190 static int handlePositionalParameter(std::function<void(const std::string&,
191 const std::string&)> addKey,
192 std::set<std::string>& seenParams,
193 std::string& errorMsg,
194 int,
195 const char** argv,
196 int paramIdx,
197 int)
198 {
199 return detail::eclPositionalParameter(addKey,
201 errorMsg,
202 argv,
203 paramIdx);
204 }
205
209 explicit FlowProblem(Simulator& simulator)
210 : ParentType(simulator)
211 , BaseType(simulator.vanguard().eclState(),
212 simulator.vanguard().schedule(),
213 simulator.vanguard().gridView())
214 , transmissibilities_(simulator.vanguard().eclState(),
215 simulator.vanguard().gridView(),
216 simulator.vanguard().cartesianIndexMapper(),
217 simulator.vanguard().grid(),
218 simulator.vanguard().cellCentroids(),
219 enableEnergy,
220 enableDiffusion,
221 enableDispersion)
222 , wellModel_(simulator)
223 , aquiferModel_(simulator)
224 , pffDofData_(simulator.gridView(), this->elementMapper())
225 , tracerModel_(simulator)
226 {
227 const auto& vanguard = simulator.vanguard();
228
229 enableDriftCompensation_ = Parameters::Get<Parameters::EnableDriftCompensation>();
230
231 enableVtkOutput_ = Parameters::Get<Parameters::EnableVtkOutput>();
232
233 this->enableTuning_ = Parameters::Get<Parameters::EnableTuning>();
234 this->initialTimeStepSize_ = Parameters::Get<Parameters::InitialTimeStepSize<Scalar>>();
235 this->maxTimeStepAfterWellEvent_ = Parameters::Get<Parameters::TimeStepAfterEventInDays<Scalar>>() * 24 * 60 * 60;
236
237 // The value N for this parameter is defined in the following order of presedence:
238 // 1. Command line value (--num-pressure-points-equil=N)
239 // 2. EQLDIMS item 2
240 // Default value is defined in opm-common/src/opm/input/eclipse/share/keywords/000_Eclipse100/E/EQLDIMS
241 if (Parameters::IsSet<Parameters::NumPressurePointsEquil>())
242 {
243 this->numPressurePointsEquil_ = Parameters::Get<Parameters::NumPressurePointsEquil>();
244 } else {
245 this->numPressurePointsEquil_ = simulator.vanguard().eclState().getTableManager().getEqldims().getNumDepthNodesP();
246 }
247
248 explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
249
250
252 relpermDiagnostics.diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper());
253 }
254
255 virtual ~FlowProblem() = default;
256
257 void prefetch(const Element& elem) const
258 { pffDofData_.prefetch(elem); }
259
271 template <class Restarter>
273 {
274 // reload the current episode/report step from the deck
275 this->beginEpisode();
276
277 // deserialize the wells
278 wellModel_.deserialize(res);
279
280 // deserialize the aquifer
281 aquiferModel_.deserialize(res);
282 }
283
290 template <class Restarter>
292 {
293 wellModel_.serialize(res);
294
295 aquiferModel_.serialize(res);
296 }
297
298 int episodeIndex() const
299 {
300 return std::max(this->simulator().episodeIndex(), 0);
301 }
302
306 virtual void beginEpisode()
307 {
309 // Proceed to the next report step
310 auto& simulator = this->simulator();
311 int episodeIdx = simulator.episodeIndex();
312 auto& eclState = simulator.vanguard().eclState();
313 const auto& schedule = simulator.vanguard().schedule();
314 const auto& events = schedule[episodeIdx].events();
315
316 if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
317 // bring the contents of the keywords to the current state of the SCHEDULE
318 // section.
319 //
320 // TODO (?): make grid topology changes possible (depending on what exactly
321 // has changed, the grid may need be re-created which has some serious
322 // implications on e.g., the solution of the simulation.)
323 const auto& miniDeck = schedule[episodeIdx].geo_keywords();
324 const auto& cc = simulator.vanguard().grid().comm();
325 eclState.apply_schedule_keywords( miniDeck );
326 eclBroadcast(cc, eclState.getTransMult() );
327
328 // Re-ordering in case of ALUGrid
329 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
330 return simulator.vanguard().gridEquilIdxToGridIdx(i);
331 };
332
333 // re-compute all quantities which may possibly be affected.
334 using TransUpdateQuantities = typename Vanguard::TransmissibilityType::TransUpdateQuantities;
335 transmissibilities_.update(true, TransUpdateQuantities::All, equilGridToGrid);
336 this->referencePorosity_[1] = this->referencePorosity_[0];
337 updateReferencePorosity_();
338 updatePffDofData_();
339 this->model().linearizer().updateDiscretizationParameters();
340 }
341
342 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
343
344 // set up the wells for the next episode.
345 wellModel_.beginEpisode();
346
347 // set up the aquifers for the next episode.
348 aquiferModel_.beginEpisode();
349
350 // set the size of the initial time step of the episode
351 Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
352 // negative value of initialTimeStepSize_ indicates no active limit from TSINIT or NEXTSTEP
353 if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
354 // allow the size of the initial time step to be set via an external parameter
355 // if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
356 dt = std::min(dt, this->initialTimeStepSize_);
357 simulator.setTimeStepSize(dt);
358 }
359
364 {
366 const int episodeIdx = this->episodeIndex();
367 const int timeStepSize = this->simulator().timeStepSize();
368
369 this->beginTimeStep_(enableExperiments,
371 this->simulator().timeStepIndex(),
372 this->simulator().startTime(),
373 this->simulator().time(),
374 timeStepSize,
375 this->simulator().endTime());
376
377 // update maximum water saturation and minimum pressure
378 // used when ROCKCOMP is activated
379 // Do not update max RS first step after a restart
380 this->updateExplicitQuantities_(episodeIdx, timeStepSize, first_step_ && (episodeIdx > 0));
381 first_step_ = false;
382
383 if (nonTrivialBoundaryConditions()) {
384 this->model().linearizer().updateBoundaryConditionData();
385 }
386
387 wellModel_.beginTimeStep();
388 aquiferModel_.beginTimeStep();
389 tracerModel_.beginTimeStep();
390
391 }
392
397 {
399 wellModel_.beginIteration();
400 aquiferModel_.beginIteration();
401 }
402
407 {
409 wellModel_.endIteration();
410 aquiferModel_.endIteration();
411 }
412
416 virtual void endTimeStep()
417 {
419
420#ifndef NDEBUG
422 // in debug mode, we don't care about performance, so we check
423 // if the model does the right thing (i.e., the mass change
424 // inside the whole reservoir must be equivalent to the fluxes
425 // over the grid's boundaries plus the source rates specified by
426 // the problem).
427 const int rank = this->simulator().gridView().comm().rank();
428 if (rank == 0) {
429 std::cout << "checking conservativeness of solution\n";
430 }
431
432 this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
433 if (rank == 0) {
434 std::cout << "solution is sufficiently conservative\n";
435 }
436 }
437#endif // NDEBUG
438
439 auto& simulator = this->simulator();
440 simulator.setTimeStepIndex(simulator.timeStepIndex()+1);
441
442 this->wellModel_.endTimeStep();
443 this->aquiferModel_.endTimeStep();
444 this->tracerModel_.endTimeStep();
445
446 // Compute flux for output
447 this->model().linearizer().updateFlowsInfo();
448
449 if (this->enableDriftCompensation_) {
451
452 const auto& residual = this->model().linearizer().residual();
453
454 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
455 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(globalDofIdx);
456 this->drift_[sfcdofIdx] = residual[sfcdofIdx] * simulator.timeStepSize();
457
459 this->drift_[sfcdofIdx] *= this->model().dofTotalVolume(sfcdofIdx);
460 }
461 }
462 }
463 }
464
468 virtual void endEpisode()
469 {
470 const int episodeIdx = this->episodeIndex();
471
472 this->wellModel_.endEpisode();
473 this->aquiferModel_.endEpisode();
474
475 const auto& schedule = this->simulator().vanguard().schedule();
476
477 // End simulation when completed.
478 if (episodeIdx + 1 >= static_cast<int>(schedule.size()) - 1) {
479 this->simulator().setFinished(true);
480 return;
481 }
482
483 // Otherwise, start next episode (report step).
484 this->simulator().startNextEpisode(schedule.stepLength(episodeIdx + 1));
485 }
486
491 void writeOutput(bool verbose = true)
492 {
494 // use the generic code to prepare the output fields and to
495 // write the desired VTK files.
496 if (Parameters::Get<Parameters::EnableWriteAllSolutions>() ||
497 this->simulator().episodeWillBeOver()) {
498 ParentType::writeOutput(verbose);
499 }
500 }
501
505 template <class Context>
506 const DimMatrix& intrinsicPermeability(const Context& context,
507 unsigned spaceIdx,
508 unsigned timeIdx) const
509 {
510 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
511 return transmissibilities_.permeability(globalSpaceIdx);
512 }
513
520 const DimMatrix& intrinsicPermeability(unsigned globalElemIdx) const
521 { return transmissibilities_.permeability(globalElemIdx); }
522
526 template <class Context>
527 Scalar transmissibility(const Context& context,
528 [[maybe_unused]] unsigned fromDofLocalIdx,
529 unsigned toDofLocalIdx) const
530 {
532 return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
533 }
534
538 Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
539 {
540 return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
541 }
542
546 template <class Context>
547 Scalar diffusivity(const Context& context,
548 [[maybe_unused]] unsigned fromDofLocalIdx,
549 unsigned toDofLocalIdx) const
550 {
552 return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
553 }
554
558 Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
559 return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
560 }
561
565 Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
566 return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
567 }
568
573 const unsigned boundaryFaceIdx) const
574 {
575 return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
576 }
577
578
579
580
584 template <class Context>
586 unsigned boundaryFaceIdx) const
587 {
588 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
589 return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
590 }
591
596 const unsigned boundaryFaceIdx) const
597 {
598 return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
599 }
600
601
606 const unsigned globalSpaceIdxOut) const
607 {
608 return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
609 }
610
614 template <class Context>
616 unsigned faceIdx,
617 unsigned timeIdx) const
618 {
619 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
620 unsigned toDofLocalIdx = face.exteriorIndex();
621 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
622 }
623
627 template <class Context>
629 unsigned faceIdx,
630 unsigned timeIdx) const
631 {
632 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
633 unsigned toDofLocalIdx = face.exteriorIndex();
634 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
635 }
636
640 template <class Context>
642 unsigned boundaryFaceIdx) const
643 {
644 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
645 return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
646 }
647
651 const typename Vanguard::TransmissibilityType& eclTransmissibilities() const
652 { return transmissibilities_; }
653
654
655 const TracerModel& tracerModel() const
656 { return tracerModel_; }
657
658 TracerModel& tracerModel()
659 { return tracerModel_; }
660
669 template <class Context>
670 Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
671 {
672 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
673 return this->porosity(globalSpaceIdx, timeIdx);
674 }
675
682 template <class Context>
683 Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
684 {
685 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
686 return this->dofCenterDepth(globalSpaceIdx);
687 }
688
695 Scalar dofCenterDepth(unsigned globalSpaceIdx) const
696 {
697 return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
698 }
699
703 template <class Context>
704 Scalar rockCompressibility(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
705 {
706 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
707 return this->rockCompressibility(globalSpaceIdx);
708 }
709
713 template <class Context>
714 Scalar rockReferencePressure(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
715 {
716 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
717 return this->rockReferencePressure(globalSpaceIdx);
718 }
719
723 template <class Context>
724 const MaterialLawParams& materialLawParams(const Context& context,
725 unsigned spaceIdx, unsigned timeIdx) const
726 {
727 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
728 return this->materialLawParams(globalSpaceIdx);
729 }
730
731 const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const
732 {
733 return materialLawManager_->materialLawParams(globalDofIdx);
734 }
735
736 const MaterialLawParams& materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const
737 {
738 return materialLawManager_->materialLawParams(globalDofIdx, facedir);
739 }
740
744 template <class Context>
745 const SolidEnergyLawParams&
747 unsigned spaceIdx,
748 unsigned timeIdx) const
749 {
750 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
751 return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
752 }
753
757 template <class Context>
758 const ThermalConductionLawParams &
759 thermalConductionLawParams(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
760 {
761 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
762 return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
763 }
764
771 std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
772 { return materialLawManager_; }
773
774 template <class FluidState>
775 void updateRelperms(
776 std::array<Evaluation,numPhases> &mobility,
777 DirectionalMobilityPtr &dirMob,
778 FluidState &fluidState,
779 unsigned globalSpaceIdx) const
780 {
781 OPM_TIMEBLOCK_LOCAL(updateRelperms);
782 {
783 // calculate relative permeabilities. note that we store the result into the
784 // mobility_ class attribute. the division by the phase viscosity happens later.
786 MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
787 Valgrind::CheckDefined(mobility);
788 }
789 if (materialLawManager_->hasDirectionalRelperms()
790 || materialLawManager_->hasDirectionalImbnum())
791 {
792 using Dir = FaceDir::DirEnum;
793 constexpr int ndim = 3;
794 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
795 Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
796 for (int i = 0; i<ndim; i++) {
798 auto& mob_array = dirMob->getArray(i);
799 MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
800 }
801 }
802 }
803
807 std::shared_ptr<EclMaterialLawManager> materialLawManager()
808 { return materialLawManager_; }
809
814 template <class Context>
815 unsigned pvtRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
816 { return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
817
822 template <class Context>
823 unsigned satnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
824 { return this->satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
825
830 template <class Context>
831 unsigned miscnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
832 { return this->miscnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
833
838 template <class Context>
839 unsigned plmixnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
840 { return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
841
842 // TODO: polymer related might need to go to the blackoil side
847 template <class Context>
848 Scalar maxPolymerAdsorption(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
849 { return this->maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); }
850
854 std::string name() const
855 { return this->simulator().vanguard().caseName(); }
856
860 template <class Context>
861 Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
862 {
863 // use the initial temperature of the DOF if temperature is not a primary
864 // variable
865 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
866 return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0);
867 }
868
869
870 Scalar temperature(unsigned globalDofIdx, unsigned /*timeIdx*/) const
871 {
872 // use the initial temperature of the DOF if temperature is not a primary
873 // variable
874 return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0);
875 }
876
877 const SolidEnergyLawParams&
879 unsigned /*timeIdx*/) const
880 {
881 return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
882 }
883 const ThermalConductionLawParams &
885 unsigned /*timeIdx*/)const
886 {
887 return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
888 }
889
899 Scalar maxOilSaturation(unsigned globalDofIdx) const
900 {
901 if (!this->vapparsActive(this->episodeIndex()))
902 return 0.0;
903
904 return this->maxOilSaturation_[globalDofIdx];
905 }
906
916 void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
917 {
918 if (!this->vapparsActive(this->episodeIndex()))
919 return;
920
921 this->maxOilSaturation_[globalDofIdx] = value;
922 }
923
928 {
929 // Calculate all intensive quantities.
930 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/0);
931
932 // We also need the intensive quantities for timeIdx == 1
933 // corresponding to the start of the current timestep, if we
934 // do not use the storage cache, or if we cannot recycle the
935 // first iteration storage.
936 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
937 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/1);
938 }
939
940 // initialize the wells. Note that this needs to be done after initializing the
941 // intrinsic permeabilities and the after applying the initial solution because
942 // the well model uses these...
943 wellModel_.init();
944
945 aquiferModel_.initialSolutionApplied();
946
947 const bool invalidateFromHyst = updateHysteresis_();
948 if (invalidateFromHyst) {
950 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
951 }
952 }
953
959 template <class Context>
960 void source(RateVector& rate,
961 const Context& context,
962 unsigned spaceIdx,
963 unsigned timeIdx) const
964 {
965 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
967 }
968
969 void source(RateVector& rate,
970 unsigned globalDofIdx,
971 unsigned timeIdx) const
972 {
974 rate = 0.0;
975
976 // Add well contribution to source here.
977 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
978
979 // convert the source term from the total mass rate of the
980 // cell to the one per unit of volume as used by the model.
981 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
982 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
983
984 Valgrind::CheckDefined(rate[eqIdx]);
986 }
987
988 // Add non-well sources.
989 addToSourceDense(rate, globalDofIdx, timeIdx);
990 }
991
992 virtual void addToSourceDense(RateVector& rate,
993 unsigned globalDofIdx,
994 unsigned timeIdx) const = 0;
995
1001 const WellModel& wellModel() const
1002 { return wellModel_; }
1003
1004 WellModel& wellModel()
1005 { return wellModel_; }
1006
1007 const AquiferModel& aquiferModel() const
1008 { return aquiferModel_; }
1009
1010 AquiferModel& mutableAquiferModel()
1011 { return aquiferModel_; }
1012
1013 bool nonTrivialBoundaryConditions() const
1014 { return nonTrivialBoundaryConditions_; }
1015
1022 Scalar nextTimeStepSize() const
1023 {
1025 // allow external code to do the timestepping
1026 if (this->nextTimeStepSize_ > 0.0)
1027 return this->nextTimeStepSize_;
1028
1029 const auto& simulator = this->simulator();
1030 int episodeIdx = simulator.episodeIndex();
1031
1032 // for the initial episode, we use a fixed time step size
1033 if (episodeIdx < 0)
1034 return this->initialTimeStepSize_;
1035
1036 // ask the newton method for a suggestion. This suggestion will be based on how
1037 // well the previous time step converged. After that, apply the runtime time
1038 // stepping constraints.
1039 const auto& newtonMethod = this->model().newtonMethod();
1040 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1041 }
1042
1048 template <class LhsEval>
1049 LhsEval rockCompPoroMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1050 {
1052 if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
1053 return 1.0;
1054
1055 unsigned tableIdx = 0;
1056 if (!this->rockTableIdx_.empty())
1057 tableIdx = this->rockTableIdx_[elementIdx];
1058
1059 const auto& fs = intQuants.fluidState();
1060 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1061 if (!this->minRefPressure_.empty())
1062 // The pore space change is irreversible
1064 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1065 this->minRefPressure_[elementIdx]);
1066
1067 if (!this->overburdenPressure_.empty())
1068 effectivePressure -= this->overburdenPressure_[elementIdx];
1069
1070
1071 if (!this->rockCompPoroMult_.empty()) {
1072 return this->rockCompPoroMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
1073 }
1074
1075 // water compaction
1076 assert(!this->rockCompPoroMultWc_.empty());
1077 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1078 LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
1079
1080 return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
1081 }
1082
1088 template <class LhsEval>
1089 LhsEval rockCompTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1090 {
1091 const bool implicit = !this->explicitRockCompaction_;
1092 return implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<LhsEval>(intQuants, elementIdx)
1093 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1094 }
1095
1096
1100 template <class LhsEval>
1101 LhsEval wellTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1102 {
1104
1105 const bool implicit = !this->explicitRockCompaction_;
1106 double trans_mult = implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<double>(intQuants, elementIdx)
1107 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1108 trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants);
1109
1110 return trans_mult;
1111 }
1112
1113 std::pair<BCType, RateVector> boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) const
1114 {
1115 OPM_TIMEBLOCK_LOCAL(boundaryCondition);
1116 if (!nonTrivialBoundaryConditions_) {
1117 return { BCType::NONE, RateVector(0.0) };
1118 }
1119 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1120 const auto& schedule = this->simulator().vanguard().schedule();
1121 if (bcindex_(dir)[globalSpaceIdx] == 0) {
1122 return { BCType::NONE, RateVector(0.0) };
1123 }
1124 if (schedule[this->episodeIndex()].bcprop.size() == 0) {
1125 return { BCType::NONE, RateVector(0.0) };
1126 }
1127 const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
1128 if (bc.bctype!=BCType::RATE) {
1129 return { bc.bctype, RateVector(0.0) };
1130 }
1131
1132 RateVector rate = 0.0;
1133 switch (bc.component) {
1134 case BCComponent::OIL:
1135 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = bc.rate;
1136 break;
1137 case BCComponent::GAS:
1138 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = bc.rate;
1139 break;
1140 case BCComponent::WATER:
1141 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate;
1142 break;
1143 case BCComponent::SOLVENT:
1144 this->handleSolventBC(bc, rate);
1145 break;
1146 case BCComponent::POLYMER:
1147 this->handlePolymerBC(bc, rate);
1148 break;
1149 case BCComponent::NONE:
1150 throw std::logic_error("you need to specify the component when RATE type is set in BC");
1151 break;
1152 }
1153 //TODO add support for enthalpy rate
1154 return {bc.bctype, rate};
1155 }
1156
1157
1158 template<class Serializer>
1159 void serializeOp(Serializer& serializer)
1160 {
1161 serializer(static_cast<BaseType&>(*this));
1162 serializer(drift_);
1163 serializer(wellModel_);
1164 serializer(aquiferModel_);
1165 serializer(tracerModel_);
1166 serializer(*materialLawManager_);
1167 }
1168
1169private:
1170 Implementation& asImp_()
1171 { return *static_cast<Implementation *>(this); }
1172
1173 const Implementation& asImp_() const
1174 { return *static_cast<const Implementation *>(this); }
1175
1176protected:
1177 template<class UpdateFunc>
1178 void updateProperty_(const std::string& failureMsg,
1180 {
1182 const auto& model = this->simulator().model();
1183 const auto& primaryVars = model.solution(/*timeIdx*/0);
1184 const auto& vanguard = this->simulator().vanguard();
1185 std::size_t numGridDof = primaryVars.size();
1186 OPM_BEGIN_PARALLEL_TRY_CATCH();
1187#ifdef _OPENMP
1188#pragma omp parallel for
1189#endif
1190 for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1191 const auto& iq = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/ 0);
1192 func(dofIdx, iq);
1193 }
1194 OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
1195 }
1196
1197 bool updateMaxOilSaturation_()
1198 {
1200 int episodeIdx = this->episodeIndex();
1201
1202 // we use VAPPARS
1203 if (this->vapparsActive(episodeIdx)) {
1204 this->updateProperty_("FlowProblem::updateMaxOilSaturation_() failed:",
1205 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1206 {
1207 this->updateMaxOilSaturation_(compressedDofIdx,iq);
1208 });
1209 return true;
1210 }
1211
1212 return false;
1213 }
1214
1215 bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1216 {
1218 const auto& fs = iq.fluidState();
1219 const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1220 auto& mos = this->maxOilSaturation_;
1221 if(mos[compressedDofIdx] < So){
1222 mos[compressedDofIdx] = So;
1223 return true;
1224 }else{
1225 return false;
1226 }
1227 }
1228
1229 bool updateMaxWaterSaturation_()
1230 {
1232 // water compaction is activated in ROCKCOMP
1233 if (this->maxWaterSaturation_.empty())
1234 return false;
1235
1236 this->maxWaterSaturation_[/*timeIdx=*/1] = this->maxWaterSaturation_[/*timeIdx=*/0];
1237 this->updateProperty_("FlowProblem::updateMaxWaterSaturation_() failed:",
1238 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1239 {
1240 this->updateMaxWaterSaturation_(compressedDofIdx,iq);
1241 });
1242 return true;
1243 }
1244
1245
1246 bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1247 {
1249 const auto& fs = iq.fluidState();
1250 const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
1251 auto& mow = this->maxWaterSaturation_;
1252 if(mow[compressedDofIdx]< Sw){
1253 mow[compressedDofIdx] = Sw;
1254 return true;
1255 }else{
1256 return false;
1257 }
1258 }
1259
1260 bool updateMinPressure_()
1261 {
1263 // IRREVERS option is used in ROCKCOMP
1264 if (this->minRefPressure_.empty())
1265 return false;
1266
1267 this->updateProperty_("FlowProblem::updateMinPressure_() failed:",
1268 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1269 {
1270 this->updateMinPressure_(compressedDofIdx,iq);
1271 });
1272 return true;
1273 }
1274
1275 bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
1277 const auto& fs = iq.fluidState();
1278 const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
1279 auto& min_pressures = this->minRefPressure_;
1282 return true;
1283 }else{
1284 return false;
1285 }
1286 }
1287
1288 // \brief Function to assign field properties of type double, on the leaf grid view.
1289 //
1290 // For CpGrid with local grid refinement, the field property of a cell on the leaf
1291 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
1292 std::function<std::vector<double>(const FieldPropsManager&, const std::string&)>
1293 fieldPropDoubleOnLeafAssigner_()
1294 {
1295 const auto& lookup = this->lookUpData_;
1296 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString)
1297 {
1298 return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
1299 };
1300 }
1301
1302 // \brief Function to assign field properties of type int, unsigned int, ..., on the leaf grid view.
1303 //
1304 // For CpGrid with local grid refinement, the field property of a cell on the leaf
1305 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
1306 template<typename IntType>
1307 std::function<std::vector<IntType>(const FieldPropsManager&, const std::string&, bool)>
1308 fieldPropIntTypeOnLeafAssigner_()
1309 {
1310 const auto& lookup = this->lookUpData_;
1311 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString, bool needsTranslation)
1312 {
1314 };
1315 }
1316
1317 void readMaterialParameters_()
1318 {
1320 const auto& simulator = this->simulator();
1321 const auto& vanguard = simulator.vanguard();
1322 const auto& eclState = vanguard.eclState();
1323
1324 // the PVT and saturation region numbers
1325 OPM_BEGIN_PARALLEL_TRY_CATCH();
1326 this->updatePvtnum_();
1327 this->updateSatnum_();
1328
1329 // the MISC region numbers (solvent model)
1330 this->updateMiscnum_();
1331 // the PLMIX region numbers (polymer model)
1332 this->updatePlmixnum_();
1333
1334 OPM_END_PARALLEL_TRY_CATCH("Invalid region numbers: ", vanguard.gridView().comm());
1336 // porosity
1337 updateReferencePorosity_();
1338 this->referencePorosity_[1] = this->referencePorosity_[0];
1340
1342 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
1343 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
1344 materialLawManager_->initFromState(eclState);
1345 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1346 this-> template fieldPropIntTypeOnLeafAssigner_<int>(),
1347 this-> lookupIdxOnLevelZeroAssigner_());
1349 }
1350
1351 void readThermalParameters_()
1352 {
1353 if constexpr (enableEnergy)
1354 {
1355 const auto& simulator = this->simulator();
1356 const auto& vanguard = simulator.vanguard();
1357 const auto& eclState = vanguard.eclState();
1358
1359 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
1360 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
1361 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1362 this-> fieldPropDoubleOnLeafAssigner_(),
1364 }
1365 }
1366
1367 void updateReferencePorosity_()
1368 {
1369 const auto& simulator = this->simulator();
1370 const auto& vanguard = simulator.vanguard();
1371 const auto& eclState = vanguard.eclState();
1372
1373 std::size_t numDof = this->model().numGridDof();
1374
1375 this->referencePorosity_[/*timeIdx=*/0].resize(numDof);
1376
1377 const auto& fp = eclState.fieldProps();
1378 const std::vector<double> porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV");
1379 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1380 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(dofIdx);
1381 Scalar poreVolume = porvData[dofIdx];
1382
1383 // we define the porosity as the accumulated pore volume divided by the
1384 // geometric volume of the element. Note that -- in pathetic cases -- it can
1385 // be larger than 1.0!
1386 Scalar dofVolume = simulator.model().dofTotalVolume(sfcdofIdx);
1387 assert(dofVolume > 0.0);
1388 this->referencePorosity_[/*timeIdx=*/0][sfcdofIdx] = poreVolume/dofVolume;
1389 }
1390 }
1391
1392 virtual void readInitialCondition_()
1393 {
1394 // TODO: whether we should move this to FlowProblemBlackoil
1395 const auto& simulator = this->simulator();
1396 const auto& vanguard = simulator.vanguard();
1397 const auto& eclState = vanguard.eclState();
1398
1399 if (eclState.getInitConfig().hasEquil())
1400 readEquilInitialCondition_();
1401 else
1402 readExplicitInitialCondition_();
1403
1404 //initialize min/max values
1405 std::size_t numElems = this->model().numGridDof();
1406 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1407 const auto& fs = asImp_().initialFluidStates()[elemIdx];
1408 if (!this->maxWaterSaturation_.empty() && waterPhaseIdx > -1)
1409 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
1410 if (!this->maxOilSaturation_.empty() && oilPhaseIdx > -1)
1411 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
1412 if (!this->minRefPressure_.empty() && refPressurePhaseIdx_() > -1)
1413 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
1414 }
1415 }
1416
1417 virtual void readEquilInitialCondition_() = 0;
1418 virtual void readExplicitInitialCondition_() = 0;
1419
1420 // update the hysteresis parameters of the material laws for the whole grid
1421 bool updateHysteresis_()
1422 {
1423 if (!materialLawManager_->enableHysteresis())
1424 return false;
1425
1426 // we need to update the hysteresis data for _all_ elements (i.e., not just the
1427 // interior ones) to avoid desynchronization of the processes in the parallel case!
1428 this->updateProperty_("FlowProblem::updateHysteresis_() failed:",
1429 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1430 {
1431 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
1432 });
1433 return true;
1434 }
1435
1436
1437 bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1438 {
1439 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
1440 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
1441 //TODO change materials to give a bool
1442 return true;
1443 }
1444
1445 Scalar getRockCompTransMultVal(std::size_t dofIdx) const
1446 {
1447 if (this->rockCompTransMultVal_.empty())
1448 return 1.0;
1449
1450 return this->rockCompTransMultVal_[dofIdx];
1451 }
1452
1453protected:
1455 {
1457 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
1460 Scalar transmissibility;
1461 };
1462
1463 // update the prefetch friendly data object
1464 void updatePffDofData_()
1465 {
1466 const auto& distFn =
1467 [this](PffDofData_& dofData,
1468 const Stencil& stencil,
1469 unsigned localDofIdx)
1470 -> void
1471 {
1472 const auto& elementMapper = this->model().elementMapper();
1473
1474 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
1475 if (localDofIdx != 0) {
1476 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(/*dofIdx=*/0));
1477 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
1478
1479 if constexpr (enableEnergy) {
1480 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
1481 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
1482 }
1483 if constexpr (enableDiffusion)
1484 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
1485 if (enableDispersion)
1486 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
1487 }
1488 };
1489
1490 pffDofData_.update(distFn);
1491 }
1492
1493 virtual void updateExplicitQuantities_(int episodeIdx, int timeStepSize, bool first_step_after_restart) = 0;
1494
1495 void readBoundaryConditions_()
1496 {
1497 const auto& simulator = this->simulator();
1498 const auto& vanguard = simulator.vanguard();
1499 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
1500 if (bcconfig.size() > 0) {
1501 nonTrivialBoundaryConditions_ = true;
1502
1503 std::size_t numCartDof = vanguard.cartesianSize();
1504 unsigned numElems = vanguard.gridView().size(/*codim=*/0);
1505 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
1506
1507 for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
1508 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
1509
1510 bcindex_.resize(numElems, 0);
1512 &vanguard](const auto& bcface,
1513 auto apply)
1514 {
1515 for (int i = bcface.i1; i <= bcface.i2; ++i) {
1516 for (int j = bcface.j1; j <= bcface.j2; ++j) {
1517 for (int k = bcface.k1; k <= bcface.k2; ++k) {
1518 std::array<int, 3> tmp = {i,j,k};
1519 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
1520 if (elemIdx >= 0)
1521 apply(elemIdx);
1522 }
1523 }
1524 }
1525 };
1526 for (const auto& bcface : bcconfig) {
1527 std::vector<int>& data = bcindex_(bcface.dir);
1528 const int index = bcface.index;
1530 [&data,index](int elemIdx)
1531 { data[elemIdx] = index; });
1532 }
1533 }
1534 }
1535
1536 // this method applies the runtime constraints specified via the deck and/or command
1537 // line parameters for the size of the next time step.
1538 Scalar limitNextTimeStepSize_(Scalar dtNext) const
1539 {
1540 if constexpr (enableExperiments) {
1541 const auto& simulator = this->simulator();
1542 const auto& schedule = simulator.vanguard().schedule();
1543 int episodeIdx = simulator.episodeIndex();
1544
1545 // first thing in the morning, limit the time step size to the maximum size
1546 Scalar maxTimeStepSize = Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60;
1547 int reportStepIdx = std::max(episodeIdx, 0);
1548 if (this->enableTuning_) {
1549 const auto& tuning = schedule[reportStepIdx].tuning();
1550 maxTimeStepSize = tuning.TSMAXZ;
1551 }
1552
1553 dtNext = std::min(dtNext, maxTimeStepSize);
1554
1555 Scalar remainingEpisodeTime =
1556 simulator.episodeStartTime() + simulator.episodeLength()
1557 - (simulator.startTime() + simulator.time());
1559
1560 // if we would have a small amount of time left over in the current episode, make
1561 // two equal time steps instead of a big and a small one
1562 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
1563 // note: limiting to the maximum time step size here is probably not strictly
1564 // necessary, but it should not hurt and is more fool-proof
1565 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
1566
1567 if (simulator.episodeStarts()) {
1568 // if a well event occurred, respect the limit for the maximum time step after
1569 // that, too
1570 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
1571 bool wellEventOccured =
1572 events.hasEvent(ScheduleEvents::NEW_WELL)
1573 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
1574 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
1575 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
1576 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
1577 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
1578 }
1579 }
1580
1581 return dtNext;
1582 }
1583
1584 int refPressurePhaseIdx_() const {
1585 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1586 return oilPhaseIdx;
1587 }
1588 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1589 return gasPhaseIdx;
1590 }
1591 else {
1592 return waterPhaseIdx;
1593 }
1594 }
1595
1596 void updateRockCompTransMultVal_()
1597 {
1598 const auto& model = this->simulator().model();
1599 std::size_t numGridDof = this->model().numGridDof();
1600 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
1601 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
1602 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, /*timeIdx=*/ 0);
1604 this->rockCompTransMultVal_[elementIdx] = trans_mult;
1605 }
1606 }
1607
1613 template <class LhsEval>
1614 LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1615 {
1617 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
1618 return 1.0;
1619
1620 unsigned tableIdx = 0;
1621 if (!this->rockTableIdx_.empty())
1622 tableIdx = this->rockTableIdx_[elementIdx];
1623
1624 const auto& fs = intQuants.fluidState();
1625 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1626
1627 if (!this->minRefPressure_.empty())
1628 // The pore space change is irreversible
1630 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1631 this->minRefPressure_[elementIdx]);
1632
1633 if (!this->overburdenPressure_.empty())
1634 effectivePressure -= this->overburdenPressure_[elementIdx];
1635
1636 if (!this->rockCompTransMult_.empty())
1637 return this->rockCompTransMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
1638
1639 // water compaction
1640 assert(!this->rockCompTransMultWc_.empty());
1641 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1642 LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
1643
1644 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
1645 }
1646
1647 typename Vanguard::TransmissibilityType transmissibilities_;
1648
1649 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
1650 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
1651
1652 bool enableDriftCompensation_;
1653 GlobalEqVector drift_;
1654
1655 WellModel wellModel_;
1656 AquiferModel aquiferModel_;
1657
1658 bool enableVtkOutput_;
1659
1660
1662 TracerModel tracerModel_;
1663
1664 template<class T>
1665 struct BCData
1666 {
1667 std::array<std::vector<T>,6> data;
1668
1669 void resize(std::size_t size, T defVal)
1670 {
1671 for (auto& d : data)
1672 d.resize(size, defVal);
1673 }
1674
1675 const std::vector<T>& operator()(FaceDir::DirEnum dir) const
1676 {
1677 if (dir == FaceDir::DirEnum::Unknown)
1678 throw std::runtime_error("Tried to access BC data for the 'Unknown' direction");
1679 int idx = 0;
1680 int div = static_cast<int>(dir);
1681 while ((div /= 2) >= 1)
1682 ++idx;
1683 assert(idx >= 0 && idx <= 5);
1684 return data[idx];
1685 }
1686
1687 std::vector<T>& operator()(FaceDir::DirEnum dir)
1688 {
1689 return const_cast<std::vector<T>&>(std::as_const(*this)(dir));
1690 }
1691 };
1692
1693 virtual void handleSolventBC(const BCProp::BCFace&, RateVector&) const = 0;
1694
1695 virtual void handlePolymerBC(const BCProp::BCFace&, RateVector&) const = 0;
1696
1697 BCData<int> bcindex_;
1698 bool nonTrivialBoundaryConditions_ = false;
1699 bool explicitRockCompaction_ = false;
1700 bool first_step_ = true;
1701
1702};
1703
1704} // namespace Opm
1705
1706#endif // OPM_FLOW_PROBLEM_HPP
Helper class for grid instantiation of ECL file-format using problems.
This is a "dummy" gradient calculator which does not do anything.
Collects necessary output values and pass it to opm-common's ECL output.
Computes the initial condition based on the EQUIL keyword from ECL.
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
A class which handles tracers as specified in by ECL.
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowGenericProblem.hpp:63
Scalar maxPolymerAdsorption(unsigned elemIdx) const
Returns the max polymer adsorption value.
Definition FlowGenericProblem_impl.hpp:709
unsigned pvtRegionIndex(unsigned elemIdx) const
Returns the index the relevant PVT region given a cell index.
Definition FlowGenericProblem_impl.hpp:668
Scalar rockReferencePressure(unsigned globalSpaceIdx) const
Direct access to rock reference pressure.
Definition FlowGenericProblem_impl.hpp:290
static std::string briefDescription()
Returns a human readable description of the problem for the help message.
Definition FlowGenericProblem_impl.hpp:97
Scalar porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
Direct indexed access to the porosity.
Definition FlowGenericProblem_impl.hpp:305
Scalar rockCompressibility(unsigned globalSpaceIdx) const
Direct access to rock compressibility.
Definition FlowGenericProblem_impl.hpp:275
unsigned miscnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant MISC region given a cell index.
Definition FlowGenericProblem_impl.hpp:688
unsigned satnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant saturation function region given a cell index.
Definition FlowGenericProblem_impl.hpp:678
unsigned plmixnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant PLMIXNUM (for polymer module) region given a cell index.
Definition FlowGenericProblem_impl.hpp:698
static std::string helpPreamble(int, const char **argv)
Returns the string that is printed before the list of command line parameters in the help message.
Definition FlowGenericProblem_impl.hpp:82
bool shouldWriteRestartFile() const
Returns true if an eWoms restart file should be written to disk.
Definition FlowGenericProblem.hpp:302
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:92
const WellModel & wellModel() const
Returns a reference to the ECL well manager used by the problem.
Definition FlowProblem.hpp:1001
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
Direct access to the transmissibility between two elements.
Definition FlowProblem.hpp:538
static int handlePositionalParameter(std::function< void(const std::string &, const std::string &)> addKey, std::set< std::string > &seenParams, std::string &errorMsg, int, const char **argv, int paramIdx, int)
Handles positional command line parameters.
Definition FlowProblem.hpp:190
const DimMatrix & intrinsicPermeability(unsigned globalElemIdx) const
This method returns the intrinsic permeability tensor given a global element index.
Definition FlowProblem.hpp:520
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:815
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:670
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition FlowProblem.hpp:396
unsigned satnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:823
Scalar thermalHalfTransmissibilityOut(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:628
Scalar transmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:585
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn, const unsigned globalSpaceIdxOut) const
Definition FlowProblem.hpp:605
Scalar rockCompressibility(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:704
Scalar maxOilSaturation(unsigned globalDofIdx) const
Returns an element's historic maximum oil phase saturation that was observed during the simulation.
Definition FlowProblem.hpp:899
std::string name() const
The problem name.
Definition FlowProblem.hpp:854
LhsEval rockCompTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1089
void endIteration()
Called by the simulator after each Newton-Raphson iteration.
Definition FlowProblem.hpp:406
FlowProblem(Simulator &simulator)
Definition FlowProblem.hpp:209
const Vanguard::TransmissibilityType & eclTransmissibilities() const
Return a reference to the object that handles the "raw" transmissibilities.
Definition FlowProblem.hpp:651
Scalar nextTimeStepSize() const
Propose the size of the next time step to the simulator.
Definition FlowProblem.hpp:1022
LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1614
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:759
Scalar transmissibility(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:527
virtual void beginEpisode()
Called by the simulator before an episode begins.
Definition FlowProblem.hpp:306
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
Direct indexed acces to the depth of an degree of freedom [m].
Definition FlowProblem.hpp:695
LhsEval wellTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Return the well transmissibility multiplier due to rock changues.
Definition FlowProblem.hpp:1101
std::shared_ptr< const EclMaterialLawManager > materialLawManager() const
Returns the ECL material law manager.
Definition FlowProblem.hpp:771
unsigned plmixnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:839
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:861
Scalar thermalHalfTransmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:641
unsigned miscnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:831
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:506
Scalar maxPolymerAdsorption(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the max polymer adsorption value.
Definition FlowProblem.hpp:848
void beginTimeStep()
Called by the simulator before each time integration.
Definition FlowProblem.hpp:363
Scalar thermalHalfTransmissibilityIn(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:615
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
Sets an element's maximum oil phase saturation observed during the simulation.
Definition FlowProblem.hpp:916
const SolidEnergyLawParams & solidEnergyLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the parameters for the energy storage law of the rock.
Definition FlowProblem.hpp:746
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the dispersivity for a face i.e.
Definition FlowProblem.hpp:565
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblem.hpp:179
virtual void endEpisode()
Called by the simulator after the end of an episode.
Definition FlowProblem.hpp:468
void source(RateVector &rate, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition FlowProblem.hpp:960
virtual void endTimeStep()
Called by the simulator after each time integration.
Definition FlowProblem.hpp:416
virtual void initialSolutionApplied()
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition FlowProblem.hpp:927
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:724
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:572
void writeOutput(bool verbose=true)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:491
Scalar diffusivity(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:547
void deserialize(Restarter &res)
This method restores the complete state of the problem and its sub-objects from disk.
Definition FlowProblem.hpp:272
std::shared_ptr< EclMaterialLawManager > materialLawManager()
Definition FlowProblem.hpp:807
Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the depth of an degree of freedom [m].
Definition FlowProblem.hpp:683
LhsEval rockCompPoroMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the porosity multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1049
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:595
void serialize(Restarter &res)
This method writes the complete state of the problem and its subobjects to disk.
Definition FlowProblem.hpp:291
Scalar rockReferencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:714
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the transmissibility for a face i.e.
Definition FlowProblem.hpp:558
A random-access container which stores data attached to a grid's degrees of freedom in a prefetch fri...
Definition pffgridvector.hh:49
This class is intend to be a relperm diagnostics, to detect wrong input of relperm table and endpoint...
Definition RelpermDiagnostics.hpp:50
void diagnosis(const EclipseState &eclState, const CartesianIndexMapper &cartesianIndexMapper)
This function is used to diagnosis relperm in eclipse data file.
Definition RelpermDiagnostics.cpp:826
A class which handles tracers as specified in by ECL.
Definition TracerModel.hpp:68
This file contains definitions related to directional mobilities.
The base class for the element-centered finite-volume discretization scheme.
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
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition propertysystem.hh:226
Definition FlowProblem.hpp:1666
Definition FlowProblem.hpp:1455