MS Wells - Fix problems of well with inactive perforations

opm/simulators/wells/MultisegmentWellGeneric.cpp

@@ -77,13 +77,17 @@ scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inle
             for (int seg = 0; seg < numberOfSegments(); ++seg) {
                 segment_rates[baseif_.numPhases() * seg + phase] *= well_phase_rate / unscaled_top_seg_rate;
             }
-        } else { // In this case, we calculated sumTw
+        } else if (baseif_.numPerfs() > 0) { // Only go here if there are actually active perforations on this process
+            // In this case, we calculated sumTw
             // only handling this specific phase
             constexpr Scalar num_single_phase = 1;
-            std::vector<Scalar> perforation_rates(num_single_phase * baseif_.numPerfs(), 0.0);
+            const int global_num_perf_this_well = ws.parallel_info.get().communication().sum(baseif_.numPerfs());
+            std::vector<Scalar> perforation_rates(num_single_phase * global_num_perf_this_well, 0.0);
+            std::cout << "perforation_rates.size() = " << perforation_rates.size() << std::endl;
             const Scalar perf_phaserate_scaled = ws.surface_rates[phase] / sumTw;
-            for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
+            for (int perf = 0; perf < global_num_perf_this_well; ++perf) { //strangely this loop does not go over the whole array...
                 perforation_rates[perf] = baseif_.wellIndex()[perf] * perf_phaserate_scaled;
+                std::cout << "perforation_rates[" << perf << "] = " << perforation_rates[perf] << std::endl;
             }
 
             std::vector<Scalar> rates;

opm/simulators/wells/MultisegmentWellSegments.cpp

@@ -80,6 +80,9 @@ MultisegmentWellSegments(const int numSegments,
 
     // initialize the segment_perforations_ and update perforation_segment_depth_diffs_
     const WellConnections& completion_set = well_.wellEcl().getConnections();
+    std::cout << "well_.numPerfs() " << well_.numPerfs() << std::endl;
+    std::cout << "well_.wellEcl().getConnections().size() " << well_.wellEcl().getConnections().size() << std::endl;
+    std::cout << "parallel_well_info.numLocalPerfs() " << parallel_well_info.numLocalPerfs() << std::endl;
     // index of the perforation within wells struct
     // there might be some perforations not active, which causes the number of the perforations in
     // well_ecl_ and wells struct different
@@ -106,10 +109,14 @@ MultisegmentWellSegments(const int numSegments,
             const Scalar segment_depth = segment_set[segment_index].depth();
             int local_perf_index = parallel_well_info.globalToLocal(i_perf_wells);
             if (local_perf_index > -1) // If local_perf_index == -1, then the perforation is not on this process
-                local_perforation_depth_diffs_[local_perf_index] = well_.perfDepth()[i_perf_wells] - segment_depth;
+                if (local_perforation_depth_diffs_.size() > static_cast<std::size_t>(local_perf_index)) {
+                    // If local_perforation_depth_diffs_.size() is not large enough, this is a SHUT connection
+                    local_perforation_depth_diffs_[local_perf_index] = well_.perfDepth()[i_perf_wells] - segment_depth;
+                }
             i_perf_wells++;
         }
     }
+    std::cout << "First loop of MultisegmentWellSegments constructor done, local_perforation_depth_diffs_.size() = " << local_perforation_depth_diffs_.size() << std::endl;
 
     // initialize the segment_inlets_
     for (const Segment& segment : segment_set) {
@@ -121,6 +128,7 @@ MultisegmentWellSegments(const int numSegments,
             inlets_[outlet_segment_index].push_back(segment_index);
         }
     }
+    std::cout << "Second loop of MultisegmentWellSegments constructor done, segment_set.size() = " << segment_set.size() << std::endl;
 
     // calculating the depth difference between the segment and its oulet_segments
     // for the top segment, we will make its zero unless we find other purpose to use this value
@@ -131,6 +139,7 @@ MultisegmentWellSegments(const int numSegments,
         const Scalar outlet_depth = outlet_segment.depth();
         depth_diffs_[seg] = segment_depth - outlet_depth;
     }
+    std::cout << "Third loop of MultisegmentWellSegments constructor done, numSegments = " << numSegments << std::endl;
 }
 
 template<class FluidSystem, class Indices>

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -350,6 +350,11 @@ namespace Opm
                 const int local_perf_index = this->pw_info_.globalToLocal(perf);
                 if (local_perf_index < 0) // then the perforation is not on this process
                     continue;
