29#ifndef OPM_INIT_STATE_EQUIL_HPP
30#define OPM_INIT_STATE_EQUIL_HPP
32#include <opm/models/utils/propertysystem.hh>
34#include <opm/material/common/Tabulated1DFunction.hpp>
35#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
48class NumericalAquifers;
60namespace Miscibility {
class RsFunction; }
67 const std::array<double,2>& span,
72 operator()(
const double x)
const;
76 std::array<double,2> span_;
77 std::vector<double> y_;
78 std::vector<double> f_;
80 double stepsize()
const;
83namespace PhasePressODE {
84template <
class Flu
idSystem>
87using TabulatedFunction = Tabulated1DFunction<double>;
89 Water(
const TabulatedFunction& tempVdTable,
90 const TabulatedFunction& saltVdTable,
91 const int pvtRegionIdx,
92 const double normGrav);
94 double operator()(
const double depth,
95 const double press)
const;
98 const TabulatedFunction& tempVdTable_;
99 const TabulatedFunction& saltVdTable_;
100 const int pvtRegionIdx_;
103 double density(
const double depth,
104 const double press)
const;
107template <
class Flu
idSystem,
class RS>
110using TabulatedFunction = Tabulated1DFunction<double>;
112 Oil(
const TabulatedFunction& tempVdTable,
114 const int pvtRegionIdx,
115 const double normGrav);
117 double operator()(
const double depth,
118 const double press)
const;
121 const TabulatedFunction& tempVdTable_;
123 const int pvtRegionIdx_;
126 double density(
const double depth,
127 const double press)
const;
130template <
class Flu
idSystem,
class RV,
class RVW>
133using TabulatedFunction = Tabulated1DFunction<double>;
135 Gas(
const TabulatedFunction& tempVdTable,
138 const int pvtRegionIdx,
139 const double normGrav);
141 double operator()(
const double depth,
142 const double press)
const;
145 const TabulatedFunction& tempVdTable_;
148 const int pvtRegionIdx_;
151 double density(
const double depth,
152 const double press)
const;
157template <
class Flu
idSystem,
class Region>
161 using VSpan = std::array<double, 2>;
172 const int samplePoints = 2000);
202 void equilibrate(
const Region& reg,
221 double oil(
const double depth)
const;
230 double gas(
const double depth)
const;
239 double water(
const double depth)
const;
243 class PressureFunction
251 explicit PressureFunction(
const ODE& ode,
256 PressureFunction(
const PressureFunction& rhs);
258 PressureFunction(PressureFunction&& rhs) =
default;
260 PressureFunction& operator=(
const PressureFunction& rhs);
262 PressureFunction& operator=(PressureFunction&& rhs);
264 double value(
const double depth)
const;
267 enum Direction : std::size_t { Up, Down, NumDir };
270 using DistrPtr = std::unique_ptr<Distribution>;
273 std::array<DistrPtr, Direction::NumDir> value_;
276 using OilPressODE = PhasePressODE::Oil<
277 FluidSystem,
typename Region::CalcDissolution
280 using GasPressODE = PhasePressODE::Gas<
281 FluidSystem,
typename Region::CalcEvaporation,
typename Region::CalcWaterEvaporation
284 using WatPressODE = PhasePressODE::Water<FluidSystem>;
286 using OPress = PressureFunction<OilPressODE>;
287 using GPress = PressureFunction<GasPressODE>;
288 using WPress = PressureFunction<WatPressODE>;
291 (const Region&, const VSpan&);
296 std::unique_ptr<OPress> oil_{};
297 std::unique_ptr<GPress> gas_{};
298 std::unique_ptr<WPress> wat_{};
300 template <
typename PressFunc>
301 void checkPtr(
const PressFunc* phasePress,
302 const std::string& phaseName)
const;
304 Strategy selectEquilibrationStrategy(
const Region& reg)
const;
306 void copyInPointers(
const PressureTable& rhs);
308 void equil_WOG(
const Region& reg,
const VSpan& span);
309 void equil_GOW(
const Region& reg,
const VSpan& span);
310 void equil_OWG(
const Region& reg,
const VSpan& span);
312 void makeOilPressure(
const typename OPress::InitCond& ic,
316 void makeGasPressure(
const typename GPress::InitCond& ic,
320 void makeWatPressure(
const typename WPress::InitCond& ic,
335 this->oil += a * rhs.oil;
336 this->gas += a * rhs.gas;
337 this->water += a * rhs.water;
353 this->oil = this->gas = this->water = 0.0;
374template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
397 const std::vector<double>& swatInit);
439 struct EvaluationPoint {
440 const Position* position{
nullptr};
441 const Region* region {
nullptr};
442 const PTable* ptable {
nullptr};
448 using FluidState = ::Opm::
449 SimpleModularFluidState<double, 3, 3,
461 using MaterialLaw =
typename MaterialLawManager::MaterialLaw;
464 using PhaseIdx = std::remove_cv_t<
465 std::remove_reference_t<
decltype(FluidSystem::oilPhaseIdx)>
469 MaterialLawManager& matLawMgr_;
472 const std::vector<double>& swatInit_;
475 PhaseQuantityValue sat_;
478 PhaseQuantityValue press_;
481 EvaluationPoint evalPt_;
484 FluidState fluidState_;
487 std::array<double, FluidSystem::numPhases> matLawCapPress_;
498 void setEvaluationPoint(
const Position& x,
505 void initializePhaseQuantities();
524 void deriveWaterSat();
528 void fixUnphysicalTransition();
532 void accountForScaledSaturations();
548 std::pair<double, bool> applySwatInit(
const double pcow);
562 std::pair<double, bool> applySwatInit(
const double pc,
const double sw);
566 void computeMaterialLawCapPress();
570 double materialLawCapPressGasOil()
const;
574 double materialLawCapPressOilWater()
const;
578 double materialLawCapPressGasWater()
const;
587 bool isConstCapPress(
const PhaseIdx phaseIdx)
const;
594 bool isOverlappingTransition()
const;
615 double fromDepthTable(
const double contactdepth,
616 const PhaseIdx phasePos,
617 const bool isincr)
const;
636 double invertCapPress(
const double pc,
637 const PhaseIdx phasePos,
638 const bool isincr)
const;
641 PhaseIdx oilPos()
const
643 return FluidSystem::oilPhaseIdx;
647 PhaseIdx gasPos()
const
649 return FluidSystem::gasPhaseIdx;
653 PhaseIdx waterPos()
const
655 return FluidSystem::waterPhaseIdx;
661template <
typename CellRange,
typename Comm>
662void verticalExtent(
const CellRange& cells,
663 const std::vector<std::pair<double, double>>& cellZMinMax,
665 std::array<double,2>& span);
667template <
class Element>
668std::pair<double,double> cellZMinMax(
const Element& element);
672namespace DeckDependent {
674template<
class FluidSystem,
678 class CartesianIndexMapper>
681 using Element =
typename GridView::template Codim<0>::Entity;
683 template<
class MaterialLawManager>
685 const EclipseState& eclipseState,
687 const GridView& gridView,
688 const CartesianIndexMapper& cartMapper,
690 const int num_pressure_points = 2000,
691 const bool applySwatInit =
true);
693 using Vec = std::vector<double>;
694 using PVec = std::vector<Vec>;
696 const Vec& temperature()
const {
return temperature_; }
697 const Vec& saltConcentration()
const {
return saltConcentration_; }
698 const Vec& saltSaturation()
const {
return saltSaturation_; }
699 const PVec& press()
const {
return pp_; }
700 const PVec& saturation()
const {
return sat_; }
701 const Vec& rs()
const {
return rs_; }
702 const Vec& rv()
const {
return rv_; }
703 const Vec& rvw()
const {
return rvw_; }
706 template <
class RMap>
707 void updateInitialTemperature_(
const EclipseState& eclState,
const RMap& reg);
709 template <
class RMap>
710 void updateInitialSaltConcentration_(
const EclipseState& eclState,
const RMap& reg);
712 template <
class RMap>
713 void updateInitialSaltSaturation_(
const EclipseState& eclState,
const RMap& reg);
715 void updateCellProps_(
const GridView& gridView,
716 const NumericalAquifers& aquifer);
718 void applyNumericalAquifers_(
const GridView& gridView,
719 const NumericalAquifers& aquifer,
720 const bool co2store_or_h2store);
723 void setRegionPvtIdx(
const EclipseState& eclState,
const RMap& reg);
725 template <
class RMap,
class MaterialLawManager,
class Comm>
726 void calcPressSatRsRv(
const RMap& reg,
727 const std::vector<EquilRecord>& rec,
728 MaterialLawManager& materialLawManager,
732 template <
class CellRange,
class EquilibrationMethod>
733 void cellLoop(
const CellRange& cells,
734 EquilibrationMethod&& eqmethod);
736 template <
class CellRange,
class PressTable,
class PhaseSat>
737 void equilibrateCellCentres(
const CellRange& cells,
739 const PressTable& ptable,
742 template <
class CellRange,
class PressTable,
class PhaseSat>
743 void equilibrateHorizontal(
const CellRange& cells,
746 const PressTable& ptable,
749 std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_;
750 std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_;
751 std::vector< std::shared_ptr<Miscibility::RsFunction> > rvwFunc_;
752 using TabulatedFunction = Tabulated1DFunction<double>;
753 std::vector<TabulatedFunction> tempVdTable_;
754 std::vector<TabulatedFunction> saltVdTable_;
755 std::vector<TabulatedFunction> saltpVdTable_;
756 std::vector<int> regionPvtIdx_;
758 Vec saltConcentration_;
765 const CartesianIndexMapper& cartesianIndexMapper_;
767 Vec cellCenterDepth_;
768 std::vector<std::pair<double,double>> cellZSpan_;
769 std::vector<std::pair<double,double>> cellZMinMax_;
770 int num_pressure_points_;
Definition InitStateEquil.hpp:680
Calculator for phase saturations.
Definition InitStateEquil.hpp:376
const PhaseQuantityValue & correctedPhasePressures() const
Retrieve saturation-corrected phase pressures.
Definition InitStateEquil.hpp:431
PressureTable< FluidSystem, Region > PTable
Convenience type alias.
Definition InitStateEquil.hpp:387
PhaseSaturations & operator=(const PhaseSaturations &)=delete
Disabled assignment operator.
const PhaseQuantityValue & deriveSaturations(const Position &x, const Region ®, const PTable &ptable)
Calculate phase saturations at particular point of the simulation model geometry.
Definition InitStateEquil_impl.hpp:574
PhaseSaturations & operator=(PhaseSaturations &&)=delete
Disabled move-assignment operator.
Definition InitStateEquil.hpp:159
PressureTable & operator=(const PressureTable &rhs)
Assignment operator.
Definition InitStateEquil_impl.hpp:990
double gas(const double depth) const
Evaluate gas phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1057
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition InitStateEquil_impl.hpp:1041
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition InitStateEquil_impl.hpp:1034
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition InitStateEquil_impl.hpp:1027
double water(const double depth) const
Evaluate water phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1067
double oil(const double depth) const
Evaluate oil phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1048
Definition InitStateEquil.hpp:64
Aggregate information base of an equilibration region.
Definition EquilibrationHelpers.hpp:600
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Simple set of per-phase (named by primary component) quantities.
Definition InitStateEquil.hpp:328
Evaluation point within a model geometry.
Definition InitStateEquil.hpp:381
Definition InitStateEquil.hpp:246