recovering well restart data for multisegment well

When multisegment wells are involved for RESTART, based on the recovered
segment-related quantities as Opm::data::Segment objects, the segment
phase rates and segment pressure are recovered as part of the
WellStateFullyImplicitBlackoil.

opm/autodiff/BlackoilWellModel.hpp

@@ -423,9 +423,14 @@ namespace Opm {
 
             // convert well data from opm-common to well state from opm-core
             void wellsToState( const data::Wells& wells,
-                               PhaseUsage phases,
+                               const PhaseUsage& phases,
-                               WellStateFullyImplicitBlackoil& state );
+                               const bool handle_ms_well,
+                               const int report_step,
+                               WellStateFullyImplicitBlackoil& state ) const;
 
+            bool anyMSWellOpen(const Wells* wells, const int report_step) const;
+
+            const Well* getWellEcl(const std::string& well_name) const;
         };
 
 

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -259,16 +259,8 @@ namespace Opm {
         well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_);
 
         // handling MS well related
-        if (param_.use_multisegment_well_) { // if we use MultisegmentWell model
+        if (param_.use_multisegment_well_&& anyMSWellOpen(wells(), timeStepIdx)) { // if we use MultisegmentWell model
-            for (const auto& well : wells_ecl_) {
+            well_state_.initWellStateMSWell(wells(), wells_ecl_, timeStepIdx, phase_usage_, &previous_well_state_);
-                // TODO: this is acutally not very accurate, because sometimes a deck just claims a MS well
-                // while keep the well shut. More accurately, we should check if the well exisits in the Wells
-                // structure here
-                if (well->isMultiSegment(timeStepIdx) ) { // there is one well is MS well
-                    well_state_.initWellStateMSWell(wells(), wells_ecl_, timeStepIdx, phase_usage_, previous_well_state_);
-                    break;
-                }
-            }
         }
 
         // update the previous well state. This is used to restart failed steps.
@@ -480,13 +472,16 @@ namespace Opm {
     initFromRestartFile(const RestartValue& restartValues)
     {
         const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames();
+
+        // The restart step value is used to identify wells present at the given
+        // time step. Wells that are added at the same time step as RESTART is initiated
+        // will not be present in a restart file. Use the previous time step to retrieve
+        // wells that have information written to the restart file.
+        const int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0);
+
         WellsManager wellsmanager(eclState(),
                                   schedule(),
-                                  // The restart step value is used to identify wells present at the given
+                                  report_step,
-                                  // time step. Wells that are added at the same time step as RESTART is initiated
-                                  // will not be present in a restart file. Use the previous time step to retrieve
-                                  // wells that have information written to the restart file.
-                                  std::max(eclState().getInitConfig().getRestartStep() - 1, 0),
                                   Opm::UgGridHelpers::numCells(grid()),
                                   Opm::UgGridHelpers::globalCell(grid()),
                                   Opm::UgGridHelpers::cartDims(grid()),
@@ -502,8 +497,9 @@ namespace Opm {
         if (nw > 0) {
             const auto phaseUsage = phaseUsageFromDeck(eclState());
             const size_t numCells = Opm::UgGridHelpers::numCells(grid());
-            well_state_.resize(wells, schedule(), numCells, phaseUsage); // Resize for restart step
+            const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpen(wells, report_step));
-            wellsToState(restartValues.wells, phaseUsage, well_state_);
+            well_state_.resize(wells, wells_ecl_, schedule(), handle_ms_well, report_step, numCells, phaseUsage); // Resize for restart step
+            wellsToState(restartValues.wells, phaseUsage, handle_ms_well, report_step, well_state_);
             previous_well_state_ = well_state_;
         }
         initial_step_ = false;
@@ -1667,8 +1663,11 @@ namespace Opm {
     void
     BlackoilWellModel<TypeTag>::
     wellsToState( const data::Wells& wells,
-                       PhaseUsage phases,
+                  const PhaseUsage& phases,
-                       WellStateFullyImplicitBlackoil& state ) {
+                  const bool handle_ms_well,
+                  const int report_step,
+                  WellStateFullyImplicitBlackoil& state ) const
+    {
 
         using rt = data::Rates::opt;
         const auto np = phases.num_phases;
@@ -1724,6 +1723,85 @@ namespace Opm {
                 }
                 ++local_comp_index;
             }
+
+            if (handle_ms_well && !well.segments.empty()) {
+                // we need the well_ecl_ information
+                const std::string& well_name = wm.first;
+                const Well* well_ecl = getWellEcl(well_name);
+                assert(well_ecl);
+
+                const WellSegments& segment_set = well_ecl->getWellSegments(report_step);
+
+                const int top_segment_index = state.topSegmentIndex(well_index);
+                const auto& segments = well.segments;
+
+                // \Note: eventually we need to hanlde the situations that some segments are shut
+                assert(segment_set.size() == segments.size());
+
+                for (const auto& segment : segments) {
+                    const int segment_index = segment_set.segmentNumberToIndex(segment.first);
+
+                    // recovering segment rates and pressure from the restart values
+                    state.segPress()[top_segment_index + segment_index] = segment.second.pressure;
+
+                    const auto& segment_rates = segment.second.rates;
+                    for (int p = 0; p < np; ++p) {
+                        state.segRates()[(top_segment_index + segment_index) * np + p] = segment_rates.get(phs[p]);
+                    }
+                }
+            }
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    bool
+    BlackoilWellModel<TypeTag>::
+    anyMSWellOpen(const Wells* wells, const int report_step) const
+    {
+        bool any_ms_well_open = false;
+
+        const int nw = wells->number_of_wells;
+        for (int w = 0; w < nw; ++w) {
+            const std::string well_name = std::string(wells->name[w]);
+
+            const Well* well_ecl = getWellEcl(well_name);
+
+            assert(well_ecl);
+
+            if (well_ecl->isMultiSegment(report_step) ) {
+                any_ms_well_open = true;
+                break;
+            }
+        }
+        return any_ms_well_open;
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    const Well*
+    BlackoilWellModel<TypeTag>::
+    getWellEcl(const std::string& well_name) const
+    {
+        // finding the location of the well in wells_ecl
+        const int nw_wells_ecl = wells_ecl_.size();
+        int index_well = 0;
+        for (; index_well < nw_wells_ecl; ++index_well) {
+            if (well_name == wells_ecl_[index_well]->name()) {
+                break;
+            }
+        }
+
+        if (index_well < nw_wells_ecl) {
+            return wells_ecl_[index_well];
+        } else {
+            return nullptr;
         }
     }
 

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -281,23 +281,28 @@ namespace Opm
                 // we need to create a trival segment related values to avoid there will be some
                 // multi-segment wells added later.
                 nseg_ = nw;
-                seg_number_.clear();
+                top_segment_index_.resize(nw);
-                top_segment_index_.reserve(nw);
+                seg_number_.resize(nw);
-                seg_number_.reserve(nw);
                 for (int w = 0; w < nw; ++w) {
-                    top_segment_index_.push_back(w);
+                    top_segment_index_[w] = w;
-                    seg_number_.push_back(1); // Top segment is segment #1
+                    seg_number_[w] = 1; // Top segment is segment #1
                 }
                 segpress_ = bhp();
                 segrates_ = wellRates();
             }
         }
 
-        void resize(const Wells* wells, const Schedule& schedule, std::size_t numCells, const PhaseUsage& pu)
+
+        void resize(const Wells* wells, const std::vector<const Well*>& wells_ecl, const Schedule& schedule,
+                    const bool handle_ms_well, const int report_step, const size_t numCells,
+                    const PhaseUsage& pu)
         {
             const std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
-            const std::vector<const Well*> wells_ecl;
             init(wells, tmp, schedule, wells_ecl, 0, nullptr, pu);
+
+            if (handle_ms_well) {
+                initWellStateMSWell(wells, wells_ecl, report_step, pu, nullptr);
+            }
         }
 
         /// Allocate and initialize if wells is non-null.  Also tries
@@ -601,9 +606,8 @@ namespace Opm
 
 
         /// init the MS well related.
-        template <typename PrevWellState>
         void initWellStateMSWell(const Wells* wells, const std::vector<const Well*>& wells_ecl,
-                                 const int time_step, const PhaseUsage& pu, const PrevWellState& prev_well_state)
+                                 const int time_step, const PhaseUsage& pu, const WellStateFullyImplicitBlackoil* prev_well_state)
         {
             // still using the order in wells
             const int nw = wells->number_of_wells;
@@ -733,13 +737,13 @@ namespace Opm
             assert(int(segpress_.size()) == nseg_);
             assert(int(segrates_.size()) == nseg_ * numPhases() );
 
-            if (!prev_well_state.wellMap().empty()) {
+            if (prev_well_state && !prev_well_state->wellMap().empty()) {
                 // copying MS well related
-                const auto& end = prev_well_state.wellMap().end();
+                const auto& end = prev_well_state->wellMap().end();
                 const int np = numPhases();
                 for (int w = 0; w < nw; ++w) {
                     const std::string name( wells->name[w] );
-                    const auto& it = prev_well_state.wellMap().find( name );
+                    const auto& it = prev_well_state->wellMap().find( name );
 
                     if (it != end) { // the well is found in the prev_well_state
                         // TODO: the well with same name can change a lot, like they might not have same number of segments
@@ -747,7 +751,7 @@ namespace Opm
                         // for now, we just copy them.
                         const int old_index_well = (*it).second[0];
                         const int new_index_well = w;
-                        const int old_top_segment_index = prev_well_state.topSegmentIndex(old_index_well);
+                        const int old_top_segment_index = prev_well_state->topSegmentIndex(old_index_well);
                         const int new_top_segmnet_index = topSegmentIndex(new_index_well);
                         int number_of_segment = 0;
                         // if it is the last well in list
@@ -758,11 +762,11 @@ namespace Opm
                         }
 
                         for (int i = 0; i < number_of_segment * np; ++i) {
-                            segrates_[new_top_segmnet_index * np + i] = prev_well_state.segRates()[old_top_segment_index * np + i];
+                            segrates_[new_top_segmnet_index * np + i] = prev_well_state->segRates()[old_top_segment_index * np + i];
                         }
 
                         for (int i = 0; i < number_of_segment; ++i) {
-                            segpress_[new_top_segmnet_index + i] = prev_well_state.segPress()[old_top_segment_index + i];
+                            segpress_[new_top_segmnet_index + i] = prev_well_state->segPress()[old_top_segment_index + i];
                         }
                     }
                 }

tests/test_wellstatefullyimplicitblackoil.cpp

@@ -64,7 +64,6 @@ namespace {
     buildWellState(const Setup& setup, const std::size_t timeStep)
     {
         auto state  = Opm::WellStateFullyImplicitBlackoil{};
-        auto state0 = Opm::WellStateFullyImplicitBlackoil{}; // Empty.
 
         const auto cpress =
             std::vector<double>(setup.grid.c_grid()->number_of_cells,
@@ -74,11 +73,13 @@ namespace {
             setup.es, setup.sched, timeStep, *setup.grid.c_grid()
         };
 
-        state.init(wmgr.c_wells(), cpress, setup.sched.getWells(timeStep),
+        state.init(wmgr.c_wells(), cpress, setup.sched,
+                   setup.sched.getWells(timeStep),
                    timeStep, nullptr, setup.pu);
 
-        state.initWellStateMSWell(wmgr.c_wells(), setup.sched.getWells(timeStep),
+        state.initWellStateMSWell(wmgr.c_wells(),
-                                  timeStep, setup.pu, state0);
+                                  setup.sched.getWells(timeStep),
+                                  timeStep, setup.pu, nullptr);
 
         return state;
     }