29#ifndef OPM_INIT_STATE_EQUIL_HPP
30#define OPM_INIT_STATE_EQUIL_HPP
34#include <opm/material/common/Tabulated1DFunction.hpp>
35#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
50class NumericalAquifers;
61template<
class Scalar>
class EquilReg;
62namespace Miscibility {
template<
class Scalar>
class RsFunction; }
65template <
class Scalar,
class RHS>
70 const std::array<Scalar,2>& span,
78 std::array<Scalar,2> span_;
79 std::vector<Scalar> y_;
80 std::vector<Scalar> f_;
82 Scalar stepsize()
const;
85namespace PhasePressODE {
86template <
class Flu
idSystem>
89 using Scalar =
typename FluidSystem::Scalar;
90 using TabulatedFunction = Tabulated1DFunction<Scalar>;
93 Water(
const TabulatedFunction& tempVdTable,
94 const TabulatedFunction& saltVdTable,
95 const int pvtRegionIdx,
96 const Scalar normGrav);
99 const Scalar press)
const;
102 const TabulatedFunction& tempVdTable_;
103 const TabulatedFunction& saltVdTable_;
104 const int pvtRegionIdx_;
107 Scalar density(
const Scalar depth,
108 const Scalar press)
const;
111template <
class Flu
idSystem,
class RS>
114 using Scalar =
typename FluidSystem::Scalar;
115 using TabulatedFunction = Tabulated1DFunction<Scalar>;
118 Oil(
const TabulatedFunction& tempVdTable,
120 const int pvtRegionIdx,
121 const Scalar normGrav);
124 const Scalar press)
const;
127 const TabulatedFunction& tempVdTable_;
129 const int pvtRegionIdx_;
132 Scalar density(
const Scalar depth,
133 const Scalar press)
const;
136template <
class Flu
idSystem,
class RV,
class RVW>
139 using Scalar =
typename FluidSystem::Scalar;
140 using TabulatedFunction = Tabulated1DFunction<Scalar>;
143 Gas(
const TabulatedFunction& tempVdTable,
146 const int pvtRegionIdx,
147 const Scalar normGrav);
150 const Scalar press)
const;
153 const TabulatedFunction& tempVdTable_;
156 const int pvtRegionIdx_;
159 Scalar density(
const Scalar depth,
160 const Scalar press)
const;
165template <
class Flu
idSystem,
class Region>
169 using Scalar =
typename FluidSystem::Scalar;
170 using VSpan = std::array<Scalar, 2>;
181 const int samplePoints = 2000);
252 class PressureFunction
260 explicit PressureFunction(
const ODE& ode,
265 PressureFunction(
const PressureFunction& rhs);
267 PressureFunction(PressureFunction&& rhs) =
default;
269 PressureFunction&
operator=(
const PressureFunction& rhs);
271 PressureFunction&
operator=(PressureFunction&& rhs);
276 enum Direction : std::size_t { Up, Down, NumDir };
279 using DistrPtr = std::unique_ptr<Distribution>;
282 std::array<DistrPtr, Direction::NumDir> value_;
286 FluidSystem,
typename Region::CalcDissolution
290 FluidSystem,
typename Region::CalcEvaporation,
typename Region::CalcWaterEvaporation
295 using OPress = PressureFunction<OilPressODE>;
296 using GPress = PressureFunction<GasPressODE>;
297 using WPress = PressureFunction<WatPressODE>;
300 (
const Region&,
const VSpan&);
305 std::unique_ptr<OPress> oil_{};
306 std::unique_ptr<GPress> gas_{};
307 std::unique_ptr<WPress> wat_{};
309 template <
typename PressFunc>
310 void checkPtr(
const PressFunc* phasePress,
311 const std::string& phaseName)
const;
313 Strategy selectEquilibrationStrategy(
const Region& reg)
const;
317 void equil_WOG(
const Region& reg,
const VSpan& span);
318 void equil_GOW(
const Region& reg,
const VSpan& span);
319 void equil_OWG(
const Region& reg,
const VSpan& span);
321 void makeOilPressure(
const typename OPress::InitCond& ic,
325 void makeGasPressure(
const typename GPress::InitCond& ic,
329 void makeWatPressure(
const typename WPress::InitCond& ic,
337template<
class Scalar>
384template <
class MaterialLawManager,
class Flu
idSystem,
class Region,
typename CellID>
388 using Scalar =
typename FluidSystem::Scalar;
408 const std::vector<Scalar>& swatInit);
450 struct EvaluationPoint {
451 const Position* position{
nullptr};
452 const Region* region {
nullptr};
453 const PTable* ptable {
nullptr};
459 using FluidState = ::Opm::
460 SimpleModularFluidState<
Scalar, 3, 3,
472 using MaterialLaw =
typename MaterialLawManager::MaterialLaw;
475 using PhaseIdx = std::remove_cv_t<
476 std::remove_reference_t<
decltype(FluidSystem::oilPhaseIdx)>
480 MaterialLawManager& matLawMgr_;
483 const std::vector<Scalar>& swatInit_;
486 PhaseQuantityValue<Scalar> sat_;
489 PhaseQuantityValue<Scalar> press_;
492 EvaluationPoint evalPt_;
495 FluidState fluidState_;
498 std::array<Scalar, FluidSystem::numPhases> matLawCapPress_;
509 void setEvaluationPoint(
const Position& x,
516 void initializePhaseQuantities();
535 void deriveWaterSat();
539 void fixUnphysicalTransition();
543 void accountForScaledSaturations();
559 std::pair<Scalar, bool> applySwatInit(
const Scalar pcow);
573 std::pair<Scalar, bool> applySwatInit(
const Scalar pc,
const Scalar sw);
577 void computeMaterialLawCapPress();
581 Scalar materialLawCapPressGasOil()
const;
585 Scalar materialLawCapPressOilWater()
const;
589 Scalar materialLawCapPressGasWater()
const;
598 bool isConstCapPress(
const PhaseIdx phaseIdx)
const;
605 bool isOverlappingTransition()
const;
627 const PhaseIdx phasePos,
628 const bool isincr)
const;
648 const PhaseIdx phasePos,
649 const bool isincr)
const;
652 PhaseIdx oilPos()
const
654 return FluidSystem::oilPhaseIdx;
658 PhaseIdx gasPos()
const
660 return FluidSystem::gasPhaseIdx;
664 PhaseIdx waterPos()
const
666 return FluidSystem::waterPhaseIdx;
672template <
typename CellRange,
class Scalar>
674 const std::vector<std::pair<Scalar, Scalar>>&
cellZMinMax,
676 std::array<Scalar,2>& span);
678template <
class Scalar,
class Element>
679std::pair<Scalar,Scalar>
cellZMinMax(
const Element& element);
683namespace DeckDependent {
685template<
class FluidSystem,
689 class CartesianIndexMapper>
692 using Element =
typename GridView::template Codim<0>::Entity;
693 using Scalar =
typename FluidSystem::Scalar;
695 template<
class MaterialLawManager>
697 const EclipseState& eclipseState,
699 const GridView& gridView,
700 const CartesianIndexMapper& cartMapper,
702 const int num_pressure_points = 2000,
703 const bool applySwatInit =
true);
705 using Vec = std::vector<Scalar>;
713 const Vec&
rs()
const {
return rs_; }
714 const Vec&
rv()
const {
return rv_; }
718 template <
class RMap>
719 void updateInitialTemperature_(
const EclipseState& eclState,
const RMap& reg);
721 template <
class RMap>
722 void updateInitialSaltConcentration_(
const EclipseState& eclState,
const RMap& reg);
724 template <
class RMap>
725 void updateInitialSaltSaturation_(
const EclipseState& eclState,
const RMap& reg);
727 void updateCellProps_(
const GridView& gridView,
728 const NumericalAquifers& aquifer);
730 void applyNumericalAquifers_(
const GridView& gridView,
731 const NumericalAquifers& aquifer,
732 const bool co2store_or_h2store);
735 void setRegionPvtIdx(
const EclipseState& eclState,
const RMap& reg);
737 template <
class RMap,
class MaterialLawManager,
class Comm>
738 void calcPressSatRsRv(
const RMap& reg,
739 const std::vector<EquilRecord>& rec,
740 MaterialLawManager& materialLawManager,
744 template <
class CellRange,
class EquilibrationMethod>
745 void cellLoop(
const CellRange& cells,
746 EquilibrationMethod&& eqmethod);
748 template <
class CellRange,
class PressTable,
class PhaseSat>
749 void equilibrateCellCentres(
const CellRange& cells,
751 const PressTable& ptable,
754 template <
class CellRange,
class PressTable,
class PhaseSat>
755 void equilibrateHorizontal(
const CellRange& cells,
758 const PressTable& ptable,
761 std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rsFunc_;
762 std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rvFunc_;
763 std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rvwFunc_;
764 using TabulatedFunction = Tabulated1DFunction<Scalar>;
765 std::vector<TabulatedFunction> tempVdTable_;
766 std::vector<TabulatedFunction> saltVdTable_;
767 std::vector<TabulatedFunction> saltpVdTable_;
768 std::vector<int> regionPvtIdx_;
770 Vec saltConcentration_;
777 const CartesianIndexMapper& cartesianIndexMapper_;
779 Vec cellCenterDepth_;
780 std::vector<std::pair<Scalar,Scalar>> cellZSpan_;
781 std::vector<std::pair<Scalar,Scalar>> cellZMinMax_;
782 int num_pressure_points_;
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:285
Definition: InitStateEquil.hpp:691
std::vector< Vec > PVec
Definition: InitStateEquil.hpp:706
const PVec & press() const
Definition: InitStateEquil.hpp:711
const Vec & rvw() const
Definition: InitStateEquil.hpp:715
InitialStateComputer(MaterialLawManager &materialLawManager, const EclipseState &eclipseState, const Grid &grid, const GridView &gridView, const CartesianIndexMapper &cartMapper, const Scalar grav, const int num_pressure_points=2000, const bool applySwatInit=true)
Definition: InitStateEquil_impl.hpp:1338
const Vec & rv() const
Definition: InitStateEquil.hpp:714
const Vec & saltSaturation() const
Definition: InitStateEquil.hpp:710
const Vec & saltConcentration() const
Definition: InitStateEquil.hpp:709
const Vec & temperature() const
Definition: InitStateEquil.hpp:708
std::vector< Scalar > Vec
Definition: InitStateEquil.hpp:705
const Vec & rs() const
Definition: InitStateEquil.hpp:713
const PVec & saturation() const
Definition: InitStateEquil.hpp:712
Definition: InitStateEquil.hpp:138
Gas(const TabulatedFunction &tempVdTable, const RV &rv, const RVW &rvw, const int pvtRegionIdx, const Scalar normGrav)
Definition: InitStateEquil_impl.hpp:338
Scalar operator()(const Scalar depth, const Scalar press) const
Definition: InitStateEquil_impl.hpp:354
Definition: InitStateEquil.hpp:113
Oil(const TabulatedFunction &tempVdTable, const RS &rs, const int pvtRegionIdx, const Scalar normGrav)
Definition: InitStateEquil_impl.hpp:290
Scalar operator()(const Scalar depth, const Scalar press) const
Definition: InitStateEquil_impl.hpp:304
Definition: InitStateEquil.hpp:88
Scalar operator()(const Scalar depth, const Scalar press) const
Definition: InitStateEquil_impl.hpp:264
Water(const TabulatedFunction &tempVdTable, const TabulatedFunction &saltVdTable, const int pvtRegionIdx, const Scalar normGrav)
Definition: InitStateEquil_impl.hpp:250
Definition: InitStateEquil.hpp:386
const PhaseQuantityValue< Scalar > & deriveSaturations(const Position &x, const Region ®, const PTable &ptable)
Definition: InitStateEquil_impl.hpp:596
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< Scalar > &swatInit)
Definition: InitStateEquil_impl.hpp:572
typename FluidSystem::Scalar Scalar
Definition: InitStateEquil.hpp:388
PressureTable< FluidSystem, Region > PTable
Convenience type alias.
Definition: InitStateEquil.hpp:398
PhaseSaturations & operator=(const PhaseSaturations &)=delete
Disabled assignment operator.
PhaseSaturations & operator=(PhaseSaturations &&)=delete
Disabled move-assignment operator.
const PhaseQuantityValue< Scalar > & correctedPhasePressures() const
Definition: InitStateEquil.hpp:442
Definition: InitStateEquil.hpp:167
PressureTable & operator=(const PressureTable &rhs)
Definition: InitStateEquil_impl.hpp:1017
Scalar water(const Scalar depth) const
Definition: InitStateEquil_impl.hpp:1097
Scalar gas(const Scalar depth) const
Definition: InitStateEquil_impl.hpp:1086
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition: InitStateEquil_impl.hpp:1068
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition: InitStateEquil_impl.hpp:1061
Scalar oil(const Scalar depth) const
Definition: InitStateEquil_impl.hpp:1076
std::array< Scalar, 2 > VSpan
Definition: InitStateEquil.hpp:170
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition: InitStateEquil_impl.hpp:1054
typename FluidSystem::Scalar Scalar
Definition: InitStateEquil.hpp:169
void equilibrate(const Region ®, const VSpan &span)
Definition: InitStateEquil_impl.hpp:1043
PressureTable(const Scalar gravity, const int samplePoints=2000)
Definition: InitStateEquil_impl.hpp:987
Definition: InitStateEquil.hpp:67
Scalar operator()(const Scalar x) const
Definition: InitStateEquil_impl.hpp:216
RK4IVP(const RHS &f, const std::array< Scalar, 2 > &span, const Scalar y0, const int N)
Definition: InitStateEquil_impl.hpp:181
Definition: EquilibrationHelpers.hpp:617
std::pair< Scalar, Scalar > cellZMinMax(const Element &element)
Definition: InitStateEquil_impl.hpp:162
void verticalExtent(const CellRange &cells, const std::vector< std::pair< Scalar, Scalar > > &cellZMinMax, const Parallel::Communication &comm, std::array< Scalar, 2 > &span)
Definition: InitStateEquil_impl.hpp:64
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilboundaryratevector.hh:37
The Opm property system, traits with inheritance.
Simple set of per-phase (named by primary component) quantities.
Definition: InitStateEquil.hpp:338
void reset()
Definition: InitStateEquil.hpp:361
Scalar gas
Definition: InitStateEquil.hpp:340
Scalar water
Definition: InitStateEquil.hpp:341
PhaseQuantityValue & operator/=(const Scalar x)
Definition: InitStateEquil.hpp:352
Scalar oil
Definition: InitStateEquil.hpp:339
PhaseQuantityValue & axpy(const PhaseQuantityValue &rhs, const Scalar a)
Definition: InitStateEquil.hpp:343
Definition: InitStateEquil.hpp:392
CellID cell
Definition: InitStateEquil.hpp:393
Scalar depth
Definition: InitStateEquil.hpp:394
Definition: InitStateEquil.hpp:255
Scalar depth
Definition: InitStateEquil.hpp:256
Scalar pressure
Definition: InitStateEquil.hpp:257