diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index 0fc574134..dca9fbcaa 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -595,6 +595,7 @@ namespace Opm {
         well_perf_data_.resize(wells_ecl_.size());
         int well_index = 0;
         for (const auto& well : wells_ecl_) {
+            std::size_t completion_index = 0;
             well_perf_data_[well_index].clear();
             well_perf_data_[well_index].reserve(well.getConnections().size());
             for (const auto& completion : well.getConnections()) {
@@ -614,6 +615,7 @@ namespace Opm {
                         pd.cell_index = active_index;
                         pd.connection_transmissibility_factor = completion.CF();
                         pd.satnum_id = completion.satTableId();
+                        pd.ecl_index = completion_index;
                         well_perf_data_[well_index].push_back(pd);
                     }
                 } else {
@@ -622,6 +624,7 @@ namespace Opm {
                                   "Completion state: " << Connection::State2String(completion.state()) << " not handled");
                     }
                 }
+                ++completion_index;
             }
             ++well_index;
         }
diff --git a/opm/simulators/wells/PerforationData.hpp b/opm/simulators/wells/PerforationData.hpp
index 47438026d..3a97bf28b 100644
--- a/opm/simulators/wells/PerforationData.hpp
+++ b/opm/simulators/wells/PerforationData.hpp
@@ -29,6 +29,8 @@ struct PerforationData
     int cell_index;
     double connection_transmissibility_factor;
     int satnum_id;
+    /// \brief The original index of the perforation in ECL Schedule
+    std::size_t ecl_index;
 };
 
 } // namespace Opm
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index 5ffb9fba8..cb9024841 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/opm/simulators/wells/WellInterface.hpp
@@ -389,7 +389,7 @@ namespace Opm
         // of states of individual well.
         int first_perf_;
 
-        std::vector<int> originalConnectionIndex_;
+        const std::vector<PerforationData>* perf_data_;
 
         std::vector<RateVector> connectionRates_;
 
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index 7c0fa6a77..52dc68621 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -48,7 +48,12 @@ namespace Opm
       , number_of_phases_(num_phases)
       , index_of_well_(index_of_well)
       , first_perf_(first_perf_index)
+      , perf_data_(&perf_data)
     {
+        assert(std::is_sorted(perf_data.begin(), perf_data.end(),
+                              [](const auto& perf1, const auto& perf2){
+                                  return perf1.ecl_index < perf2.ecl_index;
+                              }));
         if (time_step < 0) {
             OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
         }
@@ -71,16 +76,6 @@ namespace Opm
                 saturation_table_number_[perf] = pd.satnum_id;
                 ++perf;
             }
-
-            int all_perf = 0;
-            originalConnectionIndex_.reserve(perf_data.size());
-            for (const auto& connection : well.getConnections()) {
-                if (connection.state() == Connection::State::OPEN) {
-                    originalConnectionIndex_.push_back(all_perf);
-                }
-                ++all_perf;
-            }
-            assert(originalConnectionIndex_.size() == perf_data.size());
         }
 
         // initialization of the completions mapping
@@ -144,15 +139,22 @@ namespace Opm
         assert(completions_.empty() );
 
         const WellConnections& connections = well_ecl_.getConnections();
-        const int num_conns = connections.size();
+        const std::size_t num_conns = connections.size();
 
         int num_active_connections = 0;
-        for (int c = 0; c < num_conns; ++c) {
+        auto my_next_perf = perf_data_->begin();
+        for (std::size_t c = 0; c < num_conns; ++c) {
+            if (my_next_perf->ecl_index > c)
+            {
+                continue;
+            }
+            assert(my_next_perf->ecl_index == c);
             if (connections[c].state() == Connection::State::OPEN) {
                 completions_[connections[c].complnum()].push_back(num_active_connections++);
             }
+            ++my_next_perf;
         }
-        assert(num_active_connections == number_of_perforations_);
+        assert(my_next_perf == perf_data_->end());
     }
 
 
@@ -1382,7 +1384,7 @@ namespace Opm
     void
     WellInterface<TypeTag>::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger)
     {
-        const auto& connection = well_ecl_.getConnections()[originalConnectionIndex_[perfIdx]];
+        const auto& connection = well_ecl_.getConnections()[(*perf_data_)[perfIdx].ecl_index];
         if (well_ecl_.getDrainageRadius() < 0) {
             if (new_well && perfIdx == 0) {
                 deferred_logger.warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero");
diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
index 5729e021b..9764c35b5 100644
--- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
+++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
@@ -1025,25 +1025,19 @@ namespace Opm
             // Create a function that calls some function
             // for all the individual data items to simplify
             // the further code.
+            auto iterateContainer = [](auto& container, auto& func) {
+                for (auto& x : container) {
+                    func(x.second);
+                }
+            };
+
             auto forAllGroupData = [&](auto& func) {
-                for (auto& x : injection_group_rein_rates) {
-                    func(x.second);
-                }
-                for (auto& x : production_group_reduction_rates) {
-                    func(x.second);
-                }
-                for (auto& x : injection_group_reduction_rates) {
-                    func(x.second);
-                }
-                for (auto& x : injection_group_reservoir_rates) {
-                    func(x.second);
-                }
-                for (auto& x : production_group_rates) {
-                    func(x.second);
-                }
-                for (auto& x : well_rates) {
-                    func(x.second);
-                }
+                iterateContainer(injection_group_rein_rates, func);
+                iterateContainer(production_group_reduction_rates, func);
+                iterateContainer(injection_group_reduction_rates, func);
+                iterateContainer(injection_group_reservoir_rates, func);
+                iterateContainer(production_group_rates, func);
+                iterateContainer(well_rates,func);
             };
 
             // Compute the size of the data.