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>
50class NumericalAquifers;
62namespace Miscibility {
class RsFunction; }
69 const std::array<double,2>& span,
78 std::array<double,2> span_;
79 std::vector<double> y_;
80 std::vector<double> f_;
82 double stepsize()
const;
85namespace PhasePressODE {
86template <
class Flu
idSystem>
89using TabulatedFunction = Tabulated1DFunction<double>;
91 Water(
const TabulatedFunction& tempVdTable,
92 const TabulatedFunction& saltVdTable,
93 const int pvtRegionIdx,
94 const double normGrav);
97 const double press)
const;
100 const TabulatedFunction& tempVdTable_;
101 const TabulatedFunction& saltVdTable_;
102 const int pvtRegionIdx_;
105 double density(
const double depth,
106 const double press)
const;
109template <
class Flu
idSystem,
class RS>
112using TabulatedFunction = Tabulated1DFunction<double>;
114 Oil(
const TabulatedFunction& tempVdTable,
116 const int pvtRegionIdx,
117 const double normGrav);
120 const double press)
const;
123 const TabulatedFunction& tempVdTable_;
125 const int pvtRegionIdx_;
128 double density(
const double depth,
129 const double press)
const;
132template <
class Flu
idSystem,
class RV,
class RVW>
135using TabulatedFunction = Tabulated1DFunction<double>;
137 Gas(
const TabulatedFunction& tempVdTable,
140 const int pvtRegionIdx,
141 const double normGrav);
144 const double press)
const;
147 const TabulatedFunction& tempVdTable_;
150 const int pvtRegionIdx_;
153 double density(
const double depth,
154 const double press)
const;
159template <
class Flu
idSystem,
class Region>
163 using VSpan = std::array<double, 2>;
174 const int samplePoints = 2000);
223 double oil(
const double depth)
const;
232 double gas(
const double depth)
const;
241 double water(
const double depth)
const;
245 class PressureFunction
253 explicit PressureFunction(
const ODE& ode,
258 PressureFunction(
const PressureFunction& rhs);
260 PressureFunction(PressureFunction&& rhs) =
default;
262 PressureFunction&
operator=(
const PressureFunction& rhs);
264 PressureFunction&
operator=(PressureFunction&& rhs);
266 double value(
const double depth)
const;
269 enum Direction : std::size_t { Up, Down, NumDir };
272 using DistrPtr = std::unique_ptr<Distribution>;
275 std::array<DistrPtr, Direction::NumDir> value_;
279 FluidSystem,
typename Region::CalcDissolution
283 FluidSystem,
typename Region::CalcEvaporation,
typename Region::CalcWaterEvaporation
288 using OPress = PressureFunction<OilPressODE>;
289 using GPress = PressureFunction<GasPressODE>;
290 using WPress = PressureFunction<WatPressODE>;
293 (
const Region&,
const VSpan&);
298 std::unique_ptr<OPress> oil_{};
299 std::unique_ptr<GPress> gas_{};
300 std::unique_ptr<WPress> wat_{};
302 template <
typename PressFunc>
303 void checkPtr(
const PressFunc* phasePress,
304 const std::string& phaseName)
const;
306 Strategy selectEquilibrationStrategy(
const Region& reg)
const;
310 void equil_WOG(
const Region& reg,
const VSpan& span);
311 void equil_GOW(
const Region& reg,
const VSpan& span);
312 void equil_OWG(
const Region& reg,
const VSpan& span);
314 void makeOilPressure(
const typename OPress::InitCond& ic,
318 void makeGasPressure(
const typename GPress::InitCond& ic,
322 void makeWatPressure(
const typename WPress::InitCond& ic,
376template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
399 const std::vector<double>& swatInit);
441 struct EvaluationPoint {
442 const Position* position{
nullptr};
443 const Region* region {
nullptr};
444 const PTable* ptable {
nullptr};
450 using FluidState = ::Opm::
451 SimpleModularFluidState<double, 3, 3,
463 using MaterialLaw =
typename MaterialLawManager::MaterialLaw;
466 using PhaseIdx = std::remove_cv_t<
467 std::remove_reference_t<
decltype(FluidSystem::oilPhaseIdx)>
471 MaterialLawManager& matLawMgr_;
474 const std::vector<double>& swatInit_;
477 PhaseQuantityValue sat_;
480 PhaseQuantityValue press_;
483 EvaluationPoint evalPt_;
486 FluidState fluidState_;
489 std::array<double, FluidSystem::numPhases> matLawCapPress_;
500 void setEvaluationPoint(
const Position& x,
507 void initializePhaseQuantities();
526 void deriveWaterSat();
530 void fixUnphysicalTransition();
534 void accountForScaledSaturations();
550 std::pair<double, bool> applySwatInit(
const double pcow);
564 std::pair<double, bool> applySwatInit(
const double pc,
const double sw);
568 void computeMaterialLawCapPress();
572 double materialLawCapPressGasOil()
const;
576 double materialLawCapPressOilWater()
const;
580 double materialLawCapPressGasWater()
const;
589 bool isConstCapPress(
const PhaseIdx phaseIdx)
const;
596 bool isOverlappingTransition()
const;
617 double fromDepthTable(
const double contactdepth,
618 const PhaseIdx phasePos,
619 const bool isincr)
const;
638 double invertCapPress(
const double pc,
639 const PhaseIdx phasePos,
640 const bool isincr)
const;
643 PhaseIdx oilPos()
const
645 return FluidSystem::oilPhaseIdx;
649 PhaseIdx gasPos()
const
651 return FluidSystem::gasPhaseIdx;
655 PhaseIdx waterPos()
const
657 return FluidSystem::waterPhaseIdx;
663template <
typename CellRange>
665 const std::vector<std::pair<double, double>>&
cellZMinMax,
667 std::array<double,2>& span);
669template <
class Element>
670std::pair<double,double>
cellZMinMax(
const Element& element);
674namespace DeckDependent {
676template<
class FluidSystem,
680 class CartesianIndexMapper>
683 using Element =
typename GridView::template Codim<0>::Entity;
685 template<
class MaterialLawManager>
687 const EclipseState& eclipseState,
689 const GridView& gridView,
690 const CartesianIndexMapper& cartMapper,
692 const int num_pressure_points = 2000,
693 const bool applySwatInit =
true);
695 using Vec = std::vector<double>;
703 const Vec&
rs()
const {
return rs_; }
704 const Vec&
rv()
const {
return rv_; }
708 template <
class RMap>
709 void updateInitialTemperature_(
const EclipseState& eclState,
const RMap& reg);
711 template <
class RMap>
712 void updateInitialSaltConcentration_(
const EclipseState& eclState,
const RMap& reg);
714 template <
class RMap>
715 void updateInitialSaltSaturation_(
const EclipseState& eclState,
const RMap& reg);
717 void updateCellProps_(
const GridView& gridView,
718 const NumericalAquifers& aquifer);
720 void applyNumericalAquifers_(
const GridView& gridView,
721 const NumericalAquifers& aquifer,
722 const bool co2store_or_h2store);
725 void setRegionPvtIdx(
const EclipseState& eclState,
const RMap& reg);
727 template <
class RMap,
class MaterialLawManager,
class Comm>
728 void calcPressSatRsRv(
const RMap& reg,
729 const std::vector<EquilRecord>& rec,
730 MaterialLawManager& materialLawManager,
734 template <
class CellRange,
class EquilibrationMethod>
735 void cellLoop(
const CellRange& cells,
736 EquilibrationMethod&& eqmethod);
738 template <
class CellRange,
class PressTable,
class PhaseSat>
739 void equilibrateCellCentres(
const CellRange& cells,
741 const PressTable& ptable,
744 template <
class CellRange,
class PressTable,
class PhaseSat>
745 void equilibrateHorizontal(
const CellRange& cells,
748 const PressTable& ptable,
751 std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_;
752 std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_;
753 std::vector< std::shared_ptr<Miscibility::RsFunction> > rvwFunc_;
754 using TabulatedFunction = Tabulated1DFunction<double>;
755 std::vector<TabulatedFunction> tempVdTable_;
756 std::vector<TabulatedFunction> saltVdTable_;
757 std::vector<TabulatedFunction> saltpVdTable_;
758 std::vector<int> regionPvtIdx_;
760 Vec saltConcentration_;
767 const CartesianIndexMapper& cartesianIndexMapper_;
769 Vec cellCenterDepth_;
770 std::vector<std::pair<double,double>> cellZSpan_;
771 std::vector<std::pair<double,double>> cellZMinMax_;
772 int num_pressure_points_;
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:274
Definition: InitStateEquil.hpp:682
std::vector< Vec > PVec
Definition: InitStateEquil.hpp:696
const PVec & press() const
Definition: InitStateEquil.hpp:701
const Vec & rvw() const
Definition: InitStateEquil.hpp:705
const Vec & rv() const
Definition: InitStateEquil.hpp:704
const Vec & saltSaturation() const
Definition: InitStateEquil.hpp:700
const Vec & saltConcentration() const
Definition: InitStateEquil.hpp:699
const Vec & temperature() const
Definition: InitStateEquil.hpp:698
const Vec & rs() const
Definition: InitStateEquil.hpp:703
std::vector< double > Vec
Definition: InitStateEquil.hpp:695
const PVec & saturation() const
Definition: InitStateEquil.hpp:702
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:134
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:111
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:88
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:378
const PhaseQuantityValue & correctedPhasePressures() const
Definition: InitStateEquil.hpp:433
PressureTable< FluidSystem, Region > PTable
Convenience type alias.
Definition: InitStateEquil.hpp:389
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:161
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:163
double oil(const double depth) const
Definition: InitStateEquil_impl.hpp:1048
Definition: InitStateEquil.hpp:66
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 Parallel::Communication &comm, std::array< double, 2 > &span)
Definition: InitStateEquil_impl.hpp:64
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: BlackoilPhases.hpp:27
Simple set of per-phase (named by primary component) quantities.
Definition: InitStateEquil.hpp:330
double water
Definition: InitStateEquil.hpp:333
PhaseQuantityValue & operator/=(const double x)
Definition: InitStateEquil.hpp:344
PhaseQuantityValue & axpy(const PhaseQuantityValue &rhs, const double a)
Definition: InitStateEquil.hpp:335
double gas
Definition: InitStateEquil.hpp:332
double oil
Definition: InitStateEquil.hpp:331
void reset()
Definition: InitStateEquil.hpp:353
Definition: InitStateEquil.hpp:383
CellID cell
Definition: InitStateEquil.hpp:384
double depth
Definition: InitStateEquil.hpp:385
Definition: InitStateEquil.hpp:248
double pressure
Definition: InitStateEquil.hpp:250
double depth
Definition: InitStateEquil.hpp:249