Merge pull request #4512 from akva2/rateconverter_compile_unit

RateConverter: introduce a compile unit

CMakeLists_files.cmake

@@ -105,6 +105,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/MultisegmentWellSegments.cpp
   opm/simulators/wells/ParallelWellInfo.cpp
   opm/simulators/wells/PerfData.cpp
+  opm/simulators/wells/RateConverter.cpp
   opm/simulators/wells/SegmentState.cpp
   opm/simulators/wells/SingleWellState.cpp
   opm/simulators/wells/StandardWellAssemble.cpp

opm/simulators/wells/RateConverter.cpp

@@ -0,0 +1,379 @@
+/*
+  Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014, 2015 Statoil ASA.
+  Copyright 2017, IRIS
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+#include <config.h>
+#include <opm/simulators/wells/RateConverter.hpp>
+
+#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+#include <string>
+
+namespace {
+
+template <typename Rates>
+std::pair<double, double>
+dissolvedVaporisedRatio(const int    io,
+                        const int    ig,
+                        const double rs,
+                        const double rv,
+                        const Rates& surface_rates)
+{
+    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 Opm {
+namespace RateConverter {
+
+template <class FluidSystem, class Region>
+void SurfaceToReservoirVoidage<FluidSystem,Region>::
+sumRates(std::unordered_map<RegionId,Attributes>& attributes_hpv,
+         std::unordered_map<RegionId,Attributes>& attributes_pv,
+         Parallel::Communication comm)
+{
+    for (const auto& reg : rmap_.activeRegions()) {
+        // Calculating first using the hydrocarbon pore volumes
+        auto& ra = attr_.attributes(reg);
+        const auto& attri_hpv = attributes_hpv[reg];
+        ra.data = attri_hpv.data;
+        comm.sum(ra.data.data(), ra.data.size());
+        // TODO: should we have some epsilon here instead of zero?
+        // No hydrocarbon pore volume, redo but this time using full pore volumes.
+        if (ra.pv <= 0.) {
+            // using the pore volume to do the averaging
+            const auto& attri_pv = attributes_pv[reg];
+            ra.data = attri_pv.data;
+            comm.sum(ra.data.data(), ra.data.size());
+            assert(ra.pv > 0.);
+        }
+        const double pv_sum = ra.pv;
+        for (double& d : ra.data)
+            d /= pv_sum;
+        ra.pv = pv_sum;
+    }
+}
+
+template <class FluidSystem, class Region>
+template <class Coeff>
+void SurfaceToReservoirVoidage<FluidSystem,Region>::
+calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
+{
+    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 int   iw = RegionAttributeHelpers::PhasePos::water(pu);
+    const int   io = RegionAttributeHelpers::PhasePos::oil  (pu);
+    const int   ig = RegionAttributeHelpers::PhasePos::gas  (pu);
+
+    std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
+
+    if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
+        // q[w]_r = q[w]_s / bw
+
+        const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, saltConcentration);
+
+        coeff[iw] = 1.0 / bw;
+    }
+
+    if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
+        const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
+        coeff[io] += 1.0 / bo;
+    }
+
+    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+        const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, 0.0);
+        coeff[ig] += 1.0 / bg;
+    }
+}
+
+template <class FluidSystem, class Region>
+template <class Coeff>
+void SurfaceToReservoirVoidage<FluidSystem,Region>::
+calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
+{
+    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 int   iw = RegionAttributeHelpers::PhasePos::water(pu);
+    const int   io = RegionAttributeHelpers::PhasePos::oil  (pu);
+    const int   ig = RegionAttributeHelpers::PhasePos::gas  (pu);
+
+    std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
+
+    // Actual Rsw and Rvw:
+    double Rsw = ra.rsw;
+    double Rvw = ra.rvw;
+    // Determinant of 'R' matrix
+    const double detRw = 1.0 - (Rsw * Rvw);
+
+    if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
+        // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s)
+
+        const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rsw, saltConcentration);
+
+        const double den = bw * detRw;
+
+        coeff[iw] += 1.0 / den;
+
+        if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+            coeff[ig] -= ra.rvw / den;
+        }
+    }
+
+    // Actual Rs and Rv:
+    double Rs = ra.rs;
+    double Rv = ra.rv;
+
+    // Determinant of 'R' matrix
+    const double detR = 1.0 - (Rs * Rv);
+
+    // Currently we only support either gas in water or gas in oil
+    // not both
+    if (detR != 1 && detRw != 1) {
+        std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR);
+        throw std::range_error(msg);
+    }
+
+    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;
+
+        coeff[io] += 1.0 / den;
+
+        if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+            coeff[ig] -= ra.rv / den;
+        }
+    }
+
+    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, Rvw);
+        if (FluidSystem::enableDissolvedGasInWater()) {
+            const double denw = bg * detRw;
+            coeff[ig] += 1.0 / denw;
+            if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
+               coeff[iw] -= ra.rsw / denw;
+            }
+        } else {
+            const double den = bg * detR;
+            coeff[ig] += 1.0 / den;
+            if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
+                coeff[io] -= ra.rs / den;
+            }
+        }
+    }
+}
+
+template <class FluidSystem, class Region>
+template <typename SurfaceRates, typename VoidageRates>
+void SurfaceToReservoirVoidage<FluidSystem,Region>::
+calcReservoirVoidageRates(const int           pvtRegionIdx,
+                          const double        p,
+                          const double        rs,
+                          const double        rv,
+                          const double        rsw,
+                          const double        rvw,
+                          const double        T,
+                          const double        saltConcentration,
+                          const SurfaceRates& surface_rates,
+                          VoidageRates&       voidage_rates) const
+{
+    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 auto [Rs, Rv] =
+        dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates);
+
+    const auto [Rsw, Rvw] =
+        dissolvedVaporisedRatio(iw, ig, rsw, rvw, surface_rates);
+
+
+    std::fill_n(&voidage_rates[0], pu.num_phases, 0.0);
+
+
+    // Determinant of 'R' matrix
+    const auto detRw = 1.0 - (Rsw * Rvw);
+
+    if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
+        // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s)
+        voidage_rates[iw] = surface_rates[iw];
+
+        const auto bw = FluidSystem::waterPvt()
+            .inverseFormationVolumeFactor(pvtRegionIdx, T, p,
+                                          Rsw,
+                                          saltConcentration);
+
+        if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+            voidage_rates[iw] -= Rvw * surface_rates[ig];
+        }
+        voidage_rates[iw] /= bw * detRw;
+    }
+
+    // Determinant of 'R' matrix
+    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)
+        voidage_rates[io] = surface_rates[io];
+        if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+            voidage_rates[io] -= Rv * surface_rates[ig];
+        }
+
+        const auto bo = FluidSystem::oilPvt()
+            .inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
+
+        voidage_rates[io] /= bo * detR;
+    }
+
+    // we only support either gas in water
+    // or gas in oil
+    if (detR != 1 && detRw != 1) {
+        std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR);
+        throw std::range_error(msg);
+    }
+    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
+        // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
+        voidage_rates[ig] = surface_rates[ig];
+        if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
+            voidage_rates[ig] -= Rs * surface_rates[io];
+        }
+        if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
+            voidage_rates[ig] -= Rsw * surface_rates[iw];
+        }
+
+        const auto bg = FluidSystem::gasPvt()
+            .inverseFormationVolumeFactor(pvtRegionIdx, T, p,
+                                          Rv, Rvw);
+
+        // we only support either gas in water or gas in oil
+        if (detRw == 1) {
+            voidage_rates[ig] /= bg * detR;
+        } else {
+            voidage_rates[ig] /= bg * detRw;
+        }
+    }
+}
+
+template <class FluidSystem, class Region>
+template <class Rates>
+void SurfaceToReservoirVoidage<FluidSystem,Region>::
+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.rsw, ra.rvw,
+                                    ra.temperature,
+                                    ra.saltConcentration,
+                                    surface_rates,
+                                    voidage_rates);
+}
+
+template <class FluidSystem, class Region>
+template <class Rates>
+std::pair<double, double>
+SurfaceToReservoirVoidage<FluidSystem,Region>::
+inferDissolvedVaporisedRatio(const double rsMax,
+                             const double rvMax,
+                             const Rates& surface_rates) const
+{
+    const auto io = RegionAttributeHelpers::PhasePos::oil(this->phaseUsage_);
+    const auto ig = RegionAttributeHelpers::PhasePos::gas(this->phaseUsage_);
+    return dissolvedVaporisedRatio(io, ig, rsMax, rvMax, surface_rates);
+}
+
+using FS = BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>;
+template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    sumRates(std::unordered_map<int,Attributes>&,
+             std::unordered_map<int,Attributes>&,
+             Parallel::Communication);
+
+template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    calcInjCoeff<std::vector<double>>(const int, const int, std::vector<double>&) const;
+template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    calcCoeff<std::vector<double>>(const int, const int, std::vector<double>&) const;
+template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    calcReservoirVoidageRates<std::vector<double>,std::vector<double>>(const int,
+                                                                       const double,
+                                                                       const double,
+                                                                       const double,
+                                                                       const double,
+                                                                       const double,
+                                                                       const double,
+                                                                       const double,
+                                                                       const std::vector<double>&,
+                                                                       std::vector<double>&) const;
+template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    calcReservoirVoidageRates<double const*,double*>(const int,
+                                                     const double,
+                                                     const double,
+                                                     const double,
+                                                     const double,
+                                                     const double,
+                                                     const double,
+                                                     const double,
+                                                     double const* const&,
+                                                     double*&) const;
+
+template void SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    calcReservoirVoidageRates<std::vector<double>>(const int,
+                                                   const int,
+                                                   const std::vector<double>&,
+                                                   std::vector<double>&) const;
+
+template std::pair<double,double>
+SurfaceToReservoirVoidage<FS,std::vector<int>>::
+    inferDissolvedVaporisedRatio<std::vector<double>::iterator>(const double,
+                                                                const double,
+                                                                const std::vector<double>::iterator&) const;
+
+} // namespace RateConverter
+} // namespace Opm

opm/simulators/wells/RateConverter.hpp

@@ -27,16 +27,13 @@
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/wells/RegionAttributeHelpers.hpp>
 
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
 #include <dune/grid/common/gridenums.hh>
 #include <dune/grid/common/rangegenerators.hh>
 
-#include <algorithm>
 #include <array>
 #include <cassert>
-#include <cmath>
-#include <memory>
-#include <stdexcept>
-#include <type_traits>
 #include <unordered_map>
 #include <utility>
 #include <vector>
@@ -108,7 +105,6 @@ namespace Opm {
                     ra.saltConcentration = 0.0;
                     ra.rsw = 0.0;
                     ra.rvw = 0.0;
-
                 }
 
                 // quantities for pore volume average
@@ -204,26 +200,9 @@ namespace Opm {
 
                 OPM_END_PARALLEL_TRY_CATCH("SurfaceToReservoirVoidage::defineState() failed: ", simulator.vanguard().grid().comm());
 
-                for (const auto& reg : rmap_.activeRegions()) {
+                this->sumRates(attributes_hpv,
-                      // Calculating first using the hydrocarbon pore volumes
+                               attributes_pv,
-                      auto& ra = attr_.attributes(reg);
+                               comm);
-                      const auto& attri_hpv = attributes_hpv[reg];
-                      ra.data = attri_hpv.data;
-                      comm.sum(ra.data.data(), ra.data.size());
-                      // TODO: should we have some epsilon here instead of zero?
-                      // No hydrocarbon pore volume, redo but this time using full pore volumes.
-                      if (ra.pv <= 0.) {
-                          // using the pore volume to do the averaging
-                          const auto& attri_pv = attributes_pv[reg];
-                          ra.data = attri_pv.data;
-                          comm.sum(ra.data.data(), ra.data.size());
-                          assert(ra.pv > 0.);
-                      }
-                      const double pv_sum = ra.pv;
-                      for (double& d : ra.data)
-                          d /= pv_sum;
-                      ra.pv = pv_sum;
-                }
             }
 
             /**
@@ -231,7 +210,7 @@ namespace Opm {
              *
              * Integral type.
              */