+                if (this->well_cells_.size() <= static_cast<std::size_t>(local_perf_index)) {// then the perforation belongs to a shut connection
+                    std::cout << "skip computeWellRatesWithBhp, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
+                    continue;
+                }
+                std::cout << "execute computeWellRatesWithBhp, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
                 const int cell_idx = this->well_cells_[local_perf_index];
                 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
                 // flux for each perforation
@@ -887,6 +892,11 @@ namespace Opm
         const int local_perf_index = this->pw_info_.globalToLocal(perf);
         if (local_perf_index < 0) // then the perforation is not on this process
             return;
+        if (this->cell_perforation_pressure_diffs_.size() <= static_cast<std::size_t>(local_perf_index)) { // then the perforation belongs to a shut connection
+            std::cout << "skip computePerfRate, local_perf_index = " << local_perf_index << ", << cell_perforation_depth_diffs_.size() = " << this->cell_perforation_pressure_diffs_.size() << std::endl;
+            return;
+        }
+        std::cout << "execute computePerfRate, local_perf_index = " << local_perf_index << ", << cell_perforation_depth_diffs_.size() = " << this->cell_perforation_pressure_diffs_.size() << std::endl;
 
         // pressure difference between the segment and the perforation
         const Value perf_seg_press_diff = this->gravity() * segment_density *
@@ -1102,25 +1112,37 @@ namespace Opm
         // perforated cell
         // TODO: later to investigate how to handle the pvt region
         int pvt_region_index;
