Added possibility to communicate values from perforations below.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -467,7 +467,7 @@ namespace Opm {
         // Clear the communication data structures for above values.
         for(auto&& pinfo : local_parallel_well_info_)
         {
-            pinfo->clearCommunicateAbove();
+            pinfo->clearCommunicateAboveBelow();
         }
     }
 

opm/simulators/wells/ParallelWellInfo.cpp

@@ -43,13 +43,13 @@ struct CommPolicy<double*>
 
 namespace Opm
 {
-CommunicateAbove::CommunicateAbove([[maybe_unused]] const Communication& comm)
+CommunicateAboveBelow::CommunicateAboveBelow([[maybe_unused]] const Communication& comm)
 #if HAVE_MPI
     : comm_(comm), interface_(comm_)
 #endif
 {}
 
-void CommunicateAbove::clear()
+void CommunicateAboveBelow::clear()
 {
 #if HAVE_MPI
     above_indices_ = {};
@@ -60,7 +60,7 @@ void CommunicateAbove::clear()
     count_ = 0;
 }
 
-void CommunicateAbove::beginReset()
+void CommunicateAboveBelow::beginReset()
 {
     clear();
 #if HAVE_MPI
@@ -72,7 +72,7 @@ void CommunicateAbove::beginReset()
 #endif
 }
 
-void CommunicateAbove::endReset()
+void CommunicateAboveBelow::endReset()
 {
 #if HAVE_MPI
     if (comm_.size() > 1)
@@ -104,9 +104,9 @@ struct CopyGatherScatter
 };
 
 
-std::vector<double> CommunicateAbove::communicate(double first_above,
+std::vector<double> CommunicateAboveBelow::communicateAbove(double first_above,
-                                                  const double* current,
+                                                            const double* current,
-                                                  std::size_t size)
+                                                            std::size_t size)
 {
     std::vector<double> above(size, first_above);
 
@@ -131,10 +131,37 @@ std::vector<double> CommunicateAbove::communicate(double first_above,
     }
     return above;
 }
+std::vector<double> CommunicateAboveBelow::communicateBelow(double last_below,
+                                                            const double* current,
+                                                            std::size_t size)
+{
+    std::vector<double> below(size, last_below);
 
-void CommunicateAbove::pushBackEclIndex([[maybe_unused]] int above,
+#if HAVE_MPI
-                                        [[maybe_unused]] int current,
+    if (comm_.size() > 1)
-                                        [[maybe_unused]] bool isOwner)
+    {
+        auto belowData = below.data();
+        // Ugly const_cast needed as my compiler says, that
+        // passing const double*& and double* as parameter is
+        // incompatible with function decl template<Data> backward(Data&, const Data&)
+        // That would need the first argument to be double* const&
+        communicator_.backward<CopyGatherScatter>(belowData, const_cast<double*>(current));
+    }
+    else
+#endif
+    {
+        if (below.size() > 1)
+        {
+            // No comunication needed, just copy.
+            std::copy(current+1, current + below.size(), below.begin());
+        }
+    }
+    return below;
+}
+
+void CommunicateAboveBelow::pushBackEclIndex([[maybe_unused]] int above,
+                                             [[maybe_unused]] int current,
+                                             [[maybe_unused]] bool isOwner)
 {
 #if HAVE_MPI
     if (comm_.size() > 1)
@@ -175,8 +202,8 @@ ParallelWellInfo::ParallelWellInfo(const std::string& name,
     : name_(name), hasLocalCells_ (hasLocalCells),
       isOwner_(true), rankWithFirstPerf_(-1),
       comm_(new Communication(Dune::MPIHelper::getLocalCommunicator())),
-      commAbove_(new CommunicateAbove(*comm_))
+      commAboveBelow_(new CommunicateAboveBelow(*comm_))
-    {}
+{}
 
 ParallelWellInfo::ParallelWellInfo(const std::pair<std::string,bool>& well_info,
                                    [[maybe_unused]] Communication allComm)
@@ -191,7 +218,7 @@ ParallelWellInfo::ParallelWellInfo(const std::pair<std::string,bool>& well_info,
 #else
     comm_.reset(new Communication(Dune::MPIHelper::getLocalCommunicator()));
 #endif
-    commAbove_.reset(new CommunicateAbove(*comm_));
+    commAboveBelow_.reset(new CommunicateAboveBelow(*comm_));
     isOwner_ = (comm_->rank() == 0);
 }
 
@@ -207,38 +234,53 @@ void ParallelWellInfo::communicateFirstPerforation(bool hasFirst)
 
 void ParallelWellInfo::pushBackEclIndex(int above, int current)
 {
-    commAbove_->pushBackEclIndex(above, current);
+    commAboveBelow_->pushBackEclIndex(above, current);
 }
 
 void ParallelWellInfo::beginReset()
 {
-    commAbove_->beginReset();
+    commAboveBelow_->beginReset();
 }
 
 
 void ParallelWellInfo::endReset()
 {
-    commAbove_->beginReset();
+    commAboveBelow_->beginReset();
 }
 
-void ParallelWellInfo::clearCommunicateAbove()
+void ParallelWellInfo::clearCommunicateAboveBelow()
 {
-    commAbove_->clear();
+    commAboveBelow_->clear();
 }
 
 std::vector<double> ParallelWellInfo::communicateAboveValues(double zero_value,
                                                              const double* current_values,
                                                              std::size_t size) const
 {
-    return commAbove_->communicate(zero_value, current_values,
+    return commAboveBelow_->communicateAbove(zero_value, current_values,
-                                   size);
+                                             size);
 }
 
 std::vector<double> ParallelWellInfo::communicateAboveValues(double zero_value,
                                                              const std::vector<double>& current_values) const
 {
-    return commAbove_->communicate(zero_value, current_values.data(),
+    return commAboveBelow_->communicateAbove(zero_value, current_values.data(),
-                                   current_values.size());
+                                             current_values.size());
+}
+
+std::vector<double> ParallelWellInfo::communicateBelowValues(double last_value,
+                                                             const double* current_values,
+                                                             std::size_t size) const
+{
+    return commAboveBelow_->communicateBelow(last_value, current_values,
+                                             size);
+}
+
+std::vector<double> ParallelWellInfo::communicateBelowValues(double last_value,
+                                                             const std::vector<double>& current_values) const
+{
+    return commAboveBelow_->communicateBelow(last_value, current_values.data(),
+                                             current_values.size());
 }
 
 bool operator<(const ParallelWellInfo& well1, const ParallelWellInfo& well2)
@@ -294,7 +336,7 @@ bool operator!=(const ParallelWellInfo& well, const std::pair<std::string, bool>
 }
 
 CheckDistributedWellConnections::CheckDistributedWellConnections(const Well& well,
-                                                               const ParallelWellInfo& info)
+                                                                 const ParallelWellInfo& info)
     : well_(well), pwinfo_(info)
 {
     foundConnections_.resize(well.getConnections().size(), 0);

opm/simulators/wells/ParallelWellInfo.hpp

@@ -33,9 +33,9 @@
 namespace Opm
 {
 
-/// \brief Class to facilitate getting values associated with the above perforation
+/// \brief Class to facilitate getting values associated with the above/below perforation
 ///
-class CommunicateAbove
+class CommunicateAboveBelow
 {
     enum Attribute {
       owner=1, overlap=2
@@ -53,7 +53,7 @@ class CommunicateAbove
 #endif
 
 public:
-    explicit CommunicateAbove(const Communication& comm);
+    explicit CommunicateAboveBelow(const Communication& comm);
     /// \brief Adds information about original index of the perforations in ECL Schedule.
     ///
     /// \warning Theses indices need to be push in the same order as they
@@ -81,10 +81,18 @@ public:
     /// \param current C-array of the values at the perforations
     /// \param size The size of the C-array and the returned vector
     /// \return a vector containing the values for the perforation above.
-    std::vector<double> communicate(double first_value,
+    std::vector<double> communicateAbove(double first_value,
-                                    const double* current,
+                                         const double* current,
-                                    std::size_t size);
+                                         std::size_t size);
 
+    /// \brief Creates an array of values for the perforation below.
+    /// \param first_value Value to use for above of the first perforation
+    /// \param current C-array of the values at the perforations
+    /// \param size The size of the C-array and the returned vector
+    /// \return a vector containing the values for the perforation above.
+    std::vector<double> communicateBelow(double first_value,
+                                         const double* current,
+                                         std::size_t size);
 private:
 #if HAVE_MPI
     Communication comm_;
@@ -157,15 +165,31 @@ public:
     /// \param current C-array of the values at the perforations
     /// \param size The size of the C-array and the returned vector
     /// \return a vector containing the values for the perforation above.
-    std::vector<double> communicateAboveValues(double zero_value,
+    std::vector<double> communicateAboveValues(double first_value,
                                                const double* current,
                                                std::size_t size) const;
 
     /// \brief Creates an array of values for the perforation above.
     /// \param first_value Value to use for above of the first perforation
     /// \param current vector of current values
-    std::vector<double> communicateAboveValues(double zero_value,
+    std::vector<double> communicateAboveValues(double first_value,
                                                const std::vector<double>& current) const;
+
+    /// \brief Creates an array of values for the perforation below.
+    /// \param last_value Value to use for below of the last perforation
+    /// \param current C-array of the values at the perforations
+    /// \param size The size of the C-array and the returned vector
+    /// \return a vector containing the values for the perforation above.
+    std::vector<double> communicateBelowValues(double last_value,
+                                               const double* current,
+                                               std::size_t size) const;
+
+    /// \brief Creates an array of values for the perforation above.
+    /// \param last_value Value to use for below of the last perforation
+    /// \param current vector of current values
+    std::vector<double> communicateBelowValues(double last_value,
+                                               const std::vector<double>& current) const;
+
     /// \brief Adds information about the ecl indices of the perforations.
     ///
     /// \warning Theses indices need to be push in the same order as they
@@ -210,7 +234,7 @@ public:
     }
 
     /// \brief Free data of communication data structures.
-    void clearCommunicateAbove();
+    void clearCommunicateAboveBelow();
 private:
 
     /// \brief Deleter that also frees custom MPI communicators
@@ -234,7 +258,7 @@ private:
     std::unique_ptr<Communication, DestroyComm> comm_;
 
     /// \brief used to communicate the values for the perforation above.
-    std::unique_ptr<CommunicateAbove> commAbove_;
+    std::unique_ptr<CommunicateAboveBelow> commAboveBelow_;
 };
 
 /// \brief Class checking that all connections are on active cells

tests/test_parallelwellinfo.cpp

@@ -153,64 +153,76 @@ BOOST_AUTO_TEST_CASE(ParallelWellComparison)
 
 }
 
-BOOST_AUTO_TEST_CASE(CommunicateAboveSelf)
+BOOST_AUTO_TEST_CASE(CommunicateAboveBelowSelf)
 {
     auto comm = Dune::MPIHelper::getLocalCommunicator();
-    Opm::CommunicateAbove commAbove{ comm };
+    Opm::CommunicateAboveBelow commAboveBelow{ comm };
     for(std::size_t count=0; count < 2; ++count)
     {
         std::vector<int> eclIndex = {0, 1, 2, 3, 7 , 8, 10, 11};
         std::vector<double> current(eclIndex.size());
         std::transform(eclIndex.begin(), eclIndex.end(), current.begin(),
                        [](double v){ return 1+10.0*v;});
-        commAbove.beginReset();
+        commAboveBelow.beginReset();
         for (std::size_t i = 0; i < current.size(); ++i)
         {
             if (i==0)
-                commAbove.pushBackEclIndex(-1, eclIndex[i]);
+                commAboveBelow.pushBackEclIndex(-1, eclIndex[i]);
             else
-                commAbove.pushBackEclIndex(eclIndex[i-1], eclIndex[i]);
+                commAboveBelow.pushBackEclIndex(eclIndex[i-1], eclIndex[i]);
         }
-        commAbove.endReset();
+        commAboveBelow.endReset();
-        auto above = commAbove.communicate(-10.0, current.data(), current.size());
+        auto above = commAboveBelow.communicateAbove(-10.0, current.data(), current.size());
         BOOST_CHECK(above[0]==-10.0);
         BOOST_CHECK(above.size() == current.size());
         auto a = above.begin()+1;
         std::for_each(current.begin(), current.begin() + (current.size()-1),
                       [&a](double v){ BOOST_CHECK(*(a++) == v);});
+        auto below = commAboveBelow.communicateBelow(-10.0, current.data(), current.size());
+        BOOST_CHECK(below.back() == -10.0);
+        BOOST_CHECK(below.size() == current.size());
+        auto b = below.begin();
+        std::for_each(current.begin()+1, current.end(),
+                      [&b](double v){ BOOST_CHECK(*(b++) == v);});
     }
 }
 
 
-BOOST_AUTO_TEST_CASE(CommunicateAboveSelf1)
+BOOST_AUTO_TEST_CASE(CommunicateAboveBelowSelf1)
 {
     auto comm = Dune::MPIHelper::getLocalCommunicator();
-    Opm::CommunicateAbove commAbove{ comm };
+    Opm::CommunicateAboveBelow commAboveBelow{ comm };
     for(std::size_t count=0; count < 2; ++count)
     {
         std::vector<int> eclIndex = {0};
         std::vector<double> current(eclIndex.size());
         std::transform(eclIndex.begin(), eclIndex.end(), current.begin(),
                        [](double v){ return 1+10.0*v;});
-        commAbove.beginReset();
+        commAboveBelow.beginReset();
         for (std::size_t i = 0; i < current.size(); ++i)
         {
             if (i==0)
-                commAbove.pushBackEclIndex(-1, eclIndex[i]);
+                commAboveBelow.pushBackEclIndex(-1, eclIndex[i]);
             else
-                commAbove.pushBackEclIndex(eclIndex[i-1], eclIndex[i]);
+                commAboveBelow.pushBackEclIndex(eclIndex[i-1], eclIndex[i]);
         }
-        commAbove.endReset();
+        commAboveBelow.endReset();
-        auto above = commAbove.communicate(-10.0, current.data(), current.size());
+        auto above = commAboveBelow.communicateAbove(-10.0, current.data(), current.size());
         BOOST_CHECK(above[0]==-10.0);
         BOOST_CHECK(above.size() == current.size());
         auto a = above.begin()+1;
         std::for_each(current.begin(), current.begin() + (current.size()-1),
                       [&a](double v){ BOOST_CHECK(*(a++) == v);});
+        auto below = commAboveBelow.communicateBelow(-10.0, current.data(), current.size());
+        BOOST_CHECK(below.back() == -10.0);
+        BOOST_CHECK(below.size() == current.size());
+        auto b = below.begin();
+        std::for_each(current.begin()+1, current.end(),
+                      [&b](double v){ BOOST_CHECK(*(b++) == v);});
     }
 }
 
-BOOST_AUTO_TEST_CASE(CommunicateAboveParalle)
+BOOST_AUTO_TEST_CASE(CommunicateAboveBelowParallel)
 {
     using MPIComm = typename Dune::MPIHelper::MPICommunicator;
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7)
@@ -220,7 +232,7 @@ BOOST_AUTO_TEST_CASE(CommunicateAboveParalle)
 #endif
     auto comm = Communication(Dune::MPIHelper::getCommunicator());
 
-    Opm::CommunicateAbove commAbove{ comm };
+    Opm::CommunicateAboveBelow commAboveBelow{ comm };
     for(std::size_t count=0; count < 2; ++count)
     {
         std::vector<int> globalEclIndex = {0, 1, 2, 3, 7 , 8, 10, 11};
@@ -244,23 +256,23 @@ BOOST_AUTO_TEST_CASE(CommunicateAboveParalle)
 
         std::vector<double> current(3);
 
-        commAbove.beginReset();
+        commAboveBelow.beginReset();
         for (std::size_t i = 0; i < current.size(); ++i)
         {
             auto gi = comm.rank() + comm.size() * i;
 
             if (gi==0)
             {
-                commAbove.pushBackEclIndex(-1, globalEclIndex[gi]);
+                commAboveBelow.pushBackEclIndex(-1, globalEclIndex[gi]);
             }
             else
             {
-                commAbove.pushBackEclIndex(globalEclIndex[gi-1], globalEclIndex[gi]);
+                commAboveBelow.pushBackEclIndex(globalEclIndex[gi-1], globalEclIndex[gi]);
             }
             current[i] = globalCurrent[gi];
         }
-        commAbove.endReset();
+        commAboveBelow.endReset();
-        auto above = commAbove.communicate(-10.0, current.data(), current.size());
+        auto above = commAboveBelow.communicateAbove(-10.0, current.data(), current.size());
         if (comm.rank() == 0)
             BOOST_CHECK(above[0]==-10.0);
 
@@ -274,5 +286,19 @@ BOOST_AUTO_TEST_CASE(CommunicateAboveParalle)
                 BOOST_CHECK(above[i]==globalCurrent[gi-1]);
             }
         }
+        auto below = commAboveBelow.communicateBelow(-10.0, current.data(), current.size());
+        if (comm.rank() == comm.size() - 1)
+            BOOST_CHECK(below.back() == -10.0);
+
+        BOOST_CHECK(below.size() == current.size());
+
+        for (std::size_t i = 0; i < current.size(); ++i)
+        {
+            auto gi = comm.rank() + comm.size() * i;
+            if (gi < globalCurrent.size() - 1)
+            {
+                BOOST_CHECK(below[i]==globalCurrent[gi+1]);
+            }
+        }
     }
 }