-            typedef typename RegionMapping<Region>::RegionId RegionId;
+            using RegionId = typename RegionMapping<Region>::RegionId;
 
             /**
              * Compute coefficients for surface-to-reservoir voidage
@@ -262,120 +241,11 @@ namespace Opm {
              */
             template <class Coeff>
             void
-            calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
+            calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const;
-            {
-                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 int   iw = RegionAttributeHelpers::PhasePos::water(pu);
-                const int   io = RegionAttributeHelpers::PhasePos::oil  (pu);
-                const int   ig = RegionAttributeHelpers::PhasePos::gas  (pu);
-
-                std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
-
-                // Actual Rsw and Rvw:
-                double Rsw = ra.rsw;
-                double Rvw = ra.rvw;
-                // Determinant of 'R' matrix
-                const double detRw = 1.0 - (Rsw * Rvw);
-
-                if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-                    // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s)
-
-                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rsw, saltConcentration);
-
-                    const double den = bw * detRw;
-
-                    coeff[iw] += 1.0 / den;
-
-                    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-                        coeff[ig] -= ra.rvw / den;
-                    }
-                }
-
-                // Actual Rs and Rv:
-                double Rs = ra.rs;
-                double Rv = ra.rv;
-
-                // Determinant of 'R' matrix
-                const double detR = 1.0 - (Rs * Rv);
-
-                // Currently we only support either gas in water or gas in oil
-                // not both
-                if (detR != 1 && detRw != 1) {
-                    std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR);
-                    throw std::range_error(msg);
-                }
-
-                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;
-
-                    coeff[io] += 1.0 / den;
-
-                    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-                        coeff[ig] -= ra.rv / den;
-                    }
-                }
-
-                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, Rvw);
-                    if (FluidSystem::enableDissolvedGasInWater()) {
-                        const double denw = bg * detRw;
-                        coeff[ig] += 1.0 / denw;
-                        if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-                           coeff[iw] -= ra.rsw / denw;
-                        }
-                    } else {
-                        const double den = bg * detR;
-                        coeff[ig] += 1.0 / den;
-                        if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
-                            coeff[io] -= ra.rs / den;
-                        }
-                    }
-                }
-            }
 
             template <class Coeff>
             void
