Merge pull request #3050 from joakim-hove/actionx-welpi

Support combination of ACTIONX and WELPI

ebos/eclpeacemanwell.hh

@@ -1410,6 +1410,33 @@ protected:
                 / rhoWaterSurface;
     }
 
+    void
+    updateEclWell(int, int)
+    {
+        throw std::logic_error("updateEclWell() method not implemented for class eclpeacemanwell");
+    }
+
+
+    void
+    updateEclWell(int, const std::string&) {
+        throw std::logic_error("updateEclWell() method not implemented for class eclpeacemanwell");
+    }
+
+
+    double
+    wellPI(int) const
+    {
+        throw std::logic_error("wellPI() method not implemented for class eclpeacemanwell");
+    }
+
+    double
+    wellPI(const std::string& ) const
+    {
+        throw std::logic_error("wellPI() method not implemented for class eclpeacemanwell");
+    }
+
+
+
     /*!
      * \brief Compute the volumetric phase rate of the complete well given a bottom hole
      *        pressure.

ebos/eclproblem.hh

@@ -94,6 +94,7 @@
 #include <opm/parser/eclipse/EclipseState/Tables/Eqldims.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Action/ActionContext.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Action/ActionX.hpp>
 #include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RockwnodTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/OverburdTable.hpp>
@@ -1291,6 +1292,53 @@ public:
     }
 
 
+    std::unordered_map<std::string, double> fetchWellPI(int reportStep,
+                                                        const Opm::Action::ActionX& action,
+                                                        const Opm::Schedule& schedule,
+                                                        const std::vector<std::string>& matching_wells) {
+
+        auto wellpi_wells = action.wellpi_wells(WellMatcher(schedule[reportStep].well_order(),
+                                                            schedule[reportStep].wlist_manager()),
+                                                matching_wells);
+
+        if (wellpi_wells.empty())
+            return {};
+
+        const auto num_wells = schedule[reportStep].well_order().size();
+        std::vector<double> wellpi_vector(num_wells);
+        for (const auto& wname : wellpi_wells) {
+            if (this->wellModel_.hasWell(wname)) {
+                const auto& well = schedule.getWell( wname, reportStep );
+                wellpi_vector[well.seqIndex()] = this->wellModel_.wellPI(wname);
+            }
+        }
+
+        const auto& comm = this->simulator().vanguard().grid().comm();
+        if (comm.size() > 1) {
+            std::vector<double> wellpi_buffer(num_wells * comm.size());
+            comm.gather( wellpi_vector.data(), wellpi_buffer.data(), num_wells, 0 );
+            if (comm.rank() == 0) {
+                for (int rank=1; rank < comm.size(); rank++) {
+                    for (std::size_t well_index=0; well_index < num_wells; well_index++) {
+                        const auto global_index = rank*num_wells + well_index;
+                        const auto value = wellpi_buffer[global_index];
+                        if (value != 0)
+                            wellpi_vector[well_index] = value;
+                    }
+                }
+            }
+            comm.broadcast(wellpi_vector.data(), wellpi_vector.size(), 0);
+        }
+
+        std::unordered_map<std::string, double> wellpi;
+        for (const auto& wname : wellpi_wells) {
+            const auto& well = schedule.getWell( wname, reportStep );
+            wellpi[wname] = wellpi_vector[ well.seqIndex() ];
+        }
+        return wellpi;
+    }
+
+
     void applyActions(int reportStep,
                       double sim_time,
                       Opm::EclipseState& ecl_state,
@@ -1330,8 +1378,17 @@ public:
                 }
                 std::string msg = "The action: " + action->name() + " evaluated to true at " + ts + " wells: " + wells_string;
                 Opm::OpmLog::info(msg);
-                schedule.applyAction(reportStep, std::chrono::system_clock::from_time_t(simTime), *action, actionResult, {});
+
+                const auto& wellpi = this->fetchWellPI(reportStep, *action, schedule, matching_wells);
+
+                schedule.applyAction(reportStep, std::chrono::system_clock::from_time_t(simTime), *action, actionResult, wellpi);
                 actionState.add_run(*action, simTime);
+
+                for ( const auto& [wname, _] : wellpi) {
+                    (void)_;
+                    if (this->wellModel_.hasWell(wname))
+                        this->wellModel_.updateEclWell(reportStep, wname);
+                }
             } else {
                 std::string msg = "The action: " + action->name() + " evaluated to false at " + ts;
                 Opm::OpmLog::info(msg);

ebos/eclwellmanager.hh

@@ -620,6 +620,35 @@ public:
     }
 
 
+    void
+    updateEclWell(int, int)
+    {
+        throw std::logic_error("wellPI() method not implemented for class eclwellmanager");
+    }
+
+
+    void
+    updateEclWell(int, const std::string&) {
+        throw std::logic_error("wellPI() method not implemented for class eclwellmanager");
+    }
+
+
+    double
+    wellPI(int ) const
+    {
+        throw std::logic_error("wellPI() method not implemented for class eclwellmanager");
+    }
+
+    double
+    wellPI(const std::string& ) const
+    {
+        throw std::logic_error("wellPI() method not implemented for class eclwellmanager");
+    }
+
+
+
+
+
 protected:
     bool wellTopologyChanged_(const Opm::EclipseState& eclState OPM_UNUSED,
                               const Opm::Schedule& schedule,
@@ -751,6 +780,7 @@ protected:
         }
     }
 
+
     void updateWellParameters_(unsigned reportStepIdx, const WellConnectionsMap& wellConnections)
     {
         const auto& deckSchedule = simulator_.vanguard().schedule();

opm/simulators/wells/BlackoilWellModel.hpp

@@ -267,8 +267,11 @@ namespace Opm {
             /// Returns true if the well was actually found and shut.
             bool forceShutWellByNameIfPredictionMode(const std::string& wellname, const double simulation_time);
 
-            double wellPI(int well_index) const;
+            void updateEclWell(int timeStepIdx, int well_index);
+            void updateEclWell(int timeStepIdx, const std::string& wname);
             bool hasWell(const std::string& wname);
+            double wellPI(int well_index) const;
+            double wellPI(const std::string& well_name) const;
         protected:
             Simulator& ebosSimulator_;
 
@@ -482,7 +485,6 @@ namespace Opm {
 
             void setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent);
 
-            void updateEclWell(int timeStepIdx, int well_index);
             void runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger);
 
             void assignWellGuideRates(data::Wells& wsrpt) const;

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -2583,6 +2583,17 @@ namespace Opm {
         this->prod_index_calc_[well_index].reInit(well);
     }
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    updateEclWell(int timeStepIdx, const std::string& wname) {
+        auto well_iter = std::find_if( this->wells_ecl_.begin(), this->wells_ecl_.end(), [wname] (const auto& well) -> bool { return well.name() == wname;});
+        if (well_iter == this->wells_ecl_.end())
+            throw std::logic_error("Could not find well: " + wname);
+
+        auto well_index = std::distance( this->wells_ecl_.begin(), well_iter );
+        this->updateEclWell(timeStepIdx, well_index);
+    }
 
 
     template<typename TypeTag>
@@ -2619,6 +2630,19 @@ namespace Opm {
         }
     }
 
+    template<typename TypeTag>
+    double
+    BlackoilWellModel<TypeTag>::
+    wellPI(const std::string& well_name) const
+    {
+        auto well_iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(), [&well_name](const Well& well) { return well.name() == well_name; });
+        if (well_iter == this->wells_ecl_.end())
+            throw std::logic_error("Could not find well: " + well_name);
+
+        auto well_index = std::distance( this->wells_ecl_.begin(), well_iter );
+        return this->wellPI(well_index);
+    }
+
 
     template<typename TypeTag>
     void