Merge pull request #2957 from joakim-hove/wbp-collect

Wbp collect

ebos/collecttoiorank.hh

@@ -655,14 +655,70 @@ public:
 
     };
 
+    class PackUnPackWBPData : public P2PCommunicatorType::DataHandleInterface
+    {
+        const std::map<std::size_t, double>& localWBPData_;
+        std::map<std::size_t, double>& globalWBPValues_;
+
+    public:
+        PackUnPackWBPData(const std::map<std::size_t, double>& localWBPData,
+                          std::map<std::size_t, double>& globalWBPValues,
+                          bool isIORank)
+            : localWBPData_(localWBPData)
+            , globalWBPValues_(globalWBPValues)
+        {
+            if (isIORank) {
+                MessageBufferType buffer;
+                pack(0, buffer);
+
+                // pass a dummyLink to satisfy virtual class
+                int dummyLink = -1;
+                unpack(dummyLink, buffer);
+            }
+        }
+
+        // pack all data associated with link
+        void pack(int link, MessageBufferType& buffer)
+        {
+            // we should only get one link
+            if (link != 0)
+                throw std::logic_error("link in method pack is not 0 as expected");
+
+            // write all block data
+            unsigned int size = localWBPData_.size();
+            buffer.write(size);
+            for (const auto& [global_index, wbp_value] : localWBPData_) {
+                buffer.write(global_index);
+                buffer.write(wbp_value);
+            }
+        }
+
+        // unpack all data associated with link
+        void unpack(int /*link*/, MessageBufferType& buffer)
+        {
+            // read all block data
+            unsigned int size = 0;
+            buffer.read(size);
+            for (size_t i = 0; i < size; ++i) {
+                std::size_t idx;
+                double data;
+                buffer.read(idx);
+                buffer.read(data);
+                globalWBPValues_[idx] = data;
+            }
+        }
+
+    };
     // gather solution to rank 0 for EclipseWriter
     void collect(const Opm::data::Solution& localCellData,
                  const std::map<std::pair<std::string, int>, double>& localBlockData,
+                 const std::map<std::size_t, double>& localWBPData,
                  const Opm::data::Wells& localWellData,
                  const Opm::data::GroupAndNetworkValues& localGroupAndNetworkData)
     {
         globalCellData_ = {};
         globalBlockData_.clear();
+        globalWBPData_.clear();
         globalWellData_.clear();
         globalGroupAndNetworkData_.clear();
 
@@ -703,12 +759,17 @@ public:
             this->isIORank()
         };
 
+        PackUnPackWBPData packUnpackWBPData {
+            localWBPData,
+            this->globalWBPData_,
+            this->isIORank()
+        };
+
         toIORankComm_.exchange(packUnpackCellData);
         toIORankComm_.exchange(packUnpackWellData);
         toIORankComm_.exchange(packUnpackGroupAndNetworkData);
         toIORankComm_.exchange(packUnpackBlockData);
-
+        toIORankComm_.exchange(packUnpackWBPData);
-
 
 #ifndef NDEBUG
         // mkae sure every process is on the same page
@@ -716,6 +777,9 @@ public:
 #endif
     }
 
+    const std::map<std::size_t, double>& globalWBPData() const
+    { return this->globalWBPData_; }
+
     const std::map<std::pair<std::string, int>, double>& globalBlockData() const
     { return globalBlockData_; }
 
@@ -772,6 +836,7 @@ protected:
     std::vector<int> globalRanks_;
     Opm::data::Solution globalCellData_;
     std::map<std::pair<std::string, int>, double> globalBlockData_;
+    std::map<std::size_t, double> globalWBPData_;
     Opm::data::Wells globalWellData_;
     Opm::data::GroupAndNetworkValues globalGroupAndNetworkData_;
     std::vector<int> localIdxToGlobalIdx_;

ebos/ecloutputblackoilmodule.hh

@@ -195,7 +195,7 @@ class EclOutputBlackOilModule
 
 public:
     template<class CollectDataToIORankType>
-    EclOutputBlackOilModule(const Simulator& simulator, const CollectDataToIORankType& collectToIORank)
+    EclOutputBlackOilModule(const Simulator& simulator, const std::vector<std::size_t>& wbp_index_list, const CollectDataToIORankType& collectToIORank)
         : simulator_(simulator)
     {
         const Opm::SummaryConfig summaryConfig = simulator_.vanguard().summaryConfig();
@@ -225,6 +225,11 @@ public:
             }
         }
 
+        for (const auto& global_index : wbp_index_list) {
+            if (collectToIORank.isCartIdxOnThisRank(global_index - 1))
+                this->wbpData_[global_index] = 0.0;
+        }
+
         forceDisableFipOutput_ = EWOMS_GET_PARAM(TypeTag, bool, ForceDisableFluidInPlaceOutput);
     }
 
@@ -837,6 +842,8 @@ public:
             if (gasConnectionSaturations_.count(cartesianIdx) > 0) {
                 gasConnectionSaturations_[cartesianIdx] = Opm::getValue(fs.saturation(gasPhaseIdx));
             }