-        {
-            // using the first perforated cell, so we look for global index 0
+
+        Scalar fsTemperature;
+        using SaltConcType = typename std::decay<decltype(std::declval<decltype(simulator.model().intensiveQuantities(0, 0).fluidState())>().saltConcentration())>::type;
+        SaltConcType fsSaltConcentration;
+        std::cout << "computeSegmentFluidProperties, this->well_cells_.size() = " << this->well_cells_.size() << std::endl;
+        if (this->well_cells_.size() > 0) {
+            // this->well_cells_ is empty if this process does not contain active perforations
+            // using the pvt region of first perforated cell, so we look for global index 0
+            // TODO: it should be a member of the WellInterface, initialized properly
             const int cell_idx = this->well_cells_[0];
             const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
             const auto& fs = intQuants.fluidState();
 
-            // The following broadcast calls are neccessary to ensure that processes that do *not* contain
-            // the first perforation get the correct temperature, saltConcentration and pvt_region_index
-            auto fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
-            fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
-            temperature.setValue(fsTemperature);
-
-            auto fsSaltConcentration = fs.saltConcentration();
-            fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
-            saltConcentration = this->extendEval(fsSaltConcentration);
-
+            fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
+            fsSaltConcentration = fs.saltConcentration();
             pvt_region_index = fs.pvtRegionIndex();
-            pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
         }
+        std::cout << "will broadcastFirstPerforationValue now" << std::endl;
+
+        // The following broadcast calls are neccessary to ensure that processes that do *not* contain
+        // the first perforation get the correct temperature, saltConcentration and pvt_region_index
+        fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
+        temperature.setValue(fsTemperature);
+        std::cout << "-" << std::endl;
+
+        fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
+        saltConcentration = this->extendEval(fsSaltConcentration);
+        std::cout << "-" << std::endl;
+
+        pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
+        std::cout << "done" << std::endl;
 
         this->segments_.computeFluidProperties(temperature,
                                                saltConcentration,
@@ -1267,6 +1289,12 @@ namespace Opm
                 const int local_perf_index = this->pw_info_.globalToLocal(perf);
                 if (local_perf_index < 0) // then the perforation is not on this process
                     continue;
+                if (this->well_cells_.size() <= static_cast<std::size_t>(local_perf_index)) {// then the perforation belongs to a shut connection
+                    std::cout << "skip updateIPR, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
+                    continue;
+                }
+                std::cout << "execute updateIPR, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
+
                 std::vector<Scalar> mob(this->num_components_, 0.0);
 
                 // TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
@@ -1778,6 +1806,12 @@ namespace Opm
                 const int local_perf_index = this->pw_info_.globalToLocal(perf);
                 if (local_perf_index < 0) // then the perforation is not on this process
                     continue;
+                if (this->well_cells_.size() <= static_cast<std::size_t>(local_perf_index)) {// then the perforation belongs to a shut connection
+                    std::cout << "skip assembleWellEqWithoutIteration, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
+                    continue;
+                }
+                std::cout << "execute assembleWellEqWithoutIteration, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
+
                 const int cell_idx = this->well_cells_[local_perf_index];
                 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
                 std::vector<EvalWell> mob(this->num_components_, 0.0);
@@ -1928,6 +1962,11 @@ namespace Opm
                 const int local_perf_index = this->pw_info_.globalToLocal(perf);
                 if (local_perf_index < 0) // then the perforation is not on this process
                     continue;
+                if (this->well_cells_.size() <= static_cast<std::size_t>(local_perf_index)) {// then the perforation belongs to a shut connection
+                    std::cout << "skip allDrawDownWrongDirection, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
+                    continue;
+                }
+                std::cout << "execute allDrawDownWrongDirection, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
 
                 const int cell_idx = this->well_cells_[local_perf_index];
                 const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
@@ -1986,26 +2025,37 @@ namespace Opm
         EvalWell temperature;
         EvalWell saltConcentration;
         int pvt_region_index;
-        {
+
+        Scalar fsTemperature;
+        using SaltConcType = typename std::decay<decltype(std::declval<decltype(simulator.model().intensiveQuantities(0, 0).fluidState())>().saltConcentration())>::type;
+        SaltConcType fsSaltConcentration;
+
+        std::cout << "getSegmentSurfaceVolume, this->well_cells_.size() = " << this->well_cells_.size() << std::endl;
+        if (this->well_cells_.size() > 0) {
+            // this->well_cells_ is empty if this process does not contain active perforations
             // using the pvt region of first perforated cell, so we look for global index 0
             // TODO: it should be a member of the WellInterface, initialized properly
             const int cell_idx = this->well_cells_[0];
             const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
             const auto& fs = intQuants.fluidState();
 
-            // The following broadcast calls are neccessary to ensure that processes that do *not* contain
-            // the first perforation get the correct temperature, saltConcentration and pvt_region_index
-            auto fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
-            fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
-            temperature.setValue(fsTemperature);
-
-            auto fsSaltConcentration = fs.saltConcentration();
-            fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
-            saltConcentration = this->extendEval(fsSaltConcentration);
-
+            fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
+            fsSaltConcentration = fs.saltConcentration();
             pvt_region_index = fs.pvtRegionIndex();
-            pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
         }
+        std::cout << "will broadcastFirstPerforationValue now" << std::endl;
+        // The following broadcast calls are neccessary to ensure that processes that do *not* contain
+        // the first perforation get the correct temperature, saltConcentration and pvt_region_index
+        fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature);
+        temperature.setValue(fsTemperature);
+        std::cout << "-" << std::endl;
+
+        fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration);
+        saltConcentration = this->extendEval(fsSaltConcentration);
+        std::cout << "-" << std::endl;
+
+        pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index);
+        std::cout << "done" << std::endl;
 
         return this->segments_.getSurfaceVolume(temperature,
                                                 saltConcentration,
@@ -2155,6 +2205,11 @@ namespace Opm
                 const int local_perf_index = this->pw_info_.globalToLocal(perf);
                 if (local_perf_index < 0) // then the perforation is not on this process
                     continue;
+                if (this->well_cells_.size() <= static_cast<std::size_t>(local_perf_index)) {// then the perforation belongs to a shut connection
+                    std::cout << "skip maxPerfPress, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
+                    continue;
+                }
+                std::cout << "execute maxPerfPress, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl;
 
                 const int cell_idx = this->well_cells_[local_perf_index];
                 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
@@ -2187,6 +2242,11 @@ namespace Opm
                 const int local_perf_index = this->pw_info_.globalToLocal(perf);
                 if (local_perf_index < 0) // then the perforation is not on this process
                     continue;
+                if (this->well_cells_.size() <= static_cast<std::size_t>(local_perf_index)) {// then the perforation belongs to a shut connection
+                    std::cout << "skip computeCurrentWellRates for segment " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() = " << this->well_cells_.size() << std::endl;
+                    continue;
+                }
+                std::cout << "execute computeCurrentWellRates for segment " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() = " << this->well_cells_.size() << std::endl;
 
                 const int cell_idx = this->well_cells_[local_perf_index];
                 const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);

