From 54cf35e8219feb656161d0e1a586c1c05f540b40 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 9 Aug 2022 15:32:19 +0200 Subject: [PATCH] equilibrationhelpers: make templates private --- CMakeLists_files.cmake | 1 + ebos/equil/equilibrationhelpers.cc | 565 +++++++++++++++++++++++++++++ ebos/equil/equilibrationhelpers.hh | 400 +++----------------- tests/test_equil.cc | 8 +- 4 files changed, 627 insertions(+), 347 deletions(-) create mode 100644 ebos/equil/equilibrationhelpers.cc diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 7f0fbb0ef..80b16b329 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -33,6 +33,7 @@ list (APPEND MAIN_SOURCE_FILES ebos/eclgenericwriter.cc ebos/eclinterregflows.cc ebos/ecltransmissibility.cc + ebos/equil/equilibrationhelpers.cc opm/core/props/BlackoilPhases.cpp opm/core/props/phaseUsageFromDeck.cpp opm/core/props/satfunc/RelpermDiagnostics.cpp diff --git a/ebos/equil/equilibrationhelpers.cc b/ebos/equil/equilibrationhelpers.cc new file mode 100644 index 000000000..be4b914a0 --- /dev/null +++ b/ebos/equil/equilibrationhelpers.cc @@ -0,0 +1,565 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include "config.h" + +#include + +#include + +#include +#include +#include + +#include + +namespace Opm { +namespace EQUIL { + +using FluidSystemSimple = BlackOilFluidSystem; + +// Adjust oil pressure according to gas saturation and cap pressure +using SatOnlyFluidState = SimpleModularFluidState; + +namespace Miscibility { + +template +RsVD::RsVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rs) + : pvtRegionIdx_(pvtRegionIdx) + , rsVsDepth_(depth, rs) +{ +} + +template +double RsVD:: +operator()(const double depth, + const double press, + const double temp, + const double satGas) const +{ + if (satGas > 0.0) { + return satRs(press, temp); + } + else { + if (rsVsDepth_.xMin() > depth) + return rsVsDepth_.valueAt(0); + else if (rsVsDepth_.xMax() < depth) + return rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1); + return std::min(satRs(press, temp), rsVsDepth_.eval(depth, /*extrapolate=*/false)); + } +} + +template +double RsVD::satRs(const double press, const double temp) const +{ + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); +} + +template +PBVD::PBVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& pbub) + : pvtRegionIdx_(pvtRegionIdx) + , pbubVsDepth_(depth, pbub) +{ +} + +template +double PBVD:: +operator()(const double depth, + const double cellPress, + const double temp, + const double satGas) const +{ + double press = cellPress; + if (satGas <= 0.0) { + if (pbubVsDepth_.xMin() > depth) + press = pbubVsDepth_.valueAt(0); + else if (pbubVsDepth_.xMax() < depth) + press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1); + else + press = pbubVsDepth_.eval(depth, /*extrapolate=*/false); + } + return satRs(std::min(press, cellPress), temp); +} + +template +double PBVD:: +satRs(const double press, const double temp) const +{ + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); +} + +template +PDVD::PDVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& pdew) + : pvtRegionIdx_(pvtRegionIdx) + , pdewVsDepth_(depth, pdew) +{ +} + +template +double PDVD:: +operator()(const double depth, + const double cellPress, + const double temp, + const double satOil) const +{ + double press = cellPress; + if (satOil <= 0.0) { + if (pdewVsDepth_.xMin() > depth) + press = pdewVsDepth_.valueAt(0); + else if (pdewVsDepth_.xMax() < depth) + press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1); + else + press = pdewVsDepth_.eval(depth, /*extrapolate=*/false); + } + return satRv(std::min(press, cellPress), temp); +} + +template +double PDVD:: +satRv(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); +} + + +template +RvVD::RvVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rv) + : pvtRegionIdx_(pvtRegionIdx) + , rvVsDepth_(depth, rv) +{ +} + +template +double RvVD:: +operator()(const double depth, + const double press, + const double temp, + const double satOil) const +{ + if (std::abs(satOil) > 1e-16) { + return satRv(press, temp); + } + else { + if (rvVsDepth_.xMin() > depth) + return rvVsDepth_.valueAt(0); + else if (rvVsDepth_.xMax() < depth) + return rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1); + return std::min(satRv(press, temp), rvVsDepth_.eval(depth, /*extrapolate=*/false)); + } +} + +template +double RvVD:: +satRv(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); +} + + +template +RsSatAtContact:: +RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) + : pvtRegionIdx_(pvtRegionIdx) +{ + rsSatContact_ = satRs(pContact, T_contact); +} + +template +double RsSatAtContact:: +operator()(const double /* depth */, + const double press, + const double temp, + const double satGas) const +{ + if (satGas > 0.0) { + return satRs(press, temp); + } + else { + return std::min(satRs(press, temp), rsSatContact_); + } +} + +template +double RsSatAtContact:: +satRs(const double press, const double temp) const +{ + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); +} + +template +RvSatAtContact:: +RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) + : pvtRegionIdx_(pvtRegionIdx) +{ + rvSatContact_ = satRv(pContact, T_contact); +} + +template +double RvSatAtContact:: +operator()(const double /*depth*/, + const double press, + const double temp, + const double satOil) const +{ + if (satOil > 0.0) { + return satRv(press, temp); + } + else { + return std::min(satRv(press, temp), rvSatContact_); + } +} + +template +double RvSatAtContact:: +satRv(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; +} + +} // namespace Miscibility + +EquilReg::EquilReg(const EquilRecord& rec, + std::shared_ptr rs, + std::shared_ptr rv, + const TabulatedFunction& saltVdTable, + const int pvtIdx) + : rec_ (rec) + , rs_ (rs) + , rv_ (rv) + , saltVdTable_ (saltVdTable) + , pvtIdx_ (pvtIdx) +{ +} + +double EquilReg::datum() const +{ + return this->rec_.datumDepth(); +} + +double EquilReg::pressure() const +{ + return this->rec_.datumDepthPressure(); +} + +double EquilReg::zwoc() const +{ + return this->rec_.waterOilContactDepth(); +} + +double EquilReg::pcowWoc() const +{ + return this->rec_.waterOilContactCapillaryPressure(); +} + +double EquilReg::zgoc() const +{ + return this->rec_.gasOilContactDepth(); +} + +double EquilReg::pcgoGoc() const +{ + return this->rec_.gasOilContactCapillaryPressure(); +} + +int EquilReg::equilibrationAccuracy() const +{ + return this->rec_.initializationTargetAccuracy(); +} + +const EquilReg::CalcDissolution& +EquilReg::dissolutionCalculator() const +{ + return *this->rs_; +} + +const EquilReg::CalcEvaporation& +EquilReg::evaporationCalculator() const +{ + return *this->rv_; +} + +const EquilReg::TabulatedFunction& +EquilReg::saltVdTable() const +{ + return saltVdTable_; +} + +int EquilReg::pvtIdx() const +{ + return this->pvtIdx_; +} + +template +PcEq:: +PcEq(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double targetPc) + : materialLawManager_(materialLawManager), + phase_(phase), + cell_(cell), + targetPc_(targetPc) +{ +} + +template +double PcEq:: +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); + 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_]; + return pcPhase - targetPc_; +} + +template +PcEqSum:: +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) +{ +} + +template +double PcEqSum:: +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); + + 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_]; + return pc1 + pc2 - targetPc_; +} + +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; + + case FluidSystem::gasPhaseIdx: + return scaledDrainageInfo.Sgl; + + case FluidSystem::oilPhaseIdx: + throw std::runtime_error("Min saturation not implemented for oil phase."); + + default: + 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; + + case FluidSystem::gasPhaseIdx: + return scaledDrainageInfo.Sgu; + + case FluidSystem::oilPhaseIdx: + throw std::runtime_error("Max saturation not implemented for oil phase."); + + default: + throw std::runtime_error("Unknown phaseIdx ."); + } + return -1.0; +} + +template +double satFromPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double targetPc, + const bool increasing) +{ + // Find minimum and maximum saturations. + double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); + double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); + + // Create the equation f(s) = pc(s) - targetPc + const PcEq f(materialLawManager, phase, cell, targetPc); + double f0 = f(s0); + double f1 = f(s1); + if (!std::isfinite(f0 + f1)) + throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1)); + + if (f0 <= 0.0) + return s0; + else if (f1 >= 0.0) + return s1; + + const double tol = 1e-10; + // should at least converge in 2 times bisection but some safety here: + const int maxIter = -2*static_cast(std::log2(tol)) + 10; + int usedIterations = -1; + const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); + return root; +} + +template +double satFromSumOfPcs(const MaterialLawManager& materialLawManager, + const int phase1, + const int phase2, + const int cell, + const double targetPc) +{ + // Find minimum and maximum saturations. + double s0 = minSaturations(materialLawManager, phase1, cell); + double s1 = maxSaturations(materialLawManager, phase1, cell); + + // Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc + const PcEqSum f(materialLawManager, phase1, phase2, cell, targetPc); + double f0 = f(s0); + double 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; + // should at least converge in 2 times bisection but some safety here: + const int maxIter = -2*static_cast(std::log2(tol)) + 10; + int usedIterations = -1; + const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); + return root; +} + +template +double satFromDepth(const MaterialLawManager& materialLawManager, + const double cellDepth, + const double contactDepth, + const int phase, + const int cell, + const bool increasing) +{ + 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; + } +} + +template +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(); +} + +using MatLaw = EclMaterialLawManager>; +using FS = BlackOilFluidSystem; + +template struct PcEq; + +template double satFromPc(const MatLaw&,const int,const int, + const double,const bool); +template double satFromSumOfPcs(const MatLaw&,const int,const int, + const int,const double); +template double satFromDepth(const MatLaw&,const double,const double, + const int,const int,const bool); +template bool isConstPc(const MatLaw&,const int,const int); + +namespace Miscibility { + template class PBVD; + template class PDVD; + template class RsVD; + template class RsSatAtContact; + template class RvSatAtContact; + template class RvVD; +} + +} // namespace Equil +} // namespace Opm diff --git a/ebos/equil/equilibrationhelpers.hh b/ebos/equil/equilibrationhelpers.hh index e57c3724b..c22098ef8 100644 --- a/ebos/equil/equilibrationhelpers.hh +++ b/ebos/equil/equilibrationhelpers.hh @@ -29,19 +29,13 @@ #ifndef EWOMS_EQUILIBRATIONHELPERS_HH #define EWOMS_EQUILIBRATIONHELPERS_HH -#include - #include -#include -#include -#include #include -#include #include #include - +#include /* ---- synopsis of EquilibrationHelpers.hpp ---- @@ -62,16 +56,16 @@ class EquilReg; - template + template struct PcEq; - template + template double satFromPc(const MaterialLawManager& materialLawManager, const int phase, const int cell, const double targetPc, const bool increasing = false) - template + template double satFromSumOfPcs(const MaterialLawManager& materialLawManager, const int phase1, const int phase2, @@ -82,7 +76,9 @@ ---- end of synopsis of EquilibrationHelpers.hpp ---- */ + namespace Opm { + /** * Types and routines that collectively implement a basic * ECLIPSE-style equilibration-based initialisation scheme. @@ -92,23 +88,6 @@ namespace Opm { */ namespace EQUIL { - -using FluidSystemSimple = BlackOilFluidSystem; - -// Adjust oil pressure according to gas saturation and cap pressure -using SatOnlyFluidState = SimpleModularFluidState; - /** * Types and routines relating to phase mixing in * equilibration calculations. @@ -198,10 +177,7 @@ public: */ RsVD(const int pvtRegionIdx, const std::vector& depth, - const std::vector& rs) - : pvtRegionIdx_(pvtRegionIdx) - , rsVsDepth_(depth, rs) - {} + const std::vector& rs); virtual ~RsVD() = default; @@ -223,19 +199,7 @@ public: double operator()(const double depth, const double press, const double temp, - const double satGas = 0.0) const - { - if (satGas > 0.0) { - return satRs(press, temp); - } - else { - if (rsVsDepth_.xMin() > depth) - return rsVsDepth_.valueAt(0); - else if (rsVsDepth_.xMax() < depth) - return rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1); - return std::min(satRs(press, temp), rsVsDepth_.eval(depth, /*extrapolate=*/false)); - } - } + const double satGas = 0.0) const; private: using RsVsDepthFunc = Tabulated1DFunction; @@ -243,12 +207,10 @@ private: const int pvtRegionIdx_; RsVsDepthFunc rsVsDepth_; - double satRs(const double press, const double temp) const - { - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); - } + double satRs(const double press, const double temp) const; }; + /** * Type that implements "dissolved gas-oil ratio" * tabulated as a function of depth policy. Data @@ -267,10 +229,7 @@ public: */ PBVD(const int pvtRegionIdx, const std::vector& depth, - const std::vector& pbub) - : pvtRegionIdx_(pvtRegionIdx) - , pbubVsDepth_(depth, pbub) - {} + const std::vector& pbub); virtual ~PBVD() = default; @@ -291,19 +250,7 @@ public: double operator()(const double depth, const double cellPress, const double temp, - const double satGas = 0.0) const - { - double press = cellPress; - if (satGas <= 0.0) { - if (pbubVsDepth_.xMin() > depth) - press = pbubVsDepth_.valueAt(0); - else if (pbubVsDepth_.xMax() < depth) - press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1); - else - press = pbubVsDepth_.eval(depth, /*extrapolate=*/false); - } - return satRs(std::min(press, cellPress), temp); - } + const double satGas = 0.0) const; private: using PbubVsDepthFunc = Tabulated1DFunction; @@ -311,12 +258,10 @@ private: const int pvtRegionIdx_; PbubVsDepthFunc pbubVsDepth_; - double satRs(const double press, const double temp) const - { - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); - } + double satRs(const double press, const double temp) const; }; + /** * Type that implements "vaporized oil-gas ratio" * tabulated as a function of depth policy. Data @@ -335,10 +280,7 @@ public: */ PDVD(const int pvtRegionIdx, const std::vector& depth, - const std::vector& pdew) - : pvtRegionIdx_(pvtRegionIdx) - , pdewVsDepth_(depth, pdew) - {} + const std::vector& pdew); virtual ~PDVD() = default; @@ -359,19 +301,7 @@ public: double operator()(const double depth, const double cellPress, const double temp, - const double satOil = 0.0) const - { - double press = cellPress; - if (satOil <= 0.0) { - if (pdewVsDepth_.xMin() > depth) - press = pdewVsDepth_.valueAt(0); - else if (pdewVsDepth_.xMax() < depth) - press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1); - else - press = pdewVsDepth_.eval(depth, /*extrapolate=*/false); - } - return satRv(std::min(press, cellPress), temp); - } + const double satOil = 0.0) const; private: using PdewVsDepthFunc = Tabulated1DFunction; @@ -379,10 +309,7 @@ private: const int pvtRegionIdx_; PdewVsDepthFunc pdewVsDepth_; - double satRv(const double press, const double temp) const - { - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); - } + double satRv(const double press, const double temp) const; }; @@ -404,10 +331,7 @@ public: */ RvVD(const int pvtRegionIdx, const std::vector& depth, - const std::vector& rv) - : pvtRegionIdx_(pvtRegionIdx) - , rvVsDepth_(depth, rv) - {} + const std::vector& rv); /** * Function call. @@ -427,19 +351,7 @@ public: double operator()(const double depth, const double press, const double temp, - const double satOil = 0.0) const - { - if (std::abs(satOil) > 1e-16) { - return satRv(press, temp); - } - else { - if (rvVsDepth_.xMin() > depth) - return rvVsDepth_.valueAt(0); - else if (rvVsDepth_.xMax() < depth) - return rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1); - return std::min(satRv(press, temp), rvVsDepth_.eval(depth, /*extrapolate=*/false)); - } - } + const double satOil = 0.0) const; private: using RvVsDepthFunc = Tabulated1DFunction; @@ -447,10 +359,7 @@ private: const int pvtRegionIdx_; RvVsDepthFunc rvVsDepth_; - double satRv(const double press, const double temp) const - { - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); - } + double satRv(const double press, const double temp) const; }; @@ -479,11 +388,7 @@ 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) - : pvtRegionIdx_(pvtRegionIdx) - { - rsSatContact_ = satRs(pContact, T_contact); - } + RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact); /** * Function call. @@ -503,24 +408,13 @@ public: double operator()(const double /* depth */, const double press, const double temp, - const double satGas = 0.0) const - { - if (satGas > 0.0) { - return satRs(press, temp); - } - else { - return std::min(satRs(press, temp), rsSatContact_); - } - } + const double satGas = 0.0) const; private: const int pvtRegionIdx_; double rsSatContact_; - double satRs(const double press, const double temp) const - { - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); - } + double satRs(const double press, const double temp) const; }; @@ -549,11 +443,7 @@ 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) - :pvtRegionIdx_(pvtRegionIdx) - { - rvSatContact_ = satRv(pContact, T_contact); - } + RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact); /** * Function call. @@ -573,24 +463,13 @@ public: double operator()(const double /*depth*/, const double press, const double temp, - const double satOil = 0.0) const - { - if (satOil > 0.0) { - return satRv(press, temp); - } - else { - return std::min(satRv(press, temp), rvSatContact_); - } - } + const double satOil = 0.0) const; private: const int pvtRegionIdx_; double rvSatContact_; - double satRv(const double press, const double temp) const - { - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; - } + double satRv(const double press, const double temp) const; }; } // namespace Miscibility @@ -631,57 +510,51 @@ public: std::shared_ptr rs, std::shared_ptr rv, const TabulatedFunction& saltVdTable, - const int pvtIdx) - : rec_ (rec) - , rs_ (rs) - , rv_ (rv) - , saltVdTable_ (saltVdTable) - , pvtIdx_ (pvtIdx) - {} + const int pvtIdx); /** * Type of dissolved gas-oil ratio calculator. */ - typedef Miscibility::RsFunction CalcDissolution; + using CalcDissolution = Miscibility::RsFunction; /** * Type of vapourised oil-gas ratio calculator. */ - typedef Miscibility::RsFunction CalcEvaporation; + using CalcEvaporation = Miscibility::RsFunction; /** * Datum depth in current region */ - double datum() const { return this->rec_.datumDepth(); } + double datum() const; /** * Pressure at datum depth in current region. */ - double pressure() const { return this->rec_.datumDepthPressure(); } + double pressure() const; /** * Depth of water-oil contact. */ - double zwoc() const { return this->rec_.waterOilContactDepth(); } + double zwoc() const; /** * water-oil capillary pressure at water-oil contact. * * \return P_o - P_w at WOC. */ - double pcowWoc() const { return this->rec_.waterOilContactCapillaryPressure(); } + double pcowWoc() const; /** * Depth of gas-oil contact. */ - double zgoc() const { return this->rec_.gasOilContactDepth(); } + double zgoc() const; /** * Gas-oil capillary pressure at gas-oil contact. * * \return P_g - P_o at GOC. */ - double pcgoGoc() const { return this->rec_.gasOilContactCapillaryPressure(); } + double pcgoGoc() const; /** * Accuracy/strategy for initial fluid-in-place calculation. @@ -690,31 +563,25 @@ public: * horizontal subdivision method with 2*(-N) intervals, and positive * (N>0) for the tilted subdivision method with 2*N intervals. */ - int equilibrationAccuracy() const - { - return this->rec_.initializationTargetAccuracy(); - } + int equilibrationAccuracy() const; /** * Retrieve dissolved gas-oil ratio calculator of current * region. */ - const CalcDissolution& - dissolutionCalculator() const { return *this->rs_; } + const CalcDissolution& dissolutionCalculator() const; /** * Retrieve vapourised oil-gas ratio calculator of current * region. */ - const CalcEvaporation& - evaporationCalculator() const { return *this->rv_; } + const CalcEvaporation& evaporationCalculator() const; - const TabulatedFunction& - saltVdTable() const { return saltVdTable_;} + const TabulatedFunction& saltVdTable() const; /** * Retrieve pvtIdx of the region. */ - int pvtIdx() const { return this->pvtIdx_; } + int pvtIdx() const; private: EquilRecord rec_; /**< Equilibration data */ @@ -729,35 +596,16 @@ private: /// Functor for inverting capillary pressure function. /// Function represented is /// f(s) = pc(s) - targetPc -template +template struct PcEq { PcEq(const MaterialLawManager& materialLawManager, const int phase, const int cell, - const double targetPc) - : materialLawManager_(materialLawManager), - phase_(phase), - cell_(cell), - targetPc_(targetPc) - {} + const double targetPc); - 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 operator()(double s) const; - 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 - targetPc_; - } private: const MaterialLawManager& materialLawManager_; const int phase_; @@ -766,124 +614,36 @@ private: }; 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; - - case FluidSystem::gasPhaseIdx: - return scaledDrainageInfo.Sgl; - - case FluidSystem::oilPhaseIdx: - throw std::runtime_error("Min saturation not implemented for oil phase."); - - default: - throw std::runtime_error("Unknown phaseIdx ."); - } - return -1.0; -} +double minSaturations(const MaterialLawManager& materialLawManager, + const int phase, const int cell); 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; - - case FluidSystem::gasPhaseIdx: - return scaledDrainageInfo.Sgu; - - case FluidSystem::oilPhaseIdx: - throw std::runtime_error("Max saturation not implemented for oil phase."); - - default: - throw std::runtime_error("Unknown phaseIdx ."); - } - return -1.0; -} - +double maxSaturations(const MaterialLawManager& materialLawManager, + const int phase, const int cell); /// Compute saturation of some phase corresponding to a given /// capillary pressure. -template +template double satFromPc(const MaterialLawManager& materialLawManager, const int phase, const int cell, const double targetPc, - const bool increasing = false) -{ - // Find minimum and maximum saturations. - double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); - double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); - - // Create the equation f(s) = pc(s) - targetPc - const PcEq f(materialLawManager, phase, cell, targetPc); - double f0 = f(s0); - double f1 = f(s1); - if (!std::isfinite(f0 + f1)) - throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1)); - - if (f0 <= 0.0) - return s0; - else if (f1 >= 0.0) - return s1; - - const double tol = 1e-10; - // should at least converge in 2 times bisection but some safety here: - const int maxIter = -2*static_cast(std::log2(tol)) + 10; - int usedIterations = -1; - const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); - return root; -} - + const bool increasing = false); /// Functor for inverting a sum of capillary pressure functions. /// Function represented is /// f(s) = pc1(s) + pc2(1 - s) - targetPc -template +template struct PcEqSum { 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 double targetPc); - 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 operator()(double s) const; - 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 - targetPc_; - } private: const MaterialLawManager& materialLawManager_; const int phase1_; @@ -892,74 +652,30 @@ private: const double 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 +template double satFromSumOfPcs(const MaterialLawManager& materialLawManager, const int phase1, const int phase2, const int cell, - const double targetPc) -{ - // Find minimum and maximum saturations. - double s0 = minSaturations(materialLawManager, phase1, cell); - double s1 = maxSaturations(materialLawManager, phase1, cell); - - // Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc - const PcEqSum f(materialLawManager, phase1, phase2, cell, targetPc); - double f0 = f(s0); - double 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; - // should at least converge in 2 times bisection but some safety here: - const int maxIter = -2*static_cast(std::log2(tol)) + 10; - int usedIterations = -1; - const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); - return root; -} + const double targetPc); /// Compute saturation from depth. Used for constant capillary pressure function -template +template 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; - } - -} + const bool increasing = false); /// Return true if capillary pressure function is constant -template +template 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(); -} + const int cell); } // namespace Equil } // namespace Opm diff --git a/tests/test_equil.cc b/tests/test_equil.cc index 749667886..83177bad2 100644 --- a/tests/test_equil.cc +++ b/tests/test_equil.cc @@ -537,8 +537,6 @@ BOOST_AUTO_TEST_CASE(CapillaryInversion) // Test setup. using TypeTag = Opm::Properties::TTag::TestEquilTypeTag; using FluidSystem = Opm::GetPropType; - using MaterialLaw = Opm::GetPropType; - using MaterialLawManager = typename Opm::GetProp::EclMaterialLawManager; auto simulator = initSimulator("equil_capillary.DATA"); @@ -552,7 +550,7 @@ BOOST_AUTO_TEST_CASE(CapillaryInversion) const std::vector s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 }; BOOST_REQUIRE_EQUAL(pc.size(), s.size()); for (size_t i = 0; i < pc.size(); ++i) { - const double s_computed = Opm::EQUIL::satFromPc(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing); + const double s_computed = Opm::EQUIL::satFromPc(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing); BOOST_CHECK_CLOSE(s_computed, s[i], reltol); } } @@ -565,7 +563,7 @@ BOOST_AUTO_TEST_CASE(CapillaryInversion) const std::vector s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 }; BOOST_REQUIRE_EQUAL(pc.size(), s.size()); for (size_t i = 0; i < pc.size(); ++i) { - const double s_computed = Opm::EQUIL::satFromPc(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing); + const double s_computed = Opm::EQUIL::satFromPc(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing); BOOST_CHECK_CLOSE(s_computed, s[i], reltol); } } @@ -578,7 +576,7 @@ BOOST_AUTO_TEST_CASE(CapillaryInversion) const std::vector s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 }; BOOST_REQUIRE_EQUAL(pc.size(), s.size()); for (size_t i = 0; i < pc.size(); ++i) { - const double s_computed = Opm::EQUIL::satFromSumOfPcs(*simulator->problem().materialLawManager(), water, gas, cell, pc[i]); + const double s_computed = Opm::EQUIL::satFromSumOfPcs(*simulator->problem().materialLawManager(), water, gas, cell, pc[i]); BOOST_CHECK_CLOSE(s_computed, s[i], reltol); } }