Update WellState::bhp()

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1954,7 +1954,7 @@ namespace Opm {
         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.bhp()[ well_index ] = rst_well.bhp;
+            well_state.update_bhp( well_index, rst_well.bhp);
             well_state.temperature()[ well_index ] = rst_well.temperature;
 
             if (rst_well.current_control.isProducer) {

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -318,7 +318,7 @@ namespace Opm
     {
         //scale segment pressures
         const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
-        const double bhp = well_state.bhp()[index_of_well_];
+        const double bhp = well_state.bhp(index_of_well_);
         const double unscaled_top_seg_pressure = well_state.segPress()[top_segment_index];
         for (int seg = 0; seg < numberOfSegments(); ++seg) {
             const int seg_index = top_segment_index + seg;
@@ -632,7 +632,7 @@ namespace Opm
             prod_controls.bhp_limit = bhp;
             well_state_copy.currentProductionControls()[index_of_well_] = Well::ProducerCMode::BHP;
         }
-        well_state_copy.bhp()[well_copy.index_of_well_] = bhp;
+        well_state_copy.update_bhp(well_copy.index_of_well_, bhp);
         well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
 
         // initialized the well rates with the potentials i.e. the well rates based on bhp
@@ -1857,7 +1857,7 @@ namespace Opm
             rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
         }
 
-        const double bhp = well_state.bhp()[index_of_well_];
+        const double bhp = well_state.bhp(index_of_well_);
 
         well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, bhp, deferred_logger);
 
@@ -2383,7 +2383,7 @@ namespace Opm
             // update the segment pressure
             well_state.segPress()[seg + top_segment_index] = primary_variables_[seg][SPres];
             if (seg == 0) { // top segment
-                well_state.bhp()[index_of_well_] = well_state.segPress()[seg + top_segment_index];
+                well_state.update_bhp(index_of_well_, well_state.segPress()[seg + top_segment_index]);
             }
         }
         updateThp(well_state, deferred_logger);

opm/simulators/wells/StandardWell_impl.hpp

@@ -839,7 +839,7 @@ namespace Opm
 
         // Store the perforation pressure for later usage.
         auto * perf_press = &well_state.perfPress()[first_perf_];
-        perf_press[perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
+        perf_press[perf] = well_state.bhp(index_of_well_) + perf_pressure_diffs_[perf];
     }
 
 
@@ -1248,7 +1248,7 @@ namespace Opm
             F[pu.phase_pos[Gas]] += F_solvent;
         }
 
-        well_state.bhp()[index_of_well_] = primary_variables_[Bhp];
+        well_state.update_bhp(index_of_well_, 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
@@ -1319,7 +1319,7 @@ namespace Opm
             rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
         }
 
-        const double bhp = well_state.bhp()[index_of_well_];
+        const double bhp = well_state.bhp(index_of_well_);
 
         well_state.thp()[index_of_well_] = calculateThpFromBhp(well_state, rates, bhp, deferred_logger);
 
@@ -1580,7 +1580,7 @@ namespace Opm
                                    const WellState& well_state,
                                    Opm::DeferredLogger& deferred_logger)
     {
-        const double bhp = well_state.bhp()[index_of_well_];
+        const double bhp = well_state.bhp(index_of_well_);
         std::vector<double> well_rates;
         computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
 
@@ -1648,7 +1648,7 @@ namespace Opm
 
         // Compute the average pressure in each well block
         const auto * perf_press = &well_state.perfPress()[first_perf_];
-        auto p_above =  this->parallel_well_info_.communicateAboveValues(well_state.bhp()[w],
+        auto p_above =  this->parallel_well_info_.communicateAboveValues(well_state.bhp(w),
                                                                          perf_press,
                                                                          nperf);
 
@@ -2409,7 +2409,7 @@ namespace Opm
         } else {
             well_state_copy.currentProductionControls()[index_of_well_] = Well::ProducerCMode::BHP;
         }
-        well_state_copy.bhp()[index_of_well_] = bhp;
+        well_state_copy.update_bhp(index_of_well_, bhp);
 
         const double dt = ebosSimulator.timeStepSize();
         bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state_copy, group_state, deferred_logger);
@@ -2782,7 +2782,7 @@ namespace Opm
 
 
         // BHP
-        primary_variables_[Bhp] = well_state.bhp()[index_of_well_];
+        primary_variables_[Bhp] = well_state.bhp(index_of_well_);
 
         // other primary variables related to polymer injection
         if (this->has_polymermw && this->isInjector()) {

opm/simulators/wells/WellInterface_impl.hpp

@@ -1722,7 +1722,7 @@ namespace Opm
                     rates[p] = well_state.wellRates()[well_index*np + p];
                 }
                 double bhp = calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
-                well_state.bhp()[well_index] = bhp;
+                well_state.update_bhp(well_index, bhp);
 
                 // if the total rates are negative or zero
                 // we try to provide a better intial well rate
