WellState will internalize PhaseUsage member

opm/simulators/wells/BlackoilWellModel.hpp

@@ -329,7 +329,7 @@ namespace Opm {
 
             Opm::data::Wells wellData() const
             {
-                auto wsrpt = this->wellState().report(phase_usage_, Opm::UgGridHelpers::globalCell(grid()),
+                auto wsrpt = this->wellState().report(Opm::UgGridHelpers::globalCell(grid()),
                                                       [this](const int well_ndex) -> bool
                                                       {
                                                           return this->wasDynamicallyShutThisTimeStep(well_ndex);
@@ -632,6 +632,7 @@ namespace Opm {
 
              void computeWellTemperature();                       
         private:
+            BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& pu);
             /*
               The various wellState members should be accessed and modified
               through the accessor functions wellState(), prevWellState(),

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -32,12 +32,12 @@
 namespace Opm {
     template<typename TypeTag>
     BlackoilWellModel<TypeTag>::
-    BlackoilWellModel(Simulator& ebosSimulator)
+    BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& phase_usage)
         : ebosSimulator_(ebosSimulator),
-          phase_usage_(phaseUsageFromDeck(ebosSimulator_.vanguard().eclState())),
+          phase_usage_(phase_usage),
-          active_well_state_(phase_usage_.num_phases),
+          active_well_state_(phase_usage),
-          last_valid_well_state_(phase_usage_.num_phases),
+          last_valid_well_state_(phase_usage),
-          nupcol_well_state_(phase_usage_.num_phases)
+          nupcol_well_state_(phase_usage)
     {
         terminal_output_ = false;
         if (ebosSimulator.gridView().comm().rank() == 0)
@@ -74,6 +74,13 @@ namespace Opm {
         alternative_well_rate_init_ = EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
     }
 
+    template<typename TypeTag>
+    BlackoilWellModel<TypeTag>::
+    BlackoilWellModel(Simulator& ebosSimulator) :
+        BlackoilWellModel(ebosSimulator, phaseUsageFromDeck(ebosSimulator.vanguard().eclState()))
+    {}
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -268,7 +275,7 @@ namespace Opm {
 
         // handling MS well related
         if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
-            this->wellState().initWellStateMSWell(wells_ecl_, phase_usage_, &this->prevWellState());
+            this->wellState().initWellStateMSWell(wells_ecl_, &this->prevWellState());
         }
 
         const Group& fieldGroup = schedule().getGroup("FIELD", timeStepIdx);
@@ -630,7 +637,7 @@ namespace Opm {
             const auto phaseUsage = phaseUsageFromDeck(eclState());
             const size_t numCells = Opm::UgGridHelpers::numCells(grid());
             const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal());
-            this->wellState().resize(wells_ecl_, local_parallel_well_info_, schedule(), handle_ms_well, numCells, phaseUsage, well_perf_data_, summaryState, globalNumWells); // Resize for restart step
+            this->wellState().resize(wells_ecl_, local_parallel_well_info_, schedule(), handle_ms_well, numCells, well_perf_data_, summaryState, globalNumWells); // Resize for restart step
             loadRestartData(restartValues.wells, restartValues.grp_nwrk, phaseUsage, handle_ms_well, this->wellState());
         }
 
@@ -758,7 +765,7 @@ namespace Opm {
         }
 
         this->wellState().init(cellPressures, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx,
-                               &this->prevWellState(), phase_usage_, well_perf_data_,
+                               &this->prevWellState(), well_perf_data_,
                                summaryState, globalNumWells);
     }
 

opm/simulators/wells/WellState.hpp

@@ -48,8 +48,8 @@ namespace Opm
 
 
 
-        explicit WellState(int num_phases) :
+        explicit WellState(const PhaseUsage& pu) :
-            np_(num_phases)
+            phase_usage_(pu)
         {}
 
 
@@ -61,7 +61,6 @@ namespace Opm
         void init(const std::vector<double>& cellPressures,
                   const std::vector<Well>& wells_ecl,
                   const std::vector<ParallelWellInfo*>& parallel_well_info,
-                  const PhaseUsage& pu,
                   const std::vector<std::vector<PerforationData>>& well_perf_data,
                   const SummaryState& summary_state)
         {
@@ -74,9 +73,8 @@ namespace Opm
             {
                 // const int nw = wells->number_of_wells;
                 const int nw = wells_ecl.size();
+                const int np = this->phase_usage_.num_phases;
                 // const int np = wells->number_of_phases;
-                const int np = pu.num_phases;
-                np_ = np;
                 status_.assign(nw, Well::Status::OPEN);
                 bhp_.resize(nw, 0.0);
                 thp_.resize(nw, 0.0);
@@ -87,7 +85,7 @@ namespace Opm
                     const Well& well = wells_ecl[w];
               
                     // Initialize bhp(), thp(), wellRates(), temperature().
-                    initSingleWell(cellPressures, w, well, *parallel_well_info[w], pu, summary_state);
+                    initSingleWell(cellPressures, w, well, *parallel_well_info[w], summary_state);
 
                     // Setup wellname -> well index mapping.
                     const int num_perf_this_well = well_perf_data[w].size();
@@ -235,9 +233,15 @@ namespace Opm
         /// The number of phases present.
         int numPhases() const
         {
-            return np_;
+            return this->phase_usage_.num_phases;
         }
 
+        const PhaseUsage& phaseUsage() const {
+            return this->phase_usage_;
+        }
+
+
+
         void openWell(int well_index) {
             this->status_[well_index] = Well::Status::OPEN;
         }
@@ -273,12 +277,12 @@ namespace Opm
         }
 
         virtual data::Wells
-        report(const PhaseUsage& pu,
+        report(const int* globalCellIdxMap,
-               const int* globalCellIdxMap,
                const std::function<bool(const int)>& wasDynamicallyClosed) const
         {
             using rt = data::Rates::opt;
 
+            const auto& pu = this->phaseUsage();
             data::Wells dw;
             for( const auto& itr : this->wellMap_ ) {
                 const auto well_index = itr.second[ 0 ];
@@ -354,13 +358,13 @@ namespace Opm
         WellState& operator=(const WellState& rhs) = default;
 
     private:
+        PhaseUsage phase_usage_;
         std::vector<double> bhp_;
         std::vector<double> thp_;
         std::vector<double> temperature_;
         std::vector<double> wellrates_;
         std::vector<double> perfrates_;
         std::vector<double> perfpress_;
-        int np_;
     protected:
         std::vector<Well::Status> status_;
     private:
@@ -397,15 +401,16 @@ namespace Opm
                             const int w,
                             const Well& well,
                             const ParallelWellInfo& well_info,
-                            const PhaseUsage& pu,
                             const SummaryState& summary_state)
         {
             assert(well.isInjector() || well.isProducer());
 
             // Set default zero initial well rates.
             // May be overwritten below.
-            for (int p = 0; p < this->np_; ++p) {
+            const auto& pu = this->phase_usage_;
-                wellrates_[this->np_*w + p] = 0.0;
+            const int np = pu.num_phases;
+            for (int p = 0; p < np; ++p) {
+                wellrates_[np*w + p] = 0.0;
             }
 
             if ( well.isInjector() ) { 
@@ -471,15 +476,15 @@ namespace Opm
                         switch (inj_controls.injector_type) {
                         case InjectorType::WATER:
                             assert(pu.phase_used[BlackoilPhases::Aqua]);
-                            wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
+                            wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
                             break;
                         case InjectorType::GAS:
                             assert(pu.phase_used[BlackoilPhases::Vapour]);
-                            wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
+                            wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
                             break;
                         case InjectorType::OIL:
                             assert(pu.phase_used[BlackoilPhases::Liquid]);
-                            wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
+                            wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
                             break;
                         case InjectorType::MULTI:
                             // Not currently handled, keep zero init.
@@ -494,15 +499,15 @@ namespace Opm
                     switch (prod_controls.cmode) {
                     case Well::ProducerCMode::ORAT:
                         assert(pu.phase_used[BlackoilPhases::Liquid]);
-                        wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
+                        wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
                         break;
                     case Well::ProducerCMode::WRAT:
                         assert(pu.phase_used[BlackoilPhases::Aqua]);
-                        wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
+                        wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
                         break;
                     case Well::ProducerCMode::GRAT:
                         assert(pu.phase_used[BlackoilPhases::Vapour]);
-                        wellrates_[this->np_*w + pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
+                        wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
                         break;
                     default:
                         // Keep zero init.

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -69,9 +69,9 @@ namespace Opm
         using BaseType :: resetConnectionTransFactors;
         using BaseType :: updateStatus;
 
-        explicit WellStateFullyImplicitBlackoil(int num_phases) :
+        explicit WellStateFullyImplicitBlackoil(const PhaseUsage& pu) :
-            WellState(num_phases),
+            WellState(pu),
-            group_state(num_phases)
+            group_state(pu.num_phases)
         {}
 
 
@@ -84,13 +84,12 @@ namespace Opm
                   const std::vector<ParallelWellInfo*>& parallel_well_info,
                   const int report_step,
                   const WellStateFullyImplicitBlackoil* prevState,
-                  const PhaseUsage& pu,
                   const std::vector<std::vector<PerforationData>>& well_perf_data,
                   const SummaryState& summary_state,
                   const int globalNumberOfWells)
         {
             // call init on base class
-            BaseType :: init(cellPressures, wells_ecl, parallel_well_info, pu, well_perf_data, summary_state);
+            BaseType :: init(cellPressures, wells_ecl, parallel_well_info, well_perf_data, summary_state);
 
             for (const auto& winfo: parallel_well_info)
             {
@@ -105,6 +104,7 @@ namespace Opm
             if( nw == 0 ) return ;
 
             // Initialize perfphaserates_, which must be done here.
+            const auto& pu = this->phaseUsage();
             const int np = pu.num_phases;
 
             int nperf = 0;
@@ -387,16 +387,15 @@ namespace Opm
                     const Schedule& schedule,
                     const bool handle_ms_well,
                     const size_t numCells,
-                    const PhaseUsage& pu,
                     const std::vector<std::vector<PerforationData>>& well_perf_data,
                     const SummaryState& summary_state,
                     const int globalNumWells)
         {
             const std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
-            init(tmp, schedule, wells_ecl, parallel_well_info, 0, nullptr, pu, well_perf_data, summary_state, globalNumWells);
+            init(tmp, schedule, wells_ecl, parallel_well_info, 0, nullptr, well_perf_data, summary_state, globalNumWells);
 
             if (handle_ms_well) {
-                initWellStateMSWell(wells_ecl, pu, nullptr);
+                initWellStateMSWell(wells_ecl, nullptr);
             }
         }
 
@@ -530,18 +529,18 @@ namespace Opm
         }
 
         data::Wells
-        report(const PhaseUsage &pu,
+        report(const int* globalCellIdxMap,
-               const int* globalCellIdxMap,
                const std::function<bool(const int)>& wasDynamicallyClosed) const override
         {
             data::Wells res =
-                WellState::report(pu, globalCellIdxMap, wasDynamicallyClosed);
+                WellState::report(globalCellIdxMap, wasDynamicallyClosed);
 
             const int nw = this->numWells();
             if (nw == 0) {
                 return res;
             }
 
+            const auto& pu = this->phaseUsage();
             const int np = pu.num_phases;
 
             using rt = data::Rates::opt;
@@ -689,10 +688,12 @@ namespace Opm
 
         /// init the MS well related.
         void initWellStateMSWell(const std::vector<Well>& wells_ecl,
-                                 const PhaseUsage& pu, const WellStateFullyImplicitBlackoil* prev_well_state)
+                                 const WellStateFullyImplicitBlackoil* prev_well_state)
         {
             // still using the order in wells
             const int nw = wells_ecl.size();
+            const auto& pu = this->phaseUsage();
+            const int np = pu.num_phases;
             if (nw == 0) {
                 return;
             }
@@ -720,7 +721,6 @@ namespace Opm
                     nseg_ += 1;
                     seg_number_.push_back(1); // Assign single segment (top) as number 1.
                     seg_press_.push_back(bhp()[w]);
-                    const int np = numPhases();
                     for (int p = 0; p < np; ++p) {
                         seg_rates_.push_back(wellRates()[np * w + p]);
                     }
@@ -763,7 +763,6 @@ namespace Opm
 
                     // for the seg_rates_, now it becomes a recursive solution procedure.
                     {
-                        const int np = numPhases();
                         const int start_perf = connpos;
                         const int start_perf_next_well = connpos + num_perf_this_well;
 
@@ -825,7 +824,6 @@ namespace Opm
 
             if (prev_well_state && !prev_well_state->wellMap().empty()) {
                 const auto& end = prev_well_state->wellMap().end();
-                const int np = numPhases();
                 for (int w = 0; w < nw; ++w) {
                     const Well& well = wells_ecl[w];
                     if (well.getStatus() == Well::Status::SHUT)

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{setup.pu.num_phases};
+        auto state  = Opm::WellStateFullyImplicitBlackoil{setup.pu};
 
         const auto cpress =
             std::vector<double>(setup.grid.c_grid()->number_of_cells,
@@ -144,11 +144,11 @@ namespace {
 
         state.init(cpress, setup.sched,
                    wells, ppinfos,
-                   timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st,
+                   timeStep, nullptr, setup.well_perf_data, setup.st,
                    wells.size());
 
         state.initWellStateMSWell(setup.sched.getWells(timeStep),
-                                  setup.pu, nullptr);
+                                  nullptr);
 
         return state;
     }
@@ -274,7 +274,7 @@ BOOST_AUTO_TEST_CASE(Pressure)
 
     setSegPress(wells, wstate);
 
-    const auto rpt = wstate.report(setup.pu, setup.grid.c_grid()->global_cell, [](const int){return false;});
+    const auto rpt = wstate.report(setup.grid.c_grid()->global_cell, [](const int){return false;});
 
     {
         const auto& xw = rpt.at("INJE01");
@@ -321,9 +321,9 @@ BOOST_AUTO_TEST_CASE(Rates)
 
     const auto& pu = setup.pu;
 
-    setSegRates(wells, pu, wstate);
+    setSegRates(wells, pu,  wstate);
 
-    const auto rpt = wstate.report(pu, setup.grid.c_grid()->global_cell, [](const int){return false;});
+    const auto rpt = wstate.report(setup.grid.c_grid()->global_cell, [](const int){return false;});
 
     const auto wat = pu.phase_used[Opm::BlackoilPhases::Aqua];
     const auto oil = pu.phase_used[Opm::BlackoilPhases::Liquid];