Changes to well model for supporting sequential approach.

The changes are:
 - Make the WellOps struct public (needed by transport solver).
 - Make it possible to store and retrieve total reservoir volume
   perforation fluxes with getStoredWellPerforationFluxes(), controlled
   by a flag set by setStoreWellPerforationFluxesFlag(), defaulting to
   false (needed by pressure solver).

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_ps[phase] = b_perfcells[phase] * cq_p;
+            cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
+            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.
+            V& wf = const_cast<V&>(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