Merge pull request #2531 from GitPaean/msw_keywords

adding the summary output for several pressure drop values for MSW

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1664,7 +1664,8 @@ namespace Opm {
                     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 pres_idx = Opm::data::SegmentPressures::Value::Pressure;
+                    state.segPress()[top_segment_index + segment_index] = segment.second.pressures[pres_idx];
 
                     const auto& segment_rates = segment.second.rates;
                     for (int p = 0; p < np; ++p) {

opm/simulators/wells/MultisegmentWell.hpp

@@ -380,7 +380,7 @@ namespace Opm
                                const Well::ProductionControls& prod_controls,
                                Opm::DeferredLogger& deferred_logger);
 
-        void assemblePressureEq(const int seg) const;
+        void assemblePressureEq(const int seg, WellState& well_state) const;
 
         // hytrostatic pressure loss
         EvalWell getHydroPressureLoss(const int seg) const;
@@ -388,7 +388,7 @@ namespace Opm
         // frictinal pressure loss
         EvalWell getFrictionPressureLoss(const int seg) const;
 
-        void handleAccelerationPressureLoss(const int seg) const;
+        void handleAccelerationPressureLoss(const int seg, WellState& well_state) const;
 
         // handling the overshooting and undershooting of the fractions
         void processFractions(const int seg) const;
@@ -470,7 +470,7 @@ namespace Opm
 
         double maxPerfPress(const Simulator& ebos_simulator) const;
 
-        void assembleSICDPressureEq(const int seg) const;
+        void assembleSICDPressureEq(const int seg, WellState& well_state) const;
 
         // TODO: when more ICD devices join, we should have a better interface to do this
         void calculateSICDFlowScalingFactors();
@@ -478,7 +478,7 @@ namespace Opm
         EvalWell pressureDropSpiralICD(const int seg) const;
 
         // assemble the pressure equation for sub-critical valve (WSEGVALV)
-        void assembleValvePressureEq(const int seg) const;
+        void assembleValvePressureEq(const int seg, WellState& well_state) const;
 
         EvalWell pressureDropValve(const int seg) const;
     };

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1984,7 +1984,7 @@ namespace Opm
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    assemblePressureEq(const int seg) const
+    assemblePressureEq(const int seg, WellState& well_state) const
     {
         assert(seg != 0); // not top segment
 
@@ -1993,10 +1993,16 @@ namespace Opm
 
         // we need to handle the pressure difference between the two segments
         // we only consider the hydrostatic pressure loss first
-        pressure_equation -= getHydroPressureLoss(seg);
+        // TODO: we might be able to add member variables to store these values, then we update well state
+        // after converged
+        const auto hydro_pressure_drop = getHydroPressureLoss(seg);
+        well_state.segPressDropHydroStatic()[seg] = hydro_pressure_drop.value();
+        pressure_equation -= hydro_pressure_drop;
 
         if (frictionalPressureLossConsidered()) {
-            pressure_equation -= getFrictionPressureLoss(seg);
+            const auto friction_pressure_drop = getFrictionPressureLoss(seg);
+            pressure_equation -= friction_pressure_drop;
+            well_state.segPressDropFriction()[seg] = friction_pressure_drop.value();
         }
 
         resWell_[seg][SPres] = pressure_equation.value();
@@ -2014,7 +2020,7 @@ namespace Opm
         }
 
         if (accelerationalPressureLossConsidered()) {
-            handleAccelerationPressureLoss(seg);
+            handleAccelerationPressureLoss(seg, well_state);
         }
     }
 