@@ -1737,7 +1737,7 @@ namespace Opm
             }
             case Well::InjectorCMode::BHP:
             {
-                well_state.bhp()[well_index] = controls.bhp_limit;
+                well_state.update_bhp(well_index, controls.bhp_limit);
                 double total_rate = 0.0;
                 for (int p = 0; p<np; ++p) {
                     total_rate += well_state.wellRates()[well_index*np + p];
@@ -1910,7 +1910,7 @@ namespace Opm
             }
             case Well::ProducerCMode::BHP:
             {
-                well_state.bhp()[well_index] = controls.bhp_limit;
+                well_state.update_bhp(well_index, controls.bhp_limit);
                 double total_rate = 0.0;
                 for (int p = 0; p<np; ++p) {
                     total_rate -= well_state.wellRates()[well_index*np + p];
@@ -1932,7 +1932,7 @@ namespace Opm
                     rates[p] = well_state.wellRates()[well_index*np + p];
                 }
                 double bhp = calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
-                well_state.bhp()[well_index] = bhp;
+                well_state.update_bhp(well_index, bhp);
 
                 // if the total rates are negative or zero
                 // we try to provide a better intial well rate
@@ -2052,7 +2052,7 @@ namespace Opm
             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 = well_state.bhp(well_index);
                 if (bhp < current_bhp) {
                     currentControl = Well::InjectorCMode::BHP;
                     return true;
@@ -2128,7 +2128,7 @@ namespace Opm
             if (controls.hasControl(Well::ProducerCMode::BHP) && currentControl != Well::ProducerCMode::BHP )
             {
                 const double bhp = controls.bhp_limit;
-                double current_bhp = well_state.bhp()[well_index];
+                double current_bhp = well_state.bhp(well_index);
                 if (bhp > current_bhp) {
                     currentControl = Well::ProducerCMode::BHP;
                     return true;

opm/simulators/wells/WellState.cpp

@@ -192,7 +192,7 @@ data::Wells WellState::report(const int* globalCellIdxMap,
         using WellT = std::remove_reference_t<decltype(dw[ itr.first ])>;
         WellT dummyWell; // dummy if we are not owner
         auto& well = pwinfo.isOwner() ? dw[ itr.first ] : dummyWell;
-        well.bhp = this->bhp().at( well_index );
+        well.bhp = this->bhp(well_index);
         well.thp = this->thp().at( well_index );
         well.temperature = this->temperature().at( well_index );
 

opm/simulators/wells/WellState.hpp

@@ -80,8 +80,8 @@ public:
                                      const std::vector<PerforationData>& well_perf_data);
 
     /// One bhp pressure per well.
-    std::vector<double>& bhp() { return bhp_; }
+    void update_bhp(std::size_t well_index, double value) { bhp_[well_index] = value; }
-    const std::vector<double>& bhp() const { return bhp_; }
+    double bhp(std::size_t well_index) const { return bhp_[well_index]; }
 
     /// One thp pressure per well.
     std::vector<double>& thp() { return thp_; }
@@ -116,7 +116,7 @@ public:
     /// The number of wells present.
     int numWells() const
     {
-        return bhp().size();
+        return wellMap_.size();
     }
 
     /// The number of phases present.

opm/simulators/wells/WellStateFullyImplicitBlackoil.cpp

@@ -187,7 +187,7 @@ void WellStateFullyImplicitBlackoil::init(const std::vector<double>& cellPressur
                 }
 
                 // bhp
-                bhp()[ newIndex ] = prevState->bhp()[ oldIndex ];
+                this->update_bhp( newIndex, prevState->bhp( oldIndex ));
 
                 // thp
                 thp()[ newIndex ] = prevState->thp()[ oldIndex ];
@@ -329,11 +329,12 @@ void WellStateFullyImplicitBlackoil::init(const std::vector<double>& cellPressur
         nseg_ = nw;
         top_segment_index_.resize(nw);
         seg_number_.resize(nw);
+        seg_press_.resize(nw);
         for (int w = 0; w < nw; ++w) {
             top_segment_index_[w] = w;
             seg_number_[w] = 1; // Top segment is segment #1
+            this->seg_press_[w] = this->bhp(w);
         }
-        seg_press_ = bhp();
         seg_rates_ = wellRates();
 
         seg_pressdrop_.assign(nw, 0.);
@@ -562,7 +563,7 @@ void WellStateFullyImplicitBlackoil::initWellStateMSWell(const std::vector<Well>
         if ( !well_ecl.isMultiSegment() ) { // not multi-segment well
             nseg_ += 1;
             seg_number_.push_back(1); // Assign single segment (top) as number 1.
-            seg_press_.push_back(bhp()[w]);
+            seg_press_.push_back(bhp(w));
             for (int p = 0; p < np; ++p) {
                 seg_rates_.push_back(wellRates()[np * w + p]);
             }
@@ -638,7 +639,7 @@ void WellStateFullyImplicitBlackoil::initWellStateMSWell(const std::vector<Well>
             // improved during the solveWellEq process
             {
                 // top segment is always the first one, and its pressure is the well bhp
-                seg_press_.push_back(bhp()[w]);
+                seg_press_.push_back(bhp(w));
                 const int top_segment = top_segment_index_[w];
                 const int start_perf = connpos;
                 const auto * perf_press = &this->perfPress()[start_perf];