Merge pull request #382 from blattms/fix-parallel-communication-whole-system

Fix parallel communication  for whole system

opm/autodiff/CPRPreconditioner.hpp

@@ -3,6 +3,7 @@
   Copyright 2014 IRIS AS.
   Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
   Copyright 2015 NTNU
+  Copyright 2015 Statoil AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -371,28 +372,31 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
                                     parallel run
         */
         CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae,
-                           const ParallelInformation& comm=ParallelInformation())
+                           const ParallelInformation& comm=ParallelInformation(),
+                           const ParallelInformation& commAe=ParallelInformation())
             : param_( param ),
               A_(A),
               Ae_(Ae),
               de_( Ae_.N() ),
               ve_( Ae_.M() ),
               dmodified_( A_.N() ),
-              opAe_(CPRSelector<M,X,Y,P>::makeOperator(Ae_, comm)),
+              opAe_(CPRSelector<M,X,Y,P>::makeOperator(Ae_, commAe)),
               precond_(), // ilu0 preconditioner for elliptic system
               amg_(),     // amg  preconditioner for elliptic system
               pre_(), // copy A will be made be the preconditioner
               vilu_( A_.N() ),
-              comm_(comm)
+              comm_(comm),
+              commAe_(commAe)
         {
             // create appropriate preconditioner for elliptic system
-            createPreconditioner( param_.cpr_use_amg_, comm );
+            createEllipticPreconditioner( param_.cpr_use_amg_, commAe_ );
 
+            // create the preconditioner for the whole system.
             if( param_.cpr_ilu_n_ == 0 ) {
-                pre_ = createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm );
+                pre_ = createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm_ );
             }
             else {
-                pre_ = createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm );
+                pre_ = createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ );
             }
         }
 
@@ -430,6 +434,8 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
             // dmodified = d - A * vfull
             dmodified_ = d;
             A_.mmv(v, dmodified_);
+            // A is not parallel, do communication manually.
+            comm_.copyOwnerToAll(dmodified_, dmodified_);
 
             // Apply Preconditioner for whole system (relax will be applied already)
             pre_->apply( vilu_, dmodified_);
@@ -470,7 +476,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
                 ScalarProductChooser;
             // the scalar product.
             std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
-                sp(ScalarProductChooser::construct(comm_));
+                sp(ScalarProductChooser::construct(commAe_));
 
             if( amg_ )
             {
@@ -534,10 +540,13 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
         //! \brief temporary variables for ILU solve
         Y vilu_;
 
-        //! \brief The information about the parallelization
+        //! \brief The information about the parallelization of the whole system.
         const P& comm_;
+        //! \brief The information about the parallelization of the elliptic part
+        //! of the system
+        const P& commAe_;
      protected:
-        void createPreconditioner( const bool amg, const P& comm )
+        void createEllipticPreconditioner( const bool amg, const P& comm )
         {
             if( amg )
             {

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -1,5 +1,8 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2014, 2015 Statoil AS
+  Copyright 2015 NTNU
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -389,7 +392,9 @@ namespace Opm {
         ///                   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 equations for each phase.
         /// \param[in]  nc    The number of cells of the local grid.
+        /// \param[in]  nw    The number of wells on the local grid.
         /// \return The total pore volume over all cells.
         double
         convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B,
@@ -398,7 +403,9 @@ namespace Opm {
                              std::array<double,MaxNumPhases>& R_sum,
                              std::array<double,MaxNumPhases>& maxCoeff,
                              std::array<double,MaxNumPhases>& B_avg,
-                             int nc) const;
+                             std::vector<double>& maxNormWell,
+                             int nc,
+                             int nw) const;
 
         void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
                                       const int it, const double relaxRelTol,

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -1,6 +1,7 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
-  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2014, 2015 Statoil AS
   Copyright 2015 NTNU
   Copyright 2015 IRIS AS
 
@@ -1254,8 +1255,24 @@ namespace detail {
 
     namespace detail
     {
-
+        /// \brief Compute the L-infinity norm of a vector
-        double infinityNorm( const ADB& a )
+        /// \warn This function is not suitable to compute on the well equations.
+        /// \param a The container to compute the infinity norm on.
+        ///          It has to have one entry for each cell.
+        /// \param info In a parallel this holds the information about the data distribution.
+        double infinityNorm( const ADB& a, const boost::any& pinfo=boost::any() )
+        {
+#if HAVE_MPI
+            if ( pinfo.type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& real_info =
+                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
+                double result=0;
+                real_info.computeReduction(a.value(), Reduction::makeGlobalMaxFunctor<double>(), result);
+                return result;
+            }
+            else
+#endif
             {
                 if( a.value().size() > 0 ) {
                     return a.value().matrix().lpNorm<Eigen::Infinity> ();
@@ -1264,6 +1281,27 @@ namespace detail {
                     return 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.
+        /// \param info In a parallel this holds the information about the data distribution.
+        double infinityNormWell( const ADB& a, const boost::any& pinfo )
+        {
+            double result=0;
+            if( a.value().size() > 0 ) {
+                result = a.value().matrix().lpNorm<Eigen::Infinity> ();
+            }
+#if HAVE_MPI
+            if ( pinfo.type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& real_info =
+                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
+                result = real_info.communicator().max(result);
+            }
+#endif
+            return result;
+        }
 
     } // namespace detail
 
@@ -1732,7 +1770,8 @@ namespace detail {
         const std::vector<ADB>::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
 
         for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
-            const double massBalanceResid = detail::infinityNorm( (*massBalanceIt) );
+            const double massBalanceResid = detail::infinityNorm( (*massBalanceIt),
+                                                                  linsolver_.parallelInformation() );
             if (!std::isfinite(massBalanceResid)) {
                 OPM_THROW(Opm::NumericalProblem,
                           "Encountered a non-finite residual");
@@ -1741,14 +1780,16 @@ namespace detail {
         }
 
         // the following residuals are not used in the oscillation detection now
-        const double wellFluxResid = detail::infinityNorm( residual_.well_flux_eq );
+        const double wellFluxResid = detail::infinityNormWell( residual_.well_flux_eq,
+                                                               linsolver_.parallelInformation() );
         if (!std::isfinite(wellFluxResid)) {
             OPM_THROW(Opm::NumericalProblem,
                "Encountered a non-finite residual");
         }
         residualNorms.push_back(wellFluxResid);
 
-        const double wellResid = detail::infinityNorm( residual_.well_eq );
+        const double wellResid = detail::infinityNormWell( residual_.well_eq,
+                                                           linsolver_.parallelInformation() );
         if (!std::isfinite(wellResid)) {
            OPM_THROW(Opm::NumericalProblem,
                "Encountered a non-finite residual");
@@ -1834,7 +1875,9 @@ namespace detail {
                                                          std::array<double,MaxNumPhases>& R_sum,
                                                          std::array<double,MaxNumPhases>& maxCoeff,
                                                          std::array<double,MaxNumPhases>& B_avg,
-                                                         int nc) const
+                                                         std::vector<double>& maxNormWell,
+                                                         int nc,
+                                                         int nw) const
     {
         // Do the global reductions
 #if HAVE_MPI
@@ -1865,12 +1908,18 @@ namespace detail {
                     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);
+                    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]));
+                    }
                 }
                 else
                 {
-                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0;
+                    maxNormWell[idx] = R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0;
                 }
             }
+            info.communicator().max(&maxNormWell[0], MaxNumPhases);
             // Compute pore volume
             return std::get<1>(nc_and_pv);
         }
@@ -1888,6 +1937,11 @@ namespace detail {
                 {
                     R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0;
                 }
+                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 geo_.poreVolume().sum();
@@ -1920,6 +1974,7 @@ namespace detail {
         Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B(nc, cols);
         Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R(nc, cols);
         Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV(nc, cols);
+        std::vector<double> maxNormWell(MaxNumPhases);
 
         for ( int idx=0; idx<MaxNumPhases; ++idx )
         {
@@ -1932,7 +1987,8 @@ namespace detail {
             }
         }
 
-        const double pvSum = convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, nc);
+        const double pvSum = convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg,
+                                                  maxNormWell, nc, nw);
 
         bool converged_MB = true;
         bool converged_CNV = true;
@@ -1944,18 +2000,13 @@ namespace detail {
             mass_balance_residual[idx]  = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
             converged_MB                = converged_MB && (mass_balance_residual[idx] < tol_mb);
             converged_CNV               = converged_CNV && (CNV[idx] < tol_cnv);
-
+            well_flux_residual[idx] = B_avg[idx] * dt * maxNormWell[idx];
-            double maxNormWell = 0.0;
-            for ( int w=0; w<nw; ++w )
-            {
-                maxNormWell  = std::max(maxNormWell, std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
-            }
-            well_flux_residual[idx] = B_avg[idx] * dt * maxNormWell;
 
             converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
         }
 
-        const double residualWell     = detail::infinityNorm(residual_.well_eq);
+        const double residualWell     = detail::infinityNormWell(residual_.well_eq,
+                                                                 linsolver_.parallelInformation());
         converged_Well   = converged_Well && (residualWell < Opm::unit::barsa);
         const bool   converged        = converged_MB && converged_CNV && converged_Well;
 

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -2,6 +2,7 @@
   Copyright 2014 SINTEF ICT, Applied Mathematics.
   Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
   Copyright 2015 NTNU
+  Copyright 2015 Statoil AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -188,11 +189,14 @@ namespace Opm
             const ParallelISTLInformation& info =
                 boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
             Comm istlComm(info.communicator());
-            info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices());
+            Comm istlAeComm(info.communicator());
+            info.copyValuesTo(istlAeComm.indexSet(), istlAeComm.remoteIndices());
+            info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(),
+                              istlAe.N(), istlA.N()/istlAe.N());
             // Construct operator, scalar product and vectors needed.
             typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
             Operator opA(istlA, istlComm);
-            constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, istlAe, x, istlb, istlComm, result);
+            constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, istlAe, x, istlb, istlComm, istlAeComm, result);
         }
         else
 #endif