@@ -2061,7 +2067,7 @@ namespace Opm
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    handleAccelerationPressureLoss(const int seg) const
+    handleAccelerationPressureLoss(const int seg, WellState& well_state) const
     {
         // TODO: this pressure loss is not significant enough to be well tested yet.
         // handle the out velcocity head
@@ -2089,6 +2095,7 @@ namespace Opm
             const EvalWell inlet_density = segment_densities_[inlet];
             const EvalWell inlet_mass_rate = segment_mass_rates_[inlet];
             const EvalWell inlet_velocity_head = mswellhelpers::velocityHead(area, inlet_mass_rate, inlet_density);
+            well_state.segPressDropAcceleration()[seg] = inlet_velocity_head.value();
             resWell_[seg][SPres] += inlet_velocity_head.value();
             for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
                 duneD_[seg][inlet][SPres][pv_idx] += inlet_velocity_head.derivative(pv_idx + numEq);
@@ -2537,15 +2544,19 @@ namespace Opm
                 // TODO: maybe the following should go to the function assemblePressureEq()
 		switch(segmentSet()[seg].segmentType()) {
 		    case Segment::SegmentType::SICD :
-		        assembleSICDPressureEq(seg);
+		        assembleSICDPressureEq(seg, well_state);
 			break;
 		    case Segment::SegmentType::VALVE :
-		        assembleValvePressureEq(seg);
+		        assembleValvePressureEq(seg, well_state);
 			break;
 		    default :
-		        assemblePressureEq(seg);
+		        assemblePressureEq(seg, well_state);
 		}
             }
+
+            well_state.segPressDrop()[seg] = well_state.segPressDropHydroStatic()[seg] +
+                                             well_state.segPressDropFriction()[seg] +
+                                             well_state.segPressDropAcceleration()[seg];
         }
     }
 
@@ -3064,7 +3075,7 @@ namespace Opm
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    assembleSICDPressureEq(const int seg) const
+    assembleSICDPressureEq(const int seg, WellState& well_state) const
     {
         // TODO: upwinding needs to be taken care of
         // top segment can not be a spiral ICD device
@@ -3076,7 +3087,9 @@ namespace Opm
 
         EvalWell pressure_equation = getSegmentPressure(seg);
 
-        pressure_equation = pressure_equation - pressureDropSpiralICD(seg);
+        const auto sicd_pressure_drop = pressureDropSpiralICD(seg);
+        pressure_equation = pressure_equation - sicd_pressure_drop;
+        well_state.segPressDropFriction()[seg] = sicd_pressure_drop.value();
 
         resWell_[seg][SPres] = pressure_equation.value();
         for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
@@ -3100,7 +3113,7 @@ namespace Opm
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    assembleValvePressureEq(const int seg) const
+    assembleValvePressureEq(const int seg, WellState& well_state) const
     {
         // TODO: upwinding needs to be taken care of
         // top segment can not be a spiral ICD device
@@ -3116,7 +3129,9 @@ namespace Opm
 
         // const int seg_upwind = upwinding_segments_[seg];
 
-        pressure_equation = pressure_equation - pressureDropValve(seg);
+        const auto valve_pressure_drop = pressureDropValve(seg);
+        pressure_equation = pressure_equation - valve_pressure_drop;
+        well_state.segPressDropFriction()[seg] = valve_pressure_drop.value();
 
         resWell_[seg][SPres] = pressure_equation.value();
         for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -279,8 +279,13 @@ namespace Opm
                     top_segment_index_[w] = w;
                     seg_number_[w] = 1; // Top segment is segment #1
                 }
-                segpress_ = bhp();
+                seg_press_ = bhp();
-                segrates_ = wellRates();
+                seg_rates_ = wellRates();
+
+                seg_pressdrop_.assign(nw, 0.);
+                seg_pressdrop_hydorstatic_.assign(nw, 0.);
+                seg_pressdrop_friction_.assign(nw, 0.);
+                seg_pressdrop_acceleration_.assign(nw, 0.);
             }
         }
 
@@ -584,10 +589,10 @@ namespace Opm
 
             top_segment_index_.clear();
             top_segment_index_.reserve(nw);
-            segpress_.clear();
+            seg_press_.clear();
-            segpress_.reserve(nw);
+            seg_press_.reserve(nw);
-            segrates_.clear();
+            seg_rates_.clear();
-            segrates_.reserve(nw * numPhases());
+            seg_rates_.reserve(nw * numPhases());
             seg_number_.clear();
 
             nseg_ = 0;
