From f0e2a1efe0c8a2854cf3765a3378749a49567b8d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 12 Oct 2022 23:01:42 +0200 Subject: [PATCH 1/2] Make Reservoir Voidage Rate Calculator Independent This commit adds a new overload of the calcReservoirVoidageRates() member function. The new overload does not depend on the user calling defineState(), but instead takes state variables as direct input arguments. Reimplement the state-dependent overload in terms of the new function. The immediate use case is calculating additional dynamic quantities per segment in a multisegmented well model, mostly for reporting purposes. --- opm/simulators/wells/RateConverter.hpp | 208 ++++++++++++++++--------- 1 file changed, 131 insertions(+), 77 deletions(-) diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 329b5b3a0..9704d2cb7 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -24,10 +24,12 @@ #include #include -#include #include +#include + #include #include + #include #include #include @@ -36,6 +38,7 @@ #include #include #include + /** * \file * Facility for converting component rates at surface conditions to @@ -56,13 +59,13 @@ namespace Opm { * The conversion uses fluid properties evaluated at average * hydrocarbon pressure in regions or field. * - * \tparam FluidSystem Fluid system class. Expected to be a BlackOilFluidSystem + * \tparam FluidSystem Fluid system class. Expected to be a + * BlackOilFluidSystem * - * \tparam Region Type of a forward region mapping. Expected - * to provide indexed access through \code operator[]() - * \endcode as well as inner types \c value_type, \c - * size_type, and \c const_iterator. Typically \code - * std::vector \endcode. + * \tparam Region Type of a forward region mapping. Expected to + * provide indexed access through \code operator[]() \endcode as + * well as inner types \c value_type, \c size_type, and \c + * const_iterator. Typically \code std::vector \endcode. */ template class SurfaceToReservoirVoidage { @@ -70,18 +73,15 @@ namespace Opm { /** * Constructor. * - * \param[in] region Forward region mapping. Often - * corresponds to the "FIPNUM" mapping of an ECLIPSE input - * deck. + * \param[in] region Forward region mapping. Often corresponds + * to the "FIPNUM" mapping of an ECLIPSE input deck. */ SurfaceToReservoirVoidage(const PhaseUsage& phaseUsage, - const Region& region) + const Region& region) : phaseUsage_(phaseUsage) - , rmap_ (region) - , attr_ (rmap_, Attributes()) - { - } - + , rmap_ (region) + , attr_ (rmap_, Attributes()) + {} /** * Compute pore volume averaged hydrocarbon state pressure, rs and rv. @@ -94,9 +94,8 @@ namespace Opm { template void defineState(const EbosSimulator& simulator) { - - // create map from cell to region - // and set all attributes to zero + // create map from cell to region and set all attributes to + // zero for (const auto& reg : rmap_.activeRegions()) { auto& ra = attr_.attributes(reg); ra.pressure = 0.0; @@ -355,103 +354,134 @@ namespace Opm { } } + /** + * Convert surface volume flow rates to reservoir voidage flow + * rates. + * + * State dependent version. Client must call \code + * defineState() \endcode prior to invoking this member + * function. + * + * \tparam Rates Type representing contiguous collection of + * surface flow rates. Must support direct indexing through + * \code operator[]() \endcode. + * + * \param[in] r Zero based fluid-in-place region index. + * + * \param[in] pvtRegionIdx Zero based PVT region index. + * + * \param[in] surface_rates surface volume flow rates for all + * active phases. + * + * \param[out] voidage_rates reservoir volume flow rates for all + * active phases. + */ + template + void calcReservoirVoidageRates(const RegionId r, + const int pvtRegionIdx, + const Rates& surface_rates, + Rates& voidage_rates) const + { + const auto& ra = this->attr_.attributes(r); + + this->calcReservoirVoidageRates(pvtRegionIdx, + ra.pressure, ra.rs, ra.rv, + ra.temperature, + ra.saltConcentration, + surface_rates, + voidage_rates); + } /** - * Converting surface volume rates to reservoir voidage rates + * Convert surface volume flow rates to reservoir voidage flow + * rates. * - * \tparam Rates Type representing contiguous collection - * of surface-to-reservoir conversion coefficients. Must - * support direct indexing through \code operator[]() - * \endcode. + * State independent version. * + * \tparam Rates Type representing contiguous collection of + * surface flow rates. Must support direct indexing through + * \code operator[]() \endcode. * - * \param[in] r Fluid-in-place region of the well - * \param[in] pvtRegionIdx PVT region of the well - * \param[in] surface_rates surface voluem rates for - * all active phases + * \param[in] pvtRegionIdx PVT region. * - * \param[out] voidage_rates reservoir volume rates for - * all active phases + * \param[in] p Fluid pressure. + * + * \param[in] rs Dissolved gas/oil ratio. + * + * \param[in] rv Vaporised oil/gas ratio. + * + * \param[in] T Temperature. Unused in non-thermal simulation + * runs. + * + * \param[in] saltConcentration Salt concentration. Unused in + * simulation runs without salt precipitation. + * + * \param[in] surface_rates Surface volume flow rates for all + * active phases. + * + * \param[out] voidage_rates Reservoir volume flow rates for all + * active phases. */ - template - void - calcReservoirVoidageRates(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates, - Rates& voidage_rates) const + template + void calcReservoirVoidageRates(const int pvtRegionIdx, + const double p, + const double rs, + const double rv, + const double T, + const double saltConcentration, + const Rates& surface_rates, + Rates& voidage_rates) const { - assert(voidage_rates.size() == surface_rates.size()); - std::fill(voidage_rates.begin(), voidage_rates.end(), 0.0); - const auto& pu = phaseUsage_; - const auto& ra = attr_.attributes(r); - const double p = ra.pressure; - const double T = ra.temperature; - const double saltConcentration = ra.saltConcentration; + const auto& pu = this->phaseUsage_; + const auto iw = RegionAttributeHelpers::PhasePos::water(pu); + const auto io = RegionAttributeHelpers::PhasePos::oil (pu); + const auto ig = RegionAttributeHelpers::PhasePos::gas (pu); - const int iw = RegionAttributeHelpers::PhasePos::water(pu); - const int io = RegionAttributeHelpers::PhasePos::oil (pu); - const int ig = RegionAttributeHelpers::PhasePos::gas (pu); + const auto [Rs, Rv] = this-> + dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates); if (RegionAttributeHelpers::PhaseUsed::water(pu)) { // q[w]_r = q[w]_s / bw - - const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration); + const auto bw = FluidSystem::waterPvt() + .inverseFormationVolumeFactor(pvtRegionIdx, T, p, + saltConcentration); voidage_rates[iw] = surface_rates[iw] / bw; } - // Use average Rs and Rv: - auto a = ra.rs; - auto b = a; - if (io >= 0 && ig >= 0) { - b = surface_rates[ig]/(surface_rates[io]+1.0e-15); - } - - double Rs = std::min(a, b); - - a = ra.rv; - b = a; - if (io >= 0 && ig >= 0) { - b = surface_rates[io]/(surface_rates[ig]+1.0e-15); - } - - double Rv = std::min(a, b); - // Determinant of 'R' matrix - const double detR = 1.0 - (Rs * Rv); + const auto detR = 1.0 - (Rs * Rv); if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) - - const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); - const double den = bo * detR; - voidage_rates[io] = surface_rates[io]; - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { voidage_rates[io] -= Rv * surface_rates[ig]; } - voidage_rates[io] /= den; + const auto bo = FluidSystem::oilPvt() + .inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); + + voidage_rates[io] /= bo * detR; } if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) - - const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, 0.0 /*=Rvw*/); - const double den = bg * detR; - voidage_rates[ig] = surface_rates[ig]; - if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { voidage_rates[ig] -= Rs * surface_rates[io]; } - voidage_rates[ig] /= den; + const auto bg = FluidSystem::gasPvt() + .inverseFormationVolumeFactor(pvtRegionIdx, T, p, + Rv, 0.0 /*=Rvw*/); + + voidage_rates[ig] /= bg * detR; } } - /** * Compute coefficients for surface-to-reservoir voidage * conversion for solvent. @@ -510,7 +540,31 @@ namespace Opm { RegionAttributeHelpers::RegionAttributes attr_; + template + std::pair + dissolvedVaporisedRatio(const int io, + const int ig, + const double rs, + const double rv, + const Rates& surface_rates) const + { + if ((io < 0) || (ig < 0)) { + return { rs, rv }; + } + + auto eps = std::copysign(1.0e-15, surface_rates[io]); + const auto Rs = surface_rates[ig] / (surface_rates[io] + eps); + + eps = std::copysign(1.0e-15, surface_rates[ig]); + const auto Rv = surface_rates[io] / (surface_rates[ig] + eps); + + return { + std::clamp(static_cast(Rs), 0.0, rs), + std::clamp(static_cast(Rv), 0.0, rv) + }; + } }; + } // namespace RateConverter } // namespace Opm From d388ce40e587f9b070d506563a62525f2efc18e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 13 Oct 2022 16:19:33 +0200 Subject: [PATCH 2/2] Revert to Original Rs/Rv Calculation Scheme for RESV Rates Keep the alternative approach in place for reference and experimental purposes, but disabled by a macro. --- opm/simulators/wells/RateConverter.hpp | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 9704d2cb7..91ac14afe 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -552,6 +552,26 @@ namespace Opm { return { rs, rv }; } +#define BURN_RESV_BRIDGES 0 + +#if !BURN_RESV_BRIDGES + + auto b = rs; + if (io >= 0 && ig >= 0) { + b = surface_rates[ig] / (surface_rates[io] + 1.0e-15); + } + const double Rs = std::min(rs, b); + + b = rv; + if (io >= 0 && ig >= 0) { + b = surface_rates[io] / (surface_rates[ig] + 1.0e-15); + } + const double Rv = std::min(rv, b); + + return { Rs, Rv }; + +#else // BURN_RESV_BRIDGES + auto eps = std::copysign(1.0e-15, surface_rates[io]); const auto Rs = surface_rates[ig] / (surface_rates[io] + eps); @@ -562,6 +582,10 @@ namespace Opm { std::clamp(static_cast(Rs), 0.0, rs), std::clamp(static_cast(Rv), 0.0, rv) }; + +#endif // BURN_RESV_BRIDGES + +#undef BURN_RESV_BRIDGES } };