+            if (this->wbpData_.count(cartesianIdx) > 0)
+                this->wbpData_[cartesianIdx] = Opm::getValue(fs.pressure(oilPhaseIdx));
 
             // tracers
             const auto& tracerModel = simulator_.problem().tracerModel();
@@ -1852,6 +1859,10 @@ public:
         return 0;
     }
 
+    const std::map<std::size_t, double>& getWBPData() const {
+        return this->wbpData_;
+    }
+
     const std::map<std::pair<std::string, int>, double>& getBlockData()
     { return blockData_; }
 
@@ -2365,6 +2376,7 @@ private:
     ScalarBuffer hydrocarbonPoreVolume_;
     ScalarBuffer pressureTimesPoreVolume_;
     ScalarBuffer pressureTimesHydrocarbonVolume_;
+    std::map<std::size_t , double> wbpData_;
     std::map<std::pair<std::string, int>, double> blockData_;
     std::map<size_t, Scalar> oilConnectionPressures_;
     std::map<size_t, Scalar> waterConnectionSaturations_;

ebos/eclwriter.hh

@@ -191,15 +191,26 @@ public:
     EclWriter(Simulator& simulator)
         : simulator_(simulator)
         , collectToIORank_(simulator_.vanguard())
-        , eclOutputModule_(simulator, collectToIORank_)
     {
+        std::vector<std::size_t> wbp_index_list;
         if (collectToIORank_.isIORank()) {
+            const auto& schedule = simulator_.vanguard().schedule();
             eclIO_.reset(new Opm::EclipseIO(simulator_.vanguard().eclState(),
                                             Opm::UgGridHelpers::createEclipseGrid(globalGrid(), simulator_.vanguard().eclState().getInputGrid()),
-                                            simulator_.vanguard().schedule(),
+                                            schedule,
                                             simulator_.vanguard().summaryConfig()));
-        }
 
+            const auto& wbp_calculators = eclIO_->summary().wbp_calculators( schedule.size() - 1 );
+            wbp_index_list = wbp_calculators.index_list();
+        }
+        if (collectToIORank_.isParallel()) {
+            const auto& comm = simulator_.vanguard().grid().comm();
+            unsigned long size = wbp_index_list.size();
+            comm.broadcast(&size, 1, collectToIORank_.ioRank);
+            if (!collectToIORank_.isIORank())
+                wbp_index_list.resize( size );
+            comm.broadcast(wbp_index_list.data(), size, collectToIORank_.ioRank);
+        }
         // create output thread if enabled and rank is I/O rank
         // async output is enabled by default if pthread are enabled
         bool enableAsyncOutput = EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput);
@@ -207,6 +218,8 @@ public:
         if (enableAsyncOutput && collectToIORank_.isIORank())
             numWorkerThreads = 1;
         taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
+
+        this->eclOutputModule_ = std::make_unique<EclOutputBlackOilModule<TypeTag>>(simulator, wbp_index_list, this->collectToIORank_);
     }
 
     ~EclWriter()
