Merge pull request #1898 from GitPaean/refactoring_wtest_for_wellcut_completions

Fixing completion water cut checking for WTEST and some related refactoring

opm/simulators/wells/WellInterface.hpp

@@ -320,6 +320,21 @@ namespace Opm
         // well bore diameter
         std::vector<double> bore_diameters_;
 
+        /*
+         *  completions_ contains the mapping from completion id to connection indices
+         *  {
+         *      2 : [ConnectionIndex, ConnectionIndex],
+         *      1 : [ConnectionIndex, ConnectionIndex, ConnectionIndex],
+         *      5 : [ConnectionIndex],
+         *      7 : [ConnectionIndex]
+         *      ...
+         *   }
+         *   The integer IDs correspond to the COMPLETION id given by the COMPLUMP keyword.
+         *   When there is no COMPLUMP keyword used, a default completion number will be assigned
+         *   based on the order of the declaration of the connections
+         */
+        std::map<int, std::vector<int>> completions_;
+
         const PhaseUsage* phase_usage_;
 
         bool getAllowCrossFlow() const;
@@ -377,6 +392,18 @@ namespace Opm
                                              const WellState& well_state,
                                              Opm::DeferredLogger& deferred_logger) const;
 
+        template <typename RatioFunc>
+        bool checkMaxRatioLimitWell(const WellState& well_state,
+                                    const double max_ratio_limit,
+                                    const RatioFunc& ratioFunc) const;
+
+        template <typename RatioFunc>
+        void checkMaxRatioLimitCompletions(const WellState& well_state,
+                                           const double max_ratio_limit,
+                                           const RatioFunc& ratioFunc,
+                                           int& worst_offending_completion,
+                                           double& violation_extent) const;
+
         double scalingFactor(const int comp_idx) const;
 
         // whether a well is specified with a non-zero and valid VFP table number
@@ -417,6 +444,8 @@ namespace Opm
 
         void scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger);
 
+        void initCompletions();
+
         // count the number of times an output log message is created in the productivity
         // index calculations
         int well_productivity_index_logger_counter_;

opm/simulators/wells/WellInterface_impl.hpp

@@ -100,6 +100,9 @@ namespace Opm
                       saturation_table_number_.begin() );
         }
 
+        // initialization of the completions mapping
+        initCompletions();
+
         well_efficiency_factor_ = 1.0;
 
         connectionRates_.resize(number_of_perforations_);
