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/Schedule/Schedule.hpp>
41#include <opm/input/eclipse/Units/Units.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 <cstddef>
76#include <functional>
77#include <set>
78#include <stdexcept>
79#include <string>
80#include <vector>
81
82namespace Opm {
83
90template <class TypeTag>
91class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
92 , public FlowGenericProblem<GetPropType<TypeTag, Properties::GridView>,
93 GetPropType<TypeTag, Properties::FluidSystem>>
94{
95protected:
99 using Implementation = GetPropType<TypeTag, Properties::Problem>;
100
108
109 // Grid and world dimension
110 enum { dim = GridView::dimension };
111 enum { dimWorld = GridView::dimensionworld };
112
113 // copy some indices for convenience
115 enum { numPhases = FluidSystem::numPhases };
116 enum { numComponents = FluidSystem::numComponents };
117
118 enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
119 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
120 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
121 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
122 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
123 enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
124 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
125 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
126 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
127 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
128 enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
129 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
130 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
131 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
132 enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
133
134 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
135 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
136 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
137
138 // TODO: later, gasCompIdx, oilCompIdx and waterCompIdx should go to the FlowProblemBlackoil in the future
139 // we do not want them in the compositional setting
140 enum { gasCompIdx = FluidSystem::gasCompIdx };
141 enum { oilCompIdx = FluidSystem::oilCompIdx };
142 enum { waterCompIdx = FluidSystem::waterCompIdx };
143
147 using Element = typename GridView::template Codim<0>::Entity;
149 using EclMaterialLawManager = typename GetProp<TypeTag, Properties::MaterialLaw>::EclMaterialLawManager;
150 using EclThermalLawManager = typename GetProp<TypeTag, Properties::SolidEnergyLaw>::EclThermalLawManager;
151 using MaterialLawParams = typename EclMaterialLawManager::MaterialLawParams;
152 using SolidEnergyLawParams = typename EclThermalLawManager::SolidEnergyLawParams;
153 using ThermalConductionLawParams = typename EclThermalLawManager::ThermalConductionLawParams;
161
162 using Toolbox = MathToolbox<Evaluation>;
163 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
164
165
167 using DirectionalMobilityPtr = Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
168
169public:
175 using BaseType::porosity;
176
180 static void registerParameters()
181 {
182 ParentType::registerParameters();
183
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 if (! Parameters::Get<Parameters::CheckSatfuncConsistency>()) {
228 // User did not enable the "new" saturation function consistency
229 // check module. Run the original checker instead. This is a
230 // temporary measure.
232 relpermDiagnostics.diagnosis(simulator.vanguard().eclState(),
233 simulator.vanguard().levelCartesianIndexMapper());
234 }
235 }
236
237 virtual ~FlowProblem() = default;
238
239 void prefetch(const Element& elem) const
240 { this->pffDofData_.prefetch(elem); }
241
253 template <class Restarter>
255 {
256 // reload the current episode/report step from the deck
257 this->beginEpisode();
258
259 // deserialize the wells
260 wellModel_.deserialize(res);
261
262 // deserialize the aquifer
263 aquiferModel_.deserialize(res);
264 }
265
272 template <class Restarter>
274 {
275 wellModel_.serialize(res);
276
277 aquiferModel_.serialize(res);
278 }
279
280 int episodeIndex() const
281 {
282 return std::max(this->simulator().episodeIndex(), 0);
283 }
284
288 virtual void beginEpisode()
289 {
291 // Proceed to the next report step
292 auto& simulator = this->simulator();
293 int episodeIdx = simulator.episodeIndex();
294 auto& eclState = simulator.vanguard().eclState();
295 const auto& schedule = simulator.vanguard().schedule();
296 const auto& events = schedule[episodeIdx].events();
297
298 if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
299 // bring the contents of the keywords to the current state of the SCHEDULE
300 // section.
301 //
302 // TODO (?): make grid topology changes possible (depending on what exactly
303 // has changed, the grid may need be re-created which has some serious
304 // implications on e.g., the solution of the simulation.)
305 const auto& miniDeck = schedule[episodeIdx].geo_keywords();
306 const auto& cc = simulator.vanguard().grid().comm();
307 eclState.apply_schedule_keywords( miniDeck );
308 eclBroadcast(cc, eclState.getTransMult() );
309
310 // Re-ordering in case of ALUGrid
311 std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
312 return simulator.vanguard().gridEquilIdxToGridIdx(i);
313 };
314
315 // re-compute all quantities which may possibly be affected.
316 using TransUpdateQuantities = typename Vanguard::TransmissibilityType::TransUpdateQuantities;
317 transmissibilities_.update(true, TransUpdateQuantities::All, equilGridToGrid);
318 this->referencePorosity_[1] = this->referencePorosity_[0];
319 updateReferencePorosity_();
320 updatePffDofData_();
321 this->model().linearizer().updateDiscretizationParameters();
322 }
323
324 bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex());
325
326 // set up the wells for the next episode.
327 wellModel_.beginEpisode();
328
329 // set up the aquifers for the next episode.
330 aquiferModel_.beginEpisode();
331
332 // set the size of the initial time step of the episode
333 Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
334 // negative value of initialTimeStepSize_ indicates no active limit from TSINIT or NEXTSTEP
335 if ( (episodeIdx == 0 || tuningEvent) && this->initialTimeStepSize_ > 0)
336 // allow the size of the initial time step to be set via an external parameter
337 // if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
338 dt = std::min(dt, this->initialTimeStepSize_);
339 simulator.setTimeStepSize(dt);
340 }
341
346 {
348 const int episodeIdx = this->episodeIndex();
349 const int timeStepSize = this->simulator().timeStepSize();
350
351 this->beginTimeStep_(enableExperiments,
353 this->simulator().timeStepIndex(),
354 this->simulator().startTime(),
355 this->simulator().time(),
356 timeStepSize,
357 this->simulator().endTime());
358
359 // update maximum water saturation and minimum pressure
360 // used when ROCKCOMP is activated
361 // Do not update max RS first step after a restart
362 this->updateExplicitQuantities_(episodeIdx, timeStepSize, first_step_ && (episodeIdx > 0));
363 first_step_ = false;
364
365 if (nonTrivialBoundaryConditions()) {
366 this->model().linearizer().updateBoundaryConditionData();
367 }
368
369 wellModel_.beginTimeStep();
370 aquiferModel_.beginTimeStep();
371 tracerModel_.beginTimeStep();
372
373 }
374
379 {
381 wellModel_.beginIteration();
382 aquiferModel_.beginIteration();
383 }
384
389 {
391 wellModel_.endIteration();
392 aquiferModel_.endIteration();
393 }
394
398 virtual void endTimeStep()
399 {
401
402#ifndef NDEBUG
404 // in debug mode, we don't care about performance, so we check
405 // if the model does the right thing (i.e., the mass change
406 // inside the whole reservoir must be equivalent to the fluxes
407 // over the grid's boundaries plus the source rates specified by
408 // the problem).
409 const int rank = this->simulator().gridView().comm().rank();
410 if (rank == 0) {
411 std::cout << "checking conservativeness of solution\n";
412 }
413
414 this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
415 if (rank == 0) {
416 std::cout << "solution is sufficiently conservative\n";
417 }
418 }
419#endif // NDEBUG
420
421 auto& simulator = this->simulator();
422 simulator.setTimeStepIndex(simulator.timeStepIndex()+1);
423
424 this->wellModel_.endTimeStep();
425 this->aquiferModel_.endTimeStep();
426 this->tracerModel_.endTimeStep();
427
428 // Compute flux for output
429 this->model().linearizer().updateFlowsInfo();
430
431 if (this->enableDriftCompensation_) {
433
434 const auto& residual = this->model().linearizer().residual();
435
436 for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
437 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(globalDofIdx);
438 this->drift_[sfcdofIdx] = residual[sfcdofIdx] * simulator.timeStepSize();
439
441 this->drift_[sfcdofIdx] *= this->model().dofTotalVolume(sfcdofIdx);
442 }
443 }
444 }
445 }
446
450 virtual void endEpisode()
451 {
452 const int episodeIdx = this->episodeIndex();
453
454 this->wellModel_.endEpisode();
455 this->aquiferModel_.endEpisode();
456
457 const auto& schedule = this->simulator().vanguard().schedule();
458
459 // End simulation when completed.
460 if (episodeIdx + 1 >= static_cast<int>(schedule.size()) - 1) {
461 this->simulator().setFinished(true);
462 return;
463 }
464
465 // Otherwise, start next episode (report step).
466 this->simulator().startNextEpisode(schedule.stepLength(episodeIdx + 1));
467 }
468
473 virtual void writeOutput(bool verbose)
474 {
476 // use the generic code to prepare the output fields and to
477 // write the desired VTK files.
478 if (Parameters::Get<Parameters::EnableWriteAllSolutions>() ||
479 this->simulator().episodeWillBeOver()) {
480 ParentType::writeOutput(verbose);
481 }
482 }
483
487 template <class Context>
488 const DimMatrix& intrinsicPermeability(const Context& context,
489 unsigned spaceIdx,
490 unsigned timeIdx) const
491 {
492 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
493 return transmissibilities_.permeability(globalSpaceIdx);
494 }
495
502 const DimMatrix& intrinsicPermeability(unsigned globalElemIdx) const
503 { return transmissibilities_.permeability(globalElemIdx); }
504
508 template <class Context>
509 Scalar transmissibility(const Context& context,
510 [[maybe_unused]] unsigned fromDofLocalIdx,
511 unsigned toDofLocalIdx) const
512 {
514 return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
515 }
516
520 Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
521 {
522 return transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
523 }
524
528 template <class Context>
529 Scalar diffusivity(const Context& context,
530 [[maybe_unused]] unsigned fromDofLocalIdx,
531 unsigned toDofLocalIdx) const
532 {
534 return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
535 }
536
540 Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
541 return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
542 }
543
547 Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
548 return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
549 }
550
555 const unsigned boundaryFaceIdx) const
556 {
557 return transmissibilities_.thermalTransmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
558 }
559
560
561
562
566 template <class Context>
567 Scalar transmissibilityBoundary(const Context& elemCtx,
568 unsigned boundaryFaceIdx) const
569 {
570 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
571 return transmissibilities_.transmissibilityBoundary(elemIdx, boundaryFaceIdx);
572 }
573
578 const unsigned boundaryFaceIdx) const
579 {
580 return transmissibilities_.transmissibilityBoundary(globalSpaceIdx, boundaryFaceIdx);
581 }
582
583
588 const unsigned globalSpaceIdxOut) const
589 {
590 return transmissibilities_.thermalHalfTrans(globalSpaceIdxIn,globalSpaceIdxOut);
591 }
592
596 template <class Context>
597 Scalar thermalHalfTransmissibilityIn(const Context& context,
598 unsigned faceIdx,
599 unsigned timeIdx) const
600 {
601 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
602 unsigned toDofLocalIdx = face.exteriorIndex();
603 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransIn;
604 }
605
609 template <class Context>
610 Scalar thermalHalfTransmissibilityOut(const Context& context,
611 unsigned faceIdx,
612 unsigned timeIdx) const
613 {
614 const auto& face = context.stencil(timeIdx).interiorFace(faceIdx);
615 unsigned toDofLocalIdx = face.exteriorIndex();
616 return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTransOut;
617 }
618
622 template <class Context>
623 Scalar thermalHalfTransmissibilityBoundary(const Context& elemCtx,
624 unsigned boundaryFaceIdx) const
625 {
626 unsigned elemIdx = elemCtx.globalSpaceIndex(/*dofIdx=*/0, /*timeIdx=*/0);
627 return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
628 }
629
633 const typename Vanguard::TransmissibilityType& eclTransmissibilities() const
634 { return transmissibilities_; }
635
636
637 const TracerModel& tracerModel() const
638 { return tracerModel_; }
639
640 TracerModel& tracerModel()
641 { return tracerModel_; }
642
651 template <class Context>
652 Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
653 {
654 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
655 return this->porosity(globalSpaceIdx, timeIdx);
656 }
657
664 template <class Context>
665 Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
666 {
667 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
668 return this->dofCenterDepth(globalSpaceIdx);
669 }
670
677 Scalar dofCenterDepth(unsigned globalSpaceIdx) const
678 {
679 return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
680 }
681
685 template <class Context>
686 Scalar rockCompressibility(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
687 {
688 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
689 return this->rockCompressibility(globalSpaceIdx);
690 }
691
695 template <class Context>
696 Scalar rockReferencePressure(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
697 {
698 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
700 }
701
706 {
707 const auto& rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
708 if (rock_config.store()) {
709 return asImp_().initialFluidState(globalSpaceIdx).pressure(refPressurePhaseIdx_());
710 }
711 else {
712 if (this->rockParams_.empty())
713 return 1e5;
714
715 unsigned tableIdx = 0;
716 if (!this->rockTableIdx_.empty()) {
717 tableIdx = this->rockTableIdx_[globalSpaceIdx];
718 }
719 return this->rockParams_[tableIdx].referencePressure;
720 }
721 }
722
726 template <class Context>
727 const MaterialLawParams& materialLawParams(const Context& context,
728 unsigned spaceIdx, unsigned timeIdx) const
729 {
730 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
731 return this->materialLawParams(globalSpaceIdx);
732 }
733
734 const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const
735 {
736 return materialLawManager_->materialLawParams(globalDofIdx);
737 }
738
739 const MaterialLawParams& materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const
740 {
741 return materialLawManager_->materialLawParams(globalDofIdx, facedir);
742 }
743
747 template <class Context>
748 const SolidEnergyLawParams&
749 solidEnergyLawParams(const Context& context,
750 unsigned spaceIdx,
751 unsigned timeIdx) const
752 {
753 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
754 return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
755 }
756
760 template <class Context>
761 const ThermalConductionLawParams &
762 thermalConductionLawParams(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
763 {
764 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
765 return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
766 }
767
774 std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
775 { return materialLawManager_; }
776
777 template <class FluidState>
778 void updateRelperms(
779 std::array<Evaluation,numPhases> &mobility,
780 DirectionalMobilityPtr &dirMob,
781 FluidState &fluidState,
782 unsigned globalSpaceIdx) const
783 {
784 OPM_TIMEBLOCK_LOCAL(updateRelperms);
785 {
786 // calculate relative permeabilities. note that we store the result into the
787 // mobility_ class attribute. the division by the phase viscosity happens later.
789 MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
790 Valgrind::CheckDefined(mobility);
791 }
792 if (materialLawManager_->hasDirectionalRelperms()
793 || materialLawManager_->hasDirectionalImbnum())
794 {
795 using Dir = FaceDir::DirEnum;
796 constexpr int ndim = 3;
797 dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
798 Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
799 for (int i = 0; i<ndim; i++) {
801 auto& mob_array = dirMob->getArray(i);
802 MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
803 }
804 }
805 }
806
810 std::shared_ptr<EclMaterialLawManager> materialLawManager()
811 { return materialLawManager_; }
812
817 template <class Context>
818 unsigned pvtRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
819 { return pvtRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
820
825 template <class Context>
826 unsigned satnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
827 { return this->satnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
828
833 template <class Context>
834 unsigned miscnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
835 { return this->miscnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
836
841 template <class Context>
842 unsigned plmixnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
843 { return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); }
844
845 // TODO: polymer related might need to go to the blackoil side
850 template <class Context>
851 Scalar maxPolymerAdsorption(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
852 { return this->maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); }
853
857 std::string name() const
858 { return this->simulator().vanguard().caseName(); }
859
863 template <class Context>
864 Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
865 {
866 // use the initial temperature of the DOF if temperature is not a primary
867 // variable
868 unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
869 return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0);
870 }
871
872
873 Scalar temperature(unsigned globalDofIdx, unsigned /*timeIdx*/) const
874 {
875 // use the initial temperature of the DOF if temperature is not a primary
876 // variable
877 return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0);
878 }
879
880 const SolidEnergyLawParams&
882 unsigned /*timeIdx*/) const
883 {
884 return this->thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
885 }
886 const ThermalConductionLawParams &
888 unsigned /*timeIdx*/)const
889 {
890 return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
891 }
892
902 Scalar maxOilSaturation(unsigned globalDofIdx) const
903 {
904 if (!this->vapparsActive(this->episodeIndex()))
905 return 0.0;
906
907 return this->maxOilSaturation_[globalDofIdx];
908 }
909
919 void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
920 {
921 if (!this->vapparsActive(this->episodeIndex()))
922 return;
923
924 this->maxOilSaturation_[globalDofIdx] = value;
925 }
926
931 {
932 // Calculate all intensive quantities.
933 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/0);
934
935 // We also need the intensive quantities for timeIdx == 1
936 // corresponding to the start of the current timestep, if we
937 // do not use the storage cache, or if we cannot recycle the
938 // first iteration storage.
939 if (!this->model().enableStorageCache() || !this->recycleFirstIterationStorage()) {
940 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/1);
941 }
942
943 // initialize the wells. Note that this needs to be done after initializing the
944 // intrinsic permeabilities and the after applying the initial solution because
945 // the well model uses these...
946 wellModel_.init();
947
948 aquiferModel_.initialSolutionApplied();
949
950 const bool invalidateFromHyst = updateHysteresis_();
951 if (invalidateFromHyst) {
953 this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
954 }
955 }
956
962 template <class Context>
963 void source(RateVector& rate,
964 const Context& context,
965 unsigned spaceIdx,
966 unsigned timeIdx) const
967 {
968 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
969 source(rate, globalDofIdx, timeIdx);
970 }
971
972 void source(RateVector& rate,
973 unsigned globalDofIdx,
974 unsigned timeIdx) const
975 {
977 rate = 0.0;
978
979 // Add well contribution to source here.
980 wellModel_.computeTotalRatesForDof(rate, globalDofIdx);
981
982 // convert the source term from the total mass rate of the
983 // cell to the one per unit of volume as used by the model.
984 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
985 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
986
987 Valgrind::CheckDefined(rate[eqIdx]);
988 assert(isfinite(rate[eqIdx]));
989 }
990
991 // Add non-well sources.
992 addToSourceDense(rate, globalDofIdx, timeIdx);
993 }
994
995 virtual void addToSourceDense(RateVector& rate,
996 unsigned globalDofIdx,
997 unsigned timeIdx) const = 0;
998
1004 const WellModel& wellModel() const
1005 { return wellModel_; }
1006
1007 WellModel& wellModel()
1008 { return wellModel_; }
1009
1010 const AquiferModel& aquiferModel() const
1011 { return aquiferModel_; }
1012
1013 AquiferModel& mutableAquiferModel()
1014 { return aquiferModel_; }
1015
1016 bool nonTrivialBoundaryConditions() const
1017 { return nonTrivialBoundaryConditions_; }
1018
1025 Scalar nextTimeStepSize() const
1026 {
1028 // allow external code to do the timestepping
1029 if (this->nextTimeStepSize_ > 0.0)
1030 return this->nextTimeStepSize_;
1031
1032 const auto& simulator = this->simulator();
1033 int episodeIdx = simulator.episodeIndex();
1034
1035 // for the initial episode, we use a fixed time step size
1036 if (episodeIdx < 0)
1037 return this->initialTimeStepSize_;
1038
1039 // ask the newton method for a suggestion. This suggestion will be based on how
1040 // well the previous time step converged. After that, apply the runtime time
1041 // stepping constraints.
1042 const auto& newtonMethod = this->model().newtonMethod();
1043 return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
1044 }
1045
1051 template <class LhsEval>
1052 LhsEval rockCompPoroMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1053 {
1055 if (this->rockCompPoroMult_.empty() && this->rockCompPoroMultWc_.empty())
1056 return 1.0;
1057
1058 unsigned tableIdx = 0;
1059 if (!this->rockTableIdx_.empty())
1060 tableIdx = this->rockTableIdx_[elementIdx];
1061
1062 const auto& fs = intQuants.fluidState();
1063 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1064 const auto& rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
1065 if (!this->minRefPressure_.empty())
1066 // The pore space change is irreversible
1068 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1069 this->minRefPressure_[elementIdx]);
1070
1071 if (!this->overburdenPressure_.empty())
1072 effectivePressure -= this->overburdenPressure_[elementIdx];
1073
1074 if (rock_config.store()) {
1075 effectivePressure -= asImp_().initialFluidState(elementIdx).pressure(refPressurePhaseIdx_());
1076 }
1077
1078 if (!this->rockCompPoroMult_.empty()) {
1079 return this->rockCompPoroMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
1080 }
1081
1082 // water compaction
1083 assert(!this->rockCompPoroMultWc_.empty());
1084 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1085 LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
1086
1087 return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
1088 }
1089
1095 template <class LhsEval>
1096 LhsEval rockCompTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1097 {
1098 const bool implicit = !this->explicitRockCompaction_;
1099 return implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<LhsEval>(intQuants, elementIdx)
1100 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1101 }
1102
1103
1107 template <class LhsEval>
1108 LhsEval wellTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1109 {
1111
1112 const bool implicit = !this->explicitRockCompaction_;
1113 double trans_mult = implicit ? this->simulator().problem().template computeRockCompTransMultiplier_<double>(intQuants, elementIdx)
1114 : this->simulator().problem().getRockCompTransMultVal(elementIdx);
1115 trans_mult *= this->simulator().problem().template permFactTransMultiplier<double>(intQuants, elementIdx);
1116
1117 return trans_mult;
1118 }
1119
1120 std::pair<BCType, RateVector> boundaryCondition(const unsigned int globalSpaceIdx, const int directionId) const
1121 {
1122 OPM_TIMEBLOCK_LOCAL(boundaryCondition);
1123 if (!nonTrivialBoundaryConditions_) {
1124 return { BCType::NONE, RateVector(0.0) };
1125 }
1126 FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId);
1127 const auto& schedule = this->simulator().vanguard().schedule();
1128 if (bcindex_(dir)[globalSpaceIdx] == 0) {
1129 return { BCType::NONE, RateVector(0.0) };
1130 }
1131 if (schedule[this->episodeIndex()].bcprop.size() == 0) {
1132 return { BCType::NONE, RateVector(0.0) };
1133 }
1134 const auto& bc = schedule[this->episodeIndex()].bcprop[bcindex_(dir)[globalSpaceIdx]];
1135 if (bc.bctype!=BCType::RATE) {
1136 return { bc.bctype, RateVector(0.0) };
1137 }
1138
1139 RateVector rate = 0.0;
1140 switch (bc.component) {
1141 case BCComponent::OIL:
1142 rate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = bc.rate;
1143 break;
1144 case BCComponent::GAS:
1145 rate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = bc.rate;
1146 break;
1147 case BCComponent::WATER:
1148 rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate;
1149 break;
1150 case BCComponent::SOLVENT:
1151 this->handleSolventBC(bc, rate);
1152 break;
1153 case BCComponent::POLYMER:
1154 this->handlePolymerBC(bc, rate);
1155 break;
1156 case BCComponent::NONE:
1157 throw std::logic_error("you need to specify the component when RATE type is set in BC");
1158 break;
1159 }
1160 //TODO add support for enthalpy rate
1161 return {bc.bctype, rate};
1162 }
1163
1164
1165 template<class Serializer>
1166 void serializeOp(Serializer& serializer)
1167 {
1168 serializer(static_cast<BaseType&>(*this));
1169 serializer(drift_);
1170 serializer(wellModel_);
1171 serializer(aquiferModel_);
1172 serializer(tracerModel_);
1173 serializer(*materialLawManager_);
1174 }
1175
1176private:
1177 Implementation& asImp_()
1178 { return *static_cast<Implementation *>(this); }
1179
1180 const Implementation& asImp_() const
1181 { return *static_cast<const Implementation *>(this); }
1182
1183protected:
1184 template<class UpdateFunc>
1185 void updateProperty_(const std::string& failureMsg,
1187 {
1189 const auto& model = this->simulator().model();
1190 const auto& primaryVars = model.solution(/*timeIdx*/0);
1191 const auto& vanguard = this->simulator().vanguard();
1192 std::size_t numGridDof = primaryVars.size();
1193 OPM_BEGIN_PARALLEL_TRY_CATCH();
1194#ifdef _OPENMP
1195#pragma omp parallel for
1196#endif
1197 for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
1198 const auto& iq = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/ 0);
1199 func(dofIdx, iq);
1200 }
1201 OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
1202 }
1203
1204 bool updateMaxOilSaturation_()
1205 {
1207 int episodeIdx = this->episodeIndex();
1208
1209 // we use VAPPARS
1210 if (this->vapparsActive(episodeIdx)) {
1211 this->updateProperty_("FlowProblem::updateMaxOilSaturation_() failed:",
1212 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1213 {
1214 this->updateMaxOilSaturation_(compressedDofIdx,iq);
1215 });
1216 return true;
1217 }
1218
1219 return false;
1220 }
1221
1222 bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1223 {
1225 const auto& fs = iq.fluidState();
1226 const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
1227 auto& mos = this->maxOilSaturation_;
1228 if(mos[compressedDofIdx] < So){
1229 mos[compressedDofIdx] = So;
1230 return true;
1231 }else{
1232 return false;
1233 }
1234 }
1235
1236 bool updateMaxWaterSaturation_()
1237 {
1239 // water compaction is activated in ROCKCOMP
1240 if (this->maxWaterSaturation_.empty())
1241 return false;
1242
1243 this->maxWaterSaturation_[/*timeIdx=*/1] = this->maxWaterSaturation_[/*timeIdx=*/0];
1244 this->updateProperty_("FlowProblem::updateMaxWaterSaturation_() failed:",
1245 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1246 {
1247 this->updateMaxWaterSaturation_(compressedDofIdx,iq);
1248 });
1249 return true;
1250 }
1251
1252
1253 bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1254 {
1256 const auto& fs = iq.fluidState();
1257 const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
1258 auto& mow = this->maxWaterSaturation_;
1259 if(mow[compressedDofIdx]< Sw){
1260 mow[compressedDofIdx] = Sw;
1261 return true;
1262 }else{
1263 return false;
1264 }
1265 }
1266
1267 bool updateMinPressure_()
1268 {
1270 // IRREVERS option is used in ROCKCOMP
1271 if (this->minRefPressure_.empty())
1272 return false;
1273
1274 this->updateProperty_("FlowProblem::updateMinPressure_() failed:",
1275 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1276 {
1277 this->updateMinPressure_(compressedDofIdx,iq);
1278 });
1279 return true;
1280 }
1281
1282 bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
1284 const auto& fs = iq.fluidState();
1285 const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
1286 auto& min_pressures = this->minRefPressure_;
1289 return true;
1290 }else{
1291 return false;
1292 }
1293 }
1294
1295 // \brief Function to assign field properties of type double, on the leaf grid view.
1296 //
1297 // For CpGrid with local grid refinement, the field property of a cell on the leaf
1298 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
1299 std::function<std::vector<double>(const FieldPropsManager&, const std::string&)>
1300 fieldPropDoubleOnLeafAssigner_()
1301 {
1302 const auto& lookup = this->lookUpData_;
1303 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString)
1304 {
1305 return lookup.assignFieldPropsDoubleOnLeaf(fieldPropManager, propString);
1306 };
1307 }
1308
1309 // \brief Function to assign field properties of type int, unsigned int, ..., on the leaf grid view.
1310 //
1311 // For CpGrid with local grid refinement, the field property of a cell on the leaf
1312 // is inherited from its parent or equivalent (when has no parent) cell on level zero.
1313 template<typename IntType>
1314 std::function<std::vector<IntType>(const FieldPropsManager&, const std::string&, bool)>
1315 fieldPropIntTypeOnLeafAssigner_()
1316 {
1317 const auto& lookup = this->lookUpData_;
1318 return [&lookup](const FieldPropsManager& fieldPropManager, const std::string& propString, bool needsTranslation)
1319 {
1321 };
1322 }
1323
1324 void readMaterialParameters_()
1325 {
1327 const auto& simulator = this->simulator();
1328 const auto& vanguard = simulator.vanguard();
1329 const auto& eclState = vanguard.eclState();
1330
1331 // the PVT and saturation region numbers
1332 OPM_BEGIN_PARALLEL_TRY_CATCH();
1333 this->updatePvtnum_();
1334 this->updateSatnum_();
1335
1336 // the MISC region numbers (solvent model)
1337 this->updateMiscnum_();
1338 // the PLMIX region numbers (polymer model)
1339 this->updatePlmixnum_();
1340
1341 OPM_END_PARALLEL_TRY_CATCH("Invalid region numbers: ", vanguard.gridView().comm());
1343 // porosity
1344 updateReferencePorosity_();
1345 this->referencePorosity_[1] = this->referencePorosity_[0];
1347
1349 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
1350 materialLawManager_ = std::make_shared<EclMaterialLawManager>();
1351 materialLawManager_->initFromState(eclState);
1352 materialLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1353 this-> template fieldPropIntTypeOnLeafAssigner_<int>(),
1354 this-> lookupIdxOnLevelZeroAssigner_());
1356 }
1357
1358 void readThermalParameters_()
1359 {
1360 if constexpr (enableEnergy)
1361 {
1362 const auto& simulator = this->simulator();
1363 const auto& vanguard = simulator.vanguard();
1364 const auto& eclState = vanguard.eclState();
1365
1366 // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
1367 thermalLawManager_ = std::make_shared<EclThermalLawManager>();
1368 thermalLawManager_->initParamsForElements(eclState, this->model().numGridDof(),
1369 this-> fieldPropDoubleOnLeafAssigner_(),
1371 }
1372 }
1373
1374 void updateReferencePorosity_()
1375 {
1376 const auto& simulator = this->simulator();
1377 const auto& vanguard = simulator.vanguard();
1378 const auto& eclState = vanguard.eclState();
1379
1380 std::size_t numDof = this->model().numGridDof();
1381
1382 this->referencePorosity_[/*timeIdx=*/0].resize(numDof);
1383
1384 const auto& fp = eclState.fieldProps();
1385 const std::vector<double> porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV");
1386 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
1387 int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(dofIdx);
1388 Scalar poreVolume = porvData[dofIdx];
1389
1390 // we define the porosity as the accumulated pore volume divided by the
1391 // geometric volume of the element. Note that -- in pathetic cases -- it can
1392 // be larger than 1.0!
1393 Scalar dofVolume = simulator.model().dofTotalVolume(sfcdofIdx);
1394 assert(dofVolume > 0.0);
1395 this->referencePorosity_[/*timeIdx=*/0][sfcdofIdx] = poreVolume/dofVolume;
1396 }
1397 }
1398
1399 virtual void readInitialCondition_()
1400 {
1401 // TODO: whether we should move this to FlowProblemBlackoil
1402 const auto& simulator = this->simulator();
1403 const auto& vanguard = simulator.vanguard();
1404 const auto& eclState = vanguard.eclState();
1405
1406 if (eclState.getInitConfig().hasEquil())
1407 readEquilInitialCondition_();
1408 else
1409 readExplicitInitialCondition_();
1410
1411 //initialize min/max values
1412 std::size_t numElems = this->model().numGridDof();
1413 for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
1414 const auto& fs = asImp_().initialFluidStates()[elemIdx];
1415 if (!this->maxWaterSaturation_.empty() && waterPhaseIdx > -1)
1416 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
1417 if (!this->maxOilSaturation_.empty() && oilPhaseIdx > -1)
1418 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
1419 if (!this->minRefPressure_.empty() && refPressurePhaseIdx_() > -1)
1420 this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
1421 }
1422 }
1423
1424 virtual void readEquilInitialCondition_() = 0;
1425 virtual void readExplicitInitialCondition_() = 0;
1426
1427 // update the hysteresis parameters of the material laws for the whole grid
1428 bool updateHysteresis_()
1429 {
1430 if (!materialLawManager_->enableHysteresis())
1431 return false;
1432
1433 // we need to update the hysteresis data for _all_ elements (i.e., not just the
1434 // interior ones) to avoid desynchronization of the processes in the parallel case!
1435 this->updateProperty_("FlowProblem::updateHysteresis_() failed:",
1436 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
1437 {
1438 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
1439 });
1440 return true;
1441 }
1442
1443
1444 bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
1445 {
1446 OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
1447 materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
1448 //TODO change materials to give a bool
1449 return true;
1450 }
1451
1452 Scalar getRockCompTransMultVal(std::size_t dofIdx) const
1453 {
1454 if (this->rockCompTransMultVal_.empty())
1455 return 1.0;
1456
1457 return this->rockCompTransMultVal_[dofIdx];
1458 }
1459
1460protected:
1462 {
1464 ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
1467 Scalar transmissibility;
1468 };
1469
1470 // update the prefetch friendly data object
1471 void updatePffDofData_()
1472 {
1473 const auto& distFn =
1474 [this](PffDofData_& dofData,
1475 const Stencil& stencil,
1476 unsigned localDofIdx)
1477 -> void
1478 {
1479 const auto& elementMapper = this->model().elementMapper();
1480
1481 unsigned globalElemIdx = elementMapper.index(stencil.entity(localDofIdx));
1482 if (localDofIdx != 0) {
1483 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(/*dofIdx=*/0));
1484 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
1485
1486 if constexpr (enableEnergy) {
1487 *dofData.thermalHalfTransIn = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
1488 *dofData.thermalHalfTransOut = transmissibilities_.thermalHalfTrans(globalElemIdx, globalCenterElemIdx);
1489 }
1490 if constexpr (enableDiffusion)
1491 *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
1492 if (enableDispersion)
1493 dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
1494 }
1495 };
1496
1497 pffDofData_.update(distFn);
1498 }
1499
1500 virtual void updateExplicitQuantities_(int episodeIdx, int timeStepSize, bool first_step_after_restart) = 0;
1501
1502 void readBoundaryConditions_()
1503 {
1504 const auto& simulator = this->simulator();
1505 const auto& vanguard = simulator.vanguard();
1506 const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
1507 if (bcconfig.size() > 0) {
1508 nonTrivialBoundaryConditions_ = true;
1509
1510 std::size_t numCartDof = vanguard.cartesianSize();
1511 unsigned numElems = vanguard.gridView().size(/*codim=*/0);
1512 std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
1513
1514 for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
1515 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
1516
1517 bcindex_.resize(numElems, 0);
1519 &vanguard](const auto& bcface,
1520 auto apply)
1521 {
1522 for (int i = bcface.i1; i <= bcface.i2; ++i) {
1523 for (int j = bcface.j1; j <= bcface.j2; ++j) {
1524 for (int k = bcface.k1; k <= bcface.k2; ++k) {
1525 std::array<int, 3> tmp = {i,j,k};
1526 auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
1527 if (elemIdx >= 0)
1528 apply(elemIdx);
1529 }
1530 }
1531 }
1532 };
1533 for (const auto& bcface : bcconfig) {
1534 std::vector<int>& data = bcindex_(bcface.dir);
1535 const int index = bcface.index;
1537 [&data,index](int elemIdx)
1538 { data[elemIdx] = index; });
1539 }
1540 }
1541 }
1542
1543 // this method applies the runtime constraints specified via the deck and/or command
1544 // line parameters for the size of the next time step.
1545 Scalar limitNextTimeStepSize_(Scalar dtNext) const
1546 {
1547 if constexpr (enableExperiments) {
1548 const auto& simulator = this->simulator();
1549 const auto& schedule = simulator.vanguard().schedule();
1550 int episodeIdx = simulator.episodeIndex();
1551
1552 // first thing in the morning, limit the time step size to the maximum size
1553 Scalar maxTimeStepSize = Parameters::Get<Parameters::SolverMaxTimeStepInDays<Scalar>>() * 24 * 60 * 60;
1554 int reportStepIdx = std::max(episodeIdx, 0);
1555 if (this->enableTuning_) {
1556 const auto& tuning = schedule[reportStepIdx].tuning();
1557 maxTimeStepSize = tuning.TSMAXZ;
1558 }
1559
1560 dtNext = std::min(dtNext, maxTimeStepSize);
1561
1562 Scalar remainingEpisodeTime =
1563 simulator.episodeStartTime() + simulator.episodeLength()
1564 - (simulator.startTime() + simulator.time());
1566
1567 // if we would have a small amount of time left over in the current episode, make
1568 // two equal time steps instead of a big and a small one
1569 if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
1570 // note: limiting to the maximum time step size here is probably not strictly
1571 // necessary, but it should not hurt and is more fool-proof
1572 dtNext = std::min(maxTimeStepSize, remainingEpisodeTime/2.0);
1573
1574 if (simulator.episodeStarts()) {
1575 // if a well event occurred, respect the limit for the maximum time step after
1576 // that, too
1577 const auto& events = simulator.vanguard().schedule()[reportStepIdx].events();
1578 bool wellEventOccured =
1579 events.hasEvent(ScheduleEvents::NEW_WELL)
1580 || events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE)
1581 || events.hasEvent(ScheduleEvents::INJECTION_UPDATE)
1582 || events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
1583 if (episodeIdx >= 0 && wellEventOccured && this->maxTimeStepAfterWellEvent_ > 0)
1584 dtNext = std::min(dtNext, this->maxTimeStepAfterWellEvent_);
1585 }
1586 }
1587
1588 return dtNext;
1589 }
1590
1591 int refPressurePhaseIdx_() const {
1592 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1593 return oilPhaseIdx;
1594 }
1595 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1596 return gasPhaseIdx;
1597 }
1598 else {
1599 return waterPhaseIdx;
1600 }
1601 }
1602
1603 void updateRockCompTransMultVal_()
1604 {
1605 const auto& model = this->simulator().model();
1606 std::size_t numGridDof = this->model().numGridDof();
1607 this->rockCompTransMultVal_.resize(numGridDof, 1.0);
1608 for (std::size_t elementIdx = 0; elementIdx < numGridDof; ++elementIdx) {
1609 const auto& iq = *model.cachedIntensiveQuantities(elementIdx, /*timeIdx=*/ 0);
1611 this->rockCompTransMultVal_[elementIdx] = trans_mult;
1612 }
1613 }
1614
1620 template <class LhsEval>
1621 LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities& intQuants, unsigned elementIdx) const
1622 {
1624 if (this->rockCompTransMult_.empty() && this->rockCompTransMultWc_.empty())
1625 return 1.0;
1626
1627 unsigned tableIdx = 0;
1628 if (!this->rockTableIdx_.empty())
1629 tableIdx = this->rockTableIdx_[elementIdx];
1630
1631 const auto& fs = intQuants.fluidState();
1632 LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
1633 const auto& rock_config = this->simulator().vanguard().eclState().getSimulationConfig().rock_config();
1634 if (!this->minRefPressure_.empty())
1635 // The pore space change is irreversible
1637 min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
1638 this->minRefPressure_[elementIdx]);
1639
1640 if (!this->overburdenPressure_.empty())
1641 effectivePressure -= this->overburdenPressure_[elementIdx];
1642
1643 if (rock_config.store()) {
1644 effectivePressure -= asImp_().initialFluidState(elementIdx).pressure(refPressurePhaseIdx_());
1645 }
1646
1647 if (!this->rockCompTransMult_.empty())
1648 return this->rockCompTransMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
1649
1650 // water compaction
1651 assert(!this->rockCompTransMultWc_.empty());
1652 LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
1653 LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx);
1654
1655 return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
1656 }
1657
1658 typename Vanguard::TransmissibilityType transmissibilities_;
1659
1660 std::shared_ptr<EclMaterialLawManager> materialLawManager_;
1661 std::shared_ptr<EclThermalLawManager> thermalLawManager_;
1662
1663 GlobalEqVector drift_;
1664
1665 WellModel wellModel_;
1666 AquiferModel aquiferModel_;
1667
1669 TracerModel tracerModel_;
1670
1671 template<class T>
1672 struct BCData
1673 {
1674 std::array<std::vector<T>,6> data;
1675
1676 void resize(std::size_t size, T defVal)
1677 {
1678 for (auto& d : data)
1679 d.resize(size, defVal);
1680 }
1681
1682 const std::vector<T>& operator()(FaceDir::DirEnum dir) const
1683 {
1684 if (dir == FaceDir::DirEnum::Unknown)
1685 throw std::runtime_error("Tried to access BC data for the 'Unknown' direction");
1686 int idx = 0;
1687 int div = static_cast<int>(dir);
1688 while ((div /= 2) >= 1)
1689 ++idx;
1690 assert(idx >= 0 && idx <= 5);
1691 return data[idx];
1692 }
1693
1694 std::vector<T>& operator()(FaceDir::DirEnum dir)
1695 {
1696 return const_cast<std::vector<T>&>(std::as_const(*this)(dir));
1697 }
1698 };
1699
1700 virtual void handleSolventBC(const BCProp::BCFace&, RateVector&) const = 0;
1701
1702 virtual void handlePolymerBC(const BCProp::BCFace&, RateVector&) const = 0;
1703
1704 BCData<int> bcindex_;
1705 bool nonTrivialBoundaryConditions_ = false;
1706 bool first_step_ = true;
1707};
1708
1709} // namespace Opm
1710
1711#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:61
Scalar maxPolymerAdsorption(unsigned elemIdx) const
Returns the max polymer adsorption value.
Definition FlowGenericProblem_impl.hpp:731
unsigned pvtRegionIndex(unsigned elemIdx) const
Returns the index the relevant PVT region given a cell index.
Definition FlowGenericProblem_impl.hpp:690
static std::string briefDescription()
Returns a human readable description of the problem for the help message.
Definition FlowGenericProblem_impl.hpp:126
Scalar porosity(unsigned globalSpaceIdx, unsigned timeIdx) const
Direct indexed access to the porosity.
Definition FlowGenericProblem_impl.hpp:327
Scalar rockCompressibility(unsigned globalSpaceIdx) const
Direct access to rock compressibility.
Definition FlowGenericProblem_impl.hpp:312
unsigned miscnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant MISC region given a cell index.
Definition FlowGenericProblem_impl.hpp:710
unsigned satnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant saturation function region given a cell index.
Definition FlowGenericProblem_impl.hpp:700
unsigned plmixnumRegionIndex(unsigned elemIdx) const
Returns the index the relevant PLMIXNUM (for polymer module) region given a cell index.
Definition FlowGenericProblem_impl.hpp:720
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:111
bool shouldWriteRestartFile() const
Returns true if an eWoms restart file should be written to disk.
Definition FlowGenericProblem.hpp:291
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition FlowProblem.hpp:94
const WellModel & wellModel() const
Returns a reference to the ECL well manager used by the problem.
Definition FlowProblem.hpp:1004
Scalar transmissibility(unsigned globalCenterElemIdx, unsigned globalElemIdx) const
Direct access to the transmissibility between two elements.
Definition FlowProblem.hpp:520
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
virtual void writeOutput(bool verbose)
Write the requested quantities of the current solution into the output files.
Definition FlowProblem.hpp:473
const DimMatrix & intrinsicPermeability(unsigned globalElemIdx) const
This method returns the intrinsic permeability tensor given a global element index.
Definition FlowProblem.hpp:502
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
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition FlowProblem.hpp:378
unsigned satnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:826
Scalar thermalHalfTransmissibilityOut(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:610
Scalar transmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:567
Scalar thermalHalfTransmissibility(const unsigned globalSpaceIdxIn, const unsigned globalSpaceIdxOut) const
Definition FlowProblem.hpp:587
Scalar rockCompressibility(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:686
Scalar maxOilSaturation(unsigned globalDofIdx) const
Returns an element's historic maximum oil phase saturation that was observed during the simulation.
Definition FlowProblem.hpp:902
std::string name() const
The problem name.
Definition FlowProblem.hpp:857
LhsEval rockCompTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1096
void endIteration()
Called by the simulator after each Newton-Raphson iteration.
Definition FlowProblem.hpp:388
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:633
Scalar nextTimeStepSize() const
Propose the size of the next time step to the simulator.
Definition FlowProblem.hpp:1025
LhsEval computeRockCompTransMultiplier_(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the transmissibility multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1621
const ThermalConductionLawParams & thermalConductionLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:762
Scalar transmissibility(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:509
virtual void beginEpisode()
Called by the simulator before an episode begins.
Definition FlowProblem.hpp:288
Scalar dofCenterDepth(unsigned globalSpaceIdx) const
Direct indexed acces to the depth of an degree of freedom [m].
Definition FlowProblem.hpp:677
LhsEval wellTransMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Return the well transmissibility multiplier due to rock changues.
Definition FlowProblem.hpp:1108
std::shared_ptr< const EclMaterialLawManager > materialLawManager() const
Returns the ECL material law manager.
Definition FlowProblem.hpp:774
unsigned plmixnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:842
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:864
Scalar thermalHalfTransmissibilityBoundary(const Context &elemCtx, unsigned boundaryFaceIdx) const
Definition FlowProblem.hpp:623
unsigned miscnumRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition FlowProblem.hpp:834
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:488
Scalar maxPolymerAdsorption(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the max polymer adsorption value.
Definition FlowProblem.hpp:851
void beginTimeStep()
Called by the simulator before each time integration.
Definition FlowProblem.hpp:345
Scalar thermalHalfTransmissibilityIn(const Context &context, unsigned faceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:597
void setMaxOilSaturation(unsigned globalDofIdx, Scalar value)
Sets an element's maximum oil phase saturation observed during the simulation.
Definition FlowProblem.hpp:919
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:749
Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the dispersivity for a face i.e.
Definition FlowProblem.hpp:547
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition FlowProblem.hpp:180
virtual void endEpisode()
Called by the simulator after the end of an episode.
Definition FlowProblem.hpp:450
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:963
virtual void endTimeStep()
Called by the simulator after each time integration.
Definition FlowProblem.hpp:398
virtual void initialSolutionApplied()
Callback used by the model to indicate that the initial solution has been determined for all degrees ...
Definition FlowProblem.hpp:930
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:727
Scalar thermalTransmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:554
Scalar diffusivity(const Context &context, unsigned fromDofLocalIdx, unsigned toDofLocalIdx) const
Definition FlowProblem.hpp:529
Scalar rockReferencePressure(unsigned globalSpaceIdx) const
Definition FlowProblem.hpp:705
void deserialize(Restarter &res)
This method restores the complete state of the problem and its sub-objects from disk.
Definition FlowProblem.hpp:254
std::shared_ptr< EclMaterialLawManager > materialLawManager()
Definition FlowProblem.hpp:810
Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the depth of an degree of freedom [m].
Definition FlowProblem.hpp:665
LhsEval rockCompPoroMultiplier(const IntensiveQuantities &intQuants, unsigned elementIdx) const
Calculate the porosity multiplier due to water induced rock compaction.
Definition FlowProblem.hpp:1052
Scalar transmissibilityBoundary(const unsigned globalSpaceIdx, const unsigned boundaryFaceIdx) const
Direct access to a boundary transmissibility.
Definition FlowProblem.hpp:577
void serialize(Restarter &res)
This method writes the complete state of the problem and its subobjects to disk.
Definition FlowProblem.hpp:273
Scalar rockReferencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition FlowProblem.hpp:696
Scalar diffusivity(const unsigned globalCellIn, const unsigned globalCellOut) const
give the transmissibility for a face i.e.
Definition FlowProblem.hpp:540
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:51
void diagnosis(const EclipseState &eclState, const LevelCartesianIndexMapper &levelCartesianIndexMapper)
This function is used to diagnosis relperm in eclipse data file.
Definition RelpermDiagnostics.cpp:829
A class which handles tracers as specified in by ECL.
Definition TracerModel.hpp:72
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:1673
Definition FlowProblem.hpp:1462