OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5377 from akva2/equil_template_scalar

Browse Source 
InitStateEquil: template Scalar type
This commit is contained in:
Bård Skaflestad 2024-05-23 10:48:04 +02:00 committed by GitHub
commit 9841c5d21c
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
7 changed files with 773 additions and 641 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
opm/simulators/flow/equil/EquilibrationHelpers.cpp
View File

@ -27,29 +27,33 @@
namespace Opm { namespace Opm {
namespace EQUIL { namespace EQUIL {
using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>; template<class Scalar>
using FS = BlackOilFluidSystem<double>; using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<Scalar,0,1,2>>;
template<class Scalar> using FS = BlackOilFluidSystem<Scalar>;
template struct PcEq<FS,MatLaw>; #define INSTANTIATE_TYPE(T) \
template struct PcEq<FS<T>,MatLaw<T>>; \
template class EquilReg<T>; \
template T satFromPc<FS<T>,MatLaw<T>>(const MatLaw<T>&, \
const int,const int, \
const T,const bool); \
template T satFromSumOfPcs<FS<T>,MatLaw<T>>(const MatLaw<T>&, \
const int,const int, \
const int,const T); \
template T satFromDepth<FS<T>,MatLaw<T>>(const MatLaw<T>&, \
const T,const T, \
const int,const int,const bool); \
template bool isConstPc<FS<T>,MatLaw<T>>(const MatLaw<T>&,const int,const int); \
template class Miscibility::PBVD<FS<T>>; \
template class Miscibility::PDVD<FS<T>>; \
template class Miscibility::RsVD<FS<T>>; \
template class Miscibility::RsSatAtContact<FS<T>>; \
template class Miscibility::RvSatAtContact<FS<T>>; \
template class Miscibility::RvwSatAtContact<FS<T>>; \
template class Miscibility::RvVD<FS<T>>; \
template class Miscibility::RvwVD<FS<T>>;
template double satFromPc<FS,MatLaw>(const MatLaw&,const int,const int, INSTANTIATE_TYPE(double)
const double,const bool);
template double satFromSumOfPcs<FS,MatLaw>(const MatLaw&,const int,const int,
const int,const double);
template double satFromDepth<FS,MatLaw>(const MatLaw&,const double,const double,
const int,const int,const bool);
template bool isConstPc<FS,MatLaw>(const MatLaw&,const int,const int);
namespace Miscibility {
template class PBVD<FS>;
template class PDVD<FS>;
template class RsVD<FS>;
template class RsSatAtContact<FS>;
template class RvSatAtContact<FS>;
template class RvwSatAtContact<FS>;
template class RvVD<FS>;
template class RvwVD<FS>;
}
} // namespace Equil } // namespace Equil
} // namespace Opm } // namespace Opm

264
opm/simulators/flow/equil/EquilibrationHelpers.hpp
View File

@ -96,6 +96,7 @@ namespace Miscibility {
/** /**
* Base class for phase mixing functions. * Base class for phase mixing functions.
*/ */
template<class Scalar>
class RsFunction class RsFunction
{ {
public: public:
@ -116,17 +117,18 @@ public:
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
virtual double operator()(const double depth, virtual Scalar operator()(const Scalar depth,
const double press, const Scalar press,
const double temp, const Scalar temp,
const double sat = 0.0) const = 0; const Scalar sat = 0.0) const = 0;
}; };
/** /**
* Type that implements "no phase mixing" policy. * Type that implements "no phase mixing" policy.
*/ */
class NoMixing : public RsFunction template<class Scalar>
class NoMixing : public RsFunction<Scalar>
{ {
public: public:
virtual ~NoMixing() = default; virtual ~NoMixing() = default;
@ -147,11 +149,11 @@ public:
* depth and pressure @c press. In "no mixing * depth and pressure @c press. In "no mixing
* policy", this is identically zero. * policy", this is identically zero.
*/ */
double Scalar
operator()(const double /* depth */, operator()(const Scalar /* depth */,
const double /* press */, const Scalar /* press */,
const double /* temp */, const Scalar /* temp */,
const double /* sat */ = 0.0) const const Scalar /* sat */ = 0.0) const override
{ {
return 0.0; return 0.0;
} }
@ -164,9 +166,10 @@ public:
* typically taken from keyword 'RSVD'. * typically taken from keyword 'RSVD'.
*/ */
template <class FluidSystem> template <class FluidSystem>
class RsVD : public RsFunction class RsVD : public RsFunction<typename FluidSystem::Scalar>
{ {
public: public:
using Scalar = typename FluidSystem::Scalar;
/** /**
* Constructor. * Constructor.
* *
@ -175,8 +178,8 @@ public:
* \param[in] rs Dissolved gas-oil ratio at @c depth. * \param[in] rs Dissolved gas-oil ratio at @c depth.
*/ */
RsVD(const int pvtRegionIdx, RsVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& rs); const std::vector<Scalar>& rs);
virtual ~RsVD() = default; virtual ~RsVD() = default;
@ -195,18 +198,18 @@ public:
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double press, const Scalar press,
const double temp, const Scalar temp,
const double satGas = 0.0) const; const Scalar satGas = 0.0) const override;
private: private:
using RsVsDepthFunc = Tabulated1DFunction<double>; using RsVsDepthFunc = Tabulated1DFunction<Scalar>;
const int pvtRegionIdx_; const int pvtRegionIdx_;
RsVsDepthFunc rsVsDepth_; RsVsDepthFunc rsVsDepth_;
double satRs(const double press, const double temp) const; Scalar satRs(const Scalar press, const Scalar temp) const;
}; };
@ -216,9 +219,10 @@ private:
* typically from keyword 'PBVD'. * typically from keyword 'PBVD'.
*/ */
template <class FluidSystem> template <class FluidSystem>
class PBVD : public RsFunction class PBVD : public RsFunction<typename FluidSystem::Scalar>
{ {
public: public:
using Scalar = typename FluidSystem::Scalar;
/** /**
* Constructor. * Constructor.
* *
@ -227,8 +231,8 @@ public:
* \param[in] pbub Bubble-point pressure at @c depth. * \param[in] pbub Bubble-point pressure at @c depth.
*/ */
PBVD(const int pvtRegionIdx, PBVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& pbub); const std::vector<Scalar>& pbub);
virtual ~PBVD() = default; virtual ~PBVD() = default;
@ -246,18 +250,18 @@ public:
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double cellPress, const Scalar cellPress,
const double temp, const Scalar temp,
const double satGas = 0.0) const; const Scalar satGas = 0.0) const override;
private: private:
using PbubVsDepthFunc = Tabulated1DFunction<double>; using PbubVsDepthFunc = Tabulated1DFunction<Scalar>;
const int pvtRegionIdx_; const int pvtRegionIdx_;
PbubVsDepthFunc pbubVsDepth_; PbubVsDepthFunc pbubVsDepth_;
double satRs(const double press, const double temp) const; Scalar satRs(const Scalar press, const Scalar temp) const;
}; };
@ -267,8 +271,9 @@ private:
* taken from keyword 'PDVD'. * taken from keyword 'PDVD'.
*/ */
template <class FluidSystem> template <class FluidSystem>
class PDVD : public RsFunction class PDVD : public RsFunction<typename FluidSystem::Scalar>
{ {
using Scalar = typename FluidSystem::Scalar;
public: public:
/** /**
* Constructor. * Constructor.
@ -278,8 +283,8 @@ public:
* \param[in] pbub Dew-point pressure at @c depth. * \param[in] pbub Dew-point pressure at @c depth.
*/ */
PDVD(const int pvtRegionIdx, PDVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& pdew); const std::vector<Scalar>& pdew);
virtual ~PDVD() = default; virtual ~PDVD() = default;
@ -297,18 +302,18 @@ public:
* \return Vaporized oil-gas ratio (RV) at depth @c * \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double cellPress, const Scalar cellPress,
const double temp, const Scalar temp,
const double satOil = 0.0) const; const Scalar satOil = 0.0) const override;
private: private:
using PdewVsDepthFunc = Tabulated1DFunction<double>; using PdewVsDepthFunc = Tabulated1DFunction<Scalar>;
const int pvtRegionIdx_; const int pvtRegionIdx_;
PdewVsDepthFunc pdewVsDepth_; PdewVsDepthFunc pdewVsDepth_;
double satRv(const double press, const double temp) const; Scalar satRv(const Scalar press, const Scalar temp) const;
}; };
@ -318,8 +323,9 @@ private:
* typically taken from keyword 'RVVD'. * typically taken from keyword 'RVVD'.
*/ */
template <class FluidSystem> template <class FluidSystem>
class RvVD : public RsFunction class RvVD : public RsFunction<typename FluidSystem::Scalar>
{ {
using Scalar = typename FluidSystem::Scalar;
public: public:
/** /**
* Constructor. * Constructor.
@ -329,8 +335,8 @@ public:
* \param[in] rv Dissolved gas-oil ratio at @c depth. * \param[in] rv Dissolved gas-oil ratio at @c depth.
*/ */
RvVD(const int pvtRegionIdx, RvVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& rv); const std::vector<Scalar>& rv);
/** /**
* Function call. * Function call.
@ -347,18 +353,18 @@ public:
* \return Vaporized oil-gas ratio (RV) at depth @c * \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double press, const Scalar press,
const double temp, const Scalar temp,
const double satOil = 0.0) const; const Scalar satOil = 0.0) const override;
private: private:
using RvVsDepthFunc = Tabulated1DFunction<double>; using RvVsDepthFunc = Tabulated1DFunction<Scalar>;
const int pvtRegionIdx_; const int pvtRegionIdx_;
RvVsDepthFunc rvVsDepth_; RvVsDepthFunc rvVsDepth_;
double satRv(const double press, const double temp) const; Scalar satRv(const Scalar press, const Scalar temp) const;
}; };
@ -368,8 +374,9 @@ private:
* typically taken from keyword 'RVWVD'. * typically taken from keyword 'RVWVD'.
*/ */
template <class FluidSystem> template <class FluidSystem>
class RvwVD : public RsFunction class RvwVD : public RsFunction<typename FluidSystem::Scalar>
{ {
using Scalar = typename FluidSystem::Scalar;
public: public:
/** /**
* Constructor. * Constructor.
@ -379,8 +386,8 @@ public:
* \param[in] rvw Evaporized water-gasl ratio at @c depth. * \param[in] rvw Evaporized water-gasl ratio at @c depth.
*/ */
RvwVD(const int pvtRegionIdx, RvwVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& rvw); const std::vector<Scalar>& rvw);
/** /**
* Function call. * Function call.
@ -397,18 +404,18 @@ public:
* \return Vaporized water-gas ratio (RVW) at depth @c * \return Vaporized water-gas ratio (RVW) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double press, const Scalar press,
const double temp, const Scalar temp,
const double satWat = 0.0) const; const Scalar satWat = 0.0) const override;
private: private:
using RvwVsDepthFunc = Tabulated1DFunction<double>; using RvwVsDepthFunc = Tabulated1DFunction<Scalar>;
const int pvtRegionIdx_; const int pvtRegionIdx_;
RvwVsDepthFunc rvwVsDepth_; RvwVsDepthFunc rvwVsDepth_;
double satRvw(const double press, const double temp) const; Scalar satRvw(const Scalar press, const Scalar temp) const;
}; };
@ -427,8 +434,9 @@ private:
* contact, and decreasing above the contact. * contact, and decreasing above the contact.
*/ */
template <class FluidSystem> template <class FluidSystem>
class RsSatAtContact : public RsFunction class RsSatAtContact : public RsFunction<typename FluidSystem::Scalar>
{ {
using Scalar = typename FluidSystem::Scalar;
public: public:
/** /**
* Constructor. * Constructor.
@ -437,7 +445,9 @@ public:
* \param[in] pContact oil pressure at the contact * \param[in] pContact oil pressure at the contact
* \param[in] T_contact temperature at the contact * \param[in] T_contact temperature at the contact
*/ */
RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact); RsSatAtContact(const int pvtRegionIdx,
const Scalar pContact,
const Scalar T_contact);
/** /**
* Function call. * Function call.
@ -454,16 +464,16 @@ public:
* \return Dissolved gas-oil ratio (RS) at depth @c * \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double /* depth */, Scalar operator()(const Scalar /* depth */,
const double press, const Scalar press,
const double temp, const Scalar temp,
const double satGas = 0.0) const; const Scalar satGas = 0.0) const override;
private: private:
const int pvtRegionIdx_; const int pvtRegionIdx_;
double rsSatContact_; Scalar rsSatContact_;
double satRs(const double press, const double temp) const; Scalar satRs(const Scalar press, const Scalar temp) const;
}; };
@ -482,8 +492,9 @@ private:
* contact, and decreasing above the contact. * contact, and decreasing above the contact.
*/ */
template <class FluidSystem> template <class FluidSystem>
class RvSatAtContact : public RsFunction class RvSatAtContact : public RsFunction<typename FluidSystem::Scalar>
{ {
using Scalar = typename FluidSystem::Scalar;
public: public:
/** /**
* Constructor. * Constructor.
@ -492,7 +503,9 @@ public:
* \param[in] pContact oil pressure at the contact * \param[in] pContact oil pressure at the contact
* \param[in] T_contact temperature at the contact * \param[in] T_contact temperature at the contact
*/ */
RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact); RvSatAtContact(const int pvtRegionIdx,
const Scalar pContact,
const Scalar T_contact);
/** /**
* Function call. * Function call.
@ -509,16 +522,16 @@ public:
* \return Dissolved oil-gas ratio (RV) at depth @c * \return Dissolved oil-gas ratio (RV) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double /*depth*/, Scalar operator()(const Scalar /*depth*/,
const double press, const Scalar press,
const double temp, const Scalar temp,
const double satOil = 0.0) const; const Scalar satOil = 0.0) const override;
private: private:
const int pvtRegionIdx_; const int pvtRegionIdx_;
double rvSatContact_; Scalar rvSatContact_;
double satRv(const double press, const double temp) const; Scalar satRv(const Scalar press, const Scalar temp) const;
}; };
/** /**
@ -536,8 +549,9 @@ private:
* contact, and decreasing above the contact. * contact, and decreasing above the contact.
*/ */
template <class FluidSystem> template <class FluidSystem>
class RvwSatAtContact : public RsFunction class RvwSatAtContact : public RsFunction<typename FluidSystem::Scalar>
{ {
using Scalar = typename FluidSystem::Scalar;
public: public:
/** /**
* Constructor. * Constructor.
@ -546,7 +560,9 @@ public:
* \param[in] pContact oil pressure at the contact * \param[in] pContact oil pressure at the contact
* \param[in] T_contact temperature at the contact * \param[in] T_contact temperature at the contact
*/ */
RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact); RvwSatAtContact(const int pvtRegionIdx,
const Scalar pContact,
const Scalar T_contact);
/** /**
* Function call. * Function call.
@ -563,16 +579,16 @@ public:
* \return Dissolved water-gas ratio (RVW) at depth @c * \return Dissolved water-gas ratio (RVW) at depth @c
* depth and pressure @c press. * depth and pressure @c press.
*/ */
double operator()(const double /*depth*/, Scalar operator()(const Scalar /*depth*/,
const double press, const Scalar press,
const double temp, const Scalar temp,
const double satWat = 0.0) const; const Scalar satWat = 0.0) const override;
private: private:
const int pvtRegionIdx_; const int pvtRegionIdx_;
double rvwSatContact_; Scalar rvwSatContact_;
double satRvw(const double press, const double temp) const; Scalar satRvw(const Scalar press, const Scalar temp) const;
}; };
} // namespace Miscibility } // namespace Miscibility
@ -596,9 +612,10 @@ private:
* that calculates the phase densities of all phases in @c * that calculates the phase densities of all phases in @c
* svol at fluid pressure @c press. * svol at fluid pressure @c press.
*/ */
template<class Scalar>
class EquilReg class EquilReg
{ {
using TabulatedFunction = Tabulated1DFunction<double>; using TabulatedFunction = Tabulated1DFunction<Scalar>;
public: public:
/** /**
@ -611,9 +628,9 @@ public:
* \param[in] pvtRegionIdx The pvt region index * \param[in] pvtRegionIdx The pvt region index
*/ */
EquilReg(const EquilRecord& rec, EquilReg(const EquilRecord& rec,
std::shared_ptr<Miscibility::RsFunction> rs, std::shared_ptr<Miscibility::RsFunction<Scalar>> rs,
std::shared_ptr<Miscibility::RsFunction> rv, std::shared_ptr<Miscibility::RsFunction<Scalar>> rv,
std::shared_ptr<Miscibility::RsFunction> rvw, std::shared_ptr<Miscibility::RsFunction<Scalar>> rvw,
const TabulatedFunction& tempVdTable, const TabulatedFunction& tempVdTable,
const TabulatedFunction& saltVdTable, const TabulatedFunction& saltVdTable,
const int pvtIdx); const int pvtIdx);
@ -621,52 +638,52 @@ public:
/** /**
* Type of dissolved gas-oil ratio calculator. * Type of dissolved gas-oil ratio calculator.
*/ */
using CalcDissolution = Miscibility::RsFunction; using CalcDissolution = Miscibility::RsFunction<Scalar>;
/** /**
* Type of vapourised oil-gas ratio calculator. * Type of vapourised oil-gas ratio calculator.
*/ */
using CalcEvaporation = Miscibility::RsFunction; using CalcEvaporation = Miscibility::RsFunction<Scalar>;
/** /**
* Type of vapourised water-gas ratio calculator. * Type of vapourised water-gas ratio calculator.
*/ */
using CalcWaterEvaporation = Miscibility::RsFunction; using CalcWaterEvaporation = Miscibility::RsFunction<Scalar>;
/** /**
* Datum depth in current region * Datum depth in current region
*/ */
double datum() const; Scalar datum() const;
/** /**
* Pressure at datum depth in current region. * Pressure at datum depth in current region.
*/ */
double pressure() const; Scalar pressure() const;
/** /**
* Depth of water-oil contact. * Depth of water-oil contact.
*/ */
double zwoc() const; Scalar zwoc() const;
/** /**
* water-oil capillary pressure at water-oil contact. * water-oil capillary pressure at water-oil contact.
* *
* \return P_o - P_w at WOC. * \return P_o - P_w at WOC.
*/ */
double pcowWoc() const; Scalar pcowWoc() const;
/** /**
* Depth of gas-oil contact. * Depth of gas-oil contact.
*/ */
double zgoc() const; Scalar zgoc() const;
/** /**
* Gas-oil capillary pressure at gas-oil contact. * Gas-oil capillary pressure at gas-oil contact.
* *
* \return P_g - P_o at GOC. * \return P_g - P_o at GOC.
*/ */
double pcgoGoc() const; Scalar pcgoGoc() const;
/** /**
* Accuracy/strategy for initial fluid-in-place calculation. * Accuracy/strategy for initial fluid-in-place calculation.
@ -705,9 +722,9 @@ public:
private: private:
EquilRecord rec_; /**< Equilibration data */ EquilRecord rec_; /**< Equilibration data */
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */ std::shared_ptr<Miscibility::RsFunction<Scalar>> rs_; /**< RS calculator */
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */ std::shared_ptr<Miscibility::RsFunction<Scalar>> rv_; /**< RV calculator */
std::shared_ptr<Miscibility::RsFunction> rvw_; /**< RVW calculator */ std::shared_ptr<Miscibility::RsFunction<Scalar>> rvw_; /**< RVW calculator */
const TabulatedFunction& tempVdTable_; const TabulatedFunction& tempVdTable_;
const TabulatedFunction& saltVdTable_; const TabulatedFunction& saltVdTable_;
const int pvtIdx_; const int pvtIdx_;
@ -721,36 +738,40 @@ private:
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
struct PcEq struct PcEq
{ {
using Scalar = typename FluidSystem::Scalar;
PcEq(const MaterialLawManager& materialLawManager, PcEq(const MaterialLawManager& materialLawManager,
const int phase, const int phase,
const int cell, const int cell,
const double targetPc); const Scalar targetPc);
double operator()(double s) const; Scalar operator()(Scalar s) const;
private: private:
const MaterialLawManager& materialLawManager_; const MaterialLawManager& materialLawManager_;
const int phase_; const int phase_;
const int cell_; const int cell_;
const double targetPc_; const Scalar targetPc_;
}; };
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double minSaturations(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
minSaturations(const MaterialLawManager& materialLawManager,
const int phase, const int cell); const int phase, const int cell);
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double maxSaturations(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const int phase, const int cell); maxSaturations(const MaterialLawManager& materialLawManager,
const int phase, const int cell);
/// Compute saturation of some phase corresponding to a given /// Compute saturation of some phase corresponding to a given
/// capillary pressure. /// capillary pressure.
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double satFromPc(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const int phase, satFromPc(const MaterialLawManager& materialLawManager,
const int cell, const int phase,
const double targetPc, const int cell,
const bool increasing = false); const typename FluidSystem::Scalar targetPc,
const bool increasing = false);
/// Functor for inverting a sum of capillary pressure functions. /// Functor for inverting a sum of capillary pressure functions.
/// Function represented is /// Function represented is
@ -758,40 +779,43 @@ double satFromPc(const MaterialLawManager& materialLawManager,
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
struct PcEqSum struct PcEqSum
{ {
using Scalar = typename FluidSystem::Scalar;
PcEqSum(const MaterialLawManager& materialLawManager, PcEqSum(const MaterialLawManager& materialLawManager,
const int phase1, const int phase1,
const int phase2, const int phase2,
const int cell, const int cell,
const double targetPc); const Scalar targetPc);
double operator()(double s) const; Scalar operator()(Scalar s) const;
private: private:
const MaterialLawManager& materialLawManager_; const MaterialLawManager& materialLawManager_;
const int phase1_; const int phase1_;
const int phase2_; const int phase2_;
const int cell_; const int cell_;
const double targetPc_; const Scalar targetPc_;
}; };
/// Compute saturation of some phase corresponding to a given /// Compute saturation of some phase corresponding to a given
/// capillary pressure, where the capillary pressure function /// capillary pressure, where the capillary pressure function
/// is given as a sum of two other functions. /// is given as a sum of two other functions.
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double satFromSumOfPcs(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const int phase1, satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase2, const int phase1,
const int cell, const int phase2,
const double targetPc); const int cell,
const typename FluidSystem::Scalar targetPc);
/// Compute saturation from depth. Used for constant capillary pressure function /// Compute saturation from depth. Used for constant capillary pressure function
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double satFromDepth(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const double cellDepth, satFromDepth(const MaterialLawManager& materialLawManager,
const double contactDepth, const typename FluidSystem::Scalar cellDepth,
const int phase, const typename FluidSystem::Scalar contactDepth,
const int cell, const int phase,
const bool increasing = false); const int cell,
const bool increasing = false);
/// Return true if capillary pressure function is constant /// Return true if capillary pressure function is constant
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>

365
opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp
View File

@ -56,22 +56,23 @@ namespace Miscibility {
template<class FluidSystem> template<class FluidSystem>
RsVD<FluidSystem>::RsVD(const int pvtRegionIdx, RsVD<FluidSystem>::RsVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& rs) const std::vector<Scalar>& rs)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
, rsVsDepth_(depth, rs) , rsVsDepth_(depth, rs)
{ {
} }
template<class FluidSystem> template<class FluidSystem>
double RsVD<FluidSystem>:: typename RsVD<FluidSystem>::Scalar
operator()(const double depth, RsVD<FluidSystem>::
const double press, operator()(const Scalar depth,
const double temp, const Scalar press,
const double satGas) const const Scalar temp,
const Scalar satGas) const
{ {
const auto sat_rs = satRs(press, temp); const auto sat_rs = satRs(press, temp);
if (satGas > std::sqrt(std::numeric_limits<double>::epsilon())) { if (satGas > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
return sat_rs; return sat_rs;
} }
else { else {
@ -84,28 +85,31 @@ operator()(const double depth,
} }
template<class FluidSystem> template<class FluidSystem>
double RsVD<FluidSystem>::satRs(const double press, const double temp) const typename RsVD<FluidSystem>::Scalar
RsVD<FluidSystem>::
satRs(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
} }
template<class FluidSystem> template<class FluidSystem>
PBVD<FluidSystem>::PBVD(const int pvtRegionIdx, PBVD<FluidSystem>::PBVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& pbub) const std::vector<Scalar>& pbub)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
, pbubVsDepth_(depth, pbub) , pbubVsDepth_(depth, pbub)
{ {
} }
template<class FluidSystem> template<class FluidSystem>
double PBVD<FluidSystem>:: typename PBVD<FluidSystem>::Scalar
operator()(const double depth, PBVD<FluidSystem>::
const double cellPress, operator()(const Scalar depth,
const double temp, const Scalar cellPress,
const double satGas) const const Scalar temp,
const Scalar satGas) const
{ {
double press = cellPress; Scalar press = cellPress;
if (satGas <= 0.0) { if (satGas <= 0.0) {
if (pbubVsDepth_.xMin() > depth) if (pbubVsDepth_.xMin() > depth)
press = pbubVsDepth_.valueAt(0); press = pbubVsDepth_.valueAt(0);
@ -118,29 +122,31 @@ operator()(const double depth,
} }
template<class FluidSystem> template<class FluidSystem>
double PBVD<FluidSystem>:: typename PBVD<FluidSystem>::Scalar
satRs(const double press, const double temp) const PBVD<FluidSystem>::
satRs(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
} }
template<class FluidSystem> template<class FluidSystem>
PDVD<FluidSystem>::PDVD(const int pvtRegionIdx, PDVD<FluidSystem>::PDVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& pdew) const std::vector<Scalar>& pdew)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
, pdewVsDepth_(depth, pdew) , pdewVsDepth_(depth, pdew)
{ {
} }
template<class FluidSystem> template<class FluidSystem>
double PDVD<FluidSystem>:: typename PDVD<FluidSystem>::Scalar
operator()(const double depth, PDVD<FluidSystem>::
const double cellPress, operator()(const Scalar depth,
const double temp, const Scalar cellPress,
const double satOil) const const Scalar temp,
const Scalar satOil) const
{ {
double press = cellPress; Scalar press = cellPress;
if (satOil <= 0.0) { if (satOil <= 0.0) {
if (pdewVsDepth_.xMin() > depth) if (pdewVsDepth_.xMin() > depth)
press = pdewVsDepth_.valueAt(0); press = pdewVsDepth_.valueAt(0);
@ -153,36 +159,37 @@ operator()(const double depth,
} }
template<class FluidSystem> template<class FluidSystem>
double PDVD<FluidSystem>:: typename PDVD<FluidSystem>::Scalar
satRv(const double press, const double temp) const PDVD<FluidSystem>::
satRv(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
} }
template<class FluidSystem> template<class FluidSystem>
RvVD<FluidSystem>::RvVD(const int pvtRegionIdx, RvVD<FluidSystem>::RvVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& rv) const std::vector<Scalar>& rv)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
, rvVsDepth_(depth, rv) , rvVsDepth_(depth, rv)
{ {
} }
template<class FluidSystem> template<class FluidSystem>
double RvVD<FluidSystem>:: typename RvVD<FluidSystem>::Scalar
operator()(const double depth, RvVD<FluidSystem>::
const double press, operator()(const Scalar depth,
const double temp, const Scalar press,
const double satOil) const const Scalar temp,
const Scalar satOil) const
{ {
if (satOil < - std::sqrt(std::numeric_limits<double>::epsilon())) { if (satOil < - std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
throw std::logic_error { throw std::logic_error {
"Must not pass negative oil saturation" "Must not pass negative oil saturation"
}; };
} }
const auto sat_rv = satRv(press, temp); const auto sat_rv = satRv(press, temp);
if (satOil > std::sqrt(std::numeric_limits<double>::epsilon())) { if (satOil > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
return sat_rv; return sat_rv;
} }
else { else {
@ -195,28 +202,29 @@ operator()(const double depth,
} }
template<class FluidSystem> template<class FluidSystem>
double RvVD<FluidSystem>:: typename RvVD<FluidSystem>::Scalar
satRv(const double press, const double temp) const RvVD<FluidSystem>::
satRv(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
} }
template<class FluidSystem> template<class FluidSystem>
RvwVD<FluidSystem>::RvwVD(const int pvtRegionIdx, RvwVD<FluidSystem>::RvwVD(const int pvtRegionIdx,
const std::vector<double>& depth, const std::vector<Scalar>& depth,
const std::vector<double>& rvw) const std::vector<Scalar>& rvw)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
, rvwVsDepth_(depth, rvw) , rvwVsDepth_(depth, rvw)
{ {
} }
template<class FluidSystem> template<class FluidSystem>
double RvwVD<FluidSystem>:: typename RvwVD<FluidSystem>::Scalar
operator()(const double depth, RvwVD<FluidSystem>::
const double press, operator()(const Scalar depth,
const double temp, const Scalar press,
const double satWat) const const Scalar temp,
const Scalar satWat) const
{ {
if (satWat < - std::sqrt(std::numeric_limits<double>::epsilon())) { if (satWat < - std::sqrt(std::numeric_limits<double>::epsilon())) {
throw std::logic_error { throw std::logic_error {
@ -238,29 +246,30 @@ operator()(const double depth,
} }
template<class FluidSystem> template<class FluidSystem>
double RvwVD<FluidSystem>:: typename RvwVD<FluidSystem>::Scalar
satRvw(const double press, const double temp) const RvwVD<FluidSystem>::
satRvw(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press); return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);
} }
template<class FluidSystem> template<class FluidSystem>
RsSatAtContact<FluidSystem>:: RsSatAtContact<FluidSystem>::
RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) RsSatAtContact(const int pvtRegionIdx, const Scalar pContact, const Scalar T_contact)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
{ {
rsSatContact_ = satRs(pContact, T_contact); rsSatContact_ = satRs(pContact, T_contact);
} }
template<class FluidSystem> template<class FluidSystem>
double RsSatAtContact<FluidSystem>:: typename RsSatAtContact<FluidSystem>::Scalar
operator()(const double /* depth */, RsSatAtContact<FluidSystem>::
const double press, operator()(const Scalar /* depth */,
const double temp, const Scalar press,
const double satGas) const const Scalar temp,
const Scalar satGas) const
{ {
if (satGas > std::sqrt(std::numeric_limits<double>::epsilon())) { if (satGas > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
return satRs(press, temp); return satRs(press, temp);
} }
else { else {
@ -269,28 +278,30 @@ operator()(const double /* depth */,
} }
template<class FluidSystem> template<class FluidSystem>
double RsSatAtContact<FluidSystem>:: typename RsSatAtContact<FluidSystem>::Scalar
satRs(const double press, const double temp) const RsSatAtContact<FluidSystem>::
satRs(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
} }
template<class FluidSystem> template<class FluidSystem>
RvSatAtContact<FluidSystem>:: RvSatAtContact<FluidSystem>::
RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) RvSatAtContact(const int pvtRegionIdx, const Scalar pContact, const Scalar T_contact)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
{ {
rvSatContact_ = satRv(pContact, T_contact); rvSatContact_ = satRv(pContact, T_contact);
} }
template<class FluidSystem> template<class FluidSystem>
double RvSatAtContact<FluidSystem>:: typename RvSatAtContact<FluidSystem>::Scalar
operator()(const double /*depth*/, RvSatAtContact<FluidSystem>::
const double press, operator()(const Scalar /*depth*/,
const double temp, const Scalar press,
const double satOil) const const Scalar temp,
const Scalar satOil) const
{ {
if (satOil > std::sqrt(std::numeric_limits<double>::epsilon())) { if (satOil > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
return satRv(press, temp); return satRv(press, temp);
} }
else { else {
@ -299,28 +310,30 @@ operator()(const double /*depth*/,
} }
template<class FluidSystem> template<class FluidSystem>
double RvSatAtContact<FluidSystem>:: typename RvSatAtContact<FluidSystem>::Scalar
satRv(const double press, const double temp) const RvSatAtContact<FluidSystem>::
satRv(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
} }
template<class FluidSystem> template<class FluidSystem>
RvwSatAtContact<FluidSystem>:: RvwSatAtContact<FluidSystem>::
RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) RvwSatAtContact(const int pvtRegionIdx, const Scalar pContact, const Scalar T_contact)
: pvtRegionIdx_(pvtRegionIdx) : pvtRegionIdx_(pvtRegionIdx)
{ {
rvwSatContact_ = satRvw(pContact, T_contact); rvwSatContact_ = satRvw(pContact, T_contact);
} }
template<class FluidSystem> template<class FluidSystem>
double RvwSatAtContact<FluidSystem>:: typename RvwSatAtContact<FluidSystem>::Scalar
operator()(const double /*depth*/, RvwSatAtContact<FluidSystem>::
const double press, operator()(const Scalar /*depth*/,
const double temp, const Scalar press,
const double satWat) const const Scalar temp,
const Scalar satWat) const
{ {
if (satWat > std::sqrt(std::numeric_limits<double>::epsilon())) { if (satWat > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
return satRvw(press, temp); return satRvw(press, temp);
} }
else { else {
@ -329,21 +342,23 @@ operator()(const double /*depth*/,
} }
template<class FluidSystem> template<class FluidSystem>
double RvwSatAtContact<FluidSystem>:: typename RvwSatAtContact<FluidSystem>::Scalar
satRvw(const double press, const double temp) const RvwSatAtContact<FluidSystem>::
satRvw(const Scalar press, const Scalar temp) const
{ {
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);; return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);;
} }
} // namespace Miscibility } // namespace Miscibility
EquilReg::EquilReg(const EquilRecord& rec, template<class Scalar>
std::shared_ptr<Miscibility::RsFunction> rs, EquilReg<Scalar>::EquilReg(const EquilRecord& rec,
std::shared_ptr<Miscibility::RsFunction> rv, std::shared_ptr<Miscibility::RsFunction<Scalar>> rs,
std::shared_ptr<Miscibility::RsFunction> rvw, std::shared_ptr<Miscibility::RsFunction<Scalar>> rv,
const TabulatedFunction& tempVdTable, std::shared_ptr<Miscibility::RsFunction<Scalar>> rvw,
const TabulatedFunction& saltVdTable, const TabulatedFunction& tempVdTable,
const int pvtIdx) const TabulatedFunction& saltVdTable,
const int pvtIdx)
: rec_ (rec) : rec_ (rec)
, rs_ (rs) , rs_ (rs)
, rv_ (rv) , rv_ (rv)
@ -354,72 +369,85 @@ EquilReg::EquilReg(const EquilRecord& rec,
{ {
} }
double EquilReg::datum() const template<class Scalar>
Scalar EquilReg<Scalar>::datum() const
{ {
return this->rec_.datumDepth(); return this->rec_.datumDepth();
} }
double EquilReg::pressure() const template<class Scalar>
Scalar EquilReg<Scalar>::pressure() const
{ {
return this->rec_.datumDepthPressure(); return this->rec_.datumDepthPressure();
} }
double EquilReg::zwoc() const template<class Scalar>
Scalar EquilReg<Scalar>::zwoc() const
{ {
return this->rec_.waterOilContactDepth(); return this->rec_.waterOilContactDepth();
} }
double EquilReg::pcowWoc() const template<class Scalar>
Scalar EquilReg<Scalar>::pcowWoc() const
{ {
return this->rec_.waterOilContactCapillaryPressure(); return this->rec_.waterOilContactCapillaryPressure();
} }
double EquilReg::zgoc() const template<class Scalar>
Scalar EquilReg<Scalar>::zgoc() const
{ {
return this->rec_.gasOilContactDepth(); return this->rec_.gasOilContactDepth();
} }
double EquilReg::pcgoGoc() const template<class Scalar>
Scalar EquilReg<Scalar>::pcgoGoc() const
{ {
return this->rec_.gasOilContactCapillaryPressure(); return this->rec_.gasOilContactCapillaryPressure();
} }
int EquilReg::equilibrationAccuracy() const template<class Scalar>
int EquilReg<Scalar>::equilibrationAccuracy() const
{ {
return this->rec_.initializationTargetAccuracy(); return this->rec_.initializationTargetAccuracy();
} }
const EquilReg::CalcDissolution& template<class Scalar>
EquilReg::dissolutionCalculator() const const typename EquilReg<Scalar>::CalcDissolution&
EquilReg<Scalar>::dissolutionCalculator() const
{ {
return *this->rs_; return *this->rs_;
} }
const EquilReg::CalcEvaporation& template<class Scalar>
EquilReg::evaporationCalculator() const const typename EquilReg<Scalar>::CalcEvaporation&
EquilReg<Scalar>::evaporationCalculator() const
{ {
return *this->rv_; return *this->rv_;
} }
const EquilReg::CalcWaterEvaporation& template<class Scalar>
EquilReg::waterEvaporationCalculator() const const typename EquilReg<Scalar>::CalcWaterEvaporation&
EquilReg<Scalar>::waterEvaporationCalculator() const
{ {
return *this->rvw_; return *this->rvw_;
} }
const EquilReg::TabulatedFunction& template<class Scalar>
EquilReg::saltVdTable() const const typename EquilReg<Scalar>::TabulatedFunction&
EquilReg<Scalar>::saltVdTable() const
{ {
return saltVdTable_; return saltVdTable_;
} }
const EquilReg::TabulatedFunction& template<class Scalar>
EquilReg::tempVdTable() const const typename EquilReg<Scalar>::TabulatedFunction&
EquilReg<Scalar>::tempVdTable() const
{ {
return tempVdTable_; return tempVdTable_;
} }
int EquilReg::pvtIdx() const template<class Scalar>
int EquilReg<Scalar>::pvtIdx() const
{ {
return this->pvtIdx_; return this->pvtIdx_;
} }
@ -429,17 +457,18 @@ PcEq<FluidSystem,MaterialLawManager>::
PcEq(const MaterialLawManager& materialLawManager, PcEq(const MaterialLawManager& materialLawManager,
const int phase, const int phase,
const int cell, const int cell,
const double targetPc) const Scalar targetPc)
: materialLawManager_(materialLawManager), : materialLawManager_(materialLawManager)
phase_(phase), , phase_(phase)
cell_(cell), , cell_(cell)
targetPc_(targetPc) , targetPc_(targetPc)
{ {
} }
template<class FluidSystem, class MaterialLawManager> template<class FluidSystem, class MaterialLawManager>
double PcEq<FluidSystem,MaterialLawManager>:: typename PcEq<FluidSystem,MaterialLawManager>::Scalar
operator()(double s) const PcEq<FluidSystem,MaterialLawManager>::
operator()(Scalar s) const
{ {
const auto& matParams = materialLawManager_.materialLawParams(cell_); const auto& matParams = materialLawManager_.materialLawParams(cell_);
SatOnlyFluidState fluidState; SatOnlyFluidState fluidState;
@ -455,11 +484,11 @@ operator()(double s) const
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0 - s); fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0 - s);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - s); fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - s);
} }
std::array<double, FluidSystem::numPhases> pc{0.0}; std::array<Scalar, FluidSystem::numPhases> pc{0.0};
using MaterialLaw = typename MaterialLawManager::MaterialLaw; using MaterialLaw = typename MaterialLawManager::MaterialLaw;
MaterialLaw::capillaryPressures(pc, matParams, fluidState); MaterialLaw::capillaryPressures(pc, matParams, fluidState);
double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; Scalar sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_]; Scalar pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
return pcPhase - targetPc_; return pcPhase - targetPc_;
} }
@ -469,18 +498,19 @@ PcEqSum(const MaterialLawManager& materialLawManager,
const int phase1, const int phase1,
const int phase2, const int phase2,
const int cell, const int cell,
const double targetPc) const Scalar targetPc)
: materialLawManager_(materialLawManager), : materialLawManager_(materialLawManager)
phase1_(phase1), , phase1_(phase1)
phase2_(phase2), , phase2_(phase2)
cell_(cell), , cell_(cell)
targetPc_(targetPc) , targetPc_(targetPc)
{ {
} }
template<class FluidSystem, class MaterialLawManager> template<class FluidSystem, class MaterialLawManager>
double PcEqSum<FluidSystem,MaterialLawManager>:: typename PcEqSum<FluidSystem,MaterialLawManager>::Scalar
operator()(double s) const PcEqSum<FluidSystem,MaterialLawManager>::
operator()(Scalar s) const
{ {
const auto& matParams = materialLawManager_.materialLawParams(cell_); const auto& matParams = materialLawManager_.materialLawParams(cell_);
SatOnlyFluidState fluidState; SatOnlyFluidState fluidState;
@ -490,19 +520,20 @@ operator()(double s) const
fluidState.setSaturation(phase1_, s); fluidState.setSaturation(phase1_, s);
fluidState.setSaturation(phase2_, 1.0 - s); fluidState.setSaturation(phase2_, 1.0 - s);
std::array<double, FluidSystem::numPhases> pc {0.0}; std::array<Scalar, FluidSystem::numPhases> pc {0.0};
using MaterialLaw = typename MaterialLawManager::MaterialLaw; using MaterialLaw = typename MaterialLawManager::MaterialLaw;
MaterialLaw::capillaryPressures(pc, matParams, fluidState); MaterialLaw::capillaryPressures(pc, matParams, fluidState);
double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; Scalar sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_]; Scalar pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; Scalar sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_]; Scalar pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
return pc1 + pc2 - targetPc_; return pc1 + pc2 - targetPc_;
} }
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double minSaturations(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
minSaturations(const MaterialLawManager& materialLawManager,
const int phase, const int cell) const int phase, const int cell)
{ {
const auto& scaledDrainageInfo = const auto& scaledDrainageInfo =
@ -526,8 +557,9 @@ double minSaturations(const MaterialLawManager& materialLawManager,
} }
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double maxSaturations(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const int phase, const int cell) maxSaturations(const MaterialLawManager& materialLawManager,
const int phase, const int cell)
{ {
const auto& scaledDrainageInfo = const auto& scaledDrainageInfo =
materialLawManager.oilWaterScaledEpsInfoDrainage(cell); materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
@ -550,20 +582,22 @@ double maxSaturations(const MaterialLawManager& materialLawManager,
} }
template <class FluidSystem, class MaterialLawManager> template <class FluidSystem, class MaterialLawManager>
double satFromPc(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const int phase, satFromPc(const MaterialLawManager& materialLawManager,
const int cell, const int phase,
const double targetPc, const int cell,
const bool increasing) const typename FluidSystem::Scalar targetPc,
const bool increasing)
{ {
using Scalar = typename FluidSystem::Scalar;
// Find minimum and maximum saturations. // Find minimum and maximum saturations.
double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell); Scalar s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell); Scalar s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
// Create the equation f(s) = pc(s) - targetPc // Create the equation f(s) = pc(s) - targetPc
const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, targetPc); const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, targetPc);
double f0 = f(s0); Scalar f0 = f(s0);
double f1 = f(s1); Scalar f1 = f(s1);
if (!std::isfinite(f0 + f1)) if (!std::isfinite(f0 + f1))
throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1)); throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1));
@ -572,53 +606,57 @@ double satFromPc(const MaterialLawManager& materialLawManager,
else if (f1 >= 0.0) else if (f1 >= 0.0)
return s1; return s1;
const double tol = 1e-10; const Scalar tol = 1e-10;
// should at least converge in 2 times bisection but some safety here: // should at least converge in 2 times bisection but some safety here:
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10; const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
int usedIterations = -1; int usedIterations = -1;
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations); const Scalar root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
return root; return root;
} }
template<class FluidSystem, class MaterialLawManager> template<class FluidSystem, class MaterialLawManager>
double satFromSumOfPcs(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const int phase1, satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase2, const int phase1,
const int cell, const int phase2,
const double targetPc) const int cell,
const typename FluidSystem::Scalar targetPc)
{ {
using Scalar = typename FluidSystem::Scalar;
// Find minimum and maximum saturations. // Find minimum and maximum saturations.
double s0 = minSaturations<FluidSystem>(materialLawManager, phase1, cell); Scalar s0 = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
double s1 = maxSaturations<FluidSystem>(materialLawManager, phase1, cell); Scalar s1 = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
// Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc // Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc
const PcEqSum<FluidSystem, MaterialLawManager> f(materialLawManager, phase1, phase2, cell, targetPc); const PcEqSum<FluidSystem, MaterialLawManager> f(materialLawManager, phase1, phase2, cell, targetPc);
double f0 = f(s0); Scalar f0 = f(s0);
double f1 = f(s1); Scalar f1 = f(s1);
if (f0 <= 0.0) if (f0 <= 0.0)
return s0; return s0;
else if (f1 >= 0.0) else if (f1 >= 0.0)
return s1; return s1;
assert(f0 > 0.0 && f1 < 0.0); assert(f0 > 0.0 && f1 < 0.0);
const double tol = 1e-10; const Scalar tol = 1e-10;
// should at least converge in 2 times bisection but some safety here: // should at least converge in 2 times bisection but some safety here:
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10; const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
int usedIterations = -1; int usedIterations = -1;
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations); const Scalar root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
return root; return root;
} }
template<class FluidSystem, class MaterialLawManager> template<class FluidSystem, class MaterialLawManager>
double satFromDepth(const MaterialLawManager& materialLawManager, typename FluidSystem::Scalar
const double cellDepth, satFromDepth(const MaterialLawManager& materialLawManager,
const double contactDepth, const typename FluidSystem::Scalar cellDepth,
const int phase, const typename FluidSystem::Scalar contactDepth,
const int cell, const int phase,
const bool increasing) const int cell,
const bool increasing)
{ {
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell); using Scalar = typename FluidSystem::Scalar;
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell); const Scalar s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
const Scalar s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
if (cellDepth < contactDepth) { if (cellDepth < contactDepth) {
return s0; return s0;
@ -633,11 +671,12 @@ bool isConstPc(const MaterialLawManager& materialLawManager,
const int phase, const int phase,
const int cell) const int cell)
{ {
using Scalar = typename FluidSystem::Scalar;
// Create the equation f(s) = pc(s); // Create the equation f(s) = pc(s);
const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, 0); const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, 0);
const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell)); const Scalar f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell)); const Scalar f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon(); return std::abs(f0 - f1) < std::numeric_limits<Scalar>::epsilon();
} }
} }

