diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index 56cffe7e5..91a5962de 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/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::makeOperator(Ae_, comm)), + opAe_(CPRSelector::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( A_, param_.cpr_relax_, comm ); + pre_ = createILU0Ptr( A_, param_.cpr_relax_, comm_ ); } else { - pre_ = createILUnPtr( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm ); + pre_ = createILUnPtr( 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 - 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 ) { diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 5fd4dc046..f67e38b30 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/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& B, @@ -398,7 +403,9 @@ namespace Opm { std::array& R_sum, std::array& maxCoeff, std::array& B_avg, - int nc) const; + std::vector& maxNormWell, + int nc, + int nw) const; void detectNewtonOscillations(const std::vector>& residual_history, const int it, const double relaxRelTol, diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index d92b4aed0..5362d142a 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/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,15 +1255,52 @@ namespace detail { namespace detail { - - double infinityNorm( const ADB& a ) + /// \brief Compute the L-infinity norm of a vector + /// \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(pinfo); + double result=0; + real_info.computeReduction(a.value(), Reduction::makeGlobalMaxFunctor(), result); + return result; + } + else +#endif + { + if( a.value().size() > 0 ) { + return a.value().matrix().lpNorm (); + } + else { // this situation can occur when no wells are present + 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 ) { - return a.value().matrix().lpNorm (); + result = a.value().matrix().lpNorm (); } - else { // this situation can occur when no wells are present - return 0.0; +#if HAVE_MPI + if ( pinfo.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& real_info = + boost::any_cast(pinfo); + result = real_info.communicator().max(result); } +#endif + return result; } } // namespace detail @@ -1732,7 +1770,8 @@ namespace detail { const std::vector::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& R_sum, std::array& maxCoeff, std::array& B_avg, - int nc) const + std::vector& maxNormWell, + int nc, + int nw) const { // Do the global reductions #if HAVE_MPI @@ -1862,15 +1905,21 @@ namespace detail { Opm::Reduction::makeGlobalMaxFunctor(), Opm::Reduction::makeGlobalSumFunctor()); 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); + 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(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 B(nc, cols); Eigen::Array R(nc, cols); Eigen::Array tempV(nc, cols); + std::vector maxNormWell(MaxNumPhases); for ( int idx=0; idx( 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 Operator; Operator opA(istlA, istlComm); - constructPreconditionerAndSolve(opA, istlAe, x, istlb, istlComm, result); + constructPreconditionerAndSolve(opA, istlAe, x, istlb, istlComm, istlAeComm, result); } else #endif @@ -201,7 +205,7 @@ namespace Opm typedef Dune::MatrixAdapter 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 diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp index eb42ea16b..6a67ecd16 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.hpp +++ b/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 ScalarProductChooser; @@ -94,7 +98,8 @@ namespace Opm // typedef Dune::SeqILU0 Preconditioner; typedef Opm::CPRPreconditioner 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. diff --git a/opm/autodiff/RedistributeDataHandles.hpp b/opm/autodiff/RedistributeDataHandles.hpp index 29617014e..9b833b089 100644 --- a/opm/autodiff/RedistributeDataHandles.hpp +++ b/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 - 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; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 2df70b396..a2740789b 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -120,6 +120,8 @@ namespace Opm RateConverterType rateConverter_; // Threshold pressures. std::vector threshold_pressures_by_face_; + // Whether this a parallel simulation or not + bool is_parallel_run_; void computeRESV(const std::size_t step, @@ -204,7 +206,8 @@ namespace Opm const ParallelISTLInformation& info = boost::any_cast(solver_.parallelInformation()); // Only rank 0 does print to std::cout - terminal_output_= (info.communicator().rank()==0); + 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);