OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-07-07 04:53:03 -05:00
Code Issues Packages Projects Releases Wiki Activity

EquilibrationHelpers: template Scalar type

Browse Source
This commit is contained in:
Arne Morten Kvarving 2024-04-12 23:30:19 +02:00
parent 3544dfcabd
commit 4cfb7a8566
7 changed files with 444 additions and 376 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 EQUIL {
using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>;
using FS = BlackOilFluidSystem<double>;
template<class Scalar>
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,
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>;
}
INSTANTIATE_TYPE(double)
} // namespace Equil
} // namespace Opm

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

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

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

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

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

@ -84,7 +84,7 @@ INSTANCE_COMP(GridViewFem, MapperFem)
} // namespace DeckDependent
namespace Details {
template class PressureTable<BlackOilFluidSystem<double>,EquilReg>;
template class PressureTable<BlackOilFluidSystem<double>,EquilReg<double>>;
template void verticalExtent(const std::vector<int>&,
const std::vector<std::pair<double,double>>&,
const Parallel::Communication&,
@ -92,7 +92,7 @@ namespace Details {
using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>;
template class PhaseSaturations<MatLaw,BlackOilFluidSystem<double>,
EquilReg,std::size_t>;
EquilReg<double>,std::size_t>;
template std::pair<double,double> cellZMinMax(const Dune::cpgrid::Entity<0>& element);
}

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

@ -58,8 +58,8 @@ class NumericalAquifers;
*/
namespace EQUIL {
class EquilReg;
namespace Miscibility { class RsFunction; }
template<class Scalar> class EquilReg;
namespace Miscibility { template<class Scalar> class RsFunction; }
namespace Details {
template <class RHS>
@ -736,21 +736,21 @@ private:
EquilibrationMethod&& eqmethod);
template <class CellRange, class PressTable, class PhaseSat>
void equilibrateCellCentres(const CellRange& cells,
const EquilReg& eqreg,
const PressTable& ptable,
PhaseSat& psat);
void equilibrateCellCentres(const CellRange& cells,
const EquilReg<double>& eqreg,
const PressTable& ptable,
PhaseSat& psat);
template <class CellRange, class PressTable, class PhaseSat>
void equilibrateHorizontal(const CellRange& cells,
const EquilReg& eqreg,
const int acc,
const PressTable& ptable,
PhaseSat& psat);
void equilibrateHorizontal(const CellRange& cells,
const EquilReg<double>& eqreg,
const int acc,
const PressTable& ptable,
PhaseSat& psat);
std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rvwFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction<double>> > rsFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction<double>> > rvFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction<double>> > rvwFunc_;
using TabulatedFunction = Tabulated1DFunction<double>;
std::vector<TabulatedFunction> tempVdTable_;
std::vector<TabulatedFunction> saltVdTable_;

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

@ -1392,7 +1392,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
}
else {
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<double>>());
}
}
@ -1439,7 +1439,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
}
else {
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<double>>());
}
}
@ -1468,7 +1468,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
if (oilActive){
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<double>>());
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"
"Rvw is set to 0.0 in all cells. \n";
@ -1485,7 +1485,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
// 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
if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) {
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing<double>>());
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"
"Rvw is set to 0.0 in all cells. \n";
@ -1502,7 +1502,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
}
else {
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<double>>());
}
}
@ -1783,10 +1783,10 @@ calcPressSatRsRv(const RMap& reg,
const double grav)
{
using PhaseSat = Details::PhaseSaturations<
MaterialLawManager, FluidSystem, EquilReg, typename RMap::CellId
MaterialLawManager, FluidSystem, EquilReg<double>, typename RMap::CellId
>;
auto ptable = Details::PressureTable<FluidSystem, EquilReg>{ grav, this->num_pressure_points_ };
auto ptable = Details::PressureTable<FluidSystem, EquilReg<double>>{ grav, this->num_pressure_points_ };
auto psat = PhaseSat { materialLawManager, this->swatInit_ };
auto vspan = std::array<double, 2>{};
@ -1914,7 +1914,7 @@ void InitialStateComputer<FluidSystem,
ElementMapper,
CartesianIndexMapper>::
equilibrateCellCentres(const CellRange& cells,
const EquilReg& eqreg,
const EquilReg<double>& eqreg,
const PressTable& ptable,
PhaseSat& psat)
{
@ -1960,11 +1960,11 @@ void InitialStateComputer<FluidSystem,
GridView,
ElementMapper,
CartesianIndexMapper>::
equilibrateHorizontal(const CellRange& cells,
const EquilReg& eqreg,
const int acc,
const PressTable& ptable,
PhaseSat& psat)
equilibrateHorizontal(const CellRange& cells,
const EquilReg<double>& eqreg,
const int acc,
const PressTable& ptable,
PhaseSat& psat)
{
using CellPos = typename PhaseSat::Position;
using CellID = std::remove_cv_t<std::remove_reference_t<

87
tests/test_equil.cpp
View File

@ -275,11 +275,12 @@ BOOST_AUTO_TEST_CASE(PhasePressure)
auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>();
const auto region = Opm::EQUIL::EquilReg {
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const auto region = Opm::EQUIL::EquilReg<double> {
record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0
@ -296,7 +297,7 @@ BOOST_AUTO_TEST_CASE(PhasePressure)
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg
FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav };
ptable.equilibrate(region, vspan);
@ -334,37 +335,37 @@ BOOST_AUTO_TEST_CASE(CellSubset)
std::vector<double> yT = {298.15,298.15};
Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT};
const Opm::EQUIL::EquilReg region[] =
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const Opm::EQUIL::EquilReg<double> region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
@ -401,7 +402,7 @@ BOOST_AUTO_TEST_CASE(CellSubset)
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg
FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav };
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};
Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT};
const Opm::EQUIL::EquilReg region[] =
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const Opm::EQUIL::EquilReg<double> region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
@ -495,7 +496,7 @@ BOOST_AUTO_TEST_CASE(RegMapping)
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg
FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav };
std::vector<int> eqlnum(simulator->vanguard().grid().size(0));