53
opm/simulators/flow/equil/InitStateEquil.cpp
View File

@ -36,34 +36,36 @@
#endif #endif
namespace Opm { namespace Opm {
template<class Scalar>
using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<Scalar,0,1,2>>;
namespace EQUIL { namespace EQUIL {
namespace DeckDependent { namespace DeckDependent {
using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>; #define INSTANTIATE_COMP(T, GridView, Mapper) \
template class InitialStateComputer<BlackOilFluidSystem<T>, \
#define INSTANCE_COMP(GridView, Mapper) \
template class InitialStateComputer<BlackOilFluidSystem<double>, \
Dune::CpGrid, \ Dune::CpGrid, \
GridView, \ GridView, \
Mapper, \ Mapper, \
Dune::CartesianIndexMapper<Dune::CpGrid>>; \ Dune::CartesianIndexMapper<Dune::CpGrid>>; \
template InitialStateComputer<BlackOilFluidSystem<double>, \ template InitialStateComputer<BlackOilFluidSystem<T>, \
Dune::CpGrid, \ Dune::CpGrid, \
GridView, \ GridView, \
Mapper, \ Mapper, \
Dune::CartesianIndexMapper<Dune::CpGrid>>::\ Dune::CartesianIndexMapper<Dune::CpGrid>>::\
InitialStateComputer(MatLaw&, \ InitialStateComputer(MatLaw<T>&, \
const EclipseState&, \ const EclipseState&, \
const Dune::CpGrid&, \ const Dune::CpGrid&, \
const GridView&, \ const GridView&, \
const Dune::CartesianIndexMapper<Dune::CpGrid>&, \ const Dune::CartesianIndexMapper<Dune::CpGrid>&, \
const double, \ const T, \
const int, \ const int, \
const bool); const bool);
using GridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>; using GridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>;
using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>; using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
INSTANCE_COMP(GridView, Mapper) INSTANTIATE_COMP(double, GridView, Mapper)
#if HAVE_DUNE_FEM #if HAVE_DUNE_FEM
#if DUNE_VERSION_GTE(DUNE_FEM, 2, 9) #if DUNE_VERSION_GTE(DUNE_FEM, 2, 9)
@ -78,23 +80,26 @@ using GridViewFem = Dune::Fem::GridPart2GridViewImpl<
false>>; false>>;
#endif #endif
using MapperFem = Dune::MultipleCodimMultipleGeomTypeMapper<GridViewFem>; using MapperFem = Dune::MultipleCodimMultipleGeomTypeMapper<GridViewFem>;
INSTANCE_COMP(GridViewFem, MapperFem)
INSTANTIATE_COMP(double, GridViewFem, MapperFem)
#endif // HAVE_DUNE_FEM #endif // HAVE_DUNE_FEM
} // namespace DeckDependent } // namespace DeckDependent
namespace Details { namespace Details {
template class PressureTable<BlackOilFluidSystem<double>,EquilReg>; #define INSTANTIATE_TYPE(T) \
template void verticalExtent(const std::vector<int>&, template class PressureTable<BlackOilFluidSystem<T>,EquilReg<T>>; \
const std::vector<std::pair<double,double>>&, template void verticalExtent(const std::vector<int>&, \
const Parallel::Communication&, const std::vector<std::pair<T,T>>&, \
std::array<double,2>&); const Parallel::Communication&, \
std::array<T,2>&); \
template class PhaseSaturations<MatLaw<T>,BlackOilFluidSystem<T>, \
EquilReg<T>,std::size_t>; \
template std::pair<T,T> cellZMinMax<T>(const Dune::cpgrid::Entity<0>&);
using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>; INSTANTIATE_TYPE(double)
template class PhaseSaturations<MatLaw,BlackOilFluidSystem<double>,
EquilReg,std::size_t>;
template std::pair<double,double> cellZMinMax(const Dune::cpgrid::Entity<0>& element);
} }
} // namespace EQUIL } // namespace EQUIL

184
opm/simulators/flow/equil/InitStateEquil.hpp
View File

@ -58,100 +58,106 @@ class NumericalAquifers;
*/ */
namespace EQUIL { namespace EQUIL {
class EquilReg; template<class Scalar> class EquilReg;
namespace Miscibility { class RsFunction; } namespace Miscibility { template<class Scalar> class RsFunction; }
namespace Details { namespace Details {
template <class RHS> template <class Scalar, class RHS>
class RK4IVP { class RK4IVP
{
public: public:
RK4IVP(const RHS& f, RK4IVP(const RHS& f,
const std::array<double,2>& span, const std::array<Scalar,2>& span,
const double y0, const Scalar y0,
const int N); const int N);
double Scalar operator()(const Scalar x) const;
operator()(const double x) const;
private: private:
int N_; int N_;
std::array<double,2> span_; std::array<Scalar,2> span_;
std::vector<double> y_; std::vector<Scalar> y_;
std::vector<double> f_; std::vector<Scalar> f_;
double stepsize() const; Scalar stepsize() const;
}; };
namespace PhasePressODE { namespace PhasePressODE {
template <class FluidSystem> template <class FluidSystem>
class Water class Water
{ {
using TabulatedFunction = Tabulated1DFunction<double>; using Scalar = typename FluidSystem::Scalar;
using TabulatedFunction = Tabulated1DFunction<Scalar>;
public: public:
Water(const TabulatedFunction& tempVdTable, Water(const TabulatedFunction& tempVdTable,
const TabulatedFunction& saltVdTable, const TabulatedFunction& saltVdTable,
const int pvtRegionIdx, const int pvtRegionIdx,
const double normGrav); const Scalar normGrav);
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double press) const; const Scalar press) const;
private: private:
const TabulatedFunction& tempVdTable_; const TabulatedFunction& tempVdTable_;
const TabulatedFunction& saltVdTable_; const TabulatedFunction& saltVdTable_;
const int pvtRegionIdx_; const int pvtRegionIdx_;
const double g_; const Scalar g_;
double density(const double depth, Scalar density(const Scalar depth,
const double press) const; const Scalar press) const;
}; };
template <class FluidSystem, class RS> template <class FluidSystem, class RS>
class Oil class Oil
{ {
using TabulatedFunction = Tabulated1DFunction<double>; using Scalar = typename FluidSystem::Scalar;
using TabulatedFunction = Tabulated1DFunction<Scalar>;
public: public:
Oil(const TabulatedFunction& tempVdTable, Oil(const TabulatedFunction& tempVdTable,
const RS& rs, const RS& rs,
const int pvtRegionIdx, const int pvtRegionIdx,
const double normGrav); const Scalar normGrav);
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double press) const; const Scalar press) const;
private: private:
const TabulatedFunction& tempVdTable_; const TabulatedFunction& tempVdTable_;
const RS& rs_; const RS& rs_;
const int pvtRegionIdx_; const int pvtRegionIdx_;
const double g_; const Scalar g_;
double density(const double depth, Scalar density(const Scalar depth,
const double press) const; const Scalar press) const;
}; };
template <class FluidSystem, class RV, class RVW> template <class FluidSystem, class RV, class RVW>
class Gas class Gas
{ {
using TabulatedFunction = Tabulated1DFunction<double>; using Scalar = typename FluidSystem::Scalar;
using TabulatedFunction = Tabulated1DFunction<Scalar>;
public: public:
Gas(const TabulatedFunction& tempVdTable, Gas(const TabulatedFunction& tempVdTable,
const RV& rv, const RV& rv,
const RVW& rvw, const RVW& rvw,
const int pvtRegionIdx, const int pvtRegionIdx,
const double normGrav); const Scalar normGrav);
double operator()(const double depth, Scalar operator()(const Scalar depth,
const double press) const; const Scalar press) const;
private: private:
const TabulatedFunction& tempVdTable_; const TabulatedFunction& tempVdTable_;
const RV& rv_; const RV& rv_;
const RVW& rvw_; const RVW& rvw_;
const int pvtRegionIdx_; const int pvtRegionIdx_;
const double g_; const Scalar g_;
double density(const double depth, Scalar density(const Scalar depth,
const double press) const; const Scalar press) const;
}; };
} // namespace PhasePressODE } // namespace PhasePressODE
@ -160,7 +166,8 @@ template <class FluidSystem, class Region>
class PressureTable class PressureTable
{ {
public: public:
using VSpan = std::array<double, 2>; using Scalar = typename FluidSystem::Scalar;
using VSpan = std::array<Scalar, 2>;
/// Constructor /// Constructor
/// ///
@ -170,7 +177,7 @@ public:
/// ///
/// \param[in] samplePoints Number of equally spaced depth sample points /// \param[in] samplePoints Number of equally spaced depth sample points
/// in each internal phase pressure table. /// in each internal phase pressure table.
explicit PressureTable(const double gravity, explicit PressureTable(const Scalar gravity,
const int samplePoints = 2000); const int samplePoints = 2000);
/// Copy constructor /// Copy constructor
@ -220,7 +227,7 @@ public:
/// \endcode. /// \endcode.
/// ///
/// \return Oil phase pressure at specified depth. /// \return Oil phase pressure at specified depth.
double oil(const double depth) const; Scalar oil(const Scalar depth) const;
/// Evaluate gas phase pressure at specified depth. /// Evaluate gas phase pressure at specified depth.
/// ///
@ -229,7 +236,7 @@ public:
/// \endcode. /// \endcode.
/// ///
/// \return Gas phase pressure at specified depth. /// \return Gas phase pressure at specified depth.
double gas(const double depth) const; Scalar gas(const Scalar depth) const;
/// Evaluate water phase pressure at specified depth. /// Evaluate water phase pressure at specified depth.
/// ///
@ -238,7 +245,7 @@ public:
/// \endcode. /// \endcode.
/// ///
/// \return Water phase pressure at specified depth. /// \return Water phase pressure at specified depth.
double water(const double depth) const; Scalar water(const Scalar depth) const;
private: private:
template <class ODE> template <class ODE>
@ -246,8 +253,8 @@ private:
{ {
public: public:
struct InitCond { struct InitCond {
double depth; Scalar depth;
double pressure; Scalar pressure;
}; };
explicit PressureFunction(const ODE& ode, explicit PressureFunction(const ODE& ode,
@ -263,12 +270,12 @@ private:
PressureFunction& operator=(PressureFunction&& rhs); PressureFunction& operator=(PressureFunction&& rhs);
double value(const double depth) const; Scalar value(const Scalar depth) const;
private: private:
enum Direction : std::size_t { Up, Down, NumDir }; enum Direction : std::size_t { Up, Down, NumDir };
using Distribution = Details::RK4IVP<ODE>; using Distribution = Details::RK4IVP<Scalar,ODE>;
using DistrPtr = std::unique_ptr<Distribution>; using DistrPtr = std::unique_ptr<Distribution>;
InitCond initial_; InitCond initial_;
@ -292,7 +299,7 @@ private:
using Strategy = void (PressureTable::*) using Strategy = void (PressureTable::*)
(const Region&, const VSpan&); (const Region&, const VSpan&);
double gravity_; Scalar gravity_;
int nsample_; int nsample_;
std::unique_ptr<OPress> oil_{}; std::unique_ptr<OPress> oil_{};
@ -327,12 +334,13 @@ private:
// =========================================================================== // ===========================================================================
/// Simple set of per-phase (named by primary component) quantities. /// Simple set of per-phase (named by primary component) quantities.
template<class Scalar>
struct PhaseQuantityValue { struct PhaseQuantityValue {
double oil{0.0}; Scalar oil{0.0};
double gas{0.0}; Scalar gas{0.0};
double water{0.0}; Scalar water{0.0};
PhaseQuantityValue& axpy(const PhaseQuantityValue& rhs, const double a) PhaseQuantityValue& axpy(const PhaseQuantityValue& rhs, const Scalar a)
{ {
this->oil += a * rhs.oil; this->oil += a * rhs.oil;
this->gas += a * rhs.gas; this->gas += a * rhs.gas;
@ -341,7 +349,7 @@ struct PhaseQuantityValue {
return *this; return *this;
} }
PhaseQuantityValue& operator/=(const double x) PhaseQuantityValue& operator/=(const Scalar x)
{ {
this->oil /= x; this->oil /= x;
this->gas /= x; this->gas /= x;
@ -377,12 +385,13 @@ template <class MaterialLawManager, class FluidSystem, class Region, typename Ce
class PhaseSaturations class PhaseSaturations
{ {
public: public:
using Scalar = typename FluidSystem::Scalar;
/// Evaluation point within a model geometry. /// Evaluation point within a model geometry.
/// ///
/// Associates a particular depth to specific cell. /// Associates a particular depth to specific cell.
struct Position { struct Position {
CellID cell; CellID cell;
double depth; Scalar depth;
}; };
/// Convenience type alias /// Convenience type alias
@ -396,7 +405,7 @@ public:
/// \param[in] swatInit Initial water saturation array (from SWATINIT /// \param[in] swatInit Initial water saturation array (from SWATINIT
/// data). Empty if SWATINIT is not used in this simulation model. /// data). Empty if SWATINIT is not used in this simulation model.
explicit PhaseSaturations(MaterialLawManager& matLawMgr, explicit PhaseSaturations(MaterialLawManager& matLawMgr,
const std::vector<double>& swatInit); const std::vector<Scalar>& swatInit);
/// Copy constructor. /// Copy constructor.
/// ///
@ -421,7 +430,7 @@ public:
/// pertaining to the equilibration region \p reg. /// pertaining to the equilibration region \p reg.
/// ///
/// \return Set of phase saturation values defined at particular point. /// \return Set of phase saturation values defined at particular point.
const PhaseQuantityValue& const PhaseQuantityValue<Scalar>&
deriveSaturations(const Position& x, deriveSaturations(const Position& x,
const Region& reg, const Region& reg,
const PTable& ptable); const PTable& ptable);
@ -430,7 +439,7 @@ public:
/// ///
/// Values associated with evaluation point of previous call to \code /// Values associated with evaluation point of previous call to \code
/// deriveSaturations() \endcode. /// deriveSaturations() \endcode.
const PhaseQuantityValue& correctedPhasePressures() const const PhaseQuantityValue<Scalar>& correctedPhasePressures() const
{ {
return this->press_; return this->press_;
} }
@ -448,7 +457,7 @@ private:
/// information needed to calculate the capillary pressure values from /// information needed to calculate the capillary pressure values from
/// the current set of material laws. /// the current set of material laws.
using FluidState = ::Opm:: using FluidState = ::Opm::
SimpleModularFluidState<double, /*numPhases=*/3, /*numComponents=*/3, SimpleModularFluidState<Scalar, /*numPhases=*/3, /*numComponents=*/3,
FluidSystem, FluidSystem,
/*storePressure=*/false, /*storePressure=*/false,
/*storeTemperature=*/false, /*storeTemperature=*/false,
@ -471,13 +480,13 @@ private:
MaterialLawManager& matLawMgr_; MaterialLawManager& matLawMgr_;
/// Client's SWATINIT data. /// Client's SWATINIT data.
const std::vector<double>& swatInit_; const std::vector<Scalar>& swatInit_;
/// Evaluated phase saturations. /// Evaluated phase saturations.
PhaseQuantityValue sat_; PhaseQuantityValue<Scalar> sat_;
/// Saturation-corrected phase pressure values. /// Saturation-corrected phase pressure values.
PhaseQuantityValue press_; PhaseQuantityValue<Scalar> press_;
/// Current evaluation point. /// Current evaluation point.
EvaluationPoint evalPt_; EvaluationPoint evalPt_;
@ -486,7 +495,7 @@ private:
FluidState fluidState_; FluidState fluidState_;
/// Evaluated capillary pressures from current set of material laws. /// Evaluated capillary pressures from current set of material laws.
std::array<double, FluidSystem::numPhases> matLawCapPress_; std::array<Scalar, FluidSystem::numPhases> matLawCapPress_;
/// Capture the input evaluation point information in internal state. /// Capture the input evaluation point information in internal state.
/// ///
@ -547,7 +556,7 @@ private:
/// \param[in] pcow O/W capillary pressure value (Po - Pw). /// \param[in] pcow O/W capillary pressure value (Po - Pw).
/// ///
/// \return Water saturation value. /// \return Water saturation value.
std::pair<double, bool> applySwatInit(const double pcow); std::pair<Scalar, bool> applySwatInit(const Scalar pcow);
/// Derive water saturation from SWATINIT data. /// Derive water saturation from SWATINIT data.
/// ///
@ -561,7 +570,7 @@ private:
/// ///
/// \return Water saturation value. Input value, possibly mollified by /// \return Water saturation value. Input value, possibly mollified by
/// current set of material laws. /// current set of material laws.
std::pair<double, bool> applySwatInit(const double pc, const double sw); std::pair<Scalar, bool> applySwatInit(const Scalar pc, const Scalar sw);
/// Invoke material law container's capillary pressure calculator on /// Invoke material law container's capillary pressure calculator on
/// current fluid state. /// current fluid state.
@ -569,15 +578,15 @@ private:
/// Extract gas/oil capillary pressure value (Pg - Po) from current /// Extract gas/oil capillary pressure value (Pg - Po) from current
/// fluid state. /// fluid state.
double materialLawCapPressGasOil() const; Scalar materialLawCapPressGasOil() const;
/// Extract oil/water capillary pressure value (Po - Pw) from current /// Extract oil/water capillary pressure value (Po - Pw) from current
/// fluid state. /// fluid state.
double materialLawCapPressOilWater() const; Scalar materialLawCapPressOilWater() const;
/// Extract gas/water capillary pressure value (Pg - Pw) from current /// Extract gas/water capillary pressure value (Pg - Pw) from current
/// fluid state. /// fluid state.
double materialLawCapPressGasWater() const; Scalar materialLawCapPressGasWater() const;
/// Predicate for whether specific phase has constant capillary pressure /// Predicate for whether specific phase has constant capillary pressure
/// curve in current cell. /// curve in current cell.
@ -614,7 +623,7 @@ private:
/// function of phase saturation. /// function of phase saturation.
/// ///
/// \return Phase saturation. /// \return Phase saturation.
double fromDepthTable(const double contactdepth, Scalar fromDepthTable(const Scalar contactdepth,
const PhaseIdx phasePos, const PhaseIdx phasePos,
const bool isincr) const; const bool isincr) const;
@ -635,7 +644,7 @@ private:
/// ///
/// \return Phase saturation at which capillary pressure attains target /// \return Phase saturation at which capillary pressure attains target
/// value. /// value.
double invertCapPress(const double pc, Scalar invertCapPress(const Scalar pc,
const PhaseIdx phasePos, const PhaseIdx phasePos,
const bool isincr) const; const bool isincr) const;
@ -660,14 +669,14 @@ private:
// =========================================================================== // ===========================================================================
template <typename CellRange> template <typename CellRange, class Scalar>
void verticalExtent(const CellRange& cells, void verticalExtent(const CellRange& cells,
const std::vector<std::pair<double, double>>& cellZMinMax, const std::vector<std::pair<Scalar, Scalar>>& cellZMinMax,
const Parallel::Communication& comm, const Parallel::Communication& comm,
std::array<double,2>& span); std::array<Scalar,2>& span);
template <class Element> template <class Scalar, class Element>
std::pair<double,double> cellZMinMax(const Element& element); std::pair<Scalar,Scalar> cellZMinMax(const Element& element);
} // namespace Details } // namespace Details
@ -681,6 +690,7 @@ template<class FluidSystem,
class InitialStateComputer class InitialStateComputer
{ {
using Element = typename GridView::template Codim<0>::Entity; using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename FluidSystem::Scalar;
public: public:
template<class MaterialLawManager> template<class MaterialLawManager>
InitialStateComputer(MaterialLawManager& materialLawManager, InitialStateComputer(MaterialLawManager& materialLawManager,
@ -688,11 +698,11 @@ public:
const Grid& grid, const Grid& grid,
const GridView& gridView, const GridView& gridView,
const CartesianIndexMapper& cartMapper, const CartesianIndexMapper& cartMapper,
const double grav, const Scalar grav,
const int num_pressure_points = 2000, const int num_pressure_points = 2000,
const bool applySwatInit = true); const bool applySwatInit = true);
using Vec = std::vector<double>; using Vec = std::vector<Scalar>;
using PVec = std::vector<Vec>; // One per phase. using PVec = std::vector<Vec>; // One per phase.
const Vec& temperature() const { return temperature_; } const Vec& temperature() const { return temperature_; }
@ -729,29 +739,29 @@ private:
const std::vector<EquilRecord>& rec, const std::vector<EquilRecord>& rec,
MaterialLawManager& materialLawManager, MaterialLawManager& materialLawManager,
const Comm& comm, const Comm& comm,
const double grav); const Scalar grav);
template <class CellRange, class EquilibrationMethod> template <class CellRange, class EquilibrationMethod>
void cellLoop(const CellRange& cells, void cellLoop(const CellRange& cells,
EquilibrationMethod&& eqmethod); EquilibrationMethod&& eqmethod);
template <class CellRange, class PressTable, class PhaseSat> template <class CellRange, class PressTable, class PhaseSat>
void equilibrateCellCentres(const CellRange& cells, void equilibrateCellCentres(const CellRange& cells,
const EquilReg& eqreg, const EquilReg<Scalar>& eqreg,
const PressTable& ptable, const PressTable& ptable,
PhaseSat& psat); PhaseSat& psat);
template <class CellRange, class PressTable, class PhaseSat> template <class CellRange, class PressTable, class PhaseSat>
void equilibrateHorizontal(const CellRange& cells, void equilibrateHorizontal(const CellRange& cells,
const EquilReg& eqreg, const EquilReg<Scalar>& eqreg,
const int acc, const int acc,
const PressTable& ptable, const PressTable& ptable,
PhaseSat& psat); PhaseSat& psat);
std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_; std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rsFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_; std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rvFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rvwFunc_; std::vector< std::shared_ptr<Miscibility::RsFunction<Scalar>> > rvwFunc_;
using TabulatedFunction = Tabulated1DFunction<double>; using TabulatedFunction = Tabulated1DFunction<Scalar>;
std::vector<TabulatedFunction> tempVdTable_; std::vector<TabulatedFunction> tempVdTable_;
std::vector<TabulatedFunction> saltVdTable_; std::vector<TabulatedFunction> saltVdTable_;
std::vector<TabulatedFunction> saltpVdTable_; std::vector<TabulatedFunction> saltpVdTable_;
@ -767,8 +777,8 @@ private:
const CartesianIndexMapper& cartesianIndexMapper_; const CartesianIndexMapper& cartesianIndexMapper_;
Vec swatInit_; Vec swatInit_;
Vec cellCenterDepth_; Vec cellCenterDepth_;
std::vector<std::pair<double,double>> cellZSpan_; std::vector<std::pair<Scalar,Scalar>> cellZSpan_;
std::vector<std::pair<double,double>> cellZMinMax_; std::vector<std::pair<Scalar,Scalar>> cellZMinMax_;
int num_pressure_points_; int num_pressure_points_;
}; };

413
opm/simulators/flow/equil/InitStateEquil_impl.hpp
View File

@ -60,14 +60,14 @@ namespace EQUIL {
namespace Details { namespace Details {
template <typename CellRange> template <typename CellRange, class Scalar>
void verticalExtent(const CellRange& cells, void verticalExtent(const CellRange& cells,
const std::vector<std::pair<double, double>>& cellZMinMax, const std::vector<std::pair<Scalar, Scalar>>& cellZMinMax,
const Parallel::Communication& comm, const Parallel::Communication& comm,
std::array<double,2>& span) std::array<Scalar,2>& span)
{ {
span[0] = std::numeric_limits<double>::max(); span[0] = std::numeric_limits<Scalar>::max();
span[1] = std::numeric_limits<double>::lowest(); span[1] = std::numeric_limits<Scalar>::lowest();
// Define vertical span as // Define vertical span as
// //
@ -85,10 +85,11 @@ void verticalExtent(const CellRange& cells,
span[1] = comm.max(span[1]); span[1] = comm.max(span[1]);
} }
void subdivisionCentrePoints(const double left, template<class Scalar>
const double right, void subdivisionCentrePoints(const Scalar left,
const Scalar right,
const int numIntervals, const int numIntervals,
std::vector<std::pair<double, double>>& subdiv) std::vector<std::pair<Scalar, Scalar>>& subdiv)
{ {
const auto h = (right - left) / numIntervals; const auto h = (right - left) / numIntervals;
@ -101,13 +102,13 @@ void subdivisionCentrePoints(const double left,
} }
} }
template <typename CellID> template <typename CellID, typename Scalar>
std::vector<std::pair<double, double>> std::vector<std::pair<Scalar, Scalar>>
horizontalSubdivision(const CellID cell, horizontalSubdivision(const CellID cell,
const std::pair<double, double> topbot, const std::pair<Scalar, Scalar> topbot,
const int numIntervals) const int numIntervals)
{ {
auto subdiv = std::vector<std::pair<double, double>>{}; auto subdiv = std::vector<std::pair<Scalar, Scalar>>{};
subdiv.reserve(2 * numIntervals); subdiv.reserve(2 * numIntervals);
if (topbot.first > topbot.second) { if (topbot.first > topbot.second) {
@ -123,12 +124,12 @@ horizontalSubdivision(const CellID cell,
return subdiv; return subdiv;
} }
template <class Element> template <class Scalar, class Element>
double cellCenterDepth(const Element& element) Scalar cellCenterDepth(const Element& element)
{ {
typedef typename Element::Geometry Geometry; typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1; static constexpr int zCoord = Element::dimension - 1;
double zz = 0.0; Scalar zz = 0.0;
const Geometry& geometry = element.geometry(); const Geometry& geometry = element.geometry();
const int corners = geometry.corners(); const int corners = geometry.corners();
@ -138,13 +139,13 @@ double cellCenterDepth(const Element& element)
return zz/corners; return zz/corners;
} }
template <class Element> template <class Scalar, class Element>
std::pair<double,double> cellZSpan(const Element& element) std::pair<Scalar,Scalar> cellZSpan(const Element& element)
{ {
typedef typename Element::Geometry Geometry; typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1; static constexpr int zCoord = Element::dimension - 1;
double bot = 0.0; Scalar bot = 0.0;
double top = 0.0; Scalar top = 0.0;
const Geometry& geometry = element.geometry(); const Geometry& geometry = element.geometry();
const int corners = geometry.corners(); const int corners = geometry.corners();
@ -157,36 +158,36 @@ std::pair<double,double> cellZSpan(const Element& element)
return std::make_pair(bot/4, top/4); return std::make_pair(bot/4, top/4);
} }
template <class Element> template <class Scalar, class Element>
std::pair<double,double> cellZMinMax(const Element& element) std::pair<Scalar,Scalar> cellZMinMax(const Element& element)
{ {
typedef typename Element::Geometry Geometry; typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1; static constexpr int zCoord = Element::dimension - 1;
const Geometry& geometry = element.geometry(); const Geometry& geometry = element.geometry();
const int corners = geometry.corners(); const int corners = geometry.corners();
assert(corners == 8); assert(corners == 8);
auto min = std::numeric_limits<double>::max(); auto min = std::numeric_limits<Scalar>::max();
auto max = std::numeric_limits<double>::lowest(); auto max = std::numeric_limits<Scalar>::lowest();
for (int i=0; i < corners; ++i) { for (int i=0; i < corners; ++i) {
min = std::min(min, geometry.corner(i)[zCoord]); min = std::min(min, static_cast<Scalar>(geometry.corner(i)[zCoord]));
max = std::max(max, geometry.corner(i)[zCoord]); max = std::max(max, static_cast<Scalar>(geometry.corner(i)[zCoord]));
} }
return std::make_pair(min, max); return std::make_pair(min, max);
} }
template<class RHS> template<class Scalar, class RHS>
RK4IVP<RHS>::RK4IVP(const RHS& f, RK4IVP<Scalar,RHS>::RK4IVP(const RHS& f,
const std::array<double,2>& span, const std::array<Scalar,2>& span,
const double y0, const Scalar y0,
const int N) const int N)
: N_(N) : N_(N)
, span_(span) , span_(span)
{ {
const double h = stepsize(); const Scalar h = stepsize();
const double h2 = h / 2; const Scalar h2 = h / 2;
const double h6 = h / 6; const Scalar h6 = h / 6;
y_.reserve(N + 1); y_.reserve(N + 1);
f_.reserve(N + 1); f_.reserve(N + 1);
@ -195,39 +196,39 @@ RK4IVP<RHS>::RK4IVP(const RHS& f,
f_.push_back(f(span_[0], y0)); f_.push_back(f(span_[0], y0));
for (int i = 0; i < N; ++i) { for (int i = 0; i < N; ++i) {
const double x = span_[0] + i*h; const Scalar x = span_[0] + i*h;
const double y = y_.back(); const Scalar y = y_.back();
const double k1 = f_[i]; const Scalar k1 = f_[i];
const double k2 = f(x + h2, y + h2*k1); const Scalar k2 = f(x + h2, y + h2*k1);
const double k3 = f(x + h2, y + h2*k2); const Scalar k3 = f(x + h2, y + h2*k2);
const double k4 = f(x + h, y + h*k3); const Scalar k4 = f(x + h, y + h*k3);
y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4)); y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
f_.push_back(f(x + h, y_.back())); f_.push_back(f(x + h, y_.back()));
} }
assert (y_.size() == std::vector<double>::size_type(N + 1)); assert (y_.size() == typename std::vector<Scalar>::size_type(N + 1));
} }
template<class RHS> template<class Scalar, class RHS>
double RK4IVP<RHS>:: Scalar RK4IVP<Scalar,RHS>::
operator()(const double x) const operator()(const Scalar x) const
{ {
// Dense output (O(h**3)) according to Shampine // Dense output (O(h**3)) according to Shampine
// (Hermite interpolation) // (Hermite interpolation)
const double h = stepsize(); const Scalar h = stepsize();
int i = (x - span_[0]) / h; int i = (x - span_[0]) / h;
const double t = (x - (span_[0] + i*h)) / h; const Scalar t = (x - (span_[0] + i*h)) / h;
// Crude handling of evaluation point outside "span_"; // Crude handling of evaluation point outside "span_";
if (i < 0) { i = 0; } if (i < 0) { i = 0; }
if (N_ <= i) { i = N_ - 1; } if (N_ <= i) { i = N_ - 1; }
const double y0 = y_[i], y1 = y_[i + 1]; const Scalar y0 = y_[i], y1 = y_[i + 1];
const double f0 = f_[i], f1 = f_[i + 1]; const Scalar f0 = f_[i], f1 = f_[i + 1];
double u = (1 - 2*t) * (y1 - y0); Scalar u = (1 - 2*t) * (y1 - y0);
u += h * ((t - 1)*f0 + t*f1); u += h * ((t - 1)*f0 + t*f1);
u *= t * (t - 1); u *= t * (t - 1);
u += (1 - t)*y0 + t*y1; u += (1 - t)*y0 + t*y1;
@ -235,8 +236,8 @@ operator()(const double x) const
return u; return u;
} }
template<class RHS> template<class Scalar, class RHS>
double RK4IVP<RHS>:: Scalar RK4IVP<Scalar,RHS>::
stepsize() const stepsize() const
{ {
return (span_[1] - span_[0]) / N_; return (span_[1] - span_[0]) / N_;
@ -249,7 +250,7 @@ Water<FluidSystem>::
Water(const TabulatedFunction& tempVdTable, Water(const TabulatedFunction& tempVdTable,
const TabulatedFunction& saltVdTable, const TabulatedFunction& saltVdTable,
const int pvtRegionIdx, const int pvtRegionIdx,
const double normGrav) const Scalar normGrav)
: tempVdTable_(tempVdTable) : tempVdTable_(tempVdTable)
, saltVdTable_(saltVdTable) , saltVdTable_(saltVdTable)
, pvtRegionIdx_(pvtRegionIdx) , pvtRegionIdx_(pvtRegionIdx)
@ -258,22 +259,28 @@ Water(const TabulatedFunction& tempVdTable,
} }
template<class FluidSystem> template<class FluidSystem>
double Water<FluidSystem>:: typename Water<FluidSystem>::Scalar
operator()(const double depth, Water<FluidSystem>::
const double press) const operator()(const Scalar depth,
const Scalar press) const
{ {
return this->density(depth, press) * g_; return this->density(depth, press) * g_;
} }
template<class FluidSystem> template<class FluidSystem>
double Water<FluidSystem>:: typename Water<FluidSystem>::Scalar
density(const double depth, Water<FluidSystem>::
const double press) const density(const Scalar depth,
const Scalar press) const
{ {
// The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate
double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true); Scalar saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true);
double temp = tempVdTable_.eval(depth, /*extrapolate=*/true); Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true);
double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0 /*=Rsw*/, saltConcentration); Scalar rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_,
temp,
press,
Scalar{0.0} /*=Rsw*/,
saltConcentration);
rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
return rho; return rho;
} }
@ -283,7 +290,7 @@ Oil<FluidSystem,RS>::
Oil(const TabulatedFunction& tempVdTable, Oil(const TabulatedFunction& tempVdTable,
const RS& rs, const RS& rs,
const int pvtRegionIdx, const int pvtRegionIdx,
const double normGrav) const Scalar normGrav)
: tempVdTable_(tempVdTable) : tempVdTable_(tempVdTable)
, rs_(rs) , rs_(rs)
, pvtRegionIdx_(pvtRegionIdx) , pvtRegionIdx_(pvtRegionIdx)
@ -292,31 +299,33 @@ Oil(const TabulatedFunction& tempVdTable,
} }
template<class FluidSystem, class RS> template<class FluidSystem, class RS>
double Oil<FluidSystem,RS>:: typename Oil<FluidSystem,RS>::Scalar
operator()(const double depth, Oil<FluidSystem,RS>::
const double press) const operator()(const Scalar depth,
const Scalar press) const
{ {
return this->density(depth, press) * g_; return this->density(depth, press) * g_;
} }
template<class FluidSystem, class RS> template<class FluidSystem, class RS>
double Oil<FluidSystem,RS>:: typename Oil<FluidSystem,RS>::Scalar
density(const double depth, Oil<FluidSystem,RS>::
const double press) const density(const Scalar depth,
const Scalar press) const
{ {
const double temp = tempVdTable_.eval(depth, /*extrapolate=*/true); const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true);
double rs = 0.0; Scalar rs = 0.0;
if(FluidSystem::enableDissolvedGas()) if (FluidSystem::enableDissolvedGas())
rs = rs_(depth, press, temp); rs = rs_(depth, press, temp);
double bOil = 0.0; Scalar bOil = 0.0;
if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press)) { if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press)) {
bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
} }
else { else {
bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rs); bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rs);
} }
double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); Scalar rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
if (FluidSystem::enableDissolvedGas()) { if (FluidSystem::enableDissolvedGas()) {
rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
} }
@ -330,7 +339,7 @@ Gas(const TabulatedFunction& tempVdTable,
const RV& rv, const RV& rv,
const RVW& rvw, const RVW& rvw,
const int pvtRegionIdx, const int pvtRegionIdx,
const double normGrav) const Scalar normGrav)
: tempVdTable_(tempVdTable) : tempVdTable_(tempVdTable)
, rv_(rv) , rv_(rv)
, rvw_(rvw) , rvw_(rvw)
@ -340,28 +349,30 @@ Gas(const TabulatedFunction& tempVdTable,
} }
template<class FluidSystem, class RV, class RVW> template<class FluidSystem, class RV, class RVW>
double Gas<FluidSystem,RV,RVW>:: typename Gas<FluidSystem,RV,RVW>::Scalar
operator()(const double depth, Gas<FluidSystem,RV,RVW>::
const double press) const operator()(const Scalar depth,
const Scalar press) const
{ {
return this->density(depth, press) * g_; return this->density(depth, press) * g_;
} }
template<class FluidSystem, class RV, class RVW> template<class FluidSystem, class RV, class RVW>
double Gas<FluidSystem,RV,RVW>:: typename Gas<FluidSystem,RV,RVW>::Scalar
density(const double depth, Gas<FluidSystem,RV,RVW>::
const double press) const density(const Scalar depth,
const Scalar press) const
{ {
const double temp = tempVdTable_.eval(depth, /*extrapolate=*/true); const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true);
double rv = 0.0; Scalar rv = 0.0;
if (FluidSystem::enableVaporizedOil()) if (FluidSystem::enableVaporizedOil())
rv = rv_(depth, press, temp); rv = rv_(depth, press, temp);
double rvw = 0.0; Scalar rvw = 0.0;
if (FluidSystem::enableVaporizedWater()) if (FluidSystem::enableVaporizedWater())
rvw = rvw_(depth, press, temp); rvw = rvw_(depth, press, temp);
double bGas = 0.0; Scalar bGas = 0.0;
if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) { if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) {
if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press) if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)
@ -371,7 +382,7 @@ density(const double depth,
} else { } else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw); bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw);
} }
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_) rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_)
+ rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
return rho; return rho;
@ -381,9 +392,13 @@ density(const double depth,
if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) { if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) {
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
} else { } else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, 0.0/*=rvw*/); bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_,
temp,
press,
rv,
Scalar{0.0}/*=rvw*/);
} }
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
return rho; return rho;
} }
@ -393,16 +408,23 @@ density(const double depth,
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
} }
else { else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0/*=rv*/, rvw); bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_,
temp,
press,
Scalar{0.0} /*=rv*/,
rvw);
} }
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
return rho; return rho;
} }
// immiscible gas // immiscible gas
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0/*=rv*/, 0.0/*=rvw*/); bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp,
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); press,
Scalar{0.0} /*=rv*/,
Scalar{0.0} /*=rvw*/);
Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
return rho; return rho;
} }
@ -471,10 +493,10 @@ operator=(PressureFunction&& rhs)
template<class FluidSystem, class Region> template<class FluidSystem, class Region>
template<class ODE> template<class ODE>
double typename PressureTable<FluidSystem,Region>::Scalar
PressureTable<FluidSystem,Region>:: PressureTable<FluidSystem,Region>::
PressureFunction<ODE>:: PressureFunction<ODE>::
value(const double depth) const value(const Scalar depth) const
{ {
if (depth < this->initial_.depth) { if (depth < this->initial_.depth) {
// Value above initial condition depth. // Value above initial condition depth.
@ -548,7 +570,7 @@ copyInPointers(const PressureTable& rhs)
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>:: PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>::
PhaseSaturations(MaterialLawManager& matLawMgr, PhaseSaturations(MaterialLawManager& matLawMgr,
const std::vector<double>& swatInit) const std::vector<Scalar>& swatInit)
: matLawMgr_(matLawMgr) : matLawMgr_(matLawMgr)
, swatInit_ (swatInit) , swatInit_ (swatInit)
{ {
@ -569,7 +591,7 @@ PhaseSaturations(const PhaseSaturations& rhs)
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
const PhaseQuantityValue& const PhaseQuantityValue<typename FluidSystem::Scalar>&
PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>:: PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>::
deriveSaturations(const Position& x, deriveSaturations(const Position& x,
const Region& reg, const Region& reg,
@ -864,17 +886,17 @@ accountForScaledSaturations()
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
std::pair<double, bool> std::pair<typename FluidSystem::Scalar, bool>
PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>:: PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
applySwatInit(const double pcow) applySwatInit(const Scalar pcow)
{ {
return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]); return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]);
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
std::pair<double, bool> std::pair<typename FluidSystem::Scalar, bool>
PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>:: PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
applySwatInit(const double pcow, const double sw) applySwatInit(const Scalar pcow, const Scalar sw)
{ {
return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw); return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw);
} }
@ -892,7 +914,8 @@ computeMaterialLawCapPress()
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>:: typename FluidSystem::Scalar
PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
materialLawCapPressGasOil() const materialLawCapPressGasOil() const
{ {
return this->matLawCapPress_[this->oilPos()] return this->matLawCapPress_[this->oilPos()]
@ -900,7 +923,8 @@ materialLawCapPressGasOil() const
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>:: typename FluidSystem::Scalar
PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
materialLawCapPressOilWater() const materialLawCapPressOilWater() const
{ {
return this->matLawCapPress_[this->oilPos()] return this->matLawCapPress_[this->oilPos()]
@ -908,7 +932,8 @@ materialLawCapPressOilWater() const
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>:: typename FluidSystem::Scalar
PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
materialLawCapPressGasWater() const materialLawCapPressGasWater() const
{ {
return this->matLawCapPress_[this->gasPos()] return this->matLawCapPress_[this->gasPos()]
@ -933,8 +958,9 @@ isOverlappingTransition() const
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>:: typename FluidSystem::Scalar
fromDepthTable(const double contactdepth, PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
fromDepthTable(const Scalar contactdepth,
const PhaseIdx phasePos, const PhaseIdx phasePos,
const bool isincr) const const bool isincr) const
{ {
@ -945,8 +971,9 @@ fromDepthTable(const double contactdepth,
} }
template <class MaterialLawManager, class FluidSystem, class Region, typename CellID> template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>:: typename FluidSystem::Scalar
invertCapPress(const double pc, PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
invertCapPress(const Scalar pc,
const PhaseIdx phasePos, const PhaseIdx phasePos,
const bool isincr) const const bool isincr) const
{ {
@ -957,7 +984,7 @@ invertCapPress(const double pc,
template<class FluidSystem, class Region> template<class FluidSystem, class Region>
PressureTable<FluidSystem,Region>:: PressureTable<FluidSystem,Region>::
PressureTable(const double gravity, PressureTable(const Scalar gravity,
const int samplePoints) const int samplePoints)
: gravity_(gravity) : gravity_(gravity)
, nsample_(samplePoints) , nsample_(samplePoints)
@ -1044,8 +1071,9 @@ waterActive() const
} }
template <class FluidSystem, class Region> template <class FluidSystem, class Region>
double PressureTable<FluidSystem,Region>:: typename FluidSystem::Scalar
oil(const double depth) const PressureTable<FluidSystem,Region>::
oil(const Scalar depth) const
{ {
this->checkPtr(this->oil_.get(), "OIL"); this->checkPtr(this->oil_.get(), "OIL");
@ -1053,8 +1081,9 @@ oil(const double depth) const
} }
template <class FluidSystem, class Region> template <class FluidSystem, class Region>
double PressureTable<FluidSystem,Region>:: typename FluidSystem::Scalar
gas(const double depth) const PressureTable<FluidSystem,Region>::
gas(const Scalar depth) const
{ {
this->checkPtr(this->gas_.get(), "GAS"); this->checkPtr(this->gas_.get(), "GAS");
@ -1063,8 +1092,9 @@ gas(const double depth) const
template <class FluidSystem, class Region> template <class FluidSystem, class Region>
double PressureTable<FluidSystem,Region>:: typename FluidSystem::Scalar
water(const double depth) const PressureTable<FluidSystem,Region>::
water(const Scalar depth) const
{ {
this->checkPtr(this->wat_.get(), "WATER"); this->checkPtr(this->wat_.get(), "WATER");
@ -1310,16 +1340,16 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
const Grid& grid, const Grid& grid,
const GridView& gridView, const GridView& gridView,
const CartesianIndexMapper& cartMapper, const CartesianIndexMapper& cartMapper,
const double grav, const Scalar grav,
const int num_pressure_points, const int num_pressure_points,
const bool applySwatInit) const bool applySwatInit)
: temperature_(grid.size(/*codim=*/0), eclipseState.getTableManager().rtemp()), : temperature_(grid.size(/*codim=*/0), eclipseState.getTableManager().rtemp()),
saltConcentration_(grid.size(/*codim=*/0)), saltConcentration_(grid.size(/*codim=*/0)),
saltSaturation_(grid.size(/*codim=*/0)), saltSaturation_(grid.size(/*codim=*/0)),
pp_(FluidSystem::numPhases, pp_(FluidSystem::numPhases,
std::vector<double>(grid.size(/*codim=*/0))), std::vector<Scalar>(grid.size(/*codim=*/0))),
sat_(FluidSystem::numPhases, sat_(FluidSystem::numPhases,
std::vector<double>(grid.size(/*codim=*/0))), std::vector<Scalar>(grid.size(/*codim=*/0))),
rs_(grid.size(/*codim=*/0)), rs_(grid.size(/*codim=*/0)),
rv_(grid.size(/*codim=*/0)), rv_(grid.size(/*codim=*/0)),
rvw_(grid.size(/*codim=*/0)), rvw_(grid.size(/*codim=*/0)),
@ -1329,7 +1359,13 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
//Check for presence of kw SWATINIT //Check for presence of kw SWATINIT
if (applySwatInit) { if (applySwatInit) {
if (eclipseState.fieldProps().has_double("SWATINIT")) { if (eclipseState.fieldProps().has_double("SWATINIT")) {
swatInit_ = eclipseState.fieldProps().get_double("SWATINIT"); if constexpr (std::is_same_v<Scalar,double>) {
swatInit_ = eclipseState.fieldProps().get_double("SWATINIT");
} else {
const auto& input = eclipseState.fieldProps().get_double("SWATINIT");
swatInit_.resize(input.size());
std::copy(input.begin(), input.end(), swatInit_.begin());
}
} }
} }
@ -1349,6 +1385,19 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
// Create Rs functions. // Create Rs functions.
rsFunc_.reserve(rec.size()); rsFunc_.reserve(rec.size());
auto getArray = [](const std::vector<double>& input)
{
if constexpr (std::is_same_v<Scalar,double>) {
return input;
} else {
std::vector<Scalar> output;
output.resize(input.size());
std::copy(input.begin(), input.end(), output.begin());
return output;
}
};
if (FluidSystem::enableDissolvedGas()) { if (FluidSystem::enableDissolvedGas()) {
for (std::size_t i = 0; i < rec.size(); ++i) { for (std::size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty()) { if (eqlmap.cells(i).empty()) {
@ -1360,16 +1409,15 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
const TableContainer& rsvdTables = tables.getRsvdTables(); const TableContainer& rsvdTables = tables.getRsvdTables();
const TableContainer& pbvdTables = tables.getPbvdTables(); const TableContainer& pbvdTables = tables.getPbvdTables();
if (rsvdTables.size() > 0) { if (rsvdTables.size() > 0) {
const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i); const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); auto depthColumn = getArray(rsvdTable.getColumn("DEPTH").vectorCopy());
std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy(); auto rsColumn = getArray(rsvdTable.getColumn("RS").vectorCopy());
rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx, rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
depthColumn, rsColumn)); depthColumn, rsColumn));
} else if (pbvdTables.size() > 0) { } else if (pbvdTables.size() > 0) {
const PbvdTable& pbvdTable = pbvdTables.getTable<PbvdTable>(i); const PbvdTable& pbvdTable = pbvdTables.getTable<PbvdTable>(i);
std::vector<double> depthColumn = pbvdTable.getColumn("DEPTH").vectorCopy(); auto depthColumn = getArray(pbvdTable.getColumn("DEPTH").vectorCopy());
std::vector<double> pbubColumn = pbvdTable.getColumn("PBUB").vectorCopy(); auto pbubColumn = getArray(pbvdTable.getColumn("PBUB").vectorCopy());
rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx, rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx,
depthColumn, pbubColumn)); depthColumn, pbubColumn));
@ -1384,15 +1432,15 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
"datum depth must be at the gas-oil-contact. " "datum depth must be at the gas-oil-contact. "
"In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold.");
} }
const double pContact = rec[i].datumDepthPressure(); const Scalar pContact = rec[i].datumDepthPressure();
const double TContact = 273.15 + 20; // standard temperature for now const Scalar TContact = 273.15 + 20; // standard temperature for now
rsFunc_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, pContact, TContact)); rsFunc_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, pContact, TContact));
} }
} }
} }
else { else {
for (std::size_t i = 0; i < rec.size(); ++i) { for (std::size_t i = 0; i < rec.size(); ++i) {
rsFunc_.push_back(std::make_shared<Miscibility::NoMixing>()); rsFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
} }
} }
@ -1410,14 +1458,14 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
if (rvvdTables.size() > 0) { if (rvvdTables.size() > 0) {
const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i); const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); auto depthColumn = getArray(rvvdTable.getColumn("DEPTH").vectorCopy());
std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy(); auto rvColumn = getArray(rvvdTable.getColumn("RV").vectorCopy());
rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx, rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
depthColumn, rvColumn)); depthColumn, rvColumn));
} else if (pdvdTables.size() > 0) { } else if (pdvdTables.size() > 0) {
const PdvdTable& pdvdTable = pdvdTables.getTable<PdvdTable>(i); const PdvdTable& pdvdTable = pdvdTables.getTable<PdvdTable>(i);
std::vector<double> depthColumn = pdvdTable.getColumn("DEPTH").vectorCopy(); auto depthColumn = getArray(pdvdTable.getColumn("DEPTH").vectorCopy());
std::vector<double> pdewColumn = pdvdTable.getColumn("PDEW").vectorCopy(); auto pdewColumn = getArray(pdvdTable.getColumn("PDEW").vectorCopy());
rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx, rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx,
depthColumn, pdewColumn)); depthColumn, pdewColumn));
} else { } else {
@ -1431,19 +1479,19 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
"datum depth must be at the gas-oil-contact. " "datum depth must be at the gas-oil-contact. "
"In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold.");
} }
const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); const Scalar pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
const double TContact = 273.15 + 20; // standard temperature for now const Scalar TContact = 273.15 + 20; // standard temperature for now
rvFunc_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact)); rvFunc_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
} }
} }
} }
else { else {
for (std::size_t i = 0; i < rec.size(); ++i) { for (std::size_t i = 0; i < rec.size(); ++i) {
rvFunc_.push_back(std::make_shared<Miscibility::NoMixing>()); rvFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
} }
} }
rvwFunc_.reserve(rec.size()); rvwFunc_.reserve(rec.size());
if (FluidSystem::enableVaporizedWater()) { if (FluidSystem::enableVaporizedWater()) {
for (std::size_t i = 0; i < rec.size(); ++i) { for (std::size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty()) { if (eqlmap.cells(i).empty()) {
@ -1456,8 +1504,8 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
if (rvwvdTables.size() > 0) { if (rvwvdTables.size() > 0) {
const RvwvdTable& rvwvdTable = rvwvdTables.getTable<RvwvdTable>(i); const RvwvdTable& rvwvdTable = rvwvdTables.getTable<RvwvdTable>(i);
std::vector<double> depthColumn = rvwvdTable.getColumn("DEPTH").vectorCopy(); auto depthColumn = getArray(rvwvdTable.getColumn("DEPTH").vectorCopy());
std::vector<double> rvwvdColumn = rvwvdTable.getColumn("RVWVD").vectorCopy(); auto rvwvdColumn = getArray(rvwvdTable.getColumn("RVWVD").vectorCopy());
rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx, rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx,
depthColumn, rvwvdColumn)); depthColumn, rvwvdColumn));
} else { } else {
@ -1466,9 +1514,9 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
} }
else { else {
const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
if (oilActive){ if (oilActive) {
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>()); rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n" const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n"
"and datum depth is not at the gas-oil-contact. \n" "and datum depth is not at the gas-oil-contact. \n"
"Rvw is set to 0.0 in all cells. \n"; "Rvw is set to 0.0 in all cells. \n";
@ -1476,8 +1524,8 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
} else { } else {
// pg = po + Pcgo = po + (pg - po) // pg = po + Pcgo = po + (pg - po)
// for gas-condensate with initial no oil zone: water-oil contact depth (OWC) equal gas-oil contact depth (GOC) // for gas-condensate with initial no oil zone: water-oil contact depth (OWC) equal gas-oil contact depth (GOC)
const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); const Scalar pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
const double TContact = 273.15 + 20; // standard temperature for now const Scalar TContact = 273.15 + 20; // standard temperature for now
rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact)); rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
} }
} }
@ -1485,15 +1533,15 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
// two-phase gas-water sytem: water-oil contact depth is taken equal to gas-water contact depth (GWC) // two-phase gas-water sytem: water-oil contact depth is taken equal to gas-water contact depth (GWC)
// and water-oil capillary pressure (Pcwo) is taken equal to gas-water capillary pressure (Pcgw) at GWC // and water-oil capillary pressure (Pcwo) is taken equal to gas-water capillary pressure (Pcgw) at GWC
if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) { if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) {
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>()); rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n" const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n"
"and datum depth is not at the gas-water-contact. \n" "and datum depth is not at the gas-water-contact. \n"
"Rvw is set to 0.0 in all cells. \n"; "Rvw is set to 0.0 in all cells. \n";
OpmLog::warning(msg); OpmLog::warning(msg);
} else { } else {
// pg = pw + Pcgw = pw + (pg - pw) // pg = pw + Pcgw = pw + (pg - pw)
const double pContact = rec[i].datumDepthPressure() + rec[i].waterOilContactCapillaryPressure(); const Scalar pContact = rec[i].datumDepthPressure() + rec[i].waterOilContactCapillaryPressure();
const double TContact = 273.15 + 20; // standard temperature for now const Scalar TContact = 273.15 + 20; // standard temperature for now
rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact)); rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
} }
} }
@ -1502,7 +1550,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
} }
else { else {
for (std::size_t i = 0; i < rec.size(); ++i) { for (std::size_t i = 0; i < rec.size(); ++i) {
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>()); rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<Scalar>>());
} }
} }
@ -1544,8 +1592,9 @@ updateInitialTemperature_(const EclipseState& eclState, const RMap& reg)
tempVdTable_.resize(numEquilReg); tempVdTable_.resize(numEquilReg);
const auto& tables = eclState.getTableManager(); const auto& tables = eclState.getTableManager();
if (!tables.hasTables("RTEMPVD")) { if (!tables.hasTables("RTEMPVD")) {
std::vector<double> x = {0.0,1.0}; std::vector<Scalar> x = {0.0,1.0};
std::vector<double> y = {tables.rtemp(),tables.rtemp()}; std::vector<Scalar> y = {static_cast<Scalar>(tables.rtemp()),
static_cast<Scalar>(tables.rtemp())};
for (auto& table : this->tempVdTable_) { for (auto& table : this->tempVdTable_) {
table.setXYContainers(x, y); table.setXYContainers(x, y);
} }
@ -1556,7 +1605,7 @@ updateInitialTemperature_(const EclipseState& eclState, const RMap& reg)
tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn()); tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn());
const auto& cells = reg.cells(i); const auto& cells = reg.cells(i);
for (const auto& cell : cells) { for (const auto& cell : cells) {
const double depth = cellCenterDepth_[cell]; const Scalar depth = cellCenterDepth_[cell];
this->temperature_[cell] = tempVdTable_[i].eval(depth, /*extrapolate=*/true); this->temperature_[cell] = tempVdTable_[i].eval(depth, /*extrapolate=*/true);
} }
} }
@ -1583,8 +1632,8 @@ updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg)
// If no saltvd table is given, we create a trivial table for the density calculations // If no saltvd table is given, we create a trivial table for the density calculations
if (saltvdTables.empty()) { if (saltvdTables.empty()) {
std::vector<double> x = {0.0,1.0}; std::vector<Scalar> x = {0.0,1.0};
std::vector<double> y = {0.0,0.0}; std::vector<Scalar> y = {0.0,0.0};
for (auto& table : this->saltVdTable_) { for (auto& table : this->saltVdTable_) {
table.setXYContainers(x, y); table.setXYContainers(x, y);
} }
@ -1595,7 +1644,7 @@ updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg)
const auto& cells = reg.cells(i); const auto& cells = reg.cells(i);
for (const auto& cell : cells) { for (const auto& cell : cells) {
const double depth = cellCenterDepth_[cell]; const Scalar depth = cellCenterDepth_[cell];
this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true); this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true);
} }
} }
@ -1626,7 +1675,7 @@ updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg)
const auto& cells = reg.cells(i); const auto& cells = reg.cells(i);
for (const auto& cell : cells) { for (const auto& cell : cells) {
const double depth = cellCenterDepth_[cell]; const Scalar depth = cellCenterDepth_[cell];
this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true); this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true);
} }
} }
@ -1658,15 +1707,15 @@ updateCellProps_(const GridView& gridView,
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++elemIt) {
const Element& element = *elemIt; const Element& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element); const unsigned int elemIdx = elemMapper.index(element);
cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); cellCenterDepth_[elemIdx] = Details::cellCenterDepth<Scalar>(element);
const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
cellZSpan_[elemIdx] = Details::cellZSpan(element); cellZSpan_[elemIdx] = Details::cellZSpan<Scalar>(element);
cellZMinMax_[elemIdx] = Details::cellZMinMax(element); cellZMinMax_[elemIdx] = Details::cellZMinMax<Scalar>(element);
if (!num_aqu_cells.empty()) { if (!num_aqu_cells.empty()) {
const auto search = num_aqu_cells.find(cartIx); const auto search = num_aqu_cells.find(cartIx);
if (search != num_aqu_cells.end()) { if (search != num_aqu_cells.end()) {
const auto* aqu_cell = num_aqu_cells.at(cartIx); const auto* aqu_cell = num_aqu_cells.at(cartIx);
const double depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx]; const Scalar depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx];
cellCenterDepth_[elemIdx] += depth_change_num_aqu; cellCenterDepth_[elemIdx] += depth_change_num_aqu;
cellZSpan_[elemIdx].first += depth_change_num_aqu; cellZSpan_[elemIdx].first += depth_change_num_aqu;
cellZSpan_[elemIdx].second += depth_change_num_aqu; cellZSpan_[elemIdx].second += depth_change_num_aqu;
@ -1731,7 +1780,7 @@ applyNumericalAquifers_(const GridView& gridView,
// if pressure is specified for numerical aquifers, we use these pressure values // if pressure is specified for numerical aquifers, we use these pressure values
// for numerical aquifer cells // for numerical aquifer cells
if (aqu_cell->init_pressure) { if (aqu_cell->init_pressure) {
const double pres = *(aqu_cell->init_pressure); const Scalar pres = *(aqu_cell->init_pressure);
this->pp_[watPos][elemIdx] = pres; this->pp_[watPos][elemIdx] = pres;
if (FluidSystem::phaseIsActive(gasPos)) { if (FluidSystem::phaseIsActive(gasPos)) {
this->pp_[gasPos][elemIdx] = pres; this->pp_[gasPos][elemIdx] = pres;
@ -1780,15 +1829,15 @@ calcPressSatRsRv(const RMap& reg,
const std::vector<EquilRecord>& rec, const std::vector<EquilRecord>& rec,
MaterialLawManager& materialLawManager, MaterialLawManager& materialLawManager,
const Comm& comm, const Comm& comm,
const double grav) const Scalar grav)
{ {
using PhaseSat = Details::PhaseSaturations< using PhaseSat = Details::PhaseSaturations<
MaterialLawManager, FluidSystem, EquilReg, typename RMap::CellId MaterialLawManager, FluidSystem, EquilReg<Scalar>, typename RMap::CellId
>; >;
auto ptable = Details::PressureTable<FluidSystem, EquilReg>{ grav, this->num_pressure_points_ }; auto ptable = Details::PressureTable<FluidSystem, EquilReg<Scalar>>{ grav, this->num_pressure_points_ };
auto psat = PhaseSat { materialLawManager, this->swatInit_ }; auto psat = PhaseSat { materialLawManager, this->swatInit_ };
auto vspan = std::array<double, 2>{}; auto vspan = std::array<Scalar, 2>{};
std::vector<int> regionIsEmpty(rec.size(), 0); std::vector<int> regionIsEmpty(rec.size(), 0);
for (std::size_t r = 0; r < rec.size(); ++r) { for (std::size_t r = 0; r < rec.size(); ++r) {
@ -1867,11 +1916,11 @@ cellLoop(const CellRange& cells,
const auto gasActive = FluidSystem::phaseIsActive(gasPos); const auto gasActive = FluidSystem::phaseIsActive(gasPos);
const auto watActive = FluidSystem::phaseIsActive(watPos); const auto watActive = FluidSystem::phaseIsActive(watPos);
auto pressures = Details::PhaseQuantityValue{}; auto pressures = Details::PhaseQuantityValue<Scalar>{};
auto saturations = Details::PhaseQuantityValue{}; auto saturations = Details::PhaseQuantityValue<Scalar>{};
auto Rs = 0.0; Scalar Rs = 0.0;
auto Rv = 0.0; Scalar Rv = 0.0;
auto Rvw = 0.0; Scalar Rvw = 0.0;
for (const auto& cell : cells) { for (const auto& cell : cells) {
eqmethod(cell, pressures, saturations, Rs, Rv, Rvw); eqmethod(cell, pressures, saturations, Rs, Rv, Rvw);
@ -1914,7 +1963,7 @@ void InitialStateComputer<FluidSystem,
ElementMapper, ElementMapper,
CartesianIndexMapper>:: CartesianIndexMapper>::
equilibrateCellCentres(const CellRange& cells, equilibrateCellCentres(const CellRange& cells,
const EquilReg& eqreg, const EquilReg<Scalar>& eqreg,
const PressTable& ptable, const PressTable& ptable,
PhaseSat& psat) PhaseSat& psat)
{ {
@ -1923,11 +1972,11 @@ equilibrateCellCentres(const CellRange& cells,
decltype(std::declval<CellPos>().cell)>>; decltype(std::declval<CellPos>().cell)>>;
this->cellLoop(cells, [this, &eqreg, &ptable, &psat] this->cellLoop(cells, [this, &eqreg, &ptable, &psat]
(const CellID cell, (const CellID cell,
Details::PhaseQuantityValue& pressures, Details::PhaseQuantityValue<Scalar>& pressures,
Details::PhaseQuantityValue& saturations, Details::PhaseQuantityValue<Scalar>& saturations,
double& Rs, Scalar& Rs,
double& Rv, Scalar& Rv,
double& Rvw) -> void Scalar& Rvw) -> void
{ {
const auto pos = CellPos { const auto pos = CellPos {
cell, cellCenterDepth_[cell] cell, cellCenterDepth_[cell]
@ -1960,11 +2009,11 @@ void InitialStateComputer<FluidSystem,
GridView, GridView,
ElementMapper, ElementMapper,
CartesianIndexMapper>:: CartesianIndexMapper>::
equilibrateHorizontal(const CellRange& cells, equilibrateHorizontal(const CellRange& cells,
const EquilReg& eqreg, const EquilReg<Scalar>& eqreg,
const int acc, const int acc,
const PressTable& ptable, const PressTable& ptable,
PhaseSat& psat) PhaseSat& psat)
{ {
using CellPos = typename PhaseSat::Position; using CellPos = typename PhaseSat::Position;
using CellID = std::remove_cv_t<std::remove_reference_t< using CellID = std::remove_cv_t<std::remove_reference_t<
@ -1972,16 +2021,16 @@ equilibrateHorizontal(const CellRange& cells,
this->cellLoop(cells, [this, acc, &eqreg, &ptable, &psat] this->cellLoop(cells, [this, acc, &eqreg, &ptable, &psat]
(const CellID cell, (const CellID cell,
Details::PhaseQuantityValue& pressures, Details::PhaseQuantityValue<Scalar>& pressures,
Details::PhaseQuantityValue& saturations, Details::PhaseQuantityValue<Scalar>& saturations,
double& Rs, Scalar& Rs,
double& Rv, Scalar& Rv,
double& Rvw) -> void Scalar& Rvw) -> void
{ {
pressures .reset(); pressures .reset();
saturations.reset(); saturations.reset();
auto totfrac = 0.0; Scalar totfrac = 0.0;
for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) { for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) {
const auto pos = CellPos { cell, depth }; const auto pos = CellPos { cell, depth };

89
tests/test_equil.cpp
View File

@ -137,7 +137,7 @@ static std::vector<std::pair<double,double>> cellVerticalExtent(const GridView&
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++elemIt) {
const auto& element = *elemIt; const auto& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element); const unsigned int elemIdx = elemMapper.index(element);
cellZMinMax[elemIdx] = Opm::EQUIL::Details::cellZMinMax(element); cellZMinMax[elemIdx] = Opm::EQUIL::Details::cellZMinMax<double>(element);
} }
return cellZMinMax; return cellZMinMax;
} }
@ -275,11 +275,12 @@ BOOST_AUTO_TEST_CASE(PhasePressure)
auto simulator = initSimulator<TypeTag>("equil_base.DATA"); auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>(); initDefaultFluidSystem<TypeTag>();
const auto region = Opm::EQUIL::EquilReg { using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const auto region = Opm::EQUIL::EquilReg<double> {
record, record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0 0
@ -296,7 +297,7 @@ BOOST_AUTO_TEST_CASE(PhasePressure)
const auto grav = 10.0; const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable< auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav }; >{ grav };
ptable.equilibrate(region, vspan); ptable.equilibrate(region, vspan);
@ -334,37 +335,37 @@ BOOST_AUTO_TEST_CASE(CellSubset)
std::vector<double> yT = {298.15,298.15}; std::vector<double> yT = {298.15,298.15};
Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT}; Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT};
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const Opm::EQUIL::EquilReg region[] = const Opm::EQUIL::EquilReg<double> region[] =
{ {
Opm::EQUIL::EquilReg(record[0], Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
, ,
Opm::EQUIL::EquilReg(record[0], Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
, ,
Opm::EQUIL::EquilReg(record[1], Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
, ,
Opm::EQUIL::EquilReg(record[1], Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
@ -401,7 +402,7 @@ BOOST_AUTO_TEST_CASE(CellSubset)
const auto grav = 10.0; const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable< auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav }; >{ grav };
auto ppress = PPress(2, PVal(simulator->vanguard().grid().size(0), 0.0)); auto ppress = PPress(2, PVal(simulator->vanguard().grid().size(0), 0.0));
@ -448,37 +449,37 @@ BOOST_AUTO_TEST_CASE(RegMapping)
std::vector<double> yT = {298.15,298.15}; std::vector<double> yT = {298.15,298.15};
Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT}; Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT};
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const Opm::EQUIL::EquilReg region[] = const Opm::EQUIL::EquilReg<double> region[] =
{ {
Opm::EQUIL::EquilReg(record[0], Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
, ,
Opm::EQUIL::EquilReg(record[0], Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
, ,
Opm::EQUIL::EquilReg(record[1], Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
, ,
Opm::EQUIL::EquilReg(record[1], Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<NoMix>(),
trivialTempVdTable, trivialTempVdTable,
trivialSaltVdTable, trivialSaltVdTable,
0) 0)
@ -495,7 +496,7 @@ BOOST_AUTO_TEST_CASE(RegMapping)
const auto grav = 10.0; const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable< auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav }; >{ grav };
std::vector<int> eqlnum(simulator->vanguard().grid().size(0)); std::vector<int> eqlnum(simulator->vanguard().grid().size(0));