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,
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);
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);
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);
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);
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_;
277 FluidSystem,
typename Region::CalcDissolution
281 FluidSystem,
typename Region::CalcEvaporation,
typename Region::CalcWaterEvaporation
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;
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,
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>
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>;
701 const Vec&
rs()
const {
return rs_; }
702 const Vec&
rv()
const {
return rv_; }
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_;
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:270
Definition: InitStateEquil.hpp:680
std::vector< Vec > PVec
Definition: InitStateEquil.hpp:694
const PVec & press() const
Definition: InitStateEquil.hpp:699
const Vec & rvw() const
Definition: InitStateEquil.hpp:703
const Vec & rv() const
Definition: InitStateEquil.hpp:702
const Vec & saltSaturation() const
Definition: InitStateEquil.hpp:698
const Vec & saltConcentration() const
Definition: InitStateEquil.hpp:697
const Vec & temperature() const
Definition: InitStateEquil.hpp:696
const Vec & rs() const
Definition: InitStateEquil.hpp:701
std::vector< double > Vec
Definition: InitStateEquil.hpp:693
const PVec & saturation() const
Definition: InitStateEquil.hpp:700
InitialStateComputer(MaterialLawManager &materialLawManager, const EclipseState &eclipseState, const Grid &grid, const GridView &gridView, const CartesianIndexMapper &cartMapper, const double grav, const int num_pressure_points=2000, const bool applySwatInit=true)
Definition: InitStateEquil_impl.hpp:1308
Definition: InitStateEquil.hpp:132
Gas(const TabulatedFunction &tempVdTable, const RV &rv, const RVW &rvw, const int pvtRegionIdx, const double normGrav)
Definition: InitStateEquil_impl.hpp:329
double operator()(const double depth, const double press) const
Definition: InitStateEquil_impl.hpp:344
Definition: InitStateEquil.hpp:109
double operator()(const double depth, const double press) const
Definition: InitStateEquil_impl.hpp:296
Oil(const TabulatedFunction &tempVdTable, const RS &rs, const int pvtRegionIdx, const double normGrav)
Definition: InitStateEquil_impl.hpp:283
Definition: InitStateEquil.hpp:86
Water(const TabulatedFunction &tempVdTable, const TabulatedFunction &saltVdTable, const int pvtRegionIdx, const double normGrav)
Definition: InitStateEquil_impl.hpp:249
double operator()(const double depth, const double press) const
Definition: InitStateEquil_impl.hpp:262
Definition: InitStateEquil.hpp:376
const PhaseQuantityValue & correctedPhasePressures() const
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)
Definition: InitStateEquil_impl.hpp:574
PhaseSaturations & operator=(PhaseSaturations &&)=delete
Disabled move-assignment operator.
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< double > &swatInit)
Definition: InitStateEquil_impl.hpp:550
Definition: InitStateEquil.hpp:159
PressureTable & operator=(const PressureTable &rhs)
Definition: InitStateEquil_impl.hpp:990
double gas(const double depth) const
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
PressureTable(const double gravity, const int samplePoints=2000)
Definition: InitStateEquil_impl.hpp:960
void equilibrate(const Region ®, const VSpan &span)
Definition: InitStateEquil_impl.hpp:1016
double water(const double depth) const
Definition: InitStateEquil_impl.hpp:1067
std::array< double, 2 > VSpan
Definition: InitStateEquil.hpp:161
double oil(const double depth) const
Definition: InitStateEquil_impl.hpp:1048
Definition: InitStateEquil.hpp:64
double operator()(const double x) const
Definition: InitStateEquil_impl.hpp:215
RK4IVP(const RHS &f, const std::array< double, 2 > &span, const double y0, const int N)
Definition: InitStateEquil_impl.hpp:180
Definition: EquilibrationHelpers.hpp:600
std::pair< double, double > cellZMinMax(const Element &element)
Definition: InitStateEquil_impl.hpp:161
void verticalExtent(const CellRange &cells, const std::vector< std::pair< double, double > > &cellZMinMax, const Comm &comm, std::array< double, 2 > &span)
Definition: InitStateEquil_impl.hpp:64
Definition: BlackoilPhases.hpp:27
Simple set of per-phase (named by primary component) quantities.
Definition: InitStateEquil.hpp:328
double water
Definition: InitStateEquil.hpp:331
PhaseQuantityValue & operator/=(const double x)
Definition: InitStateEquil.hpp:342
PhaseQuantityValue & axpy(const PhaseQuantityValue &rhs, const double a)
Definition: InitStateEquil.hpp:333
double gas
Definition: InitStateEquil.hpp:330
double oil
Definition: InitStateEquil.hpp:329
void reset()
Definition: InitStateEquil.hpp:351
Definition: InitStateEquil.hpp:381
CellID cell
Definition: InitStateEquil.hpp:382
double depth
Definition: InitStateEquil.hpp:383
Definition: InitStateEquil.hpp:246
double pressure
Definition: InitStateEquil.hpp:248
double depth
Definition: InitStateEquil.hpp:247