Merge pull request #3620 from blattms/removeSetupCart-rebased

[refactor] Remove cartesianToCompressed Mapping from well model

ebos/eclbasevanguard.hh

@@ -296,6 +296,8 @@ public:
     /*!
      * \brief Return compressed index from cartesian index
      *
+     * \return compressed index of cell is in interior, -1 otherwise
+     *
      */
     int compressedIndex(int cartesianCellIdx) const
     {
@@ -306,6 +308,25 @@ public:
             return -1;
     }
 
+    /*!
+     * \brief Return compressed index from cartesian index only in interior
+     *
+     * \return compressed index of cell is in interior, -1 otherwise
+     *
+     */
+    int compressedIndexForInterior(int cartesianCellIdx) const
+    {
+        auto index_pair = cartesianToCompressed_.find(cartesianCellIdx);
+        if (index_pair == cartesianToCompressed_.end() ||
+            !is_interior_[index_pair->second])
+        {
+            return -1;
+        }
+        else
+        {
+            return index_pair->second;
+        }
+    }
     /*!
      * \brief Extract Cartesian index triplet (i,j,k) of an active cell.
      *
@@ -424,9 +445,22 @@ protected:
     void updateCartesianToCompressedMapping_()
     {
         size_t num_cells = asImp_().grid().leafGridView().size(0);
-        for (unsigned i = 0; i < num_cells; ++i) {
+        is_interior_.resize(num_cells);
-            unsigned cartesianCellIdx = cartesianIndex(i);
+        
-            cartesianToCompressed_[cartesianCellIdx] = i;
+        ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
+        for (const auto& element : elements(this->gridView()))
+        {
+            const auto elemIdx = elemMapper.index(element);
+            unsigned cartesianCellIdx = cartesianIndex(elemIdx);
+            cartesianToCompressed_[cartesianCellIdx] = elemIdx;
+            if (element.partitionType() == Dune::InteriorEntity)
+            {
+                is_interior_[elemIdx] = 1;
+            }
+            else
+            {
+                is_interior_[elemIdx] = 0;
+            }
         }
     }
 
@@ -515,6 +549,10 @@ protected:
     /*! \brief Cell thickness
      */
     std::vector<Scalar> cellThickness_;
+
+    /*! \brief Whether a cells is in the interior.
+     */
+    std::vector<int> is_interior_;
 };
 
 } // namespace Opm

opm/simulators/wells/BlackoilWellModel.hpp