opm/simulators/wells/ParallelWellInfo.cpp

@@ -74,7 +74,7 @@ GlobalPerfContainerFactory<Scalar>::
 GlobalPerfContainerFactory(const IndexSet& local_indices,
                            const Parallel::Communication comm,
                            const int num_local_perfs)
-    : local_indices_(local_indices), comm_(comm)
+    : local_indices_(local_indices), comm_(comm), num_local_perfs_(num_local_perfs)
 {
     if ( comm_.size() > 1 )
     {
@@ -165,6 +165,12 @@ int GlobalPerfContainerFactory<Scalar>::globalToLocal(const int globalIndex) con
     return it->second;
 }
 
+template<class Scalar>
+int GlobalPerfContainerFactory<Scalar>::numLocalPerfs() const {
+    return num_local_perfs_;
+}
+
+
 template<class Scalar>
 std::vector<Scalar> GlobalPerfContainerFactory<Scalar>::
 createGlobal(const std::vector<Scalar>& local_perf_container,
@@ -539,6 +545,14 @@ int ParallelWellInfo<Scalar>::globalToLocal(const int globalIndex) const {
         return globalIndex;
 }
 
+template<class Scalar>
+int ParallelWellInfo<Scalar>::numLocalPerfs() const {
+    if(globalPerfCont_)
+        return globalPerfCont_->numLocalPerfs();
+    else // If globalPerfCont_ is not set up, then this is a sequential run and local and the number of  indices are the same.
+        return 0;
+}
+
 template<class Scalar>
 void ParallelWellInfo<Scalar>::communicateFirstPerforation(bool hasFirst)
 {

opm/simulators/wells/ParallelWellInfo.hpp

@@ -164,6 +164,7 @@ class GlobalPerfContainerFactory
     int numGlobalPerfs() const;
     int globalToLocal(const int globalIndex) const;
     int localToGlobal(std::size_t localIndex) const;
+    int numLocalPerfs() const;
 
 private:
     void buildLocalToGlobalMap() const;
@@ -174,6 +175,7 @@ class GlobalPerfContainerFactory
     mutable bool g2l_map_built_ = false;
     const IndexSet& local_indices_;
     Parallel::Communication comm_;
+    int num_local_perfs_;
     int num_global_perfs_;
     /// \brief sizes for allgatherv
     std::vector<int> sizes_;
@@ -219,6 +221,7 @@ class ParallelWellInfo
 
     int globalToLocal(const int globalIndex) const;
     int localToGlobal(std::size_t localIndex) const;
+    int numLocalPerfs() const;
 
     /// If the well does not have any open connections the member rankWithFirstPerf
     /// is not initialized, and no broadcast is performed. In this case the argument

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -92,6 +92,8 @@ WellInterfaceGeneric(const Well& well,
 
     ref_depth_ = well.getRefDepth();
 
+    std::cout << "In WellInterfaceGeneric: this->well_ecl_.getConnections().size() = " << this->well_ecl_.getConnections().size() << std::endl;
+    std::cout << "In WellInterfaceGeneric: perf_data.size() = " << perf_data.size() << std::endl;
     // We do not want to count SHUT perforations here, so
     // it would be wrong to use wells.getConnections().size().
     // This is the number_of_local_perforations_ on this process only!