Reorganize the wellModel

- All well-related stuff is moved to StandardWellsDense.hpp
- Some unused functions are removed from BlackoilModelEbos.hpp

opm/autodiff/BlackoilDetails.hpp

@@ -188,6 +188,108 @@ namespace detail {
                 return product;
             }
         }
+        /// \brief Compute the reduction within the convergence check.
+        /// \param[in] B     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. tempV.col(i) contains the
+        ///                   values
+        ///                  for phase i.
+        /// \param[in] R     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
+        ///                   of R for the phase i.
+        /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
+        ///                   maximum of tempV for the phase i.
+        /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
+        ///                   of B for the phase i.
+        /// \param[out] maxNormWell The maximum of the well flux equations for each phase.
+        /// \param[in]  nc    The number of cells of the local grid.
+        /// \return The total pore volume over all cells.
+        inline
+        double
+        convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
+                             const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
+                             const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
+                             std::vector<double>& R_sum,
+                             std::vector<double>& maxCoeff,
+                             std::vector<double>& B_avg,
+                             std::vector<double>& maxNormWell,
+                             int nc,
+                             int np,
+                             const std::vector<double> pv,
+                             LinearisedBlackoilResidual residual)
+        {
+            const int nw = residual.well_flux_eq.size() / np;
+            assert(nw * np == int(residual.well_flux_eq.size()));
+
+            // Do the global reductions
+#if HAVE_MPI
+            if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& info =
+                    boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
+
+                // Compute the global number of cells and porevolume
+                std::vector<int> v(nc, 1);
+                auto nc_and_pv = std::tuple<int, double>(0, 0.0);
+                auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
+                                                           Opm::Reduction::makeGlobalSumFunctor<double>());
+                auto nc_and_pv_containers  = std::make_tuple(v, pv);
+                info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
+
+                for ( int idx = 0; idx < np; ++idx )
+                {
+                    auto values     = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
+                    auto containers = std::make_tuple(B.col(idx),
+                                                      tempV.col(idx),
+                                                      R.col(idx));
+                    auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalMaxFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalSumFunctor<double>());
+                    info.computeReduction(containers, operators, values);
+                    B_avg[idx]       = std::get<0>(values)/std::get<0>(nc_and_pv);
+                    maxCoeff[idx]    = std::get<1>(values);
+                    R_sum[idx]       = std::get<2>(values);
+                    assert(np >= np);
+                    if (idx < np) {
+                        maxNormWell[idx] = 0.0;
+                        for ( int w = 0; w < nw; ++w ) {
+                            maxNormWell[idx]  = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
+                        }
+                    }
+                }
+                info.communicator().max(maxNormWell.data(), np);
+                // Compute pore volume
+                return std::get<1>(nc_and_pv);
+            }
+            else
+#endif
+            {
+                B_avg.resize(np);
+                maxCoeff.resize(np);
+                R_sum.resize(np);
+                maxNormWell.resize(np);
+                for ( int idx = 0; idx < np; ++idx )
+                {
+                    B_avg[idx] = B.col(idx).sum()/nc;
+                    maxCoeff[idx] = tempV.col(idx).maxCoeff();
+                    R_sum[idx] = R.col(idx).sum();
+
+                    assert(np >= np);
+                    if (idx < np) {
+                        maxNormWell[idx] = 0.0;
+                        for ( int w = 0; w < nw; ++w ) {
+                            maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual.well_flux_eq.value()[nw*idx + w]));
+                        }
+                    }
+                }
+                // Compute total pore volume
+                return std::accumulate(pv.begin(), pv.end(), 0.0);
+            }
+        }
 
         /// \brief Compute the L-infinity norm of a vector representing a well equation.
         /// \param a The container to compute the infinity norm on.

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -33,7 +33,7 @@
 namespace Opm {
 
 class SimulatorFullyImplicitBlackoilEbos;
-class StandardWellsDense;
+//class StandardWellsDense<FluidSystem>;
 
 /// a simulator for the blackoil model
 class SimulatorFullyImplicitBlackoilEbos
@@ -42,13 +42,16 @@ public:
     typedef typename TTAG(EclFlowProblem) TypeTag;
     typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
     typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
+
     typedef WellStateFullyImplicitBlackoil WellState;
     typedef BlackoilState ReservoirState;
     typedef BlackoilOutputWriter OutputWriter;
     typedef BlackoilModelEbos Model;
     typedef BlackoilModelParameters ModelParameters;
     typedef NonlinearSolver<Model> Solver;
-    typedef StandardWellsDense WellModel;
+    typedef StandardWellsDense<FluidSystem, BlackoilIndices> WellModel;
 
 
     /// Initialise from parameters and objects to observe.
@@ -229,7 +232,8 @@ public:
             // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
 
-            const WellModel well_model(wells);
+            const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
+            const WellModel well_model(wells, model_param_, terminal_output_, pv);
 
             auto solver = createSolver(well_model);