moving updateWellControls to MultisegmentWells

opm/autodiff/BlackoilMultiSegmentModel.hpp

@@ -168,10 +168,6 @@ namespace Opm {
 
         const MultisegmentWells::MultisegmentWellOps& msWellOps() const { return msWells().wellOps(); }
 
-
-        void updateWellControls(WellState& xw) const;
-
-
         // TODO: kept for now. to be removed soon.
         void updateWellState(const V& dwells,
                              WellState& well_state);

opm/autodiff/BlackoilMultiSegmentModel_impl.hpp

@@ -444,7 +444,7 @@ namespace Opm {
 
         // Possibly switch well controls and updating well state to
         // get reasonable initial conditions for the wells
-        asImpl().updateWellControls(well_state);
+        msWells().updateWellControls(terminal_output_, well_state);
 
         // Create the primary variables.
         SolutionState state = asImpl().variableState(reservoir_state, well_state);
@@ -567,97 +567,6 @@ namespace Opm {
 
 
 
-    template <class Grid>
-    void BlackoilMultiSegmentModel<Grid>::updateWellControls(WellState& xw) const
-    {
-        if( ! wellsActive() ) return ;
-
-        std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
-        // Find, for each well, if any constraints are broken. If so,
-        // switch control to first broken constraint.
-        const int np = wellsMultiSegment()[0]->numberOfPhases();
-        const int nw = wellsMultiSegment().size();
-        for (int w = 0; w < nw; ++w) {
-            const WellControls* wc = wellsMultiSegment()[w]->wellControls();
-            // The current control in the well state overrides
-            // the current control set in the Wells struct, which
-            // is instead treated as a default.
-            int current = xw.currentControls()[w];
-            // Loop over all controls except the current one, and also
-            // skip any RESERVOIR_RATE controls, since we cannot
-            // handle those.
-            const int nwc = well_controls_get_num(wc);
-            int ctrl_index = 0;
-            for (; ctrl_index < nwc; ++ctrl_index) {
-                if (ctrl_index == current) {
-                    // This is the currently used control, so it is
-                    // used as an equation. So this is not used as an
-                    // inequality constraint, and therefore skipped.
-                    continue;
-                }
-                if (wellhelpers::constraintBroken(
-                        xw.bhp(), xw.thp(), xw.wellRates(),
-                        w, np, wellsMultiSegment()[w]->wellType(), wc, ctrl_index)) {
-                    // ctrl_index will be the index of the broken constraint after the loop.
-                    break;
-                }
-            }
-
-            if (ctrl_index != nwc) {
-                // Constraint number ctrl_index was broken, switch to it.
-                if (terminal_output_)
-                {
-                    std::cout << "Switching control mode for well " << wellsMultiSegment()[w]->name()
-                              << " from " << modestring[well_controls_iget_type(wc, current)]
-                              << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
-                }
-                xw.currentControls()[w] = ctrl_index;
-                current = xw.currentControls()[w];
-            }
-
-            // Get gravity for THP hydrostatic corrrection
-            // const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
-
-            // Updating well state and primary variables.
-            // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
-            const double target = well_controls_iget_target(wc, current);
-            const double* distr = well_controls_iget_distr(wc, current);
-            switch (well_controls_iget_type(wc, current)) {
-            case BHP:
-                xw.bhp()[w] = target;
-                xw.segPress()[xw.topSegmentLoc()[w]] = target;
-                break;
-
-            case THP: {
-                OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
-            }
-
-            case RESERVOIR_RATE:
-                // No direct change to any observable quantity at
-                // surface condition.  In this case, use existing
-                // flow rates as initial conditions as reservoir
-                // rate acts only in aggregate.
-                break;
-
-            case SURFACE_RATE:
-                for (int phase = 0; phase < np; ++phase) {
-                    if (distr[phase] > 0.0) {
-                        xw.wellRates()[np * w + phase] = target * distr[phase];
-                        // TODO: consider changing all (not just top) segment rates
-                        // to make them consistent, it could possibly improve convergence.
-                        xw.segPhaseRates()[np * xw.topSegmentLoc()[w] + phase] = target * distr[phase];
-                    }
-                }
-                break;
-            }
-
-        }
-    }
-
-
-
-
-
     template <class Grid>
     bool BlackoilMultiSegmentModel<Grid>::solveWellEq(const std::vector<ADB>& mob_perfcells,
                                                       const std::vector<ADB>& b_perfcells,
@@ -787,7 +696,7 @@ namespace Opm {
                 const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
                 assert(dx.size() == total_residual_v.size());
                 asImpl().updateWellState(dx.array(), well_state);
-                updateWellControls(well_state);
+                msWells().updateWellControls(terminal_output_, well_state);
             }
         } while (it < 15);
 

opm/autodiff/MultisegmentWells.hpp

@@ -36,6 +36,7 @@
 #include <opm/autodiff/BlackoilModelEnums.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
+#include <opm/autodiff/WellHelpers.hpp>
 
 #include <opm/autodiff/WellMultiSegment.hpp>
 
@@ -177,6 +178,11 @@ namespace Opm {
                              const std::vector<bool>& active,
                              LinearisedBlackoilResidual& residual);
 
+            template <class WellState>
+            void
+            updateWellControls(const bool terminal_output,
+                               WellState& xw) const;
+
     protected:
         // TODO: probably a wells_active_ will be required here.
         const std::vector<WellMultiSegmentConstPtr> wells_multisegment_;

opm/autodiff/MultisegmentWells_impl.hpp

@@ -738,5 +738,101 @@ namespace Opm
                             others_residual;
     }
 
+
+
+
+
+    template <class WellState>
+    void
+    MultisegmentWells::
+    updateWellControls(const bool terminal_output,
+                       WellState& xw) const
+    {
+        if( wells().empty() ) return ;
+
+        std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
+        // Find, for each well, if any constraints are broken. If so,
+        // switch control to first broken constraint.
+        const int np = wells()[0]->numberOfPhases();
+        const int nw = wells().size();
+        for (int w = 0; w < nw; ++w) {
+            const WellControls* wc = wells()[w]->wellControls();
+            // The current control in the well state overrides
+            // the current control set in the Wells struct, which
+            // is instead treated as a default.
+            int current = xw.currentControls()[w];
+            // Loop over all controls except the current one, and also
+            // skip any RESERVOIR_RATE controls, since we cannot
+            // handle those.
+            const int nwc = well_controls_get_num(wc);
+            int ctrl_index = 0;
+            for (; ctrl_index < nwc; ++ctrl_index) {
+                if (ctrl_index == current) {
+                    // This is the currently used control, so it is
+                    // used as an equation. So this is not used as an
+                    // inequality constraint, and therefore skipped.
+                    continue;
+                }
+                if (wellhelpers::constraintBroken(
+                        xw.bhp(), xw.thp(), xw.wellRates(),
+                        w, np, wells()[w]->wellType(), wc, ctrl_index)) {
+                    // ctrl_index will be the index of the broken constraint after the loop.
+                    break;
+                }
+            }
+
+            if (ctrl_index != nwc) {
+                // Constraint number ctrl_index was broken, switch to it.
+                if (terminal_output)
+                {
+                    std::cout << "Switching control mode for well " << wells()[w]->name()
+                              << " from " << modestring[well_controls_iget_type(wc, current)]
+                              << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
+                }
+                xw.currentControls()[w] = ctrl_index;
+                current = xw.currentControls()[w];
+            }
+
+            // Get gravity for THP hydrostatic corrrection
+            // const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
+
+            // Updating well state and primary variables.
+            // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
+            const double target = well_controls_iget_target(wc, current);
+            const double* distr = well_controls_iget_distr(wc, current);
+            switch (well_controls_iget_type(wc, current)) {
+            case BHP:
+                xw.bhp()[w] = target;
+                xw.segPress()[top_well_segments_[w]] = target;
+                break;
+
+            case THP: {
+                OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
+            }
+
+            case RESERVOIR_RATE:
+                // No direct change to any observable quantity at
+                // surface condition.  In this case, use existing
+                // flow rates as initial conditions as reservoir
+                // rate acts only in aggregate.
+                break;
+
+            case SURFACE_RATE:
+                for (int phase = 0; phase < np; ++phase) {
+                    if (distr[phase] > 0.0) {
+                        xw.wellRates()[np * w + phase] = target * distr[phase];
+                        // TODO: consider changing all (not just top) segment rates
+                        // to make them consistent, it could possibly improve convergence.
+                        xw.segPhaseRates()[np * xw.topSegmentLoc()[w] + phase] = target * distr[phase];
+                    }
+                }
+                break;
+            }
+
+        }
+
+    }
+
+
 }
 #endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED