Assign number of phases in WellState constructor

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -33,7 +33,11 @@ namespace Opm {
     template<typename TypeTag>
     BlackoilWellModel<TypeTag>::
     BlackoilWellModel(Simulator& ebosSimulator)
-        : ebosSimulator_(ebosSimulator)
+        : ebosSimulator_(ebosSimulator),
+          phase_usage_(phaseUsageFromDeck(ebosSimulator_.vanguard().eclState())),
+          active_well_state_(phase_usage_.num_phases),
+          last_valid_well_state_(phase_usage_.num_phases),
+          nupcol_well_state_(phase_usage_.num_phases)
     {
         terminal_output_ = false;
         if (ebosSimulator.gridView().comm().rank() == 0)
@@ -75,13 +79,9 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     init()
     {
-        const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState();
-
         extractLegacyCellPvtRegionIndex_();
         extractLegacyDepth_();
 
-        phase_usage_ = phaseUsageFromDeck(eclState);
-
         gravity_ = ebosSimulator_.problem().gravity()[2];
 
         initial_step_ = true;

opm/simulators/wells/WellState.hpp

@@ -46,6 +46,13 @@ namespace Opm
         typedef std::array< int, 3 >  mapentry_t;
         typedef std::map< std::string, mapentry_t > WellMapType;
 
+
+
+        explicit WellState(int num_phases) :
+            np_(num_phases)
+        {}
+
+
         /// Allocate and initialize if wells is non-null.
         /// Also tries to give useful initial values to the bhp() and
         /// wellRates() fields, depending on controls.  The
@@ -397,9 +404,8 @@ namespace Opm
 
             // Set default zero initial well rates.
             // May be overwritten below.
-            const int np = pu.num_phases;
+            for (int p = 0; p < this->np_; ++p) {
-            for (int p = 0; p < np; ++p) {
+                wellrates_[this->np_*w + p] = 0.0;
-                wellrates_[np*w + p] = 0.0;
             }
 
             if ( well.isInjector() ) { 
@@ -465,15 +471,15 @@ namespace Opm
                         switch (inj_controls.injector_type) {
                         case InjectorType::WATER:
                             assert(pu.phase_used[BlackoilPhases::Aqua]);
-                            wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
+                            wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
                             break;
                         case InjectorType::GAS:
                             assert(pu.phase_used[BlackoilPhases::Vapour]);
-                            wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
+                            wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
                             break;
                         case InjectorType::OIL:
                             assert(pu.phase_used[BlackoilPhases::Liquid]);
-                            wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
+                            wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
                             break;
                         case InjectorType::MULTI:
                             // Not currently handled, keep zero init.
@@ -488,15 +494,15 @@ namespace Opm
                     switch (prod_controls.cmode) {
                     case Well::ProducerCMode::ORAT:
                         assert(pu.phase_used[BlackoilPhases::Liquid]);
-                        wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
+                        wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
                         break;
                     case Well::ProducerCMode::WRAT:
                         assert(pu.phase_used[BlackoilPhases::Aqua]);
-                        wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
+                        wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
                         break;
                     case Well::ProducerCMode::GRAT:
                         assert(pu.phase_used[BlackoilPhases::Vapour]);
-                        wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
+                        wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
                         break;
                     default:
                         // Keep zero init.

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -68,6 +68,11 @@ namespace Opm
         using BaseType :: resetConnectionTransFactors;
         using BaseType :: updateStatus;
 
+        explicit WellStateFullyImplicitBlackoil(int num_phases) :
+            WellState(num_phases)
+        {}
+
+
         /// Allocate and initialize if wells is non-null.  Also tries
         /// to give useful initial values to the bhp(), wellRates()
         /// and perfPhaseRates() fields, depending on controls
@@ -105,7 +110,7 @@ namespace Opm
                 nperf += wpd.size();
             }
 
-            well_reservoir_rates_.resize(nw * np, 0.0);
+            well_reservoir_rates_.resize(nw * this->numPhases(), 0.0);
             well_dissolved_gas_rates_.resize(nw, 0.0);
             well_vaporized_oil_rates_.resize(nw, 0.0);
 
@@ -129,7 +134,7 @@ namespace Opm
 
             // Ensure that we start out with zero rates by default.
             perfphaserates_.clear();
-            perfphaserates_.resize(nperf * np, 0.0);
+            perfphaserates_.resize(nperf * this->numPhases(), 0.0);
 
             // these are only used to monitor the injectivity
             perf_water_throughput_.clear();
@@ -152,8 +157,8 @@ namespace Opm
 
                 for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) {
                     if (wells_ecl[w].getStatus() == Well::Status::OPEN) {
-                        for (int p = 0; p < np; ++p) {
+                        for (int p = 0; p < this->numPhases(); ++p) {
-                            perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(global_num_perf_this_well);
+                            perfphaserates_[this->numPhases()*perf + p] = wellRates()[this->numPhases()*w + p] / double(global_num_perf_this_well);
                         }
                     }
                     perfPress()[perf] = cellPressures[well_perf_data[w][perf-connpos].cell_index];
@@ -182,9 +187,9 @@ namespace Opm
 
             perfRateSolvent_.clear();
             perfRateSolvent_.resize(nperf, 0.0);
-            productivity_index_.resize(nw * np, 0.0);
+            productivity_index_.resize(nw * this->numPhases(), 0.0);
-            conn_productivity_index_.resize(nperf * np, 0.0);
+            conn_productivity_index_.resize(nperf * this->numPhases(), 0.0);
-            well_potentials_.resize(nw * np, 0.0);
+            well_potentials_.resize(nw * this->numPhases(), 0.0);
 
             perfRatePolymer_.clear();
             perfRatePolymer_.resize(nperf, 0.0);

tests/test_wellstatefullyimplicitblackoil.cpp

@@ -121,7 +121,7 @@ namespace {
     buildWellState(const Setup& setup, const std::size_t timeStep,
                    std::vector<Opm::ParallelWellInfo>& pinfos)
     {
-        auto state  = Opm::WellStateFullyImplicitBlackoil{};
+        auto state  = Opm::WellStateFullyImplicitBlackoil{setup.pu.num_phases};
 
         const auto cpress =
             std::vector<double>(setup.grid.c_grid()->number_of_cells,