changed: put calculation of injection gas-oil perf rates in separate method

opm/simulators/wells/StandardWell.hpp

@@ -467,6 +467,15 @@ namespace Opm
                                 const std::vector<EvalWell>& cq_s,
                                 const IntensiveQuantities& intQuants) const;
 
+        template<class Value>
+        void gasOilPerfRateInj(const std::vector<Value>& cq_s,
+                               PerforationRates& perf_rates,
+                               const Value& rv,
+                               const Value& rs,
+                               const Value& pressure,
+                               const Value& rvw,
+                               DeferredLogger& deferred_logger) const;
+
         template<class Value>
         void gasOilPerfRateProd(std::vector<Value>& cq_s,
                                 PerforationRates& perf_rates,

opm/simulators/wells/StandardWell_impl.hpp

@@ -292,42 +292,8 @@ namespace Opm
             // calculating the perforation solution gas rate and solution oil rates
             if (this->isProducer()) {
                 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                    const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
+                    gasOilPerfRateInj(cq_s, perf_rates,
-                    const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+                                      rv, rs, pressure, rvw, deferred_logger);
-                    // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
-                    // s means standard condition, r means reservoir condition
-                    // q_os = q_or * b_o + rv * q_gr * b_g
-                    // q_gs = q_gr * b_g + rs * q_or * b_o
-                    // d = 1.0 - rs * rv
-                    // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
-                    // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
-
-                    const double d = 1.0 - getValue(rv) * getValue(rs);
-
-                    if (d <= 0.0) {
-                        std::ostringstream sstr;
-                        sstr << "Problematic d value " << d << " obtained for well " << this->name()
-                             << " during computePerfRate calculations with rs " << rs
-                            << ", rv " << rv << " and pressure " << pressure
-                            << " obtaining d " << d
-                            << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
-                            << " for this connection.";
-                        deferred_logger.debug(sstr.str());
-                    } else {
-                        // vaporized oil into gas
-                        // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
-                        perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
-                        // dissolved of gas in oil
-                        // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
-                        perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
-                    }
-                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                        // q_ws = q_wr * b_w + rvw * q_gr * b_g
-                        // q_wr = 1 / b_w * (q_ws - rvw * q_gr * b_g) = 1 / b_w * (q_ws - rvw * 1 / d  (q_gs - rs * q_os))
-                        // vaporized water in gas
-                        // rvw * q_gr * b_g = q_ws -q_wr *b_w = rvw * (q_gs -rs *q_os) / d
-                        perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
-                    }
                 }
                 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                     //no oil
@@ -2401,6 +2367,60 @@ namespace Opm
         return {Base::restrictEval(cq_s_zfrac_effective), cq_s_zfrac_effective};
     }
 
+
+    template <typename TypeTag>
+    template<class Value>
+    void
+    StandardWell<TypeTag>::
+    gasOilPerfRateInj(const std::vector<Value>& cq_s,
+                      PerforationRates& perf_rates,
+                      const Value& rv,
+                      const Value& rs,
+                      const Value& pressure,
+                      const Value& rvw,
+                      DeferredLogger& deferred_logger) const
+    {
+        const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
+        const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+        // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
+        // s means standard condition, r means reservoir condition
+        // q_os = q_or * b_o + rv * q_gr * b_g
+        // q_gs = q_gr * b_g + rs * q_or * b_o
+        // d = 1.0 - rs * rv
+        // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
+        // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
+
+        const double d = 1.0 - getValue(rv) * getValue(rs);
+
+        if (d <= 0.0) {
+            std::ostringstream sstr;
+            sstr << "Problematic d value " << d << " obtained for well " << this->name()
+                 << " during computePerfRate calculations with rs " << rs
+                << ", rv " << rv << " and pressure " << pressure
+                << " obtaining d " << d
+                << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
+                << " for this connection.";
+            deferred_logger.debug(sstr.str());
+        } else {
+            // vaporized oil into gas
+            // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
+            perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
+            // dissolved of gas in oil
+            // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
+            perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
+        }
+
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            // q_ws = q_wr * b_w + rvw * q_gr * b_g
+            // q_wr = 1 / b_w * (q_ws - rvw * q_gr * b_g) = 1 / b_w * (q_ws - rvw * 1 / d  (q_gs - rs * q_os))
+            // vaporized water in gas
+            // rvw * q_gr * b_g = q_ws -q_wr *b_w = rvw * (q_gs -rs *q_os) / d
+            perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
+        }
+    }
+
+
+
     template <typename TypeTag>
     template<class Value>
     void