-            calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
+            calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const;
-            {
-                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 int   iw = RegionAttributeHelpers::PhasePos::water(pu);
-                const int   io = RegionAttributeHelpers::PhasePos::oil  (pu);
-                const int   ig = RegionAttributeHelpers::PhasePos::gas  (pu);
-
-                std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
-
-                if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-                    // q[w]_r = q[w]_s / bw
-
-                    const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, saltConcentration);
-
-                    coeff[iw] = 1.0 / bw;
-                }
-
-                if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
-                    const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
-                    coeff[io] += 1.0 / bo;
-                }
-
-                if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-                    const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, 0.0);
-                    coeff[ig] += 1.0 / bg;
-                }
-            }
 
             /**
              * Convert surface volume flow rates to reservoir voidage flow
@@ -403,18 +273,7 @@ namespace Opm {
             void calcReservoirVoidageRates(const RegionId r,
                                            const int      pvtRegionIdx,
                                            const Rates&   surface_rates,
-                                           Rates&         voidage_rates) const
+                                           Rates&         voidage_rates) const;
-            {
-                const auto& ra = this->attr_.attributes(r);
-
-                this->calcReservoirVoidageRates(pvtRegionIdx,
-                                                ra.pressure, ra.rs, ra.rv, 
-                                                ra.rsw, ra.rvw, 
-                                                ra.temperature,
-                                                ra.saltConcentration,
-                                                surface_rates,
-                                                voidage_rates);
-            }
 
             /**
              * Convert surface volume flow rates to reservoir voidage flow
@@ -460,96 +319,13 @@ namespace Opm {
                                            const double        T,
                                            const double        saltConcentration,
                                            const SurfaceRates& surface_rates,
-                                           VoidageRates&       voidage_rates) const
+                                           VoidageRates&       voidage_rates) const;
-            {
-                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 auto [Rs, Rv] = this->
-                    dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates);
-
-                const auto [Rsw, Rvw] = this->
-                    dissolvedVaporisedRatio(iw, ig, rsw, rvw, surface_rates);
-
-
-                std::fill_n(&voidage_rates[0], pu.num_phases, 0.0);
-
-
-                // Determinant of 'R' matrix
-                const auto detRw = 1.0 - (Rsw * Rvw);
-
-                if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-                    // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s)
-                    voidage_rates[iw] = surface_rates[iw];
-
-                    const auto bw = FluidSystem::waterPvt()
-                        .inverseFormationVolumeFactor(pvtRegionIdx, T, p,
-                                                      Rsw,
-                                                      saltConcentration);
-
-                    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-                        voidage_rates[iw] -= Rvw * surface_rates[ig];
-                    }
-                    voidage_rates[iw] /= bw * detRw;
-                }
-
-                // Determinant of 'R' matrix
-                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)
-                    voidage_rates[io] = surface_rates[io];
-                    if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-                        voidage_rates[io] -= Rv * surface_rates[ig];
-                    }
-
-                    const auto bo = FluidSystem::oilPvt()
-                        .inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
-
-                    voidage_rates[io] /= bo * detR;
-                }
-
-                // we only support either gas in water
-                // or gas in oil
-                if (detR != 1 && detRw != 1) {
-                    std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR);
-                    throw std::range_error(msg);
-                }
-                if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
-                    // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
-                    voidage_rates[ig] = surface_rates[ig];
-                    if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
-                        voidage_rates[ig] -= Rs * surface_rates[io];
-                    }
-                    if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
-                        voidage_rates[ig] -= Rsw * surface_rates[iw];
-                    }
-
-                    const auto bg = FluidSystem::gasPvt()
-                        .inverseFormationVolumeFactor(pvtRegionIdx, T, p,
-                                                      Rv, Rvw);
-
-                    // we only support either gas in water or gas in oil
-                    if (detRw == 1) {
-                        voidage_rates[ig] /= bg * detR;
-                    } else { 
-                        voidage_rates[ig] /= bg * detRw;
-                    }
-                }
-            }
 
             template <class Rates>
             std::pair<double, double>
             inferDissolvedVaporisedRatio(const double rsMax,
                                          const double rvMax,
-                                         const Rates& surface_rates) const
+                                         const Rates& surface_rates) const;
-            {
-                const auto io = RegionAttributeHelpers::PhasePos::oil(this->phaseUsage_);
-                const auto ig = RegionAttributeHelpers::PhasePos::gas(this->phaseUsage_);
-                return this->dissolvedVaporisedRatio(io, ig, rsMax, rvMax, surface_rates);
-            }
 
             /**
              * Compute coefficients for surface-to-reservoir voidage
@@ -625,30 +401,11 @@ namespace Opm {
                 double& saltConcentration;
             };
 
+            void sumRates(std::unordered_map<RegionId,Attributes>& attributes_hpv,
+                          std::unordered_map<RegionId,Attributes>& attributes_pv,
+                          Parallel::Communication comm);
+
             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