Extract functionality to update target rates to SingleWellState

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1239,7 +1239,7 @@ updateEclWells(const int timeStepIdx,
         }
         auto& ws = this->wellState().well(well_index);
 
-        this->wellState().updateStatus(well_index, well.getStatus());
+        ws.updateStatus( well.getStatus() );
         ws.reset_connection_factors(pd);
         this->prod_index_calc_[well_index].reInit(well);
     }

opm/simulators/wells/PerfData.cpp

@@ -24,8 +24,9 @@ namespace Opm
 {
 
 
-PerfData::PerfData(std::size_t num_perf, bool injector_, std::size_t num_phases)
+PerfData::PerfData(std::size_t num_perf, double pressure_first_connection_, bool injector_, std::size_t num_phases)
     : injector(injector_)
+    , pressure_first_connection(pressure_first_connection_)
     , pressure(num_perf)
     , rates(num_perf)
     , phase_rates(num_perf * num_phases)
@@ -61,6 +62,7 @@ bool PerfData::try_assign(const PerfData& other) {
     if (this->injector != other.injector)
         return false;
 
+    this->pressure_first_connection = other.pressure_first_connection;
     this->pressure = other.pressure;
     this->rates = other.rates;
     this->phase_rates = other.phase_rates;

opm/simulators/wells/PerfData.hpp

@@ -32,12 +32,13 @@ private:
     bool injector;
 
 public:
-    PerfData(std::size_t num_perf, bool injector_, std::size_t num_phases);
+    PerfData(std::size_t num_perf, double pressure_first_connection_, bool injector_, std::size_t num_phases);
     std::size_t size() const;
     bool empty() const;
     bool try_assign(const PerfData& other);
 
 
+    double pressure_first_connection;
     std::vector<double> pressure;
     std::vector<double> rates;
     std::vector<double> phase_rates;

opm/simulators/wells/SingleWellState.cpp

@@ -22,16 +22,23 @@
 
 namespace Opm {
 
-SingleWellState::SingleWellState(const std::string& name_, const ParallelWellInfo& pinfo, bool is_producer, const std::vector<PerforationData>& perf_input, std::size_t num_phases, double temp)
+SingleWellState::SingleWellState(const std::string& name_,
+                                 const ParallelWellInfo& pinfo,
+                                 bool is_producer,
+                                 double pressure_first_connection,
+                                 const std::vector<PerforationData>& perf_input,
+                                 const PhaseUsage& pu_,
+                                 double temp)
     : name(name_)
     , parallel_info(pinfo)
     , producer(is_producer)
+    , pu(pu_)
     , temperature(temp)
-    , well_potentials(num_phases)
-    , productivity_index(num_phases)
-    , surface_rates(num_phases)
-    , reservoir_rates(num_phases)
-    , perf_data(perf_input.size(), !is_producer, num_phases)
+    , well_potentials(pu_.num_phases)
+    , productivity_index(pu_.num_phases)
+    , surface_rates(pu_.num_phases)
+    , reservoir_rates(pu_.num_phases)
+    , perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu_.num_phases)
 {
     for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
         this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
@@ -143,7 +150,110 @@ double SingleWellState::sum_solvent_rates() const {
     return this->sum_connection_rates(this->perf_data.solvent_rates);
 }
 
+
+void SingleWellState::update_producer_targets(const Well& ecl_well, const SummaryState& st) {
+    const double bhp_safety_factor = 0.99;
+    const auto& prod_controls = ecl_well.productionControls(st);
+    if (prod_controls.hasControl(Well::ProducerCMode::THP))
+        this->thp = prod_controls.thp_limit;
+
+
+    auto cmode_is_bhp = (prod_controls.cmode == Well::ProducerCMode::BHP);
+    auto bhp_limit = prod_controls.bhp_limit;
+
+    if (ecl_well.getStatus() == Well::Status::STOP) {
+        if (cmode_is_bhp)
+            this->bhp = bhp_limit;
+        else
+            this->bhp = this->perf_data.pressure_first_connection;
+
+        return;
+    }
+
+    if (prod_controls.cmode == Well::ProducerCMode::GRUP) {
+        this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor;
+        return;
+    }
+
+    switch (prod_controls.cmode) {
+    case Well::ProducerCMode::ORAT:
+        assert(this->pu.phase_used[BlackoilPhases::Liquid]);
+        this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
+        break;
+    case Well::ProducerCMode::WRAT:
+        assert(this->pu.phase_used[BlackoilPhases::Aqua]);
+        this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
+        break;
+    case Well::ProducerCMode::GRAT:
+        assert(this->pu.phase_used[BlackoilPhases::Vapour]);
+        this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
+        break;
+    default:
+        // Keep zero init.
+        break;
+    }
+
+    if (cmode_is_bhp)
+        this->bhp = bhp_limit;
+    else
+        this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor;
+
+}
+
+void SingleWellState::update_injector_targets(const Well& ecl_well, const SummaryState& st) {
+    const double bhp_safety_factor = 1.01;
+    const auto& inj_controls = ecl_well.injectionControls(st);
+
+    if (inj_controls.hasControl(Well::InjectorCMode::THP))
+        this->thp = inj_controls.thp_limit;
+
+    auto cmode_is_bhp = (inj_controls.cmode == Well::InjectorCMode::BHP);
+    auto bhp_limit = inj_controls.bhp_limit;
+
+    if (ecl_well.getStatus() == Well::Status::STOP) {
+        if (cmode_is_bhp)
+            this->bhp = bhp_limit;
+        else
+            this->bhp = this->perf_data.pressure_first_connection;
+
+        return;
+    }
+
+
+    if (inj_controls.cmode == Well::InjectorCMode::GRUP) {
+        this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor;
+        return;
+    }
+
+    if (inj_controls.cmode == Well::InjectorCMode::RATE) {
+        auto inj_surf_rate = inj_controls.surface_rate;
+        switch (inj_controls.injector_type) {
+        case InjectorType::WATER:
+            assert(pu.phase_used[BlackoilPhases::Aqua]);
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
+            break;
+        case InjectorType::GAS:
+            assert(pu.phase_used[BlackoilPhases::Vapour]);
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
+            break;
+        case InjectorType::OIL:
+            assert(pu.phase_used[BlackoilPhases::Liquid]);
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
+            break;
+        case InjectorType::MULTI:
+            // Not currently handled, keep zero init.
+            break;
+        }
+    }
+
+    if (cmode_is_bhp)
+        this->bhp = bhp_limit;
+    else
+        this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor;
+}
 }
 
 
 
+
+

opm/simulators/wells/SingleWellState.hpp

@@ -23,11 +23,13 @@
 #include <functional>
 #include <vector>
 
-#include <opm/simulators/wells/SegmentState.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
+
+#include <opm/simulators/wells/SegmentState.hpp>
 #include <opm/simulators/wells/PerfData.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
 
 namespace Opm {
 
@@ -35,13 +37,20 @@ struct PerforationData;
 
 class SingleWellState {
 public:
-    SingleWellState(const std::string& name, const ParallelWellInfo& pinfo, bool is_producer, const std::vector<PerforationData>& perf_input, std::size_t num_phases, double temp);
+    SingleWellState(const std::string& name,
+                    const ParallelWellInfo& pinfo,
+                    bool is_producer,
+                    double presssure_first_connection,
+                    const std::vector<PerforationData>& perf_input,
+                    const PhaseUsage& pu,
+                    double temp);
 
     std::string name;
     std::reference_wrapper<const ParallelWellInfo> parallel_info;
 
     Well::Status status{Well::Status::OPEN};
     bool producer;
+    PhaseUsage pu;
     double bhp{0};
     double thp{0};
     double temperature{0};
@@ -65,6 +74,8 @@ public:
     ///    PerforationData::connection_transmissibility_factor actually
     ///    used (overwrites existing internal values).
     void reset_connection_factors(const std::vector<PerforationData>& new_perf_data);
+    void update_producer_targets(const Well& ecl_well, const SummaryState& st);
+    void update_injector_targets(const Well& ecl_well, const SummaryState& st);
     void updateStatus(Well::Status status);
     void init_timestep(const SingleWellState& other);
     void shut();

opm/simulators/wells/WellState.cpp

@@ -54,17 +54,15 @@ void WellState::base_init(const std::vector<double>& cellPressures,
     }
 }
 
-void WellState::initSingleProducer(const std::vector<double>& cellPressures,
-                                   const Well& well,
-                                   const std::vector<PerforationData>& well_perf_data,
+void WellState::initSingleProducer(const Well& well,
                                    const ParallelWellInfo& well_info,
+                                   double pressure_first_connection,
+                                   const std::vector<PerforationData>& well_perf_data,
                                    const SummaryState& summary_state) {
-    const double bhp_safety_factor = 0.99;
     const auto& pu = this->phase_usage_;
-    const int np = pu.num_phases;
     const double temp = 273.15 + 15.56;
 
-    auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, true, well_perf_data, static_cast<std::size_t>(np), temp});
+    auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, true, pressure_first_connection, well_perf_data, pu, temp});
     if ( ws.perf_data.empty())
         return;
 
@@ -72,67 +70,21 @@ void WellState::initSingleProducer(const std::vector<double>& cellPressures,
         ws.status = Well::Status::OPEN;
     }
 
-    const auto& prod_controls = well.isProducer() ? well.productionControls(summary_state) : Well::ProductionControls(0);
-    if (prod_controls.hasControl(Well::ProducerCMode::THP))
-        ws.thp = prod_controls.thp_limit;
-
-
-    auto cmode_is_bhp = (prod_controls.cmode == Well::ProducerCMode::BHP);
-    auto bhp_limit = prod_controls.bhp_limit;
-    auto pressure_first_connection = well_info.broadcastFirstPerforationValue(cellPressures[ws.perf_data.cell_index[0]]);
-
-    if (well.getStatus() == Well::Status::STOP) {
-        if (cmode_is_bhp)
-            ws.bhp = bhp_limit;
-        else
-            ws.bhp = pressure_first_connection;
-
-        return;
-    }
-
-    if (prod_controls.cmode == Well::ProducerCMode::GRUP) {
-        ws.bhp = pressure_first_connection * bhp_safety_factor;
-        return;
-    }
-
-    switch (prod_controls.cmode) {
-    case Well::ProducerCMode::ORAT:
-        assert(pu.phase_used[BlackoilPhases::Liquid]);
-        ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
-        break;
-    case Well::ProducerCMode::WRAT:
-        assert(pu.phase_used[BlackoilPhases::Aqua]);
-        ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
-        break;
-    case Well::ProducerCMode::GRAT:
-        assert(pu.phase_used[BlackoilPhases::Vapour]);
-        ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
-        break;
-    default:
-        // Keep zero init.
-        break;
-    }
-
-    if (cmode_is_bhp)
-        ws.bhp = bhp_limit;
-    else
-        ws.bhp = pressure_first_connection * bhp_safety_factor;
+    ws.update_producer_targets(well, summary_state);
 }
 
 
