From 669d3fcd2a0f8d39d3487e3c6486167074f8206d Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 2 Jan 2018 12:43:56 +0100 Subject: [PATCH] equil init: re-indent EquilibrationHelpers.hpp --- ebos/equil/EquilibrationHelpers.hpp | 1442 +++++++++++++-------------- 1 file changed, 721 insertions(+), 721 deletions(-) diff --git a/ebos/equil/EquilibrationHelpers.hpp b/ebos/equil/EquilibrationHelpers.hpp index 517300e96..4d062a9e1 100644 --- a/ebos/equil/EquilibrationHelpers.hpp +++ b/ebos/equil/EquilibrationHelpers.hpp @@ -43,762 +43,762 @@ /* ----- synopsis of EquilibrationHelpers.hpp ---- + ---- synopsis of EquilibrationHelpers.hpp ---- -namespace Opm -{ - namespace EQUIL { + namespace Opm + { + namespace EQUIL { - namespace Miscibility { - class RsFunction; - class NoMixing; - template - class RsVD; - template - class RsSatAtContact; - } + namespace Miscibility { + class RsFunction; + class NoMixing; + template + class RsVD; + template + class RsSatAtContact; + } - class EquilReg; + class EquilReg; - template - struct PcEq; + template + struct PcEq; - template - inline double satFromPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double target_pc, - const bool increasing = false) - template - inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double target_pc) - } // namespace Equil -} // namespace Opm + template + inline double satFromPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double target_pc, + const bool increasing = false) + template + inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, + const int phase1, + const int phase2, + const int cell, + const double target_pc) + } // namespace Equil + } // namespace Opm ----- end of synopsis of EquilibrationHelpers.hpp ---- + ---- end of synopsis of EquilibrationHelpers.hpp ---- */ namespace Ewoms { +/** + * Types and routines that collectively implement a basic + * ECLIPSE-style equilibration-based initialisation scheme. + * + * This namespace is intentionally nested to avoid name clashes + * with other parts of OPM. + */ +namespace EQUIL { + + +typedef Opm::FluidSystems::BlackOil FluidSystemSimple; + +// Adjust oil pressure according to gas saturation and cap pressure +typedef Opm::SimpleModularFluidState SatOnlyFluidState; + +/** + * Types and routines relating to phase mixing in + * equilibration calculations. + */ +namespace Miscibility { + +/** + * Base class for phase mixing functions. + */ +class RsFunction +{ +public: /** - * Types and routines that collectively implement a basic - * ECLIPSE-style equilibration-based initialisation scheme. + * Function call operator. * - * This namespace is intentionally nested to avoid name clashes - * with other parts of OPM. + * \param[in] depth Depth at which to calculate RS + * value. + * + * \param[in] press Pressure at which to calculate RS + * value. + * + * \param[in] temp Temperature at which to calculate RS + * value. + * + * \return Dissolved gas-oil ratio (RS) at depth @c + * depth and pressure @c press. */ - namespace EQUIL { + virtual double operator()(const double depth, + const double press, + const double temp, + const double sat = 0.0) const = 0; +}; - typedef Opm::FluidSystems::BlackOil FluidSystemSimple; - - // Adjust oil pressure according to gas saturation and cap pressure - typedef Opm::SimpleModularFluidState SatOnlyFluidState; - - /** - * Types and routines relating to phase mixing in - * equilibration calculations. - */ - namespace Miscibility { - - /** - * Base class for phase mixing functions. - */ - class RsFunction - { - public: - /** - * Function call operator. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \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; - }; +/** + * Type that implements "no phase mixing" policy. + */ +class NoMixing : public RsFunction { +public: + /** + * Function call. + * + * \param[in] depth Depth at which to calculate RS + * value. + * + * \param[in] press Pressure at which to calculate RS + * value. + * + * \param[in] temp Temperature at which to calculate RS + * value. + * + * \return Dissolved gas-oil ratio (RS) at depth @c + * 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 + { + return 0.0; + } +}; - /** - * Type that implements "no phase mixing" policy. - */ - class NoMixing : public RsFunction { - public: - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \return Dissolved gas-oil ratio (RS) at depth @c - * 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 - { - return 0.0; - } - }; +/** + * Type that implements "dissolved gas-oil ratio" + * tabulated as a function of depth policy. Data + * typically taken from keyword 'RSVD'. + */ +template +class RsVD : public RsFunction { +public: + /** + * Constructor. + * + * \param[in] pvtRegionIdx The pvt region index + * \param[in] depth Depth nodes. + * \param[in] rs Dissolved gas-oil ratio at @c depth. + */ + RsVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rs) + : pvtRegionIdx_(pvtRegionIdx) + , depth_(depth) + , rs_(rs) + { + } - - /** - * Type that implements "dissolved gas-oil ratio" - * tabulated as a function of depth policy. Data - * typically taken from keyword 'RSVD'. - */ - template - class RsVD : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] depth Depth nodes. - * \param[in] rs Dissolved gas-oil ratio at @c depth. - */ - RsVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& rs) - : pvtRegionIdx_(pvtRegionIdx) - , depth_(depth) - , rs_(rs) - { - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \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 sat_gas = 0.0) const - { - if (sat_gas > 0.0) { - return satRs(press, temp); - } else { - return std::min(satRs(press, temp), Opm::linearInterpolationNoExtrapolation(depth_, rs_, depth)); - } - } - - private: - const int pvtRegionIdx_; - std::vector depth_; /**< Depth nodes */ - std::vector rs_; /**< Dissolved gas-oil ratio */ - - double satRs(const double press, const double temp) const - { - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); - } - }; - - - /** - * Type that implements "vaporized oil-gas ratio" - * tabulated as a function of depth policy. Data - * typically taken from keyword 'RVVD'. - */ - template - class RvVD : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] depth Depth nodes. - * \param[in] rv Dissolved gas-oil ratio at @c depth. - */ - RvVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& rv) - : pvtRegionIdx_(pvtRegionIdx) - , depth_(depth) - , rv_(rv) - { - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RV - * value. - * - * \param[in] press Pressure at which to calculate RV - * value. - * - * \param[in] temp Temperature at which to calculate RV - * value. - * - * \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 sat_oil = 0.0 ) const - { - if (std::abs(sat_oil) > 1e-16) { - return satRv(press, temp); - } else { - return std::min(satRv(press, temp), Opm::linearInterpolationNoExtrapolation(depth_, rv_, depth)); - } - } - - private: - const int pvtRegionIdx_; - std::vector depth_; /**< Depth nodes */ - std::vector rv_; /**< Vaporized oil-gas ratio */ - - double satRv(const double press, const double temp) const - { - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); - } - }; - - - /** - * Class that implements "dissolved gas-oil ratio" (Rs) - * as function of depth and pressure as follows: - * - * 1. The Rs at the gas-oil contact is equal to the - * saturated Rs value, Rs_sat_contact. - * - * 2. The Rs elsewhere is equal to Rs_sat_contact, but - * constrained to the saturated value as given by the - * local pressure. - * - * This should yield Rs-values that are constant below the - * contact, and decreasing above the contact. - */ - template - class RsSatAtContact : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] p_contact oil pressure at the contact - * \param[in] T_contact temperature at the contact - */ - RsSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact) - : pvtRegionIdx_(pvtRegionIdx) - { - rs_sat_contact_ = satRs(p_contact, T_contact); - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \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 sat_gas = 0.0) const - { - if (sat_gas > 0.0) { - return satRs(press, temp); - } else { - return std::min(satRs(press, temp), rs_sat_contact_); - } - } - - private: - const int pvtRegionIdx_; - double rs_sat_contact_; - - double satRs(const double press, const double temp) const - { - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); - } - }; - - - /** - * Class that implements "vaporized oil-gas ratio" (Rv) - * as function of depth and pressure as follows: - * - * 1. The Rv at the gas-oil contact is equal to the - * saturated Rv value, Rv_sat_contact. - * - * 2. The Rv elsewhere is equal to Rv_sat_contact, but - * constrained to the saturated value as given by the - * local pressure. - * - * This should yield Rv-values that are constant below the - * contact, and decreasing above the contact. - */ - template - class RvSatAtContact : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] p_contact oil pressure at the contact - * \param[in] T_contact temperature at the contact - */ - RvSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact) - :pvtRegionIdx_(pvtRegionIdx) - { - rv_sat_contact_ = satRv(p_contact, T_contact); - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RV - * value. - * - * \param[in] press Pressure at which to calculate RV - * value. - * - * \param[in] temp Temperature at which to calculate RV - * value. - * - * \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 sat_oil = 0.0) const - { - if (sat_oil > 0.0) { - return satRv(press, temp); - } else { - return std::min(satRv(press, temp), rv_sat_contact_); - } - } - - private: - const int pvtRegionIdx_; - double rv_sat_contact_; - - double satRv(const double press, const double temp) const - { - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; - } - }; - - } // namespace Miscibility - - /** - * Aggregate information base of an equilibration region. - * - * Provides inquiry methods for retrieving depths of contacs - * and pressure values as well as a means of calculating fluid - * densities, dissolved gas-oil ratio and vapourised oil-gas - * ratios. - * - * \tparam DensCalc Type that provides access to a phase - * density calculation facility. Must implement an operator() - * declared as - * - * std::vector - * operator()(const double press, - * const std::vector& svol ) - * - * that calculates the phase densities of all phases in @c - * svol at fluid pressure @c press. - */ - class EquilReg { - public: - /** - * Constructor. - * - * \param[in] rec Equilibration data of current region. - * \param[in] rs Calculator of dissolved gas-oil ratio. - * \param[in] rv Calculator of vapourised oil-gas ratio. - * \param[in] pvtRegionIdx The pvt region index - */ - EquilReg(const Opm::EquilRecord& rec, - std::shared_ptr rs, - std::shared_ptr rv, - const int pvtIdx) - : rec_ (rec) - , rs_ (rs) - , rv_ (rv) - , pvtIdx_ (pvtIdx) - { - } - - /** - * Type of dissolved gas-oil ratio calculator. - */ - typedef Miscibility::RsFunction CalcDissolution; - - /** - * Type of vapourised oil-gas ratio calculator. - */ - typedef Miscibility::RsFunction CalcEvaporation; - - /** - * Datum depth in current region - */ - double datum() const { return this->rec_.datumDepth(); } - - /** - * Pressure at datum depth in current region. - */ - double pressure() const { return this->rec_.datumDepthPressure(); } - - /** - * Depth of water-oil contact. - */ - double zwoc() const { return this->rec_.waterOilContactDepth(); } - - /** - * water-oil capillary pressure at water-oil contact. - * - * \return P_o - P_w at WOC. - */ - double pcow_woc() const { return this->rec_.waterOilContactCapillaryPressure(); } - - /** - * Depth of gas-oil contact. - */ - double zgoc() const { return this->rec_.gasOilContactDepth(); } - - /** - * Gas-oil capillary pressure at gas-oil contact. - * - * \return P_g - P_o at GOC. - */ - double pcgo_goc() const { return this->rec_.gasOilContactCapillaryPressure(); } - - - /** - * Retrieve dissolved gas-oil ratio calculator of current - * region. - */ - const CalcDissolution& - dissolutionCalculator() const { return *this->rs_; } - - /** - * Retrieve vapourised oil-gas ratio calculator of current - * region. - */ - const CalcEvaporation& - evaporationCalculator() const { return *this->rv_; } - - /** - * Retrieve pvtIdx of the region. - */ - int pvtIdx() const { return this->pvtIdx_; } - - - private: - Opm::EquilRecord rec_; /**< Equilibration data */ - std::shared_ptr rs_; /**< RS calculator */ - std::shared_ptr rv_; /**< RV calculator */ - const int pvtIdx_; - }; - - - - /// Functor for inverting capillary pressure function. - /// Function represented is - /// f(s) = pc(s) - target_pc - template - struct PcEq - { - PcEq(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double target_pc) - : materialLawManager_(materialLawManager), - phase_(phase), - cell_(cell), - target_pc_(target_pc) - { - - } - double operator()(double s) const - { - const auto& matParams = materialLawManager_.materialLawParams(cell_); - SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); - fluidState.setSaturation(phase_, s); - - double pc[FluidSystem::numPhases]; - std::fill(pc, pc + FluidSystem::numPhases, 0.0); - MaterialLaw::capillaryPressures(pc, matParams, fluidState); - double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; - double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_]; - return pcPhase - target_pc_; - } - private: - const MaterialLawManager& materialLawManager_; - const int phase_; - const int cell_; - const double target_pc_; - }; - - template - double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); - - // Find minimum and maximum saturations. - switch(phase) { - case FluidSystem::waterPhaseIdx : - { - return scaledDrainageInfo.Swl; - break; - } - case FluidSystem::gasPhaseIdx : - { - return scaledDrainageInfo.Sgl; - break; - } - case FluidSystem::oilPhaseIdx : - { - OPM_THROW(std::runtime_error, "Min saturation not implemented for oil phase."); - break; - } - default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); - } - return -1.0; + /** + * Function call. + * + * \param[in] depth Depth at which to calculate RS + * value. + * + * \param[in] press Pressure at which to calculate RS + * value. + * + * \param[in] temp Temperature at which to calculate RS + * value. + * + * \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 sat_gas = 0.0) const + { + if (sat_gas > 0.0) { + return satRs(press, temp); + } else { + return std::min(satRs(press, temp), Opm::linearInterpolationNoExtrapolation(depth_, rs_, depth)); } + } - template - double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); +private: + const int pvtRegionIdx_; + std::vector depth_; /**< Depth nodes */ + std::vector rs_; /**< Dissolved gas-oil ratio */ - // Find minimum and maximum saturations. - switch(phase) { - case FluidSystem::waterPhaseIdx : - { - return scaledDrainageInfo.Swu; - break; - } - case FluidSystem::gasPhaseIdx : - { - return scaledDrainageInfo.Sgu; - break; - } - case FluidSystem::oilPhaseIdx : - { - OPM_THROW(std::runtime_error, "Max saturation not implemented for oil phase."); - break; - } - default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); - } - return -1.0; + double satRs(const double press, const double temp) const + { + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); + } +}; + + +/** + * Type that implements "vaporized oil-gas ratio" + * tabulated as a function of depth policy. Data + * typically taken from keyword 'RVVD'. + */ +template +class RvVD : public RsFunction { +public: + /** + * Constructor. + * + * \param[in] pvtRegionIdx The pvt region index + * \param[in] depth Depth nodes. + * \param[in] rv Dissolved gas-oil ratio at @c depth. + */ + RvVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rv) + : pvtRegionIdx_(pvtRegionIdx) + , depth_(depth) + , rv_(rv) + { + } + + /** + * Function call. + * + * \param[in] depth Depth at which to calculate RV + * value. + * + * \param[in] press Pressure at which to calculate RV + * value. + * + * \param[in] temp Temperature at which to calculate RV + * value. + * + * \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 sat_oil = 0.0 ) const + { + if (std::abs(sat_oil) > 1e-16) { + return satRv(press, temp); + } else { + return std::min(satRv(press, temp), Opm::linearInterpolationNoExtrapolation(depth_, rv_, depth)); } + } + +private: + const int pvtRegionIdx_; + std::vector depth_; /**< Depth nodes */ + std::vector rv_; /**< Vaporized oil-gas ratio */ + + double satRv(const double press, const double temp) const + { + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); + } +}; - /// Compute saturation of some phase corresponding to a given - /// capillary pressure. - template - inline double satFromPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double target_pc, - const bool increasing = false) - { - // Find minimum and maximum saturations. - const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); - const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); +/** + * Class that implements "dissolved gas-oil ratio" (Rs) + * as function of depth and pressure as follows: + * + * 1. The Rs at the gas-oil contact is equal to the + * saturated Rs value, Rs_sat_contact. + * + * 2. The Rs elsewhere is equal to Rs_sat_contact, but + * constrained to the saturated value as given by the + * local pressure. + * + * This should yield Rs-values that are constant below the + * contact, and decreasing above the contact. + */ +template +class RsSatAtContact : public RsFunction { +public: + /** + * Constructor. + * + * \param[in] pvtRegionIdx The pvt region index + * \param[in] p_contact oil pressure at the contact + * \param[in] T_contact temperature at the contact + */ + RsSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact) + : pvtRegionIdx_(pvtRegionIdx) + { + rs_sat_contact_ = satRs(p_contact, T_contact); + } - // Create the equation f(s) = pc(s) - target_pc - const PcEq f(materialLawManager, phase, cell, target_pc); - const double f0 = f(s0); - const double f1 = f(s1); - - if (f0 <= 0.0) { - return s0; - } else if (f1 > 0.0) { - return s1; - } else { - const int max_iter = 60; - const double tol = 1e-6; - int iter_used = -1; - typedef Opm::RegulaFalsi ScalarSolver; - const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used); - return sol; - } + /** + * Function call. + * + * \param[in] depth Depth at which to calculate RS + * value. + * + * \param[in] press Pressure at which to calculate RS + * value. + * + * \param[in] temp Temperature at which to calculate RS + * value. + * + * \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 sat_gas = 0.0) const + { + if (sat_gas > 0.0) { + return satRs(press, temp); + } else { + return std::min(satRs(press, temp), rs_sat_contact_); } + } + +private: + const int pvtRegionIdx_; + double rs_sat_contact_; + + double satRs(const double press, const double temp) const + { + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); + } +}; - /// Functor for inverting a sum of capillary pressure functions. - /// Function represented is - /// f(s) = pc1(s) + pc2(1 - s) - target_pc - template - struct PcEqSum - { - PcEqSum(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double target_pc) - : materialLawManager_(materialLawManager), - phase1_(phase1), - phase2_(phase2), - cell_(cell), - target_pc_(target_pc) - { - } - double operator()(double s) const - { - const auto& matParams = materialLawManager_.materialLawParams(cell_); - SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); - fluidState.setSaturation(phase1_, s); - fluidState.setSaturation(phase2_, 1.0 - s); +/** + * Class that implements "vaporized oil-gas ratio" (Rv) + * as function of depth and pressure as follows: + * + * 1. The Rv at the gas-oil contact is equal to the + * saturated Rv value, Rv_sat_contact. + * + * 2. The Rv elsewhere is equal to Rv_sat_contact, but + * constrained to the saturated value as given by the + * local pressure. + * + * This should yield Rv-values that are constant below the + * contact, and decreasing above the contact. + */ +template +class RvSatAtContact : public RsFunction { +public: + /** + * Constructor. + * + * \param[in] pvtRegionIdx The pvt region index + * \param[in] p_contact oil pressure at the contact + * \param[in] T_contact temperature at the contact + */ + RvSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact) + :pvtRegionIdx_(pvtRegionIdx) + { + rv_sat_contact_ = satRv(p_contact, T_contact); + } - double pc[FluidSystem::numPhases]; - std::fill(pc, pc + FluidSystem::numPhases, 0.0); - - 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_]; - return pc1 + pc2 - target_pc_; - } - private: - const MaterialLawManager& materialLawManager_; - const int phase1_; - const int phase2_; - const int cell_; - const double target_pc_; - }; - - - - - /// 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 - inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double target_pc) - { - // Find minimum and maximum saturations. - const double smin = minSaturations(materialLawManager, phase1, cell); - const double smax = maxSaturations(materialLawManager, phase1, cell); - - // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc - const PcEqSum f(materialLawManager, phase1, phase2, cell, target_pc); - const double f0 = f(smin); - const double f1 = f(smax); - if (f0 <= 0.0) { - return smin; - } else if (f1 > 0.0) { - return smax; - } else { - const int max_iter = 30; - const double tol = 1e-6; - int iter_used = -1; - typedef Opm::RegulaFalsi ScalarSolver; - const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used); - return sol; - } + /** + * Function call. + * + * \param[in] depth Depth at which to calculate RV + * value. + * + * \param[in] press Pressure at which to calculate RV + * value. + * + * \param[in] temp Temperature at which to calculate RV + * value. + * + * \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 sat_oil = 0.0) const + { + if (sat_oil > 0.0) { + return satRv(press, temp); + } else { + return std::min(satRv(press, temp), rv_sat_contact_); } + } - /// Compute saturation from depth. Used for constant capillary pressure function - template - inline double satFromDepth(const MaterialLawManager& materialLawManager, - const double cellDepth, - const double contactDepth, - const int phase, - const int cell, - const bool increasing = false) - { - const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); - const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); +private: + const int pvtRegionIdx_; + double rv_sat_contact_; - if (cellDepth < contactDepth){ - return s0; - } else { - return s1; - } + double satRv(const double press, const double temp) const + { + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; + } +}; - } +} // namespace Miscibility - /// Return true if capillary pressure function is constant - template - inline bool isConstPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell) - { - // Create the equation f(s) = pc(s); - const PcEq f(materialLawManager, phase, cell, 0); - const double f0 = f(minSaturations(materialLawManager, phase, cell)); - const double f1 = f(maxSaturations(materialLawManager, phase, cell)); - return std::abs(f0 - f1) < std::numeric_limits::epsilon(); - } +/** + * Aggregate information base of an equilibration region. + * + * Provides inquiry methods for retrieving depths of contacs + * and pressure values as well as a means of calculating fluid + * densities, dissolved gas-oil ratio and vapourised oil-gas + * ratios. + * + * \tparam DensCalc Type that provides access to a phase + * density calculation facility. Must implement an operator() + * declared as + * + * std::vector + * operator()(const double press, + * const std::vector& svol ) + * + * that calculates the phase densities of all phases in @c + * svol at fluid pressure @c press. + */ +class EquilReg { +public: + /** + * Constructor. + * + * \param[in] rec Equilibration data of current region. + * \param[in] rs Calculator of dissolved gas-oil ratio. + * \param[in] rv Calculator of vapourised oil-gas ratio. + * \param[in] pvtRegionIdx The pvt region index + */ + EquilReg(const Opm::EquilRecord& rec, + std::shared_ptr rs, + std::shared_ptr rv, + const int pvtIdx) + : rec_ (rec) + , rs_ (rs) + , rv_ (rv) + , pvtIdx_ (pvtIdx) + { + } - } // namespace Equil + /** + * Type of dissolved gas-oil ratio calculator. + */ + typedef Miscibility::RsFunction CalcDissolution; + + /** + * Type of vapourised oil-gas ratio calculator. + */ + typedef Miscibility::RsFunction CalcEvaporation; + + /** + * Datum depth in current region + */ + double datum() const { return this->rec_.datumDepth(); } + + /** + * Pressure at datum depth in current region. + */ + double pressure() const { return this->rec_.datumDepthPressure(); } + + /** + * Depth of water-oil contact. + */ + double zwoc() const { return this->rec_.waterOilContactDepth(); } + + /** + * water-oil capillary pressure at water-oil contact. + * + * \return P_o - P_w at WOC. + */ + double pcow_woc() const { return this->rec_.waterOilContactCapillaryPressure(); } + + /** + * Depth of gas-oil contact. + */ + double zgoc() const { return this->rec_.gasOilContactDepth(); } + + /** + * Gas-oil capillary pressure at gas-oil contact. + * + * \return P_g - P_o at GOC. + */ + double pcgo_goc() const { return this->rec_.gasOilContactCapillaryPressure(); } + + + /** + * Retrieve dissolved gas-oil ratio calculator of current + * region. + */ + const CalcDissolution& + dissolutionCalculator() const { return *this->rs_; } + + /** + * Retrieve vapourised oil-gas ratio calculator of current + * region. + */ + const CalcEvaporation& + evaporationCalculator() const { return *this->rv_; } + + /** + * Retrieve pvtIdx of the region. + */ + int pvtIdx() const { return this->pvtIdx_; } + + +private: + Opm::EquilRecord rec_; /**< Equilibration data */ + std::shared_ptr rs_; /**< RS calculator */ + std::shared_ptr rv_; /**< RV calculator */ + const int pvtIdx_; +}; + + + +/// Functor for inverting capillary pressure function. +/// Function represented is +/// f(s) = pc(s) - target_pc +template +struct PcEq +{ + PcEq(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double target_pc) + : materialLawManager_(materialLawManager), + phase_(phase), + cell_(cell), + target_pc_(target_pc) + { + + } + double operator()(double s) const + { + const auto& matParams = materialLawManager_.materialLawParams(cell_); + SatOnlyFluidState fluidState; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); + fluidState.setSaturation(phase_, s); + + double pc[FluidSystem::numPhases]; + std::fill(pc, pc + FluidSystem::numPhases, 0.0); + MaterialLaw::capillaryPressures(pc, matParams, fluidState); + double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; + double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_]; + return pcPhase - target_pc_; + } +private: + const MaterialLawManager& materialLawManager_; + const int phase_; + const int cell_; + const double target_pc_; +}; + +template +double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { + const auto& scaledDrainageInfo = + materialLawManager.oilWaterScaledEpsInfoDrainage(cell); + + // Find minimum and maximum saturations. + switch(phase) { + case FluidSystem::waterPhaseIdx : + { + return scaledDrainageInfo.Swl; + break; + } + case FluidSystem::gasPhaseIdx : + { + return scaledDrainageInfo.Sgl; + break; + } + case FluidSystem::oilPhaseIdx : + { + OPM_THROW(std::runtime_error, "Min saturation not implemented for oil phase."); + break; + } + default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); + } + return -1.0; +} + +template +double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { + const auto& scaledDrainageInfo = + materialLawManager.oilWaterScaledEpsInfoDrainage(cell); + + // Find minimum and maximum saturations. + switch(phase) { + case FluidSystem::waterPhaseIdx : + { + return scaledDrainageInfo.Swu; + break; + } + case FluidSystem::gasPhaseIdx : + { + return scaledDrainageInfo.Sgu; + break; + } + case FluidSystem::oilPhaseIdx : + { + OPM_THROW(std::runtime_error, "Max saturation not implemented for oil phase."); + break; + } + default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); + } + return -1.0; +} + + +/// Compute saturation of some phase corresponding to a given +/// capillary pressure. +template +inline double satFromPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double target_pc, + const bool increasing = false) +{ + // Find minimum and maximum saturations. + const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); + const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); + + // Create the equation f(s) = pc(s) - target_pc + const PcEq f(materialLawManager, phase, cell, target_pc); + const double f0 = f(s0); + const double f1 = f(s1); + + if (f0 <= 0.0) { + return s0; + } else if (f1 > 0.0) { + return s1; + } else { + const int max_iter = 60; + const double tol = 1e-6; + int iter_used = -1; + typedef Opm::RegulaFalsi ScalarSolver; + const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used); + return sol; + } +} + + +/// Functor for inverting a sum of capillary pressure functions. +/// Function represented is +/// f(s) = pc1(s) + pc2(1 - s) - target_pc +template +struct PcEqSum +{ + PcEqSum(const MaterialLawManager& materialLawManager, + const int phase1, + const int phase2, + const int cell, + const double target_pc) + : materialLawManager_(materialLawManager), + phase1_(phase1), + phase2_(phase2), + cell_(cell), + target_pc_(target_pc) + { + } + double operator()(double s) const + { + const auto& matParams = materialLawManager_.materialLawParams(cell_); + SatOnlyFluidState fluidState; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); + fluidState.setSaturation(phase1_, s); + fluidState.setSaturation(phase2_, 1.0 - s); + + double pc[FluidSystem::numPhases]; + std::fill(pc, pc + FluidSystem::numPhases, 0.0); + + 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_]; + return pc1 + pc2 - target_pc_; + } +private: + const MaterialLawManager& materialLawManager_; + const int phase1_; + const int phase2_; + const int cell_; + const double target_pc_; +}; + + + + +/// 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 +inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, + const int phase1, + const int phase2, + const int cell, + const double target_pc) +{ + // Find minimum and maximum saturations. + const double smin = minSaturations(materialLawManager, phase1, cell); + const double smax = maxSaturations(materialLawManager, phase1, cell); + + // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc + const PcEqSum f(materialLawManager, phase1, phase2, cell, target_pc); + const double f0 = f(smin); + const double f1 = f(smax); + if (f0 <= 0.0) { + return smin; + } else if (f1 > 0.0) { + return smax; + } else { + const int max_iter = 30; + const double tol = 1e-6; + int iter_used = -1; + typedef Opm::RegulaFalsi ScalarSolver; + const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used); + return sol; + } +} + +/// Compute saturation from depth. Used for constant capillary pressure function +template +inline double satFromDepth(const MaterialLawManager& materialLawManager, + const double cellDepth, + const double contactDepth, + const int phase, + const int cell, + const bool increasing = false) +{ + const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); + const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); + + if (cellDepth < contactDepth){ + return s0; + } else { + return s1; + } + +} + +/// Return true if capillary pressure function is constant +template +inline bool isConstPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell) +{ + // Create the equation f(s) = pc(s); + const PcEq f(materialLawManager, phase, cell, 0); + const double f0 = f(minSaturations(materialLawManager, phase, cell)); + const double f1 = f(maxSaturations(materialLawManager, phase, cell)); + return std::abs(f0 - f1) < std::numeric_limits::epsilon(); +} + +} // namespace Equil } // namespace Ewoms #endif // EWOMS_EQUILIBRATIONHELPERS_HEADER_INCLUDED