@@ -201,7 +205,7 @@ namespace Opm
             typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
             Operator opA(istlA);
             Dune::Amg::SequentialInformation info;
-            constructPreconditionerAndSolve(opA, istlAe, x, istlb, info, result);
+            constructPreconditionerAndSolve(opA, istlAe, x, istlb, info, info, result);
         }
 
         // store number of iterations

opm/autodiff/NewtonIterationBlackoilCPR.hpp

@@ -1,6 +1,9 @@
 /*
   Copyright 2014 SINTEF ICT, Applied Mathematics.
   Copyright 2015 IRIS AS
+  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+  Copyright 2015 Statoil AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -85,6 +88,7 @@ namespace Opm
         void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe,
                                              Vector& x, Vector& istlb,
                                              const P& parallelInformation,
+                                             const P& parallelInformationAe,
                                              Dune::InverseOperatorResult& result) const
         {
             typedef Dune::ScalarProductChooser<Vector,P,category> ScalarProductChooser;
@@ -94,7 +98,8 @@ namespace Opm
             // typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
            typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
             parallelInformation.copyOwnerToAll(istlb, istlb);
-            Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation);
+            Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation,
+                                   parallelInformationAe);
 
             // TODO: Revise when linear solvers interface opm-core is done
             // Construct linear solver.

opm/autodiff/RedistributeDataHandles.hpp

@@ -1,6 +1,7 @@
 /*
   Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services.
   Coypright 2015 NTNU
+  Copyright 2015 Statoil AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -71,7 +72,7 @@ public:
         }
     }
     template<class B, class T>
-    void scatter(B& buffer, const T& e, std::size_t size)
+    void scatter(B& buffer, const T& e, std::size_t /* size */)
     {
         assert( T::codimension == 0);
         double val;

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -120,6 +120,8 @@ namespace Opm
         RateConverterType rateConverter_;
         // Threshold pressures.
         std::vector<double> threshold_pressures_by_face_;
+        // Whether this a parallel simulation or not
+        bool is_parallel_run_;
 
         void
         computeRESV(const std::size_t               step,
@@ -205,6 +207,7 @@ namespace Opm
                     boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
                 // Only rank 0 does print to std::cout
                 terminal_output_ = ( info.communicator().rank() == 0 );
+                is_parallel_run_ = ( info.communicator().size() > 1 );
             }
         }
 #endif
@@ -269,7 +272,8 @@ namespace Opm
                                        Opm::UgGridHelpers::dimensions(grid_),
                                        Opm::UgGridHelpers::cell2Faces(grid_),
                                        Opm::UgGridHelpers::beginFaceCentroids(grid_),
-                                       props_.permeability());
+                                       props_.permeability(),
+                                       is_parallel_run_);
             const Wells* wells = wells_manager.c_wells();
             WellStateFullyImplicitBlackoil well_state;
             well_state.init(wells, state, prev_well_state);