diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index d52df9c96..b565ad0d1 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -66,6 +66,8 @@ struct CPRParameter double cpr_solver_tol_; int cpr_ilu_n_; MILU_VARIANT cpr_ilu_milu_; + bool cpr_ilu_redblack_; + bool cpr_ilu_reorder_sphere_; int cpr_max_ell_iter_; bool cpr_use_amg_; bool cpr_use_bicgstab_; @@ -79,13 +81,15 @@ struct CPRParameter // reset values to default reset(); - cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_); - cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_); - cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_); - cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_); - cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_); - cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_); - cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_); + cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_); + cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_); + cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_); + cpr_ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); + cpr_ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); + cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_); + cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_); + cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_); + cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_); cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_); std::string milu("ILU"); @@ -94,14 +98,16 @@ struct CPRParameter void reset() { - cpr_relax_ = 1.0; - cpr_solver_tol_ = 1e-2; - cpr_ilu_n_ = 0; - cpr_ilu_milu_ = MILU_VARIANT::ILU; - cpr_max_ell_iter_ = 25; - cpr_use_amg_ = true; - cpr_use_bicgstab_ = true; - cpr_solver_verbose_ = false; + cpr_relax_ = 1.0; + cpr_solver_tol_ = 1e-2; + cpr_ilu_n_ = 0; + cpr_ilu_milu_ = MILU_VARIANT::ILU; + cpr_ilu_redblack_ = false; + cpr_ilu_reorder_sphere_ = true; + cpr_max_ell_iter_ = 25; + cpr_use_amg_ = true; + cpr_use_bicgstab_ = true; + cpr_solver_verbose_ = false; cpr_pressure_aggregation_ = false; } }; diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index 745b9e60e..7692dbbdd 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/opm/autodiff/ISTLSolver.hpp @@ -497,7 +497,9 @@ namespace Detail const double relax = parameters_.ilu_relaxation_; const int ilu_fillin = parameters_.ilu_fillin_level_; const MILU_VARIANT ilu_milu = parameters_.ilu_milu_; - std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu)); + const bool ilu_redblack = parameters_.ilu_redblack_; + const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_; + std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres)); return precond; } @@ -520,7 +522,9 @@ namespace Detail typedef std::unique_ptr Pointer; const double relax = parameters_.ilu_relaxation_; const MILU_VARIANT ilu_milu = parameters_.ilu_milu_; - return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu)); + const bool ilu_redblack = parameters_.ilu_redblack_; + const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_; + return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres)); } #endif diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp index 668f2b4f8..683993bfc 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp @@ -46,6 +46,8 @@ namespace Opm int linear_solver_verbosity_; int ilu_fillin_level_; Opm::MILU_VARIANT ilu_milu_; + bool ilu_redblack_; + bool ilu_reorder_sphere_; bool newton_use_gmres_; bool require_full_sparsity_pattern_; bool ignoreConvergenceFailure_; @@ -71,6 +73,8 @@ namespace Opm linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ ); ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ ); ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ ); + ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_); + ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_); std::string milu("ILU"); ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu)); @@ -94,6 +98,8 @@ namespace Opm ilu_fillin_level_ = 0; ilu_relaxation_ = 0.9; ilu_milu_ = MILU_VARIANT::ILU; + ilu_redblack_ = false; + ilu_reorder_sphere_ = true; } }; diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp index 66005471d..7a92145ba 100644 --- a/opm/autodiff/ParallelOverlappingILU0.hpp +++ b/opm/autodiff/ParallelOverlappingILU0.hpp @@ -20,11 +20,13 @@ #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED +#include #include #include #include #include #include +#include #include #include @@ -396,8 +398,6 @@ class ParallelOverlappingILU0 public: - enum{ - }; //! \brief The matrix type the preconditioner is for. typedef typename std::remove_const::type matrix_type; //! \brief The domain type of the preconditioner. @@ -479,11 +479,17 @@ public: \param n ILU fill in level (for testing). This does not work in parallel. \param w The relaxation factor. \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. + \param redblack Whether to use a red-black ordering. + \param reorder_sphere If true, we start the reordering at a root node. + The vertices on each layer aound it (same distance) are + ordered consecutivly. If false, we preserver the order of + the vertices with the same color. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, const int n, const field_type w, - MILU_VARIANT milu) + MILU_VARIANT milu, bool redblack=false, + bool reorder_sphere=true) : lower_(), upper_(), inv_(), @@ -492,7 +498,8 @@ public: { // BlockMatrix is a Subclass of FieldMatrix that just adds // methods. Therefore this cast should be safe. - init( reinterpret_cast(A), n, milu ); + init( reinterpret_cast(A), n, milu, redblack, + reorder_sphere ); } /*! \brief Constructor gets all parameters to operate the prec. @@ -501,11 +508,17 @@ public: \param n ILU fill in level (for testing). This does not work in parallel. \param w The relaxation factor. \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. + \param redblack Whether to use a red-black ordering. + \param reorder_sphere If true, we start the reordering at a root node. + The vertices on each layer aound it (same distance) are + ordered consecutivly. If false, we preserver the order of + the vertices with the same color. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, const ParallelInfo& comm, const int n, const field_type w, - MILU_VARIANT milu) + MILU_VARIANT milu, bool redblack=false, + bool reorder_sphere=true) : lower_(), upper_(), inv_(), @@ -514,7 +527,8 @@ public: { // BlockMatrix is a Subclass of FieldMatrix that just adds // methods. Therefore this cast should be safe. - init( reinterpret_cast(A), n, milu ); + init( reinterpret_cast(A), n, milu, redblack, + reorder_sphere ); } /*! \brief Constructor. @@ -523,11 +537,17 @@ public: \param A The matrix to operate on. \param w The relaxation factor. \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. + \param redblack Whether to use a red-black ordering. + \param reorder_sphere If true, we start the reordering at a root node. + The vertices on each layer aound it (same distance) are + ordered consecutivly. If false, we preserver the order of + the vertices with the same color. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, - const field_type w, MILU_VARIANT milu) - : ParallelOverlappingILU0( A, 0, w, milu ) + const field_type w, MILU_VARIANT milu, bool redblack=false, + bool reorder_sphere=true) + : ParallelOverlappingILU0( A, 0, w, milu, redblack, reorder_sphere ) { } @@ -538,11 +558,17 @@ public: \param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication \param w The relaxation factor. \param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT. + \param redblack Whether to use a red-black ordering. + \param reorder_sphere If true, we start the reordering at a root node. + The vertices on each layer aound it (same distance) are + ordered consecutivly. If false, we preserver the order of + the vertices with the same color. */ template ParallelOverlappingILU0 (const Dune::BCRSMatrix& A, const ParallelInfo& comm, const field_type w, - MILU_VARIANT milu) + MILU_VARIANT milu, bool redblack=false, + bool reorder_sphere=true) : lower_(), upper_(), inv_(), @@ -551,7 +577,8 @@ public: { // BlockMatrix is a Subclass of FieldMatrix that just adds // methods. Therefore this cast should be safe. - init( reinterpret_cast(A), 0, milu ); + init( reinterpret_cast(A), 0, milu, redblack, + reorder_sphere ); } /*! @@ -572,7 +599,8 @@ public: */ virtual void apply (Domain& v, const Range& d) { - Range& md = const_cast(d); + Range& md = reorderD(d); + Domain& mv = reorderV(v); copyOwnerToAll( md ); // iterator types @@ -589,41 +617,42 @@ public: // lower triangular solve for( size_type i=0; i @@ -645,7 +674,7 @@ public: } protected: - void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu ) + void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres ) { // (For older DUNE versions the communicator might be // invalid if redistribution in AMG happened on the coarset level. @@ -669,11 +698,65 @@ protected: std::unique_ptr< Matrix > ILU; + if ( redBlack && iluIteration==0 ) + { + using Graph = Dune::Amg::MatrixGraph; + Graph graph(A); + auto colorsTuple = colorVerticesWelshPowell(graph); + const auto& colors = std::get<0>(colorsTuple); + const auto& verticesPerColor = std::get<2>(colorsTuple); + auto noColors = std::get<1>(colorsTuple); + if ( reorderSpheres ) + { + ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor, + graph, 0); + } + else + { + ordering_ = reorderVerticesPreserving(colors, noColors, verticesPerColor, + graph); + } + } + + std::vector inverseOrdering(ordering_.size()); + std::size_t index = 0; + for( auto newIndex: ordering_) + { + inverseOrdering[newIndex] = index++; + } + try { if( iluIteration == 0 ) { // create ILU-0 decomposition - ILU.reset( new Matrix( A ) ); + if ( ordering_.empty() ) + { + ILU.reset( new Matrix( A ) ); + } + else + { + ILU.reset( new Matrix(A.N(), A.M(), A.nonzeroes(), Matrix::row_wise)); + auto& newA = *ILU; + // Create sparsity pattern + for(auto iter=newA.createbegin(), iend = newA.createend(); iter != iend; ++iter) + { + const auto& row = A[inverseOrdering[iter.index()]]; + for(auto col = row.begin(), cend = row.end(); col != cend; ++col) + { + iter.insert(ordering_[col.index()]); + } + } + // Copy values + for(auto iter = A.begin(), iend = A.end(); iter != iend; ++iter) + { + auto& newRow = newA[ordering_[iter.index()]]; + for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col) + { + newRow[ordering_[col.index()]] = *col; + } + } + } + switch ( milu ) { case MILU_VARIANT::MILU_1: @@ -727,11 +810,66 @@ protected: detail::convertToCRS( *ILU, lower_, upper_, inv_ ); } + /// \brief Reorder D if needed and return a reference to it. + Range& reorderD(const Range& d) + { + if ( ordering_.empty()) + { + return const_cast(d); + } + else + { + reorderedD_.resize(d.size()); + std::size_t i = 0; + for(auto index: ordering_) + { + reorderedD_[index]=d[i++]; + } + return reorderedD_; + } + } + + /// \brief Reorder V if needed and return a reference to it. + Domain& reorderV(Domain& v) + { + if ( ordering_.empty()) + { + return v; + } + else + { + reorderedV_.resize(v.size()); + std::size_t i = 0; + for(auto index: ordering_) + { + reorderedV_[index]=v[i++]; + } + return reorderedV_; + } + } + + void reorderBack(const Range& reorderedV, Range& v) + { + if ( !ordering_.empty() ) + { + std::size_t i = 0; + for(auto index: ordering_) + { + v[i++] = reorderedV[index]; + } + } + } protected: //! \brief The ILU0 decomposition of the matrix. CRS lower_; CRS upper_; std::vector< block_type > inv_; + //! \brief the reordering of the unknowns + std::vector< std::size_t > ordering_; + //! \brief The reordered right hand side + Range reorderedD_; + //! \brief The reordered left hand side. + Domain reorderedV_; const ParallelInfo* comm_; //! \brief The relaxation factor to use.