@@ -601,10 +606,10 @@ namespace Opm
                 if ( !well_ecl.isMultiSegment() ) { // not multi-segment well
                     nseg_ += 1;
                     seg_number_.push_back(1); // Assign single segment (top) as number 1.
-                    segpress_.push_back(bhp()[w]);
+                    seg_press_.push_back(bhp()[w]);
                     const int np = numPhases();
                     for (int p = 0; p < np; ++p) {
-                        segrates_.push_back(wellRates()[np * w + p]);
+                        seg_rates_.push_back(wellRates()[np * w + p]);
                     }
                 } else { // it is a multi-segment well
                     const WellSegments& segment_set = well_ecl.getSegments();
@@ -642,7 +647,7 @@ namespace Opm
                     }
 
 
-                    // for the segrates_, now it becomes a recursive solution procedure.
+                    // for the seg_rates_, now it becomes a recursive solution procedure.
                     {
                         const int np = numPhases();
                         const int start_perf = connpos;
@@ -668,7 +673,7 @@ namespace Opm
                                                                     perfPhaseRates().begin() + np * start_perf_next_well); // the perforation rates for this well
                         std::vector<double> segment_rates;
                         calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, segment_rates);
-                        std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(segrates_));
+                        std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(seg_rates_));
                     }
 
                     // for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment
@@ -678,26 +683,32 @@ namespace Opm
                     // improved during the solveWellEq process
                     {
                         // top segment is always the first one, and its pressure is the well bhp
-                        segpress_.push_back(bhp()[w]);
+                        seg_press_.push_back(bhp()[w]);
                         const int top_segment = top_segment_index_[w];
                         const int start_perf = connpos;
                         for (int seg = 1; seg < well_nseg; ++seg) {
                             if ( !segment_perforations[seg].empty() ) {
                                 const int first_perf = segment_perforations[seg][0];
-                                segpress_.push_back(perfPress()[start_perf + first_perf]);
+                                seg_press_.push_back(perfPress()[start_perf + first_perf]);
                             } else {
-                                // segpress_.push_back(bhp); // may not be a good decision
+                                // seg_press_.push_back(bhp); // may not be a good decision
                                 // using the outlet segment pressure // it needs the ordering is correct
                                 const int outlet_seg = segment_set[seg].outletSegment();
-                                segpress_.push_back(segpress_[top_segment + segment_set.segmentNumberToIndex(outlet_seg)]);
+                                seg_press_.push_back(
+                                    seg_press_[top_segment + segment_set.segmentNumberToIndex(outlet_seg)]);
                             }
                         }
                     }
                 }
                 connpos += num_perf_this_well;
             }
-            assert(int(segpress_.size()) == nseg_);
+            assert(int(seg_press_.size()) == nseg_);
-            assert(int(segrates_.size()) == nseg_ * numPhases() );
+            assert(int(seg_rates_.size()) == nseg_ * numPhases() );
+
+            seg_pressdrop_.assign(nseg_, 0.);
+            seg_pressdrop_hydorstatic_.assign(nseg_, 0.);
+            seg_pressdrop_friction_.assign(nseg_, 0.);
+            seg_pressdrop_acceleration_.assign(nseg_, 0.);
 
             if (prev_well_state && !prev_well_state->wellMap().empty()) {
                 // copying MS well related
@@ -723,11 +734,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];
+                            seg_rates_[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];
+                            seg_press_[new_top_segmnet_index + i] = prev_well_state->segPress()[old_top_segment_index + i];
                         }
                     }
                 }
@@ -810,22 +821,62 @@ namespace Opm
 
         const std::vector<double>& segRates() const
         {
-            return segrates_;
+            return seg_rates_;
         }
 
         std::vector<double>& segRates()
         {
-            return segrates_;
+            return seg_rates_;
         }
 
         const std::vector<double>& segPress() const
         {
-            return segpress_;
+            return seg_press_;
+        }
+
+        std::vector<double>& segPressDrop()
+        {
+            return seg_pressdrop_;
+        }
+
+        const std::vector<double>& segPressDrop() const
+        {
+            return seg_pressdrop_;
+        }
+
+        std::vector<double>& segPressDropFriction()
+        {
+            return seg_pressdrop_friction_;
+        }
+
+        const std::vector<double>& segPressDropFriction() const
+        {
+            return seg_pressdrop_friction_;
+        }
+
+        std::vector<double>& segPressDropHydroStatic()
+        {
+            return seg_pressdrop_hydorstatic_;
+        }
+
+        const std::vector<double>& segPressDropHydroStatic() const
+        {
+            return seg_pressdrop_hydorstatic_;
+        }
+
+        std::vector<double>& segPressDropAcceleration()
+        {
+            return seg_pressdrop_acceleration_;
+        }
+
+        const std::vector<double>& segPressDropAcceleration() const
+        {
+            return seg_pressdrop_acceleration_;
         }
 
         std::vector<double>& segPress()
         {
-            return segpress_;
+            return seg_press_;
         }
 
         int numSegment() const