-void WellState::initSingleInjector(const std::vector<double>& cellPressures,
-                                   const Well& well,
-                                   const std::vector<PerforationData>& well_perf_data,
+void WellState::initSingleInjector(const Well& well,
                                    const ParallelWellInfo& well_info,
+                                   double pressure_first_connection,
+                                   const std::vector<PerforationData>& well_perf_data,
                                    const SummaryState& summary_state) {
 
-    const double bhp_safety_factor = 1.01;
     const auto& pu = this->phase_usage_;
-    const int np = pu.num_phases;
     const auto& inj_controls = well.injectionControls(summary_state);
     const double temp = inj_controls.temperature;
 
-    auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, well.isProducer(), well_perf_data, static_cast<std::size_t>(np), temp});
+    auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, false, pressure_first_connection, well_perf_data, pu, temp});
     if ( ws.perf_data.empty())
         return;
 
@@ -140,53 +92,7 @@ void WellState::initSingleInjector(const std::vector<double>& cellPressures,
         ws.status = Well::Status::OPEN;
     }
 
-    if (inj_controls.hasControl(Well::InjectorCMode::THP))
-        ws.thp = inj_controls.thp_limit;
-
-    auto cmode_is_bhp = (inj_controls.cmode == Well::InjectorCMode::BHP);
-    auto bhp_limit = inj_controls.bhp_limit;
-    auto pressure_first_connection = well_info.broadcastFirstPerforationValue(cellPressures[ws.perf_data.cell_index[0]]);
-
-    if (well.getStatus() == Well::Status::STOP) {
-        if (cmode_is_bhp)
-            ws.bhp = bhp_limit;
-        else
-            ws.bhp = pressure_first_connection;
-
-        return;
-    }
-
-
-    if (inj_controls.cmode == Well::InjectorCMode::GRUP) {
-        ws.bhp = pressure_first_connection * bhp_safety_factor;
-        return;
-    }
-
-    if (inj_controls.cmode == Well::InjectorCMode::RATE) {
-        auto inj_surf_rate = inj_controls.surface_rate;
-        switch (inj_controls.injector_type) {
-        case InjectorType::WATER:
-            assert(pu.phase_used[BlackoilPhases::Aqua]);
-            ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
-            break;
-        case InjectorType::GAS:
-            assert(pu.phase_used[BlackoilPhases::Vapour]);
-            ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
-            break;
-        case InjectorType::OIL:
-            assert(pu.phase_used[BlackoilPhases::Liquid]);
-            ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
-            break;
-        case InjectorType::MULTI:
-            // Not currently handled, keep zero init.
-            break;
-        }
-    }
-
-    if (cmode_is_bhp)
-        ws.bhp = bhp_limit;
-    else
-        ws.bhp = pressure_first_connection * bhp_safety_factor;
+    ws.update_injector_targets(well, summary_state);
 }
 
 
