Merge remote-tracking branch 'remotes/opm/master' into format-step-output

opm/autodiff/MultisegmentWells.cpp

@@ -141,7 +141,7 @@ namespace Opm {
 
     MultisegmentWells::
     MultisegmentWells(const Wells* wells_arg,
-                      const std::vector<WellConstPtr>& wells_ecl,
+                      const std::vector< const Well* >& wells_ecl,
                       const int time_step)
       : wells_multisegment_( createMSWellVector(wells_arg, wells_ecl, time_step) )
       , wops_ms_(wells_multisegment_)
@@ -179,7 +179,7 @@ namespace Opm {
 
     std::vector<WellMultiSegmentConstPtr>
     MultisegmentWells::createMSWellVector(const Wells* wells_arg,
-                                          const std::vector<WellConstPtr>& wells_ecl,
+                                          const std::vector< const Well* >& wells_ecl,
                                           const int time_step)
     {
         std::vector<WellMultiSegmentConstPtr> wells_multisegment;

opm/autodiff/MultisegmentWells.hpp

@@ -83,11 +83,11 @@ namespace Opm {
             // TODO: using a vector of WellMultiSegmentConstPtr for now
             // TODO: it should use const Wells or something else later.
             MultisegmentWells(const Wells* wells_arg,
-                              const std::vector<WellConstPtr>& wells_ecl,
+                              const std::vector< const Well* >& wells_ecl,
                               const int time_step);
 
             std::vector<WellMultiSegmentConstPtr> createMSWellVector(const Wells* wells_arg,
-                                                                     const std::vector<WellConstPtr>& wells_ecl,
+                                                                     const std::vector< const Well* >& wells_ecl,
                                                                      const int time_step);
 
             void init(const BlackoilPropsAdInterface* fluid_arg,

opm/autodiff/SimulatorBase_impl.hpp

@@ -264,14 +264,14 @@ namespace Opm
     }
 
     namespace SimFIBODetails {
-        typedef std::unordered_map<std::string, WellConstPtr> WellMap;
+        typedef std::unordered_map<std::string, const Well* > WellMap;
 
         inline WellMap
-        mapWells(const std::vector<WellConstPtr>& wells)
+        mapWells(const std::vector< const Well* >& wells)
         {
             WellMap wmap;
 
-            for (std::vector<WellConstPtr>::const_iterator
+            for (std::vector< const Well* >::const_iterator
                      w = wells.begin(), e = wells.end();
                  w != e; ++w)
             {
@@ -313,7 +313,7 @@ namespace Opm
             WellMap::const_iterator i = wmap.find(name);
 
             if (i != wmap.end()) {
-                WellConstPtr wp = i->second;
+                const Well* wp = i->second;
 
                 match = (wp->isProducer(step) &&
                          wp->getProductionProperties(step)
@@ -435,7 +435,7 @@ namespace Opm
     {
         typedef SimFIBODetails::WellMap WellMap;
 
-        const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
+        const auto w_ecl = eclipse_state_->getSchedule()->getWells(step);
         const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
 
         const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
@@ -513,7 +513,7 @@ namespace Opm
                     WellMap::const_iterator i = wmap.find(wells->name[*rp]);
 
                     if (i != wmap.end()) {
-                        WellConstPtr wp = i->second;
+                        const auto* wp = i->second;
 
                         const WellProductionProperties& p =
                             wp->getProductionProperties(step);
@@ -570,7 +570,7 @@ namespace Opm
                 if (!is_producer && wells->name[w] != 0) {
                     WellMap::const_iterator i = wmap.find(wells->name[w]);
                     if (i != wmap.end()) {
-                        WellConstPtr wp = i->second;
+                        const auto* wp = i->second;
                         const WellInjectionProperties& injector = wp->getInjectionProperties(step);
                         if (!injector.predictionMode) {
                             //History matching WCONINJEH

opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp

@@ -114,7 +114,7 @@ namespace Opm
             WellState well_state;
             // well_state.init(wells, state, prev_well_state);
 
-            const std::vector<WellConstPtr>& wells_ecl = eclipse_state_->getSchedule()->getWells(timer.currentStepNum());
+            const auto wells_ecl = eclipse_state_->getSchedule()->getWells(timer.currentStepNum());
             const int current_time_step = timer.currentStepNum();
 
             const WellModel well_model(wells, wells_ecl, current_time_step);

opm/autodiff/StandardWells.hpp

@@ -40,7 +40,7 @@ namespace Opm {
 
         /// Class for handling the standard well model.
         class StandardWells {
-        protected:
+        public:
             struct WellOps {
                 explicit WellOps(const Wells* wells);
                 Eigen::SparseMatrix<double> w2p;              // well -> perf (scatter)
@@ -48,7 +48,6 @@ namespace Opm {
                 std::vector<int> well_cells;                  // the set of perforated cells
             };
 
-        public:
             // ---------      Types      ---------
             using ADB = AutoDiffBlock<double>;
             using Vector = ADB::V;
@@ -60,7 +59,7 @@ namespace Opm {
                                             Eigen::Dynamic,
                                             Eigen::RowMajor>;
             // ---------  Public methods  ---------
-            StandardWells(const Wells* wells_arg);
+            explicit StandardWells(const Wells* wells_arg);
 
             void init(const BlackoilPropsAdInterface* fluid_arg,
                       const std::vector<bool>* active_arg,
@@ -161,6 +160,15 @@ namespace Opm {
             variableWellStateInitials(const WellState& xw,
                                       std::vector<Vector>& vars0) const;
 
+            /// If set, computeWellFlux() will additionally store the
+            /// total reservoir volume perforation fluxes.
+            void setStoreWellPerforationFluxesFlag(const bool store_fluxes);
+
+            /// Retrieves the stored fluxes. It is an error to call this
+            /// unless setStoreWellPerforationFluxesFlag(true) has been
+            /// called.
+            const Vector& getStoredWellPerforationFluxes() const;
+
         protected:
             bool wells_active_;
             const Wells*   wells_;
@@ -179,6 +187,9 @@ namespace Opm {
             Vector well_perforation_densities_;
             Vector well_perforation_pressure_diffs_;
 
+            bool store_well_perforation_fluxes_;
+            Vector well_perforation_fluxes_;
+
             // protected methods
             template <class SolutionState, class WellState>
             void computePropertiesForWellConnectionPressures(const SolutionState& state,

opm/autodiff/StandardWells_impl.hpp

@@ -81,6 +81,7 @@ namespace Opm
       , vfp_properties_(nullptr)
       , well_perforation_densities_(Vector())
       , well_perforation_pressure_diffs_(Vector())
+      , store_well_perforation_fluxes_(false)
     {
     }
 
@@ -441,10 +442,11 @@ namespace Opm
 
         // HANDLE FLOW INTO WELLBORE
         // compute phase volumetric rates at standard conditions
+        std::vector<ADB> cq_p(np, ADB::null());
         std::vector<ADB> cq_ps(np, ADB::null());
         for (int phase = 0; phase < np; ++phase) {
-            const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
+            cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
-            cq_ps[phase] = b_perfcells[phase] * cq_p;
+            cq_ps[phase] = b_perfcells[phase] * cq_p[phase];
         }
         const Opm::PhaseUsage& pu = fluid_->phaseUsage();
         if ((*active_)[Oil] && (*active_)[Gas]) {
@@ -467,6 +469,16 @@ namespace Opm
         // injection perforations total volume rates
         const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
 
+        // Store well perforation total fluxes (reservor volumes) if requested.
+        if (store_well_perforation_fluxes_) {
+            // Ugly const-cast, but unappealing alternatives.
+            Vector& wf = const_cast<Vector&>(well_perforation_fluxes_);
+            wf = cqt_i.value();
+            for (int phase = 0; phase < np; ++phase) {
+                wf += cq_p[phase].value();
+            }
+        }
+
         // compute wellbore mixture for injecting perforations
         // The wellbore mixture depends on the inflow from the reservoar
         // and the well injection rates.
@@ -1206,4 +1218,26 @@ namespace Opm
         }
     }
 
-}
+
+
+
+
+    void
+    StandardWells::setStoreWellPerforationFluxesFlag(const bool store_fluxes)
+    {
+        store_well_perforation_fluxes_ = store_fluxes;
+    }
+
+
+
+
+
+    const StandardWells::Vector&
+    StandardWells::getStoredWellPerforationFluxes() const
+    {
+        assert(store_well_perforation_fluxes_);
+        return well_perforation_fluxes_;
+    }
+
+
+} // namespace Opm

opm/autodiff/WellMultiSegment.cpp

@@ -26,7 +26,7 @@
 namespace Opm
 {
 
-    WellMultiSegment::WellMultiSegment(WellConstPtr well, size_t time_step, const Wells* wells) {
+    WellMultiSegment::WellMultiSegment(const Well* well, size_t time_step, const Wells* wells) {
         m_well_name_ = well->name();
         if (well->isMultiSegment(time_step)) {
             initMultiSegmentWell(well, time_step, wells);
@@ -36,7 +36,7 @@ namespace Opm
         updateWellOps();
     }
 
-    void WellMultiSegment::initMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells) {
+    void WellMultiSegment::initMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells) {
 
         CompletionSetConstPtr completion_set = well->getCompletions(time_step);
 
@@ -169,7 +169,7 @@ namespace Opm
 
     }
 
-    void WellMultiSegment::initNonMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells) {
+    void WellMultiSegment::initNonMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells) {
 
         CompletionSetConstPtr completion_set = well->getCompletions(time_step);
 

opm/autodiff/WellMultiSegment.hpp

@@ -54,7 +54,7 @@ namespace Opm
         /// \param[in]       well information from EclipseState
         /// \param[in]       current time step
         /// \param[in[       pointer to Wells structure, to be removed eventually
-        WellMultiSegment(WellConstPtr well, size_t time_step, const Wells* wells);
+        WellMultiSegment(const Well* well, size_t time_step, const Wells* wells);
 
         /// Well name.
         const std::string& name() const;
@@ -146,8 +146,8 @@ namespace Opm
     private:
         // for the moment, we use the information from wells.
         // TODO: remove the dependency on wells from opm-core.
-        void initMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells);
+        void initMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells);
-        void initNonMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells);
+        void initNonMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells);
         void updateWellOps();
 
     private:

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp

@@ -186,12 +186,12 @@ namespace Opm
                                    cartesian_to_compressed);
 
         ScheduleConstPtr schedule = eclipseState->getSchedule();
-        std::vector<WellConstPtr> wells    = schedule->getWells(timeStep);
+        auto wells       = schedule->getWells(timeStep);
 
         int well_index = 0;
 
         for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
-             WellConstPtr well = (*wellIter);
+             const auto* well = (*wellIter);
 
              if (well->getStatus(timeStep) == WellCommon::SHUT) {
                  continue;

tests/test_multisegmentwells.cpp

@@ -104,7 +104,7 @@ struct SetupMSW {
                                         false);
 
         const Wells* wells = wells_manager.c_wells();
-        const std::vector<Opm::WellConstPtr>& wells_ecl = ecl_state->getSchedule()->getWells(current_timestep);
+        const auto wells_ecl = ecl_state->getSchedule()->getWells(current_timestep);
 
         ms_wells.reset(new Opm::MultisegmentWells(wells, wells_ecl, current_timestep));
     };