@@ -1026,8 +1077,17 @@ namespace Opm
 
         // MS well related
         // for StandardWell, the number of segments will be one
-        std::vector<double> segrates_;
+        std::vector<double> seg_rates_;
-        std::vector<double> segpress_;
+        std::vector<double> seg_press_;
+        // The following data are only recorded for output
+        // pressure drop
+        std::vector<double> seg_pressdrop_;
+        // frictional pressure drop
+        std::vector<double> seg_pressdrop_friction_;
+        // hydrostatic pressure drop
+        std::vector<double> seg_pressdrop_hydorstatic_;
+        // accelerational pressure drop
+        std::vector<double> seg_pressdrop_acceleration_;
         // the index of the top segments, which is used to locate the
         // multisegment well related information in WellState
         std::vector<int> top_segment_index_;
@@ -1062,7 +1122,15 @@ namespace Opm
             const auto* rate =
                 &this->segRates()[seg_dof * this->numPhases()];
 
-            seg_res.pressure = this->segPress()[seg_dof];
+            {
+                using Value =::Opm::data::SegmentPressures::Value;
+                auto& segpress = seg_res.pressures;
+                segpress[Value::Pressure] = this->segPress()[seg_dof];
+                segpress[Value::PDrop] = this->segPressDrop()[seg_dof];
+                segpress[Value::PDropHydrostatic] = this->segPressDropHydroStatic()[seg_dof];
+                segpress[Value::PDropFriction] = this->segPressDropFriction()[seg_dof];
+                segpress[Value::PDropAccel] = this->segPressDropAcceleration()[seg_dof];
+            }
 
             if (pu.phase_used[Water]) {
                 seg_res.rates.set(data::Rates::opt::wat,

tests/test_ParallelRestart.cpp

@@ -173,7 +173,8 @@ Opm::data::Segment getSegment()
     Opm::data::Segment seg1;
     seg1.rates = getRates();
     seg1.segNumber = 1;
-    seg1.pressure = 2.0;
+    const auto pres_idx = Opm::data::SegmentPressures::Value::Pressure;
+    seg1.pressures[pres_idx] = 2.0;
     return seg1;
 }
 

tests/test_wellstatefullyimplicitblackoil.cpp

@@ -261,7 +261,8 @@ BOOST_AUTO_TEST_CASE(Pressure)
         const auto& xseg = xw.segments.at(1);
 
         BOOST_CHECK_EQUAL(xseg.segNumber, 1);
-        BOOST_CHECK_CLOSE(xseg.pressure, prod01_first ? 100.0 : 0.0, 1.0e-10);
+        const auto pres_idx = Opm::data::SegmentPressures::Value::Pressure;
+        BOOST_CHECK_CLOSE(xseg.pressures[pres_idx], prod01_first ? 100.0 : 0.0, 1.0e-10);
     }
 
     {
@@ -276,7 +277,8 @@ BOOST_AUTO_TEST_CASE(Pressure)
             const auto& xseg = xw.segments.at(segID + 1);
 
             BOOST_CHECK_EQUAL(xseg.segNumber, segID + 1);
-            BOOST_CHECK_CLOSE(xseg.pressure, pressTop + 1.0*segID, 1.0e-10);
+            const auto pres_idx = Opm::data::SegmentPressures::Value::Pressure;
+            BOOST_CHECK_CLOSE(xseg.pressures[pres_idx], pressTop + 1.0*segID, 1.0e-10);
         }
     }
 }