@@ -241,7 +254,6 @@ public:
 
     void evalSummaryState(bool isSubStep)
     {
-        PAvgCalculatorCollection wbp_calculators;
         const int reportStepNum = simulator_.episodeIndex() + 1;
         /*
           The summary data is not evaluated for timestep 0, that is
@@ -279,23 +291,35 @@ public:
         this->prepareLocalCellData(isSubStep, reportStepNum);
 
         if (collectToIORank_.isParallel())
-            collectToIORank_.collect({}, eclOutputModule_.getBlockData(),
+            collectToIORank_.collect({},
-                                     localWellData, localGroupAndNetworkData);
+                                     eclOutputModule_->getBlockData(),
+                                     eclOutputModule_->getWBPData(),
+                                     localWellData,
+                                     localGroupAndNetworkData);
 
         std::map<std::string, double> miscSummaryData;
         std::map<std::string, std::vector<double>> regionData;
-        auto inplace = eclOutputModule_.outputFipLog(miscSummaryData, regionData, isSubStep);
+        auto inplace = eclOutputModule_->outputFipLog(miscSummaryData, regionData, isSubStep);
 
         bool forceDisableProdOutput = false;
         bool forceDisableInjOutput = false;
         bool forceDisableCumOutput = false;
-        eclOutputModule_.outputProdLog(reportStepNum, isSubStep, forceDisableProdOutput);
+        eclOutputModule_->outputProdLog(reportStepNum, isSubStep, forceDisableProdOutput);
-        eclOutputModule_.outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput);
+        eclOutputModule_->outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput);
-        eclOutputModule_.outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput);
+        eclOutputModule_->outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput);
+
 
         std::vector<char> buffer;
         if (this->collectToIORank_.isIORank()) {
             const auto& summary = eclIO_->summary();
+            auto wbp_calculators = summary.wbp_calculators(reportStepNum);
+            const auto& wbpData
+                = this->collectToIORank_.isParallel()
+                ? this->collectToIORank_.globalWBPData()
+                : this->eclOutputModule_->getWBPData();
+
+            for (const auto& [global_index, pressure] : wbpData)
+                wbp_calculators.add_pressure( global_index, pressure );
 
             // Add TCPU
             if (totalCpuTime != 0.0) {
@@ -316,7 +340,7 @@ public:
             const auto& blockData
                 = this->collectToIORank_.isParallel()
                 ? this->collectToIORank_.globalBlockData()
-                : this->eclOutputModule_.getBlockData();
+                : this->eclOutputModule_->getBlockData();
 
             summary.eval(summaryState(),
                          reportStepNum,
@@ -324,7 +348,7 @@ public:
                          wellData,
                          groupAndNetworkData,
                          miscSummaryData,
-                         eclOutputModule_.initialInplace(),
+                         eclOutputModule_->initialInplace(),
                          inplace,
                          wbp_calculators,
                          regionData,
@@ -366,7 +390,7 @@ public:
         const int reportStepNum = simulator_.episodeIndex() + 1;
 
         this->prepareLocalCellData(isSubStep, reportStepNum);
-        this->eclOutputModule_.outputErrorLog();
+        this->eclOutputModule_->outputErrorLog();
 
         // output using eclWriter if enabled
         auto localWellData = simulator_.problem().wellModel().wellData();
@@ -375,15 +399,18 @@ public:
 
         Opm::data::Solution localCellData = {};
         if (! isSubStep) {
-            this->eclOutputModule_.assignToSolution(localCellData);
+            this->eclOutputModule_->assignToSolution(localCellData);
 
             // add cell data to perforations for Rft output
-            this->eclOutputModule_.addRftDataToWells(localWellData, reportStepNum);
+            this->eclOutputModule_->addRftDataToWells(localWellData, reportStepNum);
         }
 
         if (this->collectToIORank_.isParallel()) {
-            collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(),
+            collectToIORank_.collect(localCellData,
-                                     localWellData, localGroupAndNetworkData);
+                                     eclOutputModule_->getBlockData(),
+                                     eclOutputModule_->getWBPData(),
+                                     localWellData,
+                                     localGroupAndNetworkData);
         }
 
         if (this->collectToIORank_.isIORank()) {
@@ -427,7 +454,7 @@ public:
 
         const auto& gridView = simulator_.vanguard().gridView();
         unsigned numElements = gridView.size(/*codim=*/0);
-        eclOutputModule_.allocBuffers(numElements, restartStepIdx, /*isSubStep=*/false, /*log=*/false, /*isRestart*/ true);
+        eclOutputModule_->allocBuffers(numElements, restartStepIdx, /*isSubStep=*/false, /*log=*/false, /*isRestart*/ true);
 
         {
             Opm::SummaryState& summaryState = simulator_.vanguard().summaryState();
@@ -436,7 +463,7 @@ public:
                                                      gridView.grid().comm());
             for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
                 unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
-                eclOutputModule_.setRestart(restartValues.solution, elemIdx, globalIdx);
+                eclOutputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
             }
 
             if (inputThpres.active()) {
@@ -459,7 +486,7 @@ public:
     {}
 
     const EclOutputBlackOilModule<TypeTag>& eclOutputModule() const
-    { return eclOutputModule_; }
+    { return *eclOutputModule_; }
 
     Scalar restartTimeStepSize() const
     { return restartTimeStepSize_; }
@@ -724,7 +751,7 @@ private:
         const int numElements = gridView.size(/*codim=*/0);
         const bool log = collectToIORank_.isIORank();
 
-        eclOutputModule_.allocBuffers(numElements, reportStepNum,
+        eclOutputModule_->allocBuffers(numElements, reportStepNum,
                                       isSubStep, log, /*isRestart*/ false);
 
         ElementContext elemCtx(simulator_);
@@ -737,7 +764,7 @@ private:
             elemCtx.updatePrimaryStencil(elem);
             elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
 
-            eclOutputModule_.processElement(elemCtx);
+            eclOutputModule_->processElement(elemCtx);
         }
     }
 
@@ -791,7 +818,7 @@ private:
 
     Simulator& simulator_;
     CollectDataToIORankType collectToIORank_;
-    EclOutputBlackOilModule<TypeTag> eclOutputModule_;
+    std::unique_ptr<EclOutputBlackOilModule<TypeTag>> eclOutputModule_;
     std::unique_ptr<Opm::EclipseIO> eclIO_;
     std::unique_ptr<TaskletRunner> taskletRunner_;
     Scalar restartTimeStepSize_;

tests/test_ecl_output.cc

@@ -149,7 +149,7 @@ BOOST_AUTO_TEST_CASE(Summary)
     using Vanguard = Opm::GetPropType<TypeTag, Opm::Properties::Vanguard>;
     typedef Opm::CollectDataToIORank< Vanguard > CollectDataToIORankType;
     CollectDataToIORankType collectToIORank(simulator->vanguard());
-    Opm::EclOutputBlackOilModule<TypeTag> eclOutputModule(*simulator, collectToIORank);
+    Opm::EclOutputBlackOilModule<TypeTag> eclOutputModule(*simulator, {}, collectToIORank);
 
     typedef Opm::EclWriter<TypeTag> EclWriterType;
     // create the actual ECL writer