@@ -200,10 +106,15 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     if (well.isInjector() == well.isProducer())
         throw std::logic_error(fmt::format("Well must be either producer or injector - logic error for well: {}", well.name()));
 
+    double pressure_first_connection = -1;
+    if (!well_perf_data.empty())
+        pressure_first_connection = cellPressures[well_perf_data[0].cell_index];
+    pressure_first_connection = well_info.broadcastFirstPerforationValue(pressure_first_connection);
+
     if (well.isInjector())
-        this->initSingleInjector(cellPressures, well, well_perf_data, well_info, summary_state);
+        this->initSingleInjector(well, well_info, pressure_first_connection, well_perf_data, summary_state);
     else
-        this->initSingleProducer(cellPressures, well, well_perf_data, well_info, summary_state);
+        this->initSingleProducer(well, well_info, pressure_first_connection, well_perf_data, summary_state);
 }
 
 

opm/simulators/wells/WellState.hpp

@@ -296,16 +296,16 @@ private:
                         const ParallelWellInfo& well_info,
                         const SummaryState& summary_state);
 
-    void initSingleProducer(const std::vector<double>& cellPressures,
-                            const Well& well,
-                            const std::vector<PerforationData>& well_perf_data,
+    void initSingleProducer(const Well& well,
                             const ParallelWellInfo& well_info,
+                            double pressure_first_connection,
+                            const std::vector<PerforationData>& well_perf_data,
                             const SummaryState& summary_state);
 
-    void initSingleInjector(const std::vector<double>& cellPressures,
-                            const Well& well,
-                            const std::vector<PerforationData>& well_perf_data,
+    void initSingleInjector(const Well& well,
                             const ParallelWellInfo& well_info,
+                            double pressure_first_connection,
+                            const std::vector<PerforationData>& well_perf_data,
                             const SummaryState& summary_state);
 
 };

tests/test_wellstate.cpp

@@ -546,10 +546,10 @@ BOOST_AUTO_TEST_CASE(TESTSegmentState2) {
 
 
 BOOST_AUTO_TEST_CASE(TESTPerfData) {
-    Opm::PerfData pd1(3, true, 3);
-    Opm::PerfData pd2(3, true, 3);
-    Opm::PerfData pd3(2, true, 3);
-    Opm::PerfData pd4(3, false, 3);
+    Opm::PerfData pd1(3, 100, true, 3);
+    Opm::PerfData pd2(3, 100, true, 3);
+    Opm::PerfData pd3(2, 100, true, 3);
+    Opm::PerfData pd4(3, 100, false, 3);
 
 
     for (std::size_t i = 0; i < 3; i++) {
@@ -575,9 +575,14 @@ BOOST_AUTO_TEST_CASE(TESTPerfData) {
 BOOST_AUTO_TEST_CASE(TestSingleWellState) {
     Opm::ParallelWellInfo pinfo;
     std::vector<Opm::PerforationData> connections = {{0,1,1,0},{1,1,1,1},{2,1,1,2}};
-    Opm::SingleWellState ws1("W1", pinfo, true,  connections, 3, 1);
-    Opm::SingleWellState ws2("W2", pinfo, true,  connections, 3, 2);
-    Opm::SingleWellState ws3("W3", pinfo, false, connections, 3, 3);
+    Opm::PhaseUsage pu;
+
+    // This is totally bonkers, but the pu needs a complete deck to initialize properly
+    pu.num_phases = 3;
+
+    Opm::SingleWellState ws1("W1", pinfo, true,  100, connections, pu, 1);
+    Opm::SingleWellState ws2("W2", pinfo, true,  100, connections, pu, 2);
+    Opm::SingleWellState ws3("W3", pinfo, false, 100, connections, pu, 3);
 
     ws1.bhp = 100;
     ws1.thp = 200;