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

@@ -24,10 +24,12 @@
 
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/grid/utility/RegionMapping.hpp>
-#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
+#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
+
 #include <dune/grid/common/gridenums.hh>
 #include <dune/grid/common/rangegenerators.hh>
+
 #include <algorithm>
 #include <cmath>
 #include <memory>
@@ -36,6 +38,7 @@
 #include <unordered_map>
 #include <utility>
 #include <vector>
+
 /**
  * \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
+         * \tparam Region Type of a forward region mapping.  Expected to
-         * to provide indexed access through \code operator[]()
+         *    provide indexed access through \code operator[]() \endcode as
-         * \endcode as well as inner types \c value_type, \c
+         *    well as inner types \c value_type, \c size_type, and \c
-         * size_type, and \c const_iterator.  Typically \code
+         *    const_iterator.  Typically \code std::vector<int> \endcode.
-         * std::vector<int> \endcode.
          */
         template <class FluidSystem, class Region>
         class SurfaceToReservoirVoidage {
@@ -70,18 +73,15 @@ namespace Opm {
             /**
              * Constructor.
              *
-             * \param[in] region Forward region mapping.  Often
+             * \param[in] region Forward region mapping.  Often corresponds
-             * corresponds to the "FIPNUM" mapping of an ECLIPSE input
+             * to the "FIPNUM" mapping of an ECLIPSE input deck.
-             * deck.
              */
             SurfaceToReservoirVoidage(const PhaseUsage& phaseUsage,
-                                      const Region&   region)
+                                      const Region&     region)
                 : phaseUsage_(phaseUsage)
-                , rmap_ (region)
+                , rmap_      (region)
-                , attr_ (rmap_, Attributes())
+                , attr_      (rmap_, Attributes())
-            {
+            {}
-            }
-
 
             /**
              * Compute pore volume averaged hydrocarbon state pressure, rs and rv.
@@ -94,9 +94,8 @@ namespace Opm {
             template <typename ElementContext, class EbosSimulator>
             void defineState(const EbosSimulator& simulator)
             {
-
+                // create map from cell to region and set all attributes to
-                // create map from cell to region
+                // zero
-                // 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 <class Rates>
+            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
+             * State independent version.
-             * of surface-to-reservoir conversion coefficients.  Must
-             * support direct indexing through \code operator[]()
-             * \endcode.
              *
+             * \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.
-             * \param[in] pvtRegionIdx PVT region of the well
-             * \param[in] surface_rates surface voluem rates for
-             * all active phases
              *
-             * \param[out] voidage_rates reservoir volume rates for
+             * \param[in] p Fluid pressure.
-             * all active phases
+             *
+             * \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 <class Rates >
+            template <class Rates>
-            void
+            void calcReservoirVoidageRates(const int    pvtRegionIdx,
-            calcReservoirVoidageRates(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates,
+                                           const double p,
-                                      Rates& voidage_rates) const
+                                           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& pu = this->phaseUsage_;
-                const auto& ra = attr_.attributes(r);
+                const auto  iw = RegionAttributeHelpers::PhasePos::water(pu);
-                const double p = ra.pressure;
+                const auto  io = RegionAttributeHelpers::PhasePos::oil  (pu);
-                const double T = ra.temperature;
+                const auto  ig = RegionAttributeHelpers::PhasePos::gas  (pu);
-                const double saltConcentration = ra.saltConcentration;
 
-                const int   iw = RegionAttributeHelpers::PhasePos::water(pu);
+                const auto [Rs, Rv] = this->
-                const int   io = RegionAttributeHelpers::PhasePos::oil  (pu);
+                    dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates);
-                const int   ig = RegionAttributeHelpers::PhasePos::gas  (pu);
 
                 if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
                     // q[w]_r = q[w]_s / bw
-
+                    const auto bw = FluidSystem::waterPvt()
-                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
+                        .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<RegionId, Attributes> attr_;
 
+            template <typename Rates>
+            std::pair<double, double>
+            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<double>(Rs), 0.0, rs),
+                    std::clamp(static_cast<double>(Rv), 0.0, rv)
+                };
+            }
         };
+
     } // namespace RateConverter
 } // namespace Opm