@@ -370,8 +370,6 @@ namespace Opm {
             // setting the well_solutions_ based on well_state.
             void updatePrimaryVariables(DeferredLogger& deferred_logger);
 
-            void setupCartesianToCompressed_(const int* global_cell, int local_num__cells);
-
             void updateAverageFormationFactor();
 
             void computePotentials(const std::size_t widx,
@@ -421,6 +419,10 @@ namespace Opm {
 
             void assignWellTracerRates(data::Wells& wsrpt) const;
 
+            int compressedIndexForInterior(int cartesian_cell_idx) const override {
+                return ebosSimulator_.vanguard().compressedIndexForInterior(cartesian_cell_idx);
+            }
+
         private:
             BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& pu);
         };

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -498,8 +498,7 @@ initializeWellPerfData()
         parallelWellInfo.beginReset();
 
         for (const auto& completion : well.getConnections()) {
-            const int active_index =
+            const int active_index = compressedIndexForInterior(completion.global_index());
-                cartesian_to_compressed_[completion.global_index()];
             if (completion.state() == Connection::State::OPEN) {
                 if (active_index >= 0) {
                     if (firstOpenCompletion)
@@ -1874,7 +1873,7 @@ BlackoilWellModelGeneric::
 setRepRadiusPerfLength()
 {
     for (const auto& well : well_container_generic_) {
-        well->setRepRadiusPerfLength(cartesian_to_compressed_);
+        well->setRepRadiusPerfLength();
     }
 }
 

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -405,6 +405,10 @@ protected:
     void runWellPIScaling(const int timeStepIdx,
                           DeferredLogger& local_deferredLogger);
 
+    /// \brief get compressed index for interior cells (-1, otherwise
+    virtual int compressedIndexForInterior(int cartesian_cell_idx) const = 0;
+
+
     Schedule& schedule_;
     const SummaryState& summaryState_;
     const EclipseState& eclState_;
@@ -432,11 +436,6 @@ protected:
 
     std::vector<WellProdIndexCalculator> prod_index_calc_;
 
-    // Map from logically cartesian cell indices to compressed ones.
-    // Cells not in the interior are not mapped. This deactivates
-    // these for distributed wells and makes the distribution non-overlapping.
-    std::vector<int> cartesian_to_compressed_;
-
     std::vector<int> pvt_region_idx_;
 
     mutable std::unordered_set<std::string> closed_this_step_;

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -51,18 +51,14 @@ namespace Opm {
         // Number of cells the global grid view
         global_num_cells_ = ebosSimulator_.vanguard().globalNumCells();
 
-        // Set up cartesian mapping.
+        // Set up parallel wells
+        auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
+
+        this->parallel_well_info_.reserve(parallel_wells.size());
+        for( const auto& name_bool: parallel_wells)
         {
-            const auto& grid = this->ebosSimulator_.vanguard().grid();
+            this->parallel_well_info_.emplace_back(name_bool,
-            const auto& cartDims = UgGridHelpers::cartDims(grid);
+                                                   ebosSimulator_.gridView().comm());
-            setupCartesianToCompressed_(UgGridHelpers::globalCell(grid),
-                                        cartDims[0] * cartDims[1] * cartDims[2]);
-
-            auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
-
-            for (const auto& wellinfo : parallel_wells) {                   
-                this->parallel_well_info_.emplace_back(wellinfo, grid.comm());         
-            }
         }
 
         this->alternative_well_rate_init_ =
@@ -129,8 +125,7 @@ namespace Opm {
             for ( size_t c=0; c < connectionSet.size(); c++ )
             {
                 const auto& connection = connectionSet.get(c);
-                int compressed_idx = cartesian_to_compressed_
+                int compressed_idx = compressedIndexForInterior(connection.global_index());
-                    .at(connection.global_index());
 
                 if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells.
                     wellCells.push_back(compressed_idx);
@@ -1431,31 +1426,6 @@ namespace Opm {
 
 
 
-    template<typename TypeTag>
-    void
-    BlackoilWellModel<TypeTag>::
-    setupCartesianToCompressed_(const int* global_cell, int number_of_cartesian_cells)
-    {
-        cartesian_to_compressed_.resize(number_of_cartesian_cells, -1);
-        if (global_cell) {
-            auto elemIt = ebosSimulator_.gridView().template begin</*codim=*/ 0>();
-            for (unsigned i = 0; i < local_num_cells_; ++i) {
-                // Skip perforations in the overlap/ghost for distributed wells.
-                if (elemIt->partitionType() == Dune::InteriorEntity)
-                {
-                    assert(ebosSimulator_.gridView().indexSet().index(*elemIt) == static_cast<int>(i));
-                    cartesian_to_compressed_[global_cell[i]] = i;
-                }
-                ++elemIt;
-            }
-        }
-        else {
-            for (unsigned i = 0; i < local_num_cells_; ++i) {
-                cartesian_to_compressed_[i] = i;
-            }
-        }
-
-    }
 
     template<typename TypeTag>
     void

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -263,7 +263,7 @@ void WellInterfaceGeneric::setWellEfficiencyFactor(const double efficiency_facto
     well_efficiency_factor_ = efficiency_factor;
 }
 
-void WellInterfaceGeneric::setRepRadiusPerfLength(const std::vector<int>& cartesian_to_compressed)
+void WellInterfaceGeneric::setRepRadiusPerfLength()
 {
     const int nperf = number_of_perforations_;
 
@@ -275,31 +275,35 @@ void WellInterfaceGeneric::setRepRadiusPerfLength(const std::vector<int>& cartes
     perf_length_.reserve(nperf);
     bore_diameters_.reserve(nperf);
 
-    // COMPDAT handling
+    const WellConnections& connections = well_ecl_.getConnections();
-    const auto& connectionSet = well_ecl_.getConnections();
+    const std::size_t num_conns = connections.size();
-    CheckDistributedWellConnections checker(well_ecl_, parallel_well_info_);
+    int num_active_connections = 0;
-    for (size_t c=0; c<connectionSet.size(); c++) {
+    auto my_next_perf = perf_data_->begin();
-        const auto& connection = connectionSet.get(c);
+    for (std::size_t c = 0; c < num_conns; ++c) {
-        const int cell =
+        if (my_next_perf == perf_data_->end())
-                cartesian_to_compressed[connection.global_index()];
-        if (connection.state() != Connection::State::OPEN || cell >= 0)
         {
-            checker.connectionFound(c);
+            break;
         }
+        if (my_next_perf->ecl_index > c)
+        {
+            continue;
+        }
+        assert(my_next_perf->ecl_index == c);
+        const auto& connection = connections[c];
         if (connection.state() == Connection::State::OPEN) {
-
+            double radius = connection.rw();
-            if (cell >= 0) {
+            double re = connection.re(); // area equivalent radius of the grid block
-                double radius = connection.rw();
+            double perf_length = connection.connectionLength(); // the length of the well perforation
-                double re = connection.re(); // area equivalent radius of the grid block
+            const double repR = std::sqrt(re * radius);
-                double perf_length = connection.connectionLength(); // the length of the well perforation
+            perf_rep_radius_.push_back(repR);
-                const double repR = std::sqrt(re * radius);
+            perf_length_.push_back(perf_length);
-                perf_rep_radius_.push_back(repR);
+            bore_diameters_.push_back(2. * radius);
-                perf_length_.push_back(perf_length);
+            num_active_connections++;
-                bore_diameters_.push_back(2. * radius);
-            }
         }
+        ++my_next_perf;
     }
-    checker.checkAllConnectionsFound();
+    assert(my_next_perf == perf_data_->end());
+    assert(num_active_connections == nperf);
 }
 
 void WellInterfaceGeneric::setWsolvent(const double wsolvent)

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -93,7 +93,7 @@ public:
     void setVFPProperties(const VFPProperties* vfp_properties_arg);
     void setGuideRate(const GuideRate* guide_rate_arg);
     void setWellEfficiencyFactor(const double efficiency_factor);
-    void setRepRadiusPerfLength(const std::vector<int>& cartesian_to_compressed);
+    void setRepRadiusPerfLength();
     void setWsolvent(const double wsolvent);
     void setDynamicThpLimit(const double thp_limit);
     void updatePerforatedCell(std::vector<bool>& is_cell_perforated);