@@ -137,6 +140,28 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    initCompletions()
+    {
+        assert(completions_.empty() );
+
+        const WellConnections& connections = well_ecl_.getConnections();
+        const int num_conns = connections.size();
+
+        assert(num_conns == number_of_perforations_);
+
+        for (int c = 0; c < num_conns; c++) {
+            completions_[connections[c].complnum()].push_back(c);
+        }
+    }
+
+
+
+
+
+
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
@@ -583,77 +608,39 @@ namespace Opm
     checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
                           const WellState& well_state) const
     {
-        bool water_cut_limit_violated = false;
-        int worst_offending_completion = INVALIDCOMPLETION;
-        double violation_extent = -1.0;
-
-        const int np = number_of_phases_;
-        const Opm::PhaseUsage& pu = phaseUsage();
-        const int well_number = index_of_well_;
 
         assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
         assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
 
-        const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
+        // function to calculate water cut based on rates
-        const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
+        auto waterCut = [](const std::vector<double>& rates,
-        const double liquid_rate = oil_rate + water_rate;
+                           const PhaseUsage& pu) {
-        double water_cut;
+
-        if (std::abs(liquid_rate) != 0.) {
+            const double oil_rate = rates[pu.phase_pos[Oil]];
-            water_cut = water_rate / liquid_rate;
+            const double water_rate = rates[pu.phase_pos[Water]];
-        } else {
+
-            water_cut = 0.0;
+            // both rate should be in the same direction
-        }
+            assert(oil_rate * water_rate >= 0.);
+
+            const double liquid_rate = oil_rate + water_rate;
+            if (liquid_rate != 0.) {
+                return (water_rate / liquid_rate);
+            } else {
+                return 0.;
+            }
+        };
 
         const double max_water_cut_limit = econ_production_limits.maxWaterCut();
-        if (water_cut > max_water_cut_limit) {
+        assert(max_water_cut_limit > 0.);
-            water_cut_limit_violated = true;
+
-        }
+        const bool water_cut_limit_violated = checkMaxRatioLimitWell(well_state, max_water_cut_limit, waterCut);
+
+        int worst_offending_completion = INVALIDCOMPLETION;
+        double violation_extent = -1.;
 
         if (water_cut_limit_violated) {
-            // need to handle the worst_offending_connection
+            checkMaxRatioLimitCompletions(well_state, max_water_cut_limit, waterCut,
-            const int perf_start = first_perf_;
+                                          worst_offending_completion, violation_extent);
-            const int perf_number = number_of_perforations_;
-
-            std::vector<double> water_cut_perf(perf_number);
-            for (int perf = 0; perf < perf_number; ++perf) {
-                const int i_perf = perf_start + perf;
-                const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
-                const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
-                const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
-                if (std::abs(liquid_perf_rate) != 0.) {
-                    water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
-                } else {
-                    water_cut_perf[perf] = 0.;
-                }
-            }
-            const auto& completions = well_ecl_.getCompletions();
-            const auto& connections = well_ecl_.getConnections();
-
-            int complnumIdx = 0;
-            std::vector<double> water_cut_in_completions(completions.size(), 0.0);
-            for (const auto& completion : completions) {
-                int complnum = completion.first;
-                for (int perf = 0; perf < perf_number; ++perf) {
-                    if (complnum == connections.get ( perf ).complnum()) {
-                        water_cut_in_completions[complnumIdx] +=  water_cut_perf[perf];
-                    }
-                }
-                complnumIdx++;
-            }
-
-            double max_water_cut_perf = 0.;
-            complnumIdx = 0;
-            for (const auto& completion : completions) {
-                if (water_cut_in_completions[complnumIdx] > max_water_cut_perf) {
-                    worst_offending_completion = completion.first;
-                    max_water_cut_perf = water_cut_in_completions[complnumIdx];
-                }
-                complnumIdx++;
-            }
-
-            assert(max_water_cut_limit != 0.);
-            assert(worst_offending_completion != INVALIDCOMPLETION);
-            violation_extent = max_water_cut_perf / max_water_cut_limit;
         }
 
         return std::make_tuple(water_cut_limit_violated, worst_offending_completion, violation_extent);
@@ -718,6 +705,80 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    template<typename RatioFunc>
+    bool
+    WellInterface<TypeTag>::
+    checkMaxRatioLimitWell(const WellState& well_state,
+                           const double max_ratio_limit,
+                           const RatioFunc& ratioFunc) const
+    {
+        const int np = number_of_phases_;
+
+        std::vector<double> well_rates(np, 0.0);
+
+        for (int p = 0; p < np; ++p) {
+            well_rates[p] = well_state.wellRates()[index_of_well_ * np + p];
+        }
+
+        const double well_ratio = ratioFunc(well_rates, phaseUsage());
+
+        return (well_ratio > max_ratio_limit);
+    }
+
+
+
+
+    template<typename TypeTag>
+    template<typename RatioFunc>
+    void
+    WellInterface<TypeTag>::
+    checkMaxRatioLimitCompletions(const WellState& well_state,
+                                  const double max_ratio_limit,
+                                  const RatioFunc& ratioFunc,
+                                  int& worst_offending_completion,
+                                  double& violation_extent) const
+    {
+        worst_offending_completion = INVALIDCOMPLETION;
+        violation_extent = -1.0;
+
+        // the maximum water cut value of the completions
+        // it is used to identify the most offending completion
+        double max_ratio_completion = 0;
+
+        // look for the worst_offending_completion
+        for (const auto& completion : completions_) {
+
+            const int np = number_of_phases_;
+            std::vector<double> completion_rates(np, 0.0);
+
+            // looping through the connections associated with the completion
+            const std::vector<int>& conns = completion.second;
+            for (const int c : conns) {
+                const int index_con = c + first_perf_;
+
+                for (int p = 0; p < np; ++p) {
+                    const double connection_rate = well_state.perfPhaseRates()[index_con * np + p];
+                    completion_rates[p] += connection_rate;
+                }
+            } // end of for (const int c : conns)
+
+            const double ratio_completion = ratioFunc(completion_rates, phaseUsage());
+
+            if (ratio_completion > max_ratio_completion) {
+                worst_offending_completion = completion.first;
+                max_ratio_completion = ratio_completion;
+            }
+        } // end of for (const auto& completion : completions_)
+
+        assert(max_ratio_completion > max_ratio_limit);
+        assert(worst_offending_completion != INVALIDCOMPLETION);
+        violation_extent = max_ratio_completion / max_ratio_limit;
+    }
+
+
+
+
 
     template<typename TypeTag>
     void