From a54daf75e84acbb6fac67aceb1e1abefb1eec0a7 Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Tue, 3 Aug 2021 20:05:14 +0200
Subject: [PATCH 1/7] Add trivial class SingleWellState with bhp and thp

---
 CMakeLists_files.cmake                        |  2 +
 .../wells/BlackoilWellModelGeneric.cpp        |  5 +-
 opm/simulators/wells/MultisegmentWellEval.cpp | 12 +++--
 .../wells/MultisegmentWellGeneric.cpp         |  3 +-
 .../wells/MultisegmentWell_impl.hpp           |  3 +-
 opm/simulators/wells/SingleWellState.cpp      | 43 +++++++++++++++++
 opm/simulators/wells/SingleWellState.hpp      | 42 +++++++++++++++++
 opm/simulators/wells/StandardWellEval.cpp     | 21 +++++----
 opm/simulators/wells/StandardWell_impl.hpp    | 11 +++--
 opm/simulators/wells/WellContainer.hpp        |  6 ++-
 .../wells/WellInterfaceFluidSystem.cpp        | 10 ++--
 opm/simulators/wells/WellInterface_impl.hpp   | 11 +++--
 opm/simulators/wells/WellState.cpp            | 47 ++++++++++---------
 opm/simulators/wells/WellState.hpp            | 29 ++++++++----
 tests/test_wellstate.cpp                      | 16 +++++++
 15 files changed, 194 insertions(+), 67 deletions(-)
 create mode 100644 opm/simulators/wells/SingleWellState.cpp
 create mode 100644 opm/simulators/wells/SingleWellState.hpp

diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index d72a0b420..9730a915f 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -70,6 +70,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/ParallelWellInfo.cpp
   opm/simulators/wells/PerfData.cpp
   opm/simulators/wells/SegmentState.cpp
+  opm/simulators/wells/SingleWellState.cpp
   opm/simulators/wells/StandardWellEval.cpp
   opm/simulators/wells/StandardWellGeneric.cpp
   opm/simulators/wells/TargetCalculator.cpp
@@ -307,6 +308,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/PerforationData.hpp
   opm/simulators/wells/RateConverter.hpp
   opm/simulators/utils/readDeck.hpp
+  opm/simulators/wells/SingleWellState.hpp
   opm/simulators/wells/TargetCalculator.hpp
   opm/simulators/wells/WellConnectionAuxiliaryModule.hpp
   opm/simulators/wells/WellState.hpp
diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
index 32cc5666b..232764e20 100644
--- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp
+++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
@@ -172,8 +172,9 @@ loadRestartData(const data::Wells& rst_wells,
     for( const auto& wm : well_state.wellMap() ) {
         const auto well_index = wm.second[ 0 ];
         const auto& rst_well = rst_wells.at( wm.first );
-        well_state.update_thp(well_index, rst_well.thp);
-        well_state.update_bhp(well_index, rst_well.bhp);
+        auto& ws = well_state.well(well_index);
+        ws.bhp = rst_well.bhp;
+        ws.thp = rst_well.thp;
         well_state.update_temperature(well_index,  rst_well.temperature);
 
         if (rst_well.current_control.isProducer) {
diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp
index 2e4a91094..f158de2f0 100644
--- a/opm/simulators/wells/MultisegmentWellEval.cpp
+++ b/opm/simulators/wells/MultisegmentWellEval.cpp
@@ -1474,7 +1474,8 @@ updateThp(WellState& well_state,
 
     // When there is no vaild VFP table provided, we set the thp to be zero.
     if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
-        well_state.update_thp(baseif_.indexOfWell(), 0.);
+        auto& well = well_state.well(baseif_.indexOfWell());
+        well.thp = 0;
         return;
     }
 
@@ -1492,9 +1493,9 @@ updateThp(WellState& well_state,
         rates[ Gas ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Gas ] ];
     }
 
-    const double bhp = well_state.bhp(baseif_.indexOfWell());
-
-    well_state.update_thp(baseif_.indexOfWell(), this->calculateThpFromBhp(rates, bhp, rho, deferred_logger));
+    auto& well = well_state.well(baseif_.indexOfWell());
+    const double bhp = well.bhp;
+    well.thp = this->calculateThpFromBhp(rates, bhp, rho, deferred_logger);
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -1614,6 +1615,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
     assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
     const int oil_pos = pu.phase_pos[Oil];
 
+    auto& well = well_state.well(baseif_.indexOfWell());
     auto& segments = well_state.segments(baseif_.indexOfWell());
     auto& segment_rates = segments.rates;
     auto& segment_pressure = segments.pressure;
@@ -1658,7 +1660,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
         // update the segment pressure
         segment_pressure[seg] = primary_variables_[seg][SPres];
         if (seg == 0) { // top segment
-            well_state.update_bhp(baseif_.indexOfWell(), segment_pressure[seg]);
+            well.bhp = segment_pressure[seg];
         }
     }
     updateThp(well_state, rho, deferred_logger);
diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp
index 246806ae3..6ccb36638 100644
--- a/opm/simulators/wells/MultisegmentWellGeneric.cpp
+++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp
@@ -136,8 +136,9 @@ void
 MultisegmentWellGeneric<Scalar>::
 scaleSegmentPressuresWithBhp(WellState& well_state) const
 {
+    auto& well = well_state.well(baseif_.indexOfWell());
     auto& segments = well_state.segments(baseif_.indexOfWell());
-    auto bhp = well_state.bhp(baseif_.indexOfWell());
+    const auto bhp = well.bhp;
     segments.scale_pressure(bhp);
 }
 
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index de29fc53e..1c2dcb2ef 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -340,6 +340,7 @@ namespace Opm
         // store a copy of the well state, we don't want to update the real well state
         WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
         const auto& group_state = ebosSimulator.problem().wellModel().groupState();
+        auto& ws = well_state_copy.well(this->index_of_well_);
 
         // Get the current controls.
         const auto& summary_state = ebosSimulator.vanguard().summaryState();
@@ -358,7 +359,7 @@ namespace Opm
             prod_controls.bhp_limit = bhp;
             well_state_copy.currentProductionControl(index_of_well_, Well::ProducerCMode::BHP);
         }
-        well_state_copy.update_bhp(well_copy.index_of_well_, bhp);
+        ws.bhp = bhp;
         well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
 
         // initialized the well rates with the potentials i.e. the well rates based on bhp
diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp
new file mode 100644
index 000000000..3f8f8a6ca
--- /dev/null
+++ b/opm/simulators/wells/SingleWellState.cpp
@@ -0,0 +1,43 @@
+/*
+  Copyright 2021 Equinor ASA
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/simulators/wells/SingleWellState.hpp>
+
+namespace Opm {
+
+
+void SingleWellState::init_timestep(const SingleWellState& other) {
+    this->bhp = other.bhp;
+    this->thp = other.thp;
+}
+
+
+void SingleWellState::shut() {
+    this->bhp = 0;
+    this->thp = 0;
+}
+
+void SingleWellState::stop() {
+    this->thp = 0;
+}
+
+}
+
+
+
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
new file mode 100644
index 000000000..d07caae94
--- /dev/null
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -0,0 +1,42 @@
+/*
+  Copyright 2021 Equinor ASA
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SINGLE_WELL_STATE_HEADER_INCLUDED
+#define OPM_SINGLE_WELL_STATE_HEADER_INCLUDED
+
+
+namespace Opm {
+
+class SingleWellState {
+public:
+    double bhp{0};
+    double thp{0};
+
+
+    void init_timestep(const SingleWellState& other);
+    void shut();
+    void stop();
+};
+
+
+}
+
+
+
+#endif
diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp
index cdc20a7be..41607275b 100644
--- a/opm/simulators/wells/StandardWellEval.cpp
+++ b/opm/simulators/wells/StandardWellEval.cpp
@@ -255,7 +255,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
     const int well_index = baseif_.indexOfWell();
     const int np = baseif_.numPhases();
     const auto& pu = baseif_.phaseUsage();
-
+    const auto& well = well_state.well(well_index);
     // the weighted total well rate
     double total_well_rate = 0.0;
     for (int p = 0; p < np; ++p) {
@@ -339,7 +339,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
 
 
     // BHP
-    primary_variables_[Bhp] = well_state.bhp(baseif_.indexOfWell());
+    primary_variables_[Bhp] = well.bhp;
 }
 
 template<class FluidSystem, class Indices, class Scalar>
@@ -553,10 +553,11 @@ updateThp(WellState& well_state,
     static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
     static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
     static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
+    auto& well = well_state.well(baseif_.indexOfWell());
 
     // When there is no vaild VFP table provided, we set the thp to be zero.
     if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
-        well_state.update_thp(baseif_.indexOfWell(), 0);
+        well.thp = 0;
         return;
     }
 
@@ -574,13 +575,12 @@ updateThp(WellState& well_state,
         rates[ Gas ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Gas ] ];
     }
 
-    const double bhp = well_state.bhp(baseif_.indexOfWell());
+    const double bhp = well.bhp;
 
-    well_state.update_thp(baseif_.indexOfWell(),
-                          this->calculateThpFromBhp(well_state,
-                                                    rates,
-                                                    bhp,
-                                                    deferred_logger));
+    well.thp = this->calculateThpFromBhp(well_state,
+                                         rates,
+                                         bhp,
+                                         deferred_logger);
 }
 
 template<class FluidSystem, class Indices, class Scalar>
@@ -652,7 +652,8 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
         F[pu.phase_pos[Gas]] += F_solvent;
     }
 
-    well_state.update_bhp(baseif_.indexOfWell(), primary_variables_[Bhp]);
+    auto& well = well_state.well(baseif_.indexOfWell());
+    well.bhp = primary_variables_[Bhp];
 
     // calculate the phase rates based on the primary variables
     // for producers, this is not a problem, while not sure for injectors here
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index 8ff47c1e3..eb5da692c 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -547,6 +547,7 @@ namespace Opm
                             cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
 
         auto& perf_data = well_state.perfData(this->index_of_well_);
+        auto& ws = well_state.well(this->index_of_well_);
         if constexpr (has_polymer && Base::has_polymermw) {
             if (this->isInjector()) {
                 // Store the original water flux computed from the reservoir quantities.
@@ -694,7 +695,7 @@ namespace Opm
         }
 
         // Store the perforation pressure for later usage.
-        perf_data.pressure[perf] = well_state.bhp(this->index_of_well_) + this->perf_pressure_diffs_[perf];
+        perf_data.pressure[perf] = ws.bhp + this->perf_pressure_diffs_[perf];
     }
 
 
@@ -1156,7 +1157,7 @@ namespace Opm
                                    const WellState& well_state,
                                    DeferredLogger& deferred_logger)
     {
-        const double bhp = well_state.bhp(index_of_well_);
+        const double bhp = well_state.well(index_of_well_).bhp;
         std::vector<double> well_rates;
         computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
 
@@ -1211,6 +1212,7 @@ namespace Opm
         b_perf.resize(nperf * num_components_);
         surf_dens_perf.resize(nperf * num_components_);
         const int w = index_of_well_;
+        const auto& ws = well_state.well(this->index_of_well_);
 
         const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
         const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
@@ -1224,7 +1226,7 @@ namespace Opm
 
         // Compute the average pressure in each well block
         const auto& perf_press = well_state.perfData(w).pressure;
-        auto p_above =  this->parallel_well_info_.communicateAboveValues(well_state.bhp(w),
+        auto p_above =  this->parallel_well_info_.communicateAboveValues(ws.bhp,
                                                                          perf_press.data(),
                                                                          nperf);
 
@@ -1681,6 +1683,7 @@ namespace Opm
         // to replace the original one
         WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
         const auto& group_state  = ebosSimulator.problem().wellModel().groupState();
+        auto& ws = well_state_copy.well(this->index_of_well_);
 
         //  Set current control to bhp, and bhp value in state, modify bhp limit in control object.
         if (well_ecl_.isInjector()) {
@@ -1688,7 +1691,7 @@ namespace Opm
         } else {
             well_state_copy.currentProductionControl(index_of_well_, Well::ProducerCMode::BHP);
         }
-        well_state_copy.update_bhp(index_of_well_, bhp);
+        ws.bhp = bhp;
 
         const double dt = ebosSimulator.timeStepSize();
         bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state_copy, group_state, deferred_logger);
diff --git a/opm/simulators/wells/WellContainer.hpp b/opm/simulators/wells/WellContainer.hpp
index e790b54f4..2ac62f50b 100644
--- a/opm/simulators/wells/WellContainer.hpp
+++ b/opm/simulators/wells/WellContainer.hpp
@@ -61,20 +61,22 @@ public:
         return this->m_data.size();
     }
 
-    void add(const std::string& name, T&& value) {
+    T& add(const std::string& name, T&& value) {
         if (index_map.count(name) != 0)
             throw std::logic_error("An object with name: " + name + " already exists in container");
 
         this->index_map.emplace(name, this->m_data.size());
         this->m_data.push_back(std::forward<T>(value));
+        return this->m_data.back();
     }
 
-    void add(const std::string& name, const T& value) {
+    T& add(const std::string& name, const T& value) {
         if (index_map.count(name) != 0)
             throw std::logic_error("An object with name: " + name + " already exists in container");
 
         this->index_map.emplace(name, this->m_data.size());
         this->m_data.push_back(value);
+        return this->m_data.back();
     }
 
     bool has(const std::string& name) const {
diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
index a7b22b90e..e9f440288 100644
--- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp
+++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
@@ -91,10 +91,11 @@ activeProductionConstraint(const WellState& well_state,
     const int well_index = this->index_of_well_;
     const auto controls = this->well_ecl_.productionControls(summaryState);
     auto currentControl = well_state.currentProductionControl(well_index);
+    auto& ws = well_state.well(well_index);
 
     if (controls.hasControl(Well::ProducerCMode::BHP) && currentControl != Well::ProducerCMode::BHP) {
         const double bhp_limit = controls.bhp_limit;
-        double current_bhp = well_state.bhp(well_index);
+        double current_bhp = ws.bhp;
         if (bhp_limit > current_bhp)
             return Well::ProducerCMode::BHP;
     }
@@ -164,7 +165,7 @@ activeProductionConstraint(const WellState& well_state,
 
     if (controls.hasControl(Well::ProducerCMode::THP) && currentControl != Well::ProducerCMode::THP) {
         const auto& thp = getTHPConstraint(summaryState);
-        double current_thp = well_state.thp(well_index);
+        double current_thp = ws.thp;
         if (thp > current_thp)
             return Well::ProducerCMode::THP;
     }
@@ -181,6 +182,7 @@ activeInjectionConstraint(const WellState& well_state,
 {
     const PhaseUsage& pu = this->phaseUsage();
     const int well_index = this->index_of_well_;
+    const auto& ws = well_state.well(well_index);
 
     const auto controls = this->well_ecl_.injectionControls(summaryState);
     auto currentControl = well_state.currentInjectionControl(well_index);
@@ -188,7 +190,7 @@ activeInjectionConstraint(const WellState& well_state,
     if (controls.hasControl(Well::InjectorCMode::BHP) && currentControl != Well::InjectorCMode::BHP)
     {
         const auto& bhp = controls.bhp_limit;
-        double current_bhp = well_state.bhp(well_index);
+        double current_bhp = ws.bhp;
         if (bhp < current_bhp)
             return Well::InjectorCMode::BHP;
     }
@@ -241,7 +243,7 @@ activeInjectionConstraint(const WellState& well_state,
     if (controls.hasControl(Well::InjectorCMode::THP) && currentControl != Well::InjectorCMode::THP)
     {
         const auto& thp = getTHPConstraint(summaryState);
-        double current_thp = well_state.thp(well_index);
+        double current_thp = ws.thp;
         if (thp < current_thp)
             return Well::InjectorCMode::THP;
     }
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index 2fe63bc97..63e39c47d 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -599,6 +599,7 @@ namespace Opm
         // only bhp and wellRates are used to initilize the primaryvariables for standard wells
         const auto& well = this->well_ecl_;
         const int well_index = this->index_of_well_;
+        auto& ws = well_state.well(well_index);
         const auto& pu = this->phaseUsage();
         const int np = well_state.numPhases();
         const auto& summaryState = ebos_simulator.vanguard().summaryState();
@@ -608,7 +609,7 @@ namespace Opm
             for (int p = 0; p<np; ++p) {
                 well_state.wellRates(well_index)[p] = 0.0;
             }
-            well_state.update_thp(well_index, 0.0);
+            ws.thp = 0;
             return;
         }
 
@@ -663,7 +664,7 @@ namespace Opm
                     rates[p] = well_state.wellRates(well_index)[p];
                 }
                 double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
-                well_state.update_bhp(well_index, bhp);
+                ws.bhp = bhp;
 
                 // if the total rates are negative or zero
                 // we try to provide a better intial well rate
@@ -678,7 +679,7 @@ namespace Opm
             }
             case Well::InjectorCMode::BHP:
             {
-                well_state.update_bhp(well_index, controls.bhp_limit);
+                ws.bhp = controls.bhp_limit;
                 double total_rate = 0.0;
                 for (int p = 0; p<np; ++p) {
                     total_rate += well_state.wellRates(well_index)[p];
@@ -864,7 +865,7 @@ namespace Opm
             }
             case Well::ProducerCMode::BHP:
             {
-                well_state.update_bhp(well_index, controls.bhp_limit);
+                ws.bhp = controls.bhp_limit;
                 double total_rate = 0.0;
                 for (int p = 0; p<np; ++p) {
                     total_rate -= well_state.wellRates(well_index)[p];
@@ -886,7 +887,7 @@ namespace Opm
                     rates[p] = well_state.wellRates(well_index)[p];
                 }
                 double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
-                well_state.update_bhp(well_index, bhp);
+                ws.bhp = bhp;
 
                 // if the total rates are negative or zero
                 // we try to provide a better intial well rate
diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp
index 97861d957..37693beaf 100644
--- a/opm/simulators/wells/WellState.cpp
+++ b/opm/simulators/wells/WellState.cpp
@@ -44,12 +44,11 @@ void WellState::base_init(const std::vector<double>& cellPressures,
     this->status_.clear();
     this->parallel_well_info_.clear();
     this->wellrates_.clear();
-    this->bhp_.clear();
-    this->thp_.clear();
     this->temperature_.clear();
     this->segment_state.clear();
     this->well_potentials_.clear();
     this->productivity_index_.clear();
+    this->wells_.clear();
     {
         // const int nw = wells->number_of_wells;
         const int nw = wells_ecl.size();
@@ -92,6 +91,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     const auto& pu = this->phase_usage_;
     const int np = pu.num_phases;
 
+    auto& ws = this->wells_.add(well.name(), SingleWellState{});
     this->status_.add(well.name(), Well::Status::OPEN);
     this->parallel_well_info_.add(well.name(), well_info);
     this->wellrates_.add(well.name(), std::vector<double>(np, 0));
@@ -99,8 +99,6 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     const int num_perf_this_well = well_info->communication().sum(well_perf_data.size());
     this->segment_state.add(well.name(), SegmentState{});
     this->perfdata.add(well.name(), PerfData{static_cast<std::size_t>(num_perf_this_well), well.isInjector(), this->phase_usage_});
-    this->bhp_.add(well.name(), 0.0);
-    this->thp_.add(well.name(), 0.0);
     this->productivity_index_.add(well.name(), std::vector<double>(np, 0));
     if ( well.isInjector() )
         this->temperature_.add(well.name(), well.injectionControls(summary_state).temperature);
@@ -136,9 +134,9 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
         //    applicable, otherwise assign equal to
         //    first perforation cell pressure.
         if (is_bhp) {
-            bhp_[w] = bhp_limit;
+            ws.bhp = bhp_limit;
         } else {
-            bhp_[w] = global_pressure;
+            ws.bhp = global_pressure;
         }
     } else if (is_grup) {
         // Well under group control.
@@ -148,7 +146,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
         //    the well is an injector or producer)
         //    pressure in first perforation cell.
         const double safety_factor = well.isInjector() ? 1.01 : 0.99;
-        bhp_[w] = safety_factor * global_pressure;
+        ws.bhp = safety_factor * global_pressure;
     } else {
         // Open well, under own control:
         // 1. Rates: initialize well rates to match
@@ -206,10 +204,10 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
         //    the well is an injector or producer)
         //    pressure in first perforation cell.
         if (is_bhp) {
-            bhp_[w] = bhp_limit;
+            ws.bhp = bhp_limit;
         } else {
             const double safety_factor = well.isInjector() ? 1.01 : 0.99;
-            bhp_[w] = safety_factor * global_pressure;
+            ws.bhp = safety_factor * global_pressure;
         }
     }
 
@@ -219,7 +217,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
                                            : prod_controls.hasControl(Well::ProducerCMode::THP);
     const double thp_limit = well.isInjector() ? inj_controls.thp_limit : prod_controls.thp_limit;
     if (has_thp) {
-        thp_[w] = thp_limit;
+        ws.thp = thp_limit;
     }
 }
 
@@ -356,11 +354,17 @@ void WellState::init(const std::vector<double>& cellPressures,
                 continue;
             }
             const auto& wname = well.name();
+            auto& new_well = this->well(w);
             auto it = prevState->wellMap().find(well.name());
             if (it != end)
             {
                 const int newIndex = w;
                 const int oldIndex = it->second[ 0 ];
+
+                const auto& prev_well = prevState->well(oldIndex);
+                new_well.init_timestep(prev_well);
+
+
                 if (prevState->status_[oldIndex] == Well::Status::SHUT) {
                     // Well was shut in previous state, do not use its values.
                     continue;
@@ -371,12 +375,6 @@ void WellState::init(const std::vector<double>& cellPressures,
                     continue;
                 }
 
-                // bhp
-                this->update_bhp( newIndex, prevState->bhp( oldIndex ));
-
-                // thp
-                this->update_thp( newIndex, prevState->thp( oldIndex ));
-
                 // If new target is set using WCONPROD, WCONINJE etc. we use the new control
                 if (!this->events_[w].hasEvent(WellState::event_mask)) {
                     current_injection_controls_[ newIndex ] = prevState->currentInjectionControl(oldIndex);
@@ -436,7 +434,7 @@ void WellState::init(const std::vector<double>& cellPressures,
                 : well.productionControls(summary_state).hasControl(Well::ProducerCMode::THP);
 
             if (!has_thp) {
-                this->update_thp(w, 0.0);
+                new_well.thp = 0;
             }
         }
     }
@@ -512,6 +510,7 @@ WellState::report(const int* globalCellIdxMap,
         {
             continue;
         }
+        const auto& ws = this->well(well_index);
 
         const auto& reservoir_rates = this->well_reservoir_rates_[well_index];
         const auto& well_potentials = this->well_potentials_[well_index];
@@ -519,8 +518,8 @@ WellState::report(const int* globalCellIdxMap,
         const auto& wv = this->wellRates(well_index);
 
         data::Well well;
-        well.bhp = this->bhp(well_index);
-        well.thp = this->thp( well_index );
+        well.bhp = ws.bhp;
+        well.thp = ws.thp;
         well.temperature = this->temperature( well_index );
 
         if( pu.phase_used[BlackoilPhases::Aqua] ) {
@@ -680,6 +679,7 @@ void WellState::initWellStateMSWell(const std::vector<Well>& wells_ecl,
     // what we do here, is to set the segment rates and perforation rates
     for (int w = 0; w < nw; ++w) {
         const auto& well_ecl = wells_ecl[w];
+        auto& ws = this->well(w);
 
         if ( well_ecl.isMultiSegment() ) {
             const WellSegments& segment_set = well_ecl.getSegments();
@@ -748,7 +748,7 @@ void WellState::initWellStateMSWell(const std::vector<Well>& wells_ecl,
             {
                 // top segment is always the first one, and its pressure is the well bhp
                 auto& segment_pressure = this->segments(w).pressure;
-                segment_pressure[0] = bhp(w);
+                segment_pressure[0] = ws.bhp;
                 const auto& perf_press = perf_data.pressure;
                 for (int seg = 1; seg < well_nseg; ++seg) {
                     if ( !segment_perforations[seg].empty() ) {
@@ -842,15 +842,16 @@ double WellState::brineWellRate(const int w) const
 
 void WellState::stopWell(int well_index)
 {
+    auto& ws = this->well(well_index);
+    ws.stop();
     this->status_[well_index] = Well::Status::STOP;
-    this->thp_[well_index] = 0;
 }
 
 void WellState::shutWell(int well_index)
 {
+    auto& ws = this->well(well_index);
+    ws.shut();
     this->status_[well_index] = Well::Status::SHUT;
-    this->thp_[well_index] = 0;
-    this->bhp_[well_index] = 0;
     const int np = numPhases();
     this->wellrates_[well_index].assign(np, 0);
 
diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp
index d8ac6b40b..997d99364 100644
--- a/opm/simulators/wells/WellState.hpp
+++ b/opm/simulators/wells/WellState.hpp
@@ -22,6 +22,7 @@
 #define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
 
 #include <opm/simulators/wells/ALQState.hpp>
+#include <opm/simulators/wells/SingleWellState.hpp>
 #include <opm/simulators/wells/GlobalWellInfo.hpp>
 #include <opm/simulators/wells/SegmentState.hpp>
 #include <opm/simulators/wells/WellContainer.hpp>
@@ -309,14 +310,6 @@ public:
         return this->phase_usage_;
     }
 
-    /// One bhp pressure per well.
-    void update_bhp(std::size_t well_index, double value) { bhp_[well_index] = value; }
-    double bhp(std::size_t well_index) const { return bhp_[well_index]; }
-
-    /// One thp pressure per well.
-    void update_thp(std::size_t well_index, double value) { thp_[well_index] = value; }
-    double thp(std::size_t well_index) const { return thp_[well_index]; }
-
     /// One temperature per well.
     void update_temperature(std::size_t well_index, double value) { temperature_[well_index] = value; }
     double temperature(std::size_t well_index) const { return temperature_[well_index]; }
@@ -352,6 +345,23 @@ public:
         return this->is_producer_[well_index];
     }
 
+    const SingleWellState& well(std::size_t well_index) const {
+        return this->wells_[well_index];
+    }
+
+    const SingleWellState& well(const std::string& well_name) const {
+        return this->wells_[well_name];
+    }
+
+    SingleWellState& well(std::size_t well_index) {
+        return this->wells_[well_index];
+    }
+
+    SingleWellState& well(const std::string& well_name) {
+        return this->wells_[well_name];
+    }
+
+
 
 private:
     WellMapType wellMap_;
@@ -362,10 +372,9 @@ private:
     ALQState alq_state;
     bool do_glift_optimization_;
 
+    WellContainer<SingleWellState> wells_;
     WellContainer<Well::Status> status_;
     WellContainer<const ParallelWellInfo*> parallel_well_info_;
-    WellContainer<double> bhp_;
-    WellContainer<double> thp_;
     WellContainer<double> temperature_;
     WellContainer<std::vector<double>> wellrates_;
     PhaseUsage phase_usage_;
diff --git a/tests/test_wellstate.cpp b/tests/test_wellstate.cpp
index 7975cd65c..12ab328b7 100644
--- a/tests/test_wellstate.cpp
+++ b/tests/test_wellstate.cpp
@@ -26,6 +26,7 @@
 #include <opm/simulators/wells/GlobalWellInfo.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
 #include <opm/simulators/wells/WellState.hpp>
+#include <opm/simulators/wells/SingleWellState.hpp>
 #include <opm/simulators/wells/SegmentState.hpp>
 #include <opm/simulators/wells/WellContainer.hpp>
 #include <opm/simulators/wells/PerfData.hpp>
@@ -575,5 +576,20 @@ GAS
 }
 
 
+BOOST_AUTO_TEST_CASE(TestSingleWellState) {
+    Opm::SingleWellState ws1;
+    Opm::SingleWellState ws2;
+
+    ws1.bhp = 100;
+    ws1.thp = 200;
+
+    ws2.init_timestep(ws1);
+    BOOST_CHECK_EQUAL(ws2.bhp, ws1.bhp);
+    BOOST_CHECK_EQUAL(ws2.thp, ws1.thp);
+}
+
+
+
+
 
 BOOST_AUTO_TEST_SUITE_END()

From 1dd9b91ad3c95c555c60bd3dd696b06ce31bc32f Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Tue, 3 Aug 2021 21:11:09 +0200
Subject: [PATCH 2/7] Add member "bool producer" to SingleWellState

---
 opm/simulators/wells/SingleWellState.cpp |  7 +++++++
 opm/simulators/wells/SingleWellState.hpp |  5 +++++
 opm/simulators/wells/WellState.cpp       | 15 ++++-----------
 opm/simulators/wells/WellState.hpp       |  6 ------
 tests/test_wellstate.cpp                 | 11 +++++++++--
 5 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp
index 3f8f8a6ca..251c58075 100644
--- a/opm/simulators/wells/SingleWellState.cpp
+++ b/opm/simulators/wells/SingleWellState.cpp
@@ -21,8 +21,15 @@
 
 namespace Opm {
 
+SingleWellState::SingleWellState(bool is_producer)
+    : producer(is_producer)
+{}
+
 
 void SingleWellState::init_timestep(const SingleWellState& other) {
+    if (this->producer != other.producer)
+        return;
+
     this->bhp = other.bhp;
     this->thp = other.thp;
 }
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
index d07caae94..6ccded3b7 100644
--- a/opm/simulators/wells/SingleWellState.hpp
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -25,6 +25,11 @@ namespace Opm {
 
 class SingleWellState {
 public:
+    explicit SingleWellState(bool is_producer);
+
+
+
+    bool producer;
     double bhp{0};
     double thp{0};
 
diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp
index 37693beaf..08d6e60f7 100644
--- a/opm/simulators/wells/WellState.cpp
+++ b/opm/simulators/wells/WellState.cpp
@@ -91,7 +91,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     const auto& pu = this->phase_usage_;
     const int np = pu.num_phases;
 
-    auto& ws = this->wells_.add(well.name(), SingleWellState{});
+    auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer()});
     this->status_.add(well.name(), Well::Status::OPEN);
     this->parallel_well_info_.add(well.name(), well_info);
     this->wellrates_.add(well.name(), std::vector<double>(np, 0));
@@ -305,12 +305,6 @@ void WellState::init(const std::vector<double>& cellPressures,
         this->well_vaporized_oil_rates_.add(wname, 0);
     }
 
-    is_producer_.clear();
-    for (int w = 0; w < nw; ++w) {
-        const auto& ecl_well = wells_ecl[w];
-        this->is_producer_.add( ecl_well.name(), ecl_well.isProducer());
-    }
-
     current_injection_controls_.clear();
     current_production_controls_.clear();
     for (int w = 0; w < nw; ++w) {
@@ -370,10 +364,9 @@ void WellState::init(const std::vector<double>& cellPressures,
                     continue;
                 }
 
-                if (is_producer_[newIndex] != prevState->is_producer_[oldIndex]) {
+                if (new_well.producer != prev_well.producer)
                     // Well changed to/from injector from/to producer, do not use its privious values.
                     continue;
-                }
 
                 // If new target is set using WCONPROD, WCONINJE etc. we use the new control
                 if (!this->events_[w].hasEvent(WellState::event_mask)) {
@@ -555,7 +548,7 @@ WellState::report(const int* globalCellIdxMap,
             well.rates.set(rt::brine, brineWellRate(well_index));
         }
 
-        if (is_producer_[well_index]) {
+        if (ws.producer) {
             well.rates.set(rt::alq, getALQ(wname));
         }
         else {
@@ -568,7 +561,7 @@ WellState::report(const int* globalCellIdxMap,
         {
             auto& curr = well.current_control;
 
-            curr.isProducer = this->is_producer_[well_index];
+            curr.isProducer = ws.producer;
             curr.prod = this->currentProductionControl(well_index);
             curr.inj  = this->currentInjectionControl(well_index);
         }
diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp
index 997d99364..e5d138bfc 100644
--- a/opm/simulators/wells/WellState.hpp
+++ b/opm/simulators/wells/WellState.hpp
@@ -341,10 +341,6 @@ public:
         return this->status_.well_name(well_index);
     }
 
-    bool producer(std::size_t well_index) const {
-        return this->is_producer_[well_index];
-    }
-
     const SingleWellState& well(std::size_t well_index) const {
         return this->wells_[well_index];
     }
@@ -380,8 +376,6 @@ private:
     PhaseUsage phase_usage_;
     WellContainer<PerfData> perfdata;
 
-    WellContainer<int> is_producer_; // Size equal to number of local wells.
-
     // vector with size number of wells +1.
     // iterate over all perforations of a given well
     // for (int perf = first_perf_index_[well_index]; perf < first_perf_index_[well_index] + num_perf_[well_index]; ++perf)
diff --git a/tests/test_wellstate.cpp b/tests/test_wellstate.cpp
index 12ab328b7..e5de1468a 100644
--- a/tests/test_wellstate.cpp
+++ b/tests/test_wellstate.cpp
@@ -577,8 +577,9 @@ GAS
 
 
 BOOST_AUTO_TEST_CASE(TestSingleWellState) {
-    Opm::SingleWellState ws1;
-    Opm::SingleWellState ws2;
+    Opm::SingleWellState ws1(true);
+    Opm::SingleWellState ws2(true);
+    Opm::SingleWellState ws3(false);
 
     ws1.bhp = 100;
     ws1.thp = 200;
@@ -586,6 +587,12 @@ BOOST_AUTO_TEST_CASE(TestSingleWellState) {
     ws2.init_timestep(ws1);
     BOOST_CHECK_EQUAL(ws2.bhp, ws1.bhp);
     BOOST_CHECK_EQUAL(ws2.thp, ws1.thp);
+
+    ws3.bhp = ws1.bhp * 2;
+    ws3.thp = ws1.thp * 2;
+    ws3.init_timestep(ws1);
+    BOOST_CHECK(ws3.bhp != ws1.bhp);
+    BOOST_CHECK(ws3.thp != ws1.thp);
 }
 
 

From d5ab308fbd95147f581ba53e36b3e274f014ab6b Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Tue, 3 Aug 2021 21:25:03 +0200
Subject: [PATCH 3/7] Add temperature to SingleWellState

---
 opm/simulators/wells/BlackoilWellModelGeneric.cpp |  2 +-
 opm/simulators/wells/BlackoilWellModel_impl.hpp   |  2 +-
 opm/simulators/wells/SingleWellState.cpp          |  4 +++-
 opm/simulators/wells/SingleWellState.hpp          |  5 ++---
 opm/simulators/wells/WellState.cpp                | 10 +++-------
 opm/simulators/wells/WellState.hpp                |  5 -----
 tests/test_wellstate.cpp                          |  6 +++---
 7 files changed, 13 insertions(+), 21 deletions(-)

diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
index 232764e20..fdf4d0d68 100644
--- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp
+++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
@@ -175,7 +175,7 @@ loadRestartData(const data::Wells& rst_wells,
         auto& ws = well_state.well(well_index);
         ws.bhp = rst_well.bhp;
         ws.thp = rst_well.thp;
-        well_state.update_temperature(well_index,  rst_well.temperature);
+        ws.temperature = rst_well.temperature;
 
         if (rst_well.current_control.isProducer) {
             well_state.currentProductionControl(well_index, rst_well.current_control.prod);
diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index 631956129..5f24060ef 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -1732,7 +1732,7 @@ namespace Opm {
             }
             weighted_temperature = well_info.communication().sum(weighted_temperature);
             total_weight = well_info.communication().sum(total_weight);
-            this->wellState().update_temperature(wellID, weighted_temperature/total_weight);
+            this->wellState().well(wellID).temperature = weighted_temperature/total_weight;
         }
     }
 
diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp
index 251c58075..2784bdebc 100644
--- a/opm/simulators/wells/SingleWellState.cpp
+++ b/opm/simulators/wells/SingleWellState.cpp
@@ -21,8 +21,9 @@
 
 namespace Opm {
 
-SingleWellState::SingleWellState(bool is_producer)
+SingleWellState::SingleWellState(bool is_producer, double temp)
     : producer(is_producer)
+    , temperature(temp)
 {}
 
 
@@ -32,6 +33,7 @@ void SingleWellState::init_timestep(const SingleWellState& other) {
 
     this->bhp = other.bhp;
     this->thp = other.thp;
+    this->temperature = other.temperature;
 }
 
 
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
index 6ccded3b7..0a1eb1ac6 100644
--- a/opm/simulators/wells/SingleWellState.hpp
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -25,14 +25,13 @@ namespace Opm {
 
 class SingleWellState {
 public:
-    explicit SingleWellState(bool is_producer);
-
+    SingleWellState(bool is_producer, double temp);
 
 
     bool producer;
     double bhp{0};
     double thp{0};
-
+    double temperature{};
 
     void init_timestep(const SingleWellState& other);
     void shut();
diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp
index 08d6e60f7..cca31b5c8 100644
--- a/opm/simulators/wells/WellState.cpp
+++ b/opm/simulators/wells/WellState.cpp
@@ -44,7 +44,6 @@ void WellState::base_init(const std::vector<double>& cellPressures,
     this->status_.clear();
     this->parallel_well_info_.clear();
     this->wellrates_.clear();
-    this->temperature_.clear();
     this->segment_state.clear();
     this->well_potentials_.clear();
     this->productivity_index_.clear();
@@ -90,8 +89,9 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     // May be overwritten below.
     const auto& pu = this->phase_usage_;
     const int np = pu.num_phases;
+    double temp = well.isInjector() ? well.injectionControls(summary_state).temperature : 273.15 + 15.56;
 
-    auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer()});
+    auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), temp});
     this->status_.add(well.name(), Well::Status::OPEN);
     this->parallel_well_info_.add(well.name(), well_info);
     this->wellrates_.add(well.name(), std::vector<double>(np, 0));
@@ -100,10 +100,6 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     this->segment_state.add(well.name(), SegmentState{});
     this->perfdata.add(well.name(), PerfData{static_cast<std::size_t>(num_perf_this_well), well.isInjector(), this->phase_usage_});
     this->productivity_index_.add(well.name(), std::vector<double>(np, 0));
-    if ( well.isInjector() )
-        this->temperature_.add(well.name(), well.injectionControls(summary_state).temperature);
-    else
-        this->temperature_.add(well.name(), 273.15 + 15.56); // standard condition temperature
 
     if ( num_perf_this_well == 0 )
         return;
@@ -513,7 +509,7 @@ WellState::report(const int* globalCellIdxMap,
         data::Well well;
         well.bhp = ws.bhp;
         well.thp = ws.thp;
-        well.temperature = this->temperature( well_index );
+        well.temperature = ws.temperature;
 
         if( pu.phase_used[BlackoilPhases::Aqua] ) {
             well.rates.set(rt::wat, wv[ pu.phase_pos[BlackoilPhases::Aqua] ] );
diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp
index e5d138bfc..884b3b324 100644
--- a/opm/simulators/wells/WellState.hpp
+++ b/opm/simulators/wells/WellState.hpp
@@ -310,10 +310,6 @@ public:
         return this->phase_usage_;
     }
 
-    /// One temperature per well.
-    void update_temperature(std::size_t well_index, double value) { temperature_[well_index] = value; }
-    double temperature(std::size_t well_index) const { return temperature_[well_index]; }
-
     /// One rate per well and phase.
     const WellContainer<std::vector<double>>& wellRates() const { return wellrates_; }
     std::vector<double>& wellRates(std::size_t well_index) { return wellrates_[well_index]; }
@@ -371,7 +367,6 @@ private:
     WellContainer<SingleWellState> wells_;
     WellContainer<Well::Status> status_;
     WellContainer<const ParallelWellInfo*> parallel_well_info_;
-    WellContainer<double> temperature_;
     WellContainer<std::vector<double>> wellrates_;
     PhaseUsage phase_usage_;
     WellContainer<PerfData> perfdata;
diff --git a/tests/test_wellstate.cpp b/tests/test_wellstate.cpp
index e5de1468a..c806d9b91 100644
--- a/tests/test_wellstate.cpp
+++ b/tests/test_wellstate.cpp
@@ -577,9 +577,9 @@ GAS
 
 
 BOOST_AUTO_TEST_CASE(TestSingleWellState) {
-    Opm::SingleWellState ws1(true);
-    Opm::SingleWellState ws2(true);
-    Opm::SingleWellState ws3(false);
+    Opm::SingleWellState ws1(true, 1);
+    Opm::SingleWellState ws2(true, 2);
+    Opm::SingleWellState ws3(false, 3);
 
     ws1.bhp = 100;
     ws1.thp = 200;

From 75eb65efb7e926794c8954af47f6b04f4e5ad7fd Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Tue, 3 Aug 2021 21:29:47 +0200
Subject: [PATCH 4/7] Use WellState::wells_ to infer well name

---
 opm/simulators/wells/WellState.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp
index 884b3b324..80563a61a 100644
--- a/opm/simulators/wells/WellState.hpp
+++ b/opm/simulators/wells/WellState.hpp
@@ -334,7 +334,7 @@ public:
     }
 
     const std::string& name(std::size_t well_index) const {
-        return this->status_.well_name(well_index);
+        return this->wells_.well_name(well_index);
     }
 
     const SingleWellState& well(std::size_t well_index) const {

From 379e938ee1253829262e96e95d46d448c37eae81 Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Wed, 4 Aug 2021 13:50:08 +0200
Subject: [PATCH 5/7] Add status to SingleWellState

---
 opm/simulators/wells/GlobalWellInfo.cpp  | 28 ++++++++++--------------
 opm/simulators/wells/GlobalWellInfo.hpp  |  5 +++--
 opm/simulators/wells/SingleWellState.cpp |  8 +++++++
 opm/simulators/wells/SingleWellState.hpp |  3 ++-
 opm/simulators/wells/WellState.cpp       | 25 +++++++++++----------
 opm/simulators/wells/WellState.hpp       |  3 +--
 6 files changed, 39 insertions(+), 33 deletions(-)

diff --git a/opm/simulators/wells/GlobalWellInfo.cpp b/opm/simulators/wells/GlobalWellInfo.cpp
index 89e091343..2a7de8228 100644
--- a/opm/simulators/wells/GlobalWellInfo.cpp
+++ b/opm/simulators/wells/GlobalWellInfo.cpp
@@ -37,10 +37,8 @@ GlobalWellInfo::GlobalWellInfo(const Schedule& sched, std::size_t report_step, c
         this->name_map.emplace( well.name(), global_well_index );
     }
 
-    for (const auto& well : local_wells) {
+    for (const auto& well : local_wells)
         this->local_map.push_back( well.seqIndex() );
-        this->is_injector.push_back( well.isInjector() );
-    }
 }
 
 
@@ -56,23 +54,19 @@ bool GlobalWellInfo::in_producing_group(const std::string& wname) const {
 }
 
 
-void GlobalWellInfo::update_group(const std::vector<Well::Status>& well_status, const std::vector<Well::InjectorCMode>& injection_cmode, const std::vector<Well::ProducerCMode>& production_cmode) {
-    if (well_status.size() != this->local_map.size())
-        throw std::logic_error("Size mismatch");
+void GlobalWellInfo::update_injector(std::size_t well_index, Well::Status well_status, Well::InjectorCMode injection_cmode) {
+    if (well_status == Well::Status::OPEN && injection_cmode == Well::InjectorCMode::GRUP)
+        this->m_in_injecting_group[this->local_map[well_index]] = 1;
+}
 
+void GlobalWellInfo::update_producer(std::size_t well_index, Well::Status well_status, Well::ProducerCMode production_cmode) {
+    if (well_status == Well::Status::OPEN && production_cmode == Well::ProducerCMode::GRUP)
+        this->m_in_producing_group[this->local_map[well_index]] = 1;
+}
+
+void GlobalWellInfo::clear() {
     this->m_in_injecting_group.assign(this->name_map.size(), 0);
     this->m_in_producing_group.assign(this->name_map.size(), 0);
-    for (std::size_t well_index = 0; well_index < well_status.size(); well_index++) {
-        if (well_status[well_index] == Well::Status::OPEN) {
-            if (this->is_injector[well_index]) {
-                if (injection_cmode[well_index] == Well::InjectorCMode::GRUP)
-                    this->m_in_injecting_group[this->local_map[well_index]] = 1;
-            } else {
-                if (production_cmode[well_index] == Well::ProducerCMode::GRUP)
-                    this->m_in_producing_group[this->local_map[well_index]] = 1;
-            }
-        }
-    }
 }
 
 
diff --git a/opm/simulators/wells/GlobalWellInfo.hpp b/opm/simulators/wells/GlobalWellInfo.hpp
index c03b9b175..ca90fd8d1 100644
--- a/opm/simulators/wells/GlobalWellInfo.hpp
+++ b/opm/simulators/wells/GlobalWellInfo.hpp
@@ -69,13 +69,14 @@ public:
     GlobalWellInfo(const Schedule& sched, std::size_t report_step, const std::vector<Well>& local_wells);
     bool in_producing_group(const std::string& wname) const;
     bool in_injecting_group(const std::string& wname) const;
-    void update_group(const std::vector<Well::Status>& well_status, const std::vector<Well::InjectorCMode>& injection_cmode, const std::vector<Well::ProducerCMode>& production_cmode);
     std::size_t well_index(const std::string& wname) const;
     const std::string& well_name(std::size_t well_index) const;
+    void update_injector(std::size_t well_index, Well::Status well_status, Well::InjectorCMode injection_cmode);
+    void update_producer(std::size_t well_index, Well::Status well_status, Well::ProducerCMode production_cmode);
+    void clear();
 
 private:
     std::vector<std::size_t> local_map;    // local_index -> global_index
-    std::vector<bool> is_injector;         // local_index -> bool
 
     std::map<std::string, std::size_t> name_map; // string -> global_index
     std::vector<int> m_in_injecting_group;       // global_index -> int/bool
diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp
index 2784bdebc..614b6c508 100644
--- a/opm/simulators/wells/SingleWellState.cpp
+++ b/opm/simulators/wells/SingleWellState.cpp
@@ -31,6 +31,12 @@ void SingleWellState::init_timestep(const SingleWellState& other) {
     if (this->producer != other.producer)
         return;
 
+    if (this->status == Well::Status::SHUT)
+        return;
+
+    if (other.status == Well::Status::SHUT)
+        return;
+
     this->bhp = other.bhp;
     this->thp = other.thp;
     this->temperature = other.temperature;
@@ -40,10 +46,12 @@ void SingleWellState::init_timestep(const SingleWellState& other) {
 void SingleWellState::shut() {
     this->bhp = 0;
     this->thp = 0;
+    this->status = Well::Status::SHUT;
 }
 
 void SingleWellState::stop() {
     this->thp = 0;
+    this->status = Well::Status::STOP;
 }
 
 }
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
index 0a1eb1ac6..f45175703 100644
--- a/opm/simulators/wells/SingleWellState.hpp
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -20,6 +20,7 @@
 #ifndef OPM_SINGLE_WELL_STATE_HEADER_INCLUDED
 #define OPM_SINGLE_WELL_STATE_HEADER_INCLUDED
 
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
 
 namespace Opm {
 
@@ -27,7 +28,7 @@ class SingleWellState {
 public:
     SingleWellState(bool is_producer, double temp);
 
-
+    Well::Status status{Well::Status::OPEN};
     bool producer;
     double bhp{0};
     double thp{0};
diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp
index cca31b5c8..29c91463e 100644
--- a/opm/simulators/wells/WellState.cpp
+++ b/opm/simulators/wells/WellState.cpp
@@ -41,7 +41,6 @@ void WellState::base_init(const std::vector<double>& cellPressures,
     // clear old name mapping
     this->wellMap_.clear();
     this->perfdata.clear();
-    this->status_.clear();
     this->parallel_well_info_.clear();
     this->wellrates_.clear();
     this->segment_state.clear();
@@ -92,7 +91,6 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     double temp = well.isInjector() ? well.injectionControls(summary_state).temperature : 273.15 + 15.56;
 
     auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), temp});
-    this->status_.add(well.name(), Well::Status::OPEN);
     this->parallel_well_info_.add(well.name(), well_info);
     this->wellrates_.add(well.name(), std::vector<double>(np, 0));
     this->well_potentials_.add(well.name(), std::vector<double>(np, 0));
@@ -120,7 +118,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
     const double global_pressure = well_info->broadcastFirstPerforationValue(local_pressure);
 
     if (well.getStatus() == Well::Status::OPEN) {
-        this->status_[w] = Well::Status::OPEN;
+        ws.status = Well::Status::OPEN;
     }
 
     if (well.getStatus() == Well::Status::STOP) {
@@ -355,7 +353,7 @@ void WellState::init(const std::vector<double>& cellPressures,
                 new_well.init_timestep(prev_well);
 
 
-                if (prevState->status_[oldIndex] == Well::Status::SHUT) {
+                if (prev_well.status == Well::Status::SHUT) {
                     // Well was shut in previous state, do not use its values.
                     continue;
                 }
@@ -494,12 +492,11 @@ WellState::report(const int* globalCellIdxMap,
     data::Wells res;
     for( const auto& [wname, winfo]: this->wellMap() ) {
         const auto well_index = winfo[ 0 ];
-        if ((this->status_[well_index] == Well::Status::SHUT) &&
-            ! wasDynamicallyClosed(well_index))
+        const auto& ws = this->well(well_index);
+        if ((ws.status == Well::Status::SHUT) && !wasDynamicallyClosed(well_index))
         {
             continue;
         }
-        const auto& ws = this->well(well_index);
 
         const auto& reservoir_rates = this->well_reservoir_rates_[well_index];
         const auto& well_potentials = this->well_potentials_[well_index];
@@ -766,7 +763,8 @@ void WellState::initWellStateMSWell(const std::vector<Well>& wells_ecl,
 
             const auto& wname = well.name();
             if (prev_well_state->segment_state.has(wname)) {
-                if (prev_well_state->status_[wname] == Well::Status::SHUT) {
+                const auto& prev_ws = prev_well_state->well(wname);
+                if (prev_ws.status == Well::Status::SHUT) {
                     continue;
                 }
 
@@ -833,14 +831,12 @@ void WellState::stopWell(int well_index)
 {
     auto& ws = this->well(well_index);
     ws.stop();
-    this->status_[well_index] = Well::Status::STOP;
 }
 
 void WellState::shutWell(int well_index)
 {
     auto& ws = this->well(well_index);
     ws.shut();
-    this->status_[well_index] = Well::Status::SHUT;
     const int np = numPhases();
     this->wellrates_[well_index].assign(np, 0);
 
@@ -925,7 +921,14 @@ void WellState::communicateGroupRates(const Comm& comm)
 template<class Comm>
 void WellState::updateGlobalIsGrup(const Comm& comm)
 {
-    this->global_well_info.value().update_group(this->status_.data(), this->current_injection_controls_.data(), this->current_production_controls_.data());
+    this->global_well_info.value().clear();
+    for (std::size_t well_index = 0; well_index < this->size(); well_index++) {
+        const auto& ws = this->well(well_index);
+        if (ws.producer)
+            this->global_well_info.value().update_producer(well_index, ws.status, this->current_production_controls_[well_index]);
+        else
+            this->global_well_info.value().update_injector(well_index, ws.status, this->current_injection_controls_[well_index]);
+    }
     this->global_well_info.value().communicate(comm);
 }
 
diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp
index 80563a61a..3a9bb8dc7 100644
--- a/opm/simulators/wells/WellState.hpp
+++ b/opm/simulators/wells/WellState.hpp
@@ -294,7 +294,7 @@ public:
     void updateStatus(int well_index, Well::Status status);
 
     void openWell(int well_index) {
-        this->status_[well_index] = Well::Status::OPEN;
+        this->wells_[well_index].status = Well::Status::OPEN;
     }
 
     void shutWell(int well_index);
@@ -365,7 +365,6 @@ private:
     bool do_glift_optimization_;
 
     WellContainer<SingleWellState> wells_;
-    WellContainer<Well::Status> status_;
     WellContainer<const ParallelWellInfo*> parallel_well_info_;
     WellContainer<std::vector<double>> wellrates_;
     PhaseUsage phase_usage_;

From c38145980441e3bd18c0e96f14b7f9f0c86112e1 Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Wed, 4 Aug 2021 10:08:31 +0200
Subject: [PATCH 6/7] Add Events member to SingleWellState

---
 opm/simulators/wells/BlackoilWellModel_impl.hpp | 4 ++--
 opm/simulators/wells/SingleWellState.hpp        | 2 ++
 opm/simulators/wells/WellState.cpp              | 7 ++-----
 opm/simulators/wells/WellState.hpp              | 9 ---------
 4 files changed, 6 insertions(+), 16 deletions(-)

diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index 5f24060ef..80d58456f 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -612,7 +612,7 @@ namespace Opm {
                 if (this->wellTestState_.hasWellClosed(well_name)) {
                     // TODO: more checking here, to make sure this standard more specific and complete
                     // maybe there is some WCON keywords will not open the well
-                    auto& events = this->wellState().events(w);
+                    auto& events = this->wellState().well(w).events;
                     if (events.hasEvent(WellState::event_mask)) {
                         if (wellTestState_.lastTestTime(well_name) == ebosSimulator_.time()) {
                             // The well was shut this timestep, we are most likely retrying
@@ -1424,7 +1424,7 @@ namespace Opm {
 
                 if (!well->isOperable() ) continue;
 
-                auto& events = this->wellState().events(well->indexOfWell());
+                auto& events = this->wellState().well(well->indexOfWell()).events;
                 if (events.hasEvent(WellState::event_mask)) {
                     well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
                     // There is no new well control change input within a report step,
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
index f45175703..76a023216 100644
--- a/opm/simulators/wells/SingleWellState.hpp
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -21,6 +21,7 @@
 #define OPM_SINGLE_WELL_STATE_HEADER_INCLUDED
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
 
 namespace Opm {
 
@@ -33,6 +34,7 @@ public:
     double bhp{0};
     double thp{0};
     double temperature{};
+    Events events;
 
     void init_timestep(const SingleWellState& other);
     void shut();
diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp
index 29c91463e..030e73834 100644
--- a/opm/simulators/wells/WellState.cpp
+++ b/opm/simulators/wells/WellState.cpp
@@ -257,15 +257,12 @@ void WellState::init(const std::vector<double>& cellPressures,
     well_dissolved_gas_rates_.clear();
     well_vaporized_oil_rates_.clear();
 
-    this->events_.clear();
     {
         const auto& wg_events = schedule[report_step].wellgroup_events();
         for (const auto& ecl_well : wells_ecl) {
             const auto& wname = ecl_well.name();
             if (wg_events.has(wname))
-                this->events_.add( wname, wg_events.at(wname) );
-            else
-                this->events_.add( wname, Events() );
+                this->well(wname).events = wg_events.at(wname);
         }
     }
 
@@ -363,7 +360,7 @@ void WellState::init(const std::vector<double>& cellPressures,
                     continue;
 
                 // If new target is set using WCONPROD, WCONINJE etc. we use the new control
-                if (!this->events_[w].hasEvent(WellState::event_mask)) {
+                if (!new_well.events.hasEvent(WellState::event_mask)) {
                     current_injection_controls_[ newIndex ] = prevState->currentInjectionControl(oldIndex);
                     current_production_controls_[ newIndex ] = prevState->currentProductionControl(oldIndex);
                 }
diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp
index 3a9bb8dc7..0b9c83099 100644
--- a/opm/simulators/wells/WellState.hpp
+++ b/opm/simulators/wells/WellState.hpp
@@ -148,10 +148,6 @@ public:
     static void calculateSegmentRates(const std::vector<std::vector<int>>& segment_inlets, const std::vector<std::vector<int>>&segment_perforations,
                                       const std::vector<double>& perforation_rates, const int np, const int segment, std::vector<double>& segment_rates);
 
-    Events& events(std::size_t well_index) {
-        return this->events_[well_index];
-    }
-
     /// One rate pr well
     double solventWellRate(const int w) const;
 
@@ -393,11 +389,6 @@ private:
     // should be zero for injection wells
     WellContainer<double> well_vaporized_oil_rates_;
 
-    // some events happens to the well, like this well is a new well
-    // or new well control keywords happens
-    // \Note: for now, only WCON* keywords, and well status change is considered
-    WellContainer<Events> events_;
-
     WellContainer<SegmentState> segment_state;
 
     // Productivity Index

From 581f571c9d24cc183824ee5b1b8bf09ded0cb10f Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Wed, 4 Aug 2021 12:03:36 +0200
Subject: [PATCH 7/7] Use SingleWellState for production and injection control

---
 .../wells/BlackoilWellModelGeneric.cpp        |  4 +--
 opm/simulators/wells/GasLiftGroupInfo.cpp     |  6 ++---
 .../wells/GasLiftSingleWellGeneric.cpp        |  5 ++--
 opm/simulators/wells/MultisegmentWellEval.cpp | 10 ++++---
 .../wells/MultisegmentWell_impl.hpp           |  9 ++++---
 opm/simulators/wells/SingleWellState.hpp      |  3 +++
 opm/simulators/wells/StandardWellGeneric.cpp  |  7 ++---
 opm/simulators/wells/StandardWell_impl.hpp    |  9 ++++---
 opm/simulators/wells/WellGroupHelpers.cpp     | 10 ++++---
 opm/simulators/wells/WellInterfaceEval.cpp    |  4 +--
 .../wells/WellInterfaceFluidSystem.cpp        | 26 ++++++++++---------
 opm/simulators/wells/WellInterface_impl.hpp   | 15 ++++++-----
 opm/simulators/wells/WellState.cpp            | 22 +++++++---------
 opm/simulators/wells/WellState.hpp            | 14 ----------
 14 files changed, 69 insertions(+), 75 deletions(-)

diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
index fdf4d0d68..7bc469dbe 100644
--- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp
+++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
@@ -178,10 +178,10 @@ loadRestartData(const data::Wells& rst_wells,
         ws.temperature = rst_well.temperature;
 
         if (rst_well.current_control.isProducer) {
-            well_state.currentProductionControl(well_index, rst_well.current_control.prod);
+            ws.production_cmode = rst_well.current_control.prod;
         }
         else {
-            well_state.currentInjectionControl(well_index, rst_well.current_control.inj);
+            ws.injection_cmode = rst_well.current_control.inj;
         }
 
         for( size_t i = 0; i < phs.size(); ++i ) {
diff --git a/opm/simulators/wells/GasLiftGroupInfo.cpp b/opm/simulators/wells/GasLiftGroupInfo.cpp
index ac1d18f42..4c4a98221 100644
--- a/opm/simulators/wells/GasLiftGroupInfo.cpp
+++ b/opm/simulators/wells/GasLiftGroupInfo.cpp
@@ -202,9 +202,9 @@ checkDoGasLiftOptimization_(const std::string &well_name)
         if (itr != this->ecl_wells_.end()) {
             //const Well *well = (itr->second).first;
             //assert(well); // Should never be nullptr
-            const int index = (itr->second).second;
-            const Well::ProducerCMode& control_mode
-                = this->well_state_.currentProductionControl(index);
+            const int well_index = (itr->second).second;
+            const auto& ws = this->well_state_.well(well_index);
+            const Well::ProducerCMode& control_mode = ws.production_cmode;
             if (control_mode != Well::ProducerCMode::THP ) {
                 displayDebugMessage_("Not THP control. Skipping.", well_name);
                 return false;
diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
index f7b9ddbc3..aadb2f0c2 100644
--- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
+++ b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
@@ -1235,9 +1235,8 @@ bool
 GasLiftSingleWellGeneric::OptimizeState::
 checkThpControl()
 {
-    const int index = this->parent.well_state_.wellIndex(this->parent.well_name_);
-    const Well::ProducerCMode& control_mode
-        = this->parent.well_state_.currentProductionControl(index);
+    const int well_index = this->parent.well_state_.wellIndex(this->parent.well_name_);
+    const Well::ProducerCMode& control_mode = this->parent.well_state_.well(well_index).production_cmode;
     return control_mode == Well::ProducerCMode::THP;
 }
 
diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp
index f158de2f0..5529fbe42 100644
--- a/opm/simulators/wells/MultisegmentWellEval.cpp
+++ b/opm/simulators/wells/MultisegmentWellEval.cpp
@@ -155,9 +155,10 @@ checkConvergenceControlEq(const WellState& well_state,
     CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
 
     const int well_index = baseif_.indexOfWell();
+    const auto& ws = well_state.well(well_index);
     if (baseif_.isInjector() )
     {
-        auto current = well_state.currentInjectionControl(well_index);
+        auto current = ws.injection_cmode;
         switch(current) {
         case Well::InjectorCMode::THP:
             ctrltype = CR::WellFailure::Type::ControlTHP;
@@ -183,7 +184,7 @@ checkConvergenceControlEq(const WellState& well_state,
 
     if (baseif_.isProducer() )
     {
-        auto current = well_state.currentProductionControl(well_index);
+        auto current = ws.production_cmode;
         switch(current) {
         case Well::ProducerCMode::THP:
             ctrltype = CR::WellFailure::Type::ControlTHP;
@@ -1797,9 +1798,10 @@ getControlTolerance(const WellState& well_state,
     double control_tolerance = 0.;
 
     const int well_index = baseif_.indexOfWell();
+    const auto& ws = well_state.well(well_index);
     if (baseif_.isInjector() )
     {
-        auto current = well_state.currentInjectionControl(well_index);
+        auto current = ws.injection_cmode;
         switch(current) {
         case Well::InjectorCMode::THP:
             control_tolerance = tolerance_pressure_ms_wells;
@@ -1821,7 +1823,7 @@ getControlTolerance(const WellState& well_state,
 
     if (baseif_.isProducer() )
     {
-        auto current = well_state.currentProductionControl(well_index);
+        auto current = ws.production_cmode;
         switch(current) {
         case Well::ProducerCMode::THP:
             control_tolerance = tolerance_pressure_ms_wells; // 0.1 bar
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index 1c2dcb2ef..2fbd2057a 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -257,8 +257,9 @@ namespace Opm
         // If the well is pressure controlled the potential equals the rate.
         bool thp_controlled_well = false;
         bool bhp_controlled_well = false;
+        const auto& ws = well_state.well(this->index_of_well_);
         if (this->isInjector()) {
-            const Well::InjectorCMode& current = well_state.currentInjectionControl(index_of_well_);
+            const Well::InjectorCMode& current = ws.injection_cmode;
             if (current == Well::InjectorCMode::THP) {
                 thp_controlled_well = true;
             }
@@ -266,7 +267,7 @@ namespace Opm
                 bhp_controlled_well = true;
             }
         } else {
-            const Well::ProducerCMode& current = well_state.currentProductionControl(index_of_well_);
+            const Well::ProducerCMode& current = ws.production_cmode;
             if (current == Well::ProducerCMode::THP) {
                 thp_controlled_well = true;
             }
@@ -354,10 +355,10 @@ namespace Opm
         //  Set current control to bhp, and bhp value in state, modify bhp limit in control object.
         if (well_copy.well_ecl_.isInjector()) {
             inj_controls.bhp_limit = bhp;
-            well_state_copy.currentInjectionControl(index_of_well_, Well::InjectorCMode::BHP);
+            ws.injection_cmode = Well::InjectorCMode::BHP;
         } else {
             prod_controls.bhp_limit = bhp;
-            well_state_copy.currentProductionControl(index_of_well_, Well::ProducerCMode::BHP);
+            ws.production_cmode = Well::ProducerCMode::BHP;
         }
         ws.bhp = bhp;
         well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
index 76a023216..f99ad13d6 100644
--- a/opm/simulators/wells/SingleWellState.hpp
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -35,6 +35,9 @@ public:
     double thp{0};
     double temperature{};
     Events events;
+    Well::InjectorCMode injection_cmode{Well::InjectorCMode::CMODE_UNDEFINED};
+    Well::ProducerCMode production_cmode{Well::ProducerCMode::CMODE_UNDEFINED};
+
 
     void init_timestep(const SingleWellState& other);
     void shut();
diff --git a/opm/simulators/wells/StandardWellGeneric.cpp b/opm/simulators/wells/StandardWellGeneric.cpp
index 639f42bfc..3cffc9e03 100644
--- a/opm/simulators/wells/StandardWellGeneric.cpp
+++ b/opm/simulators/wells/StandardWellGeneric.cpp
@@ -265,7 +265,7 @@ doGasLiftOptimize(const WellState &well_state,
     }
     if (glift_optimize_only_thp_wells) {
         const int well_index = baseif_.indexOfWell();
-        auto control_mode = well_state.currentProductionControl(well_index);
+        auto control_mode = well_state.well(well_index).production_cmode;
         if (control_mode != Well::ProducerCMode::THP ) {
             gliftDebug("Not THP control", deferred_logger);
             return false;
@@ -698,13 +698,14 @@ checkConvergenceControlEq(const WellState& well_state,
     CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
 
     const int well_index = baseif_.indexOfWell();
+    const auto& ws = well_state.well(well_index);
     if (baseif_.wellIsStopped()) {
         ctrltype = CR::WellFailure::Type::ControlRate;
         control_tolerance = 1.e-6; // use smaller tolerance for zero control?
     }
     else if (baseif_.isInjector() )
     {
-        auto current = well_state.currentInjectionControl(well_index);
+        auto current = ws.injection_cmode;
         switch(current) {
         case Well::InjectorCMode::THP:
             ctrltype = CR::WellFailure::Type::ControlTHP;
@@ -729,7 +730,7 @@ checkConvergenceControlEq(const WellState& well_state,
     }
     else if (baseif_.isProducer() )
     {
-        auto current = well_state.currentProductionControl(well_index);
+        auto current = ws.production_cmode;
         switch(current) {
         case Well::ProducerCMode::THP:
             ctrltype = CR::WellFailure::Type::ControlTHP;
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index eb5da692c..06c06c62a 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -1687,9 +1687,9 @@ namespace Opm
 
         //  Set current control to bhp, and bhp value in state, modify bhp limit in control object.
         if (well_ecl_.isInjector()) {
-            well_state_copy.currentInjectionControl(index_of_well_, Well::InjectorCMode::BHP);
+            ws.injection_cmode = Well::InjectorCMode::BHP;
         } else {
-            well_state_copy.currentProductionControl(index_of_well_, Well::ProducerCMode::BHP);
+            ws.production_cmode = Well::ProducerCMode::BHP;
         }
         ws.bhp = bhp;
 
@@ -1839,8 +1839,9 @@ namespace Opm
         // If the well is pressure controlled the potential equals the rate.
         bool thp_controlled_well = false;
         bool bhp_controlled_well = false;
+        const auto& ws = well_state.well(this->index_of_well_);
         if (this->isInjector()) {
-            const Well::InjectorCMode& current = well_state.currentInjectionControl(index_of_well_);
+            const Well::InjectorCMode& current = ws.injection_cmode;
             if (current == Well::InjectorCMode::THP) {
                 thp_controlled_well = true;
             }
@@ -1848,7 +1849,7 @@ namespace Opm
                 bhp_controlled_well = true;
             }
         } else {
-            const Well::ProducerCMode& current = well_state.currentProductionControl(index_of_well_);
+            const Well::ProducerCMode& current = ws.production_cmode;
             if (current == Well::ProducerCMode::THP) {
                 thp_controlled_well = true;
             }
diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp
index fbf7c06f5..192485043 100644
--- a/opm/simulators/wells/WellGroupHelpers.cpp
+++ b/opm/simulators/wells/WellGroupHelpers.cpp
@@ -404,13 +404,14 @@ namespace WellGroupHelpers
 
             const double efficiency = wellTmp.getEfficiencyFactor();
             // add contributino from wells not under group control
+            const auto& ws = wellState.well(well_index);
             if (isInjector) {
-                if (wellState.currentInjectionControl(well_index) != Well::InjectorCMode::GRUP)
+                if (ws.injection_cmode != Well::InjectorCMode::GRUP)
                     for (int phase = 0; phase < np; phase++) {
                         groupTargetReduction[phase] += wellStateNupcol.wellRates(well_index)[phase] * efficiency;
                     }
             } else {
-                if (wellState.currentProductionControl(well_index) != Well::ProducerCMode::GRUP)
+                if (ws.production_cmode != Well::ProducerCMode::GRUP)
                     for (int phase = 0; phase < np; phase++) {
                         groupTargetReduction[phase] -= wellStateNupcol.wellRates(well_index)[phase] * efficiency;
                     }
@@ -480,13 +481,14 @@ namespace WellGroupHelpers
             }
 
             // scale rates
+            const auto& ws = wellState.well(well_index);
             if (isInjector) {
-                if (wellState.currentInjectionControl(well_index) == Well::InjectorCMode::GRUP)
+                if (ws.injection_cmode == Well::InjectorCMode::GRUP)
                     for (int phase = 0; phase < np; phase++) {
                         wellState.wellRates(well_index)[phase] *= scale;
                     }
             } else {
-                if (wellState.currentProductionControl(well_index) == Well::ProducerCMode::GRUP)
+                if (ws.production_cmode == Well::ProducerCMode::GRUP)
                     for (int phase = 0; phase < np; phase++) {
                         wellState.wellRates(well_index)[phase] *= scale;
                     }
diff --git a/opm/simulators/wells/WellInterfaceEval.cpp b/opm/simulators/wells/WellInterfaceEval.cpp
index 180147371..5ac1cbea5 100644
--- a/opm/simulators/wells/WellInterfaceEval.cpp
+++ b/opm/simulators/wells/WellInterfaceEval.cpp
@@ -278,7 +278,7 @@ assembleControlEqProd_(const WellState& well_state,
                        EvalWell& control_eq,
                        DeferredLogger& deferred_logger) const
 {
-    auto current = well_state.currentProductionControl(baseif_.indexOfWell());
+    const auto current = well_state.well(baseif_.indexOfWell()).production_cmode;
     const auto& pu = baseif_.phaseUsage();
     const double efficiencyFactor = baseif_.wellEcl().getEfficiencyFactor();
 
@@ -392,7 +392,7 @@ assembleControlEqInj_(const WellState& well_state,
                       EvalWell& control_eq,
                       DeferredLogger& deferred_logger) const
 {
-    auto current = well_state.currentInjectionControl(baseif_.indexOfWell());
+    auto current = well_state.well(baseif_.indexOfWell()).injection_cmode;
     const InjectorType injectorType = controls.injector_type;
     const auto& pu = baseif_.phaseUsage();
     const double efficiencyFactor = baseif_.wellEcl().getEfficiencyFactor();
diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
index e9f440288..95bd904b7 100644
--- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp
+++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
@@ -90,8 +90,8 @@ activeProductionConstraint(const WellState& well_state,
     const PhaseUsage& pu = this->phaseUsage();
     const int well_index = this->index_of_well_;
     const auto controls = this->well_ecl_.productionControls(summaryState);
-    auto currentControl = well_state.currentProductionControl(well_index);
     auto& ws = well_state.well(well_index);
+    const auto currentControl = ws.production_cmode;
 
     if (controls.hasControl(Well::ProducerCMode::BHP) && currentControl != Well::ProducerCMode::BHP) {
         const double bhp_limit = controls.bhp_limit;
@@ -170,7 +170,7 @@ activeProductionConstraint(const WellState& well_state,
             return Well::ProducerCMode::THP;
     }
 
-    return well_state.currentProductionControl(well_index);
+    return currentControl;
 }
 
 
@@ -185,7 +185,7 @@ activeInjectionConstraint(const WellState& well_state,
     const auto& ws = well_state.well(well_index);
 
     const auto controls = this->well_ecl_.injectionControls(summaryState);
-    auto currentControl = well_state.currentInjectionControl(well_index);
+    const auto currentControl = ws.injection_cmode;
 
     if (controls.hasControl(Well::InjectorCMode::BHP) && currentControl != Well::InjectorCMode::BHP)
     {
@@ -248,7 +248,7 @@ activeInjectionConstraint(const WellState& well_state,
             return Well::InjectorCMode::THP;
     }
 
-    return well_state.currentInjectionControl(well_index);
+    return currentControl;
 }
 
 template <typename FluidSystem>
@@ -258,18 +258,19 @@ checkIndividualConstraints(WellState& well_state,
                            const SummaryState& summaryState) const
 {
     const int well_index = this->index_of_well_;
+    auto& ws = well_state.well(well_index);
     if (this->well_ecl_.isProducer()) {
         auto new_cmode = this->activeProductionConstraint(well_state, summaryState);
-        if (new_cmode != well_state.currentProductionControl(well_index)) {
-            well_state.currentProductionControl(well_index, new_cmode);
+        if (new_cmode != ws.production_cmode) {
+            ws.production_cmode = new_cmode;
             return true;
         }
     }
 
     if (this->well_ecl_.isInjector()) {
         auto new_cmode = this->activeInjectionConstraint(well_state, summaryState);
-        if (new_cmode != well_state.currentInjectionControl(well_index)) {
-            well_state.currentInjectionControl(well_index, new_cmode);
+        if (new_cmode != ws.injection_cmode) {
+            ws.injection_cmode = new_cmode;
             return true;
         }
     }
@@ -376,9 +377,10 @@ checkGroupConstraints(WellState& well_state,
 {
     const auto& well = well_ecl_;
     const int well_index = index_of_well_;
+    auto& ws = well_state.well(well_index);
 
     if (well.isInjector()) {
-        auto currentControl = well_state.currentInjectionControl(well_index);
+        const auto currentControl = ws.injection_cmode;
 
         if (currentControl != Well::InjectorCMode::GRUP) {
             // This checks only the first encountered group limit,
@@ -395,7 +397,7 @@ checkGroupConstraints(WellState& well_state,
             // If a group constraint was broken, we set the current well control to
             // be GRUP.
             if (group_constraint.first) {
-                well_state.currentInjectionControl(index_of_well_, Well::InjectorCMode::GRUP);
+                ws.injection_cmode = Well::InjectorCMode::GRUP;
                 const int np = well_state.numPhases();
                 for (int p = 0; p<np; ++p) {
                     well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
@@ -406,7 +408,7 @@ checkGroupConstraints(WellState& well_state,
     }
 
     if (well.isProducer( )) {
-        auto currentControl = well_state.currentProductionControl(well_index);
+        const auto currentControl = ws.production_cmode;
 
         if (currentControl != Well::ProducerCMode::GRUP) {
             // This checks only the first encountered group limit,
@@ -423,7 +425,7 @@ checkGroupConstraints(WellState& well_state,
             // If a group constraint was broken, we set the current well control to
             // be GRUP.
             if (group_constraint.first) {
-                well_state.currentProductionControl(index_of_well_, Well::ProducerCMode::GRUP);
+                ws.production_cmode = Well::ProducerCMode::GRUP;
                 const int np = well_state.numPhases();
                 for (int p = 0; p<np; ++p) {
                     well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index 63e39c47d..579863936 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -165,11 +165,12 @@ namespace Opm
         const auto& summaryState = ebos_simulator.vanguard().summaryState();
         const auto& schedule = ebos_simulator.vanguard().schedule();
         const auto& well = this->well_ecl_;
+        auto& ws = well_state.well(this->index_of_well_);
         std::string from;
         if (well.isInjector()) {
-            from = Well::InjectorCMode2String(well_state.currentInjectionControl(this->index_of_well_));
+            from = Well::InjectorCMode2String(ws.injection_cmode);
         } else {
-            from = Well::ProducerCMode2String(well_state.currentProductionControl(this->index_of_well_));
+            from = Well::ProducerCMode2String(ws.production_cmode);
         }
 
         bool changed = false;
@@ -188,9 +189,9 @@ namespace Opm
         if (changed) {
             std::string to;
             if (well.isInjector()) {
-                to = Well::InjectorCMode2String(well_state.currentInjectionControl(this->index_of_well_));
+                to = Well::InjectorCMode2String(ws.injection_cmode);
             } else {
-                to = Well::ProducerCMode2String(well_state.currentProductionControl(this->index_of_well_));
+                to = Well::ProducerCMode2String(ws.production_cmode);
             }
             std::ostringstream ss;
             ss << "    Switching control mode for well " << this->name()
@@ -573,7 +574,7 @@ namespace Opm
     {
         this->operability_status_.reset();
 
-        auto current_control = well_state.currentProductionControl(this->index_of_well_);
+        auto current_control = well_state.well(this->index_of_well_).production_cmode;
         // Operability checking is not free
         // Only check wells under BHP and THP control
         if(current_control == Well::ProducerCMode::BHP || current_control == Well::ProducerCMode::THP) {
@@ -639,7 +640,7 @@ namespace Opm
                 OPM_DEFLOG_THROW(std::runtime_error, "Expected WATER, OIL or GAS as type for injectors "  + this->name(), deferred_logger );
             }
 
-            auto current = well_state.currentInjectionControl(well_index);
+            const auto current = ws.injection_cmode;
 
             switch(current) {
             case Well::InjectorCMode::RATE:
@@ -722,7 +723,7 @@ namespace Opm
         //Producer
         else
         {
-            auto current = well_state.currentProductionControl(well_index);
+            const auto current = ws.production_cmode;
             const auto& controls = well.productionControls(summaryState);
             switch (current) {
             case Well::ProducerCMode::ORAT:
diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp
index 030e73834..87b79e1d8 100644
--- a/opm/simulators/wells/WellState.cpp
+++ b/opm/simulators/wells/WellState.cpp
@@ -296,19 +296,15 @@ void WellState::init(const std::vector<double>& cellPressures,
         this->well_vaporized_oil_rates_.add(wname, 0);
     }
 
-    current_injection_controls_.clear();
-    current_production_controls_.clear();
     for (int w = 0; w < nw; ++w) {
-        const auto& wname = wells_ecl[w].name();
-        current_production_controls_.add(wname, Well::ProducerCMode::CMODE_UNDEFINED);
-        current_injection_controls_.add(wname, Well::InjectorCMode::CMODE_UNDEFINED);
+        auto& ws = this->well(w);
         if (wells_ecl[w].isProducer()) {
             const auto controls = wells_ecl[w].productionControls(summary_state);
-            currentProductionControl(w, controls.cmode);
+            ws.production_cmode = controls.cmode;
         }
         else {
             const auto controls = wells_ecl[w].injectionControls(summary_state);
-            currentInjectionControl(w, controls.cmode);
+            ws.injection_cmode = controls.cmode;
         }
     }
 
@@ -361,8 +357,8 @@ void WellState::init(const std::vector<double>& cellPressures,
 
                 // If new target is set using WCONPROD, WCONINJE etc. we use the new control
                 if (!new_well.events.hasEvent(WellState::event_mask)) {
-                    current_injection_controls_[ newIndex ] = prevState->currentInjectionControl(oldIndex);
-                    current_production_controls_[ newIndex ] = prevState->currentProductionControl(oldIndex);
+                    new_well.injection_cmode = prev_well.injection_cmode;
+                    new_well.production_cmode = prev_well.production_cmode;
                 }
 
                 wellRates(w) = prevState->wellRates(oldIndex);
@@ -552,8 +548,8 @@ WellState::report(const int* globalCellIdxMap,
             auto& curr = well.current_control;
 
             curr.isProducer = ws.producer;
-            curr.prod = this->currentProductionControl(well_index);
-            curr.inj  = this->currentInjectionControl(well_index);
+            curr.prod = ws.production_cmode;
+            curr.inj  = ws.injection_cmode;
         }
 
         const auto& pwinfo = *this->parallel_well_info_[well_index];
@@ -922,9 +918,9 @@ void WellState::updateGlobalIsGrup(const Comm& comm)
     for (std::size_t well_index = 0; well_index < this->size(); well_index++) {
         const auto& ws = this->well(well_index);
         if (ws.producer)
-            this->global_well_info.value().update_producer(well_index, ws.status, this->current_production_controls_[well_index]);
+            this->global_well_info.value().update_producer(well_index, ws.status, ws.production_cmode);
         else
-            this->global_well_info.value().update_injector(well_index, ws.status, this->current_injection_controls_[well_index]);
+            this->global_well_info.value().update_injector(well_index, ws.status, ws.injection_cmode);
     }
     this->global_well_info.value().communicate(comm);
 }
diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp
index 0b9c83099..248ac7c3d 100644
--- a/opm/simulators/wells/WellState.hpp
+++ b/opm/simulators/wells/WellState.hpp
@@ -108,14 +108,6 @@ public:
                 const std::vector<std::vector<PerforationData>>& well_perf_data,
                 const SummaryState& summary_state);
 
-    /// One current control per injecting well.
-    Well::InjectorCMode currentInjectionControl(std::size_t well_index) const { return current_injection_controls_[well_index]; }
-    void currentInjectionControl(std::size_t well_index, Well::InjectorCMode cmode) { current_injection_controls_[well_index] = cmode; }
-
-    /// One current control per producing well.
-    Well::ProducerCMode currentProductionControl(std::size_t well_index) const { return current_production_controls_[well_index]; }
-    void currentProductionControl(std::size_t well_index, Well::ProducerCMode cmode) { current_production_controls_[well_index] = cmode; }
-
     void setCurrentWellRates(const std::string& wellName, const std::vector<double>& new_rates ) {
         auto& [owner, rates] = this->well_rates.at(wellName);
         if (owner)
@@ -366,12 +358,6 @@ private:
     PhaseUsage phase_usage_;
     WellContainer<PerfData> perfdata;
 
-    // vector with size number of wells +1.
-    // iterate over all perforations of a given well
-    // for (int perf = first_perf_index_[well_index]; perf < first_perf_index_[well_index] + num_perf_[well_index]; ++perf)
-    WellContainer<Opm::Well::InjectorCMode> current_injection_controls_;
-    WellContainer<Well::ProducerCMode> current_production_controls_;
-
     // The well_rates variable is defined for all wells on all processors. The
     // bool in the value pair is whether the current process owns the well or
     // not.