From 5e9b28d7a9fabac87bdbc0f9a7b9cc7fdad62de3 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 18 Mar 2013 10:04:03 +0100 Subject: [PATCH 1/2] Added a fast amg version of AMG (with one step of Gaus-Seidel smoothing) and AMG with Krylov-cycle. The former is only available when using the inofficial 2.2.1 cmake release. The latter is currently not optimized. --- opm/core/linalg/LinearSolverIstl.cpp | 198 ++++++++++++++++++++++----- opm/core/linalg/LinearSolverIstl.hpp | 8 +- 2 files changed, 173 insertions(+), 33 deletions(-) diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index a3f9583ee..69353b10e 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -39,6 +39,7 @@ #include #include #include +#include #include @@ -46,14 +47,13 @@ namespace Opm { - using namespace Dune; // While not great, it's okay in a cpp file like this. namespace { - typedef FieldVector VectorBlockType; - typedef FieldMatrix MatrixBlockType; - typedef BCRSMatrix Mat; - typedef BlockVector Vector; - typedef MatrixAdapter Operator; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector Vector; + typedef Dune::MatrixAdapter Operator; LinearSolverInterface::LinearSolverReport solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); @@ -61,6 +61,14 @@ namespace Opm LinearSolverInterface::LinearSolverReport solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); + LinearSolverInterface::LinearSolverReport + solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); + +#ifdef HAS_DUNE_FAST_AMG + LinearSolverInterface::LinearSolverReport + solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); +#endif + LinearSolverInterface::LinearSolverReport solveBiCGStab_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); } // anonymous namespace @@ -73,7 +81,9 @@ namespace Opm linsolver_verbosity_(0), linsolver_type_(CG_AMG), linsolver_save_system_(false), - linsolver_max_iterations_(0) + linsolver_max_iterations_(0), + linsolver_smooth_steps_(2), + linsolver_prolongate_factor_(1.6) { } @@ -85,7 +95,9 @@ namespace Opm linsolver_verbosity_(0), linsolver_type_(CG_AMG), linsolver_save_system_(false), - linsolver_max_iterations_(0) + linsolver_max_iterations_(0), + linsolver_smooth_steps_(2), + linsolver_prolongate_factor_(1.6) { linsolver_residual_tolerance_ = param.getDefault("linsolver_residual_tolerance", linsolver_residual_tolerance_); linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_); @@ -95,6 +107,9 @@ namespace Opm linsolver_save_filename_ = param.getDefault("linsolver_save_filename", std::string("linsys")); } linsolver_max_iterations_ = param.getDefault("linsolver_max_iterations", linsolver_max_iterations_); + linsolver_smooth_steps_ = param.getDefault("linsolver_smooth_steps", linsolver_smooth_steps_); + linsolver_prolongate_factor_ = param.getDefault("linsolver_prolongate_factor", linsolver_prolongate_factor_); + } @@ -151,10 +166,10 @@ namespace Opm int maxit = linsolver_max_iterations_; if (maxit == 0) { - maxit = A.N(); + maxit = 5000; } - LinearSolverReport res; + LinearSolverReport res; switch (linsolver_type_) { case CG_ILU0: res = solveCG_ILU0(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); @@ -162,6 +177,17 @@ namespace Opm case CG_AMG: res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); break; + case KAMG: + res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); + break; + case FastAMG: +#ifdef HAS_DUNE_FAST_AMG + res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); +#else + #warning "Fast AMG is not available; falling back to CG preconditioned with the normal one" + res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); +#endif + break; case BiCGStab_ILU0: res = solveBiCGStab_ILU0(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); break; @@ -194,13 +220,13 @@ namespace Opm Operator opA(A); // Construct preconditioner. - SeqILU0 precond(A, 1.0); + Dune::SeqILU0 precond(A, 1.0); // Construct linear solver. - CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); // Solve system. - InverseOperatorResult result; + Dune::InverseOperatorResult result; linsolve.apply(x, b, result); // Output results. @@ -213,35 +239,35 @@ namespace Opm +#define FIRST_DIAGONAL 1 +#define SYMMETRIC 1 +#define SMOOTHER_ILU 0 +#define ANISOTROPIC_3D 0 LinearSolverInterface::LinearSolverReport solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) { // Solve with AMG solver. -#define FIRST_DIAGONAL 1 -#define SYMMETRIC 1 -#define SMOOTHER_ILU 1 -#define ANISOTROPIC_3D 0 #if FIRST_DIAGONAL - typedef Amg::FirstDiagonal CouplingMetric; + typedef Dune::Amg::FirstDiagonal CouplingMetric; #else - typedef Amg::RowSum CouplingMetric; + typedef Dune::Amg::RowSum CouplingMetric; #endif #if SYMMETRIC - typedef Amg::SymmetricCriterion CriterionBase; + typedef Dune::Amg::SymmetricCriterion CriterionBase; #else - typedef Amg::UnSymmetricCriterion CriterionBase; + typedef Dune::Amg::UnSymmetricCriterion CriterionBase; #endif #if SMOOTHER_ILU - typedef SeqILU0 Smoother; + typedef Dune::SeqILU0 Smoother; #else - typedef SeqSSOR Smoother; + typedef Dune::SeqSOR Smoother; #endif - typedef Amg::CoarsenCriterion Criterion; - typedef Amg::AMG Precond; + typedef Dune::Amg::CoarsenCriterion Criterion; + typedef Dune::Amg::AMG Precond; Operator opA(A); @@ -254,13 +280,13 @@ namespace Opm #if ANISOTROPIC_3D criterion.setDefaultValuesAnisotropic(3, 2); #endif - Precond precond(opA, criterion, smootherArgs); - + Precond precond(opA, criterion, smootherArgs, 1,1,1); + // Construct linear solver. - CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); // Solve system. - InverseOperatorResult result; + Dune::InverseOperatorResult result; linsolve.apply(x, b, result); // Output results. @@ -272,6 +298,114 @@ namespace Opm } + LinearSolverInterface::LinearSolverReport + solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) + { + // Solve with AMG solver. + +#if FIRST_DIAGONAL + typedef Dune::Amg::FirstDiagonal CouplingMetric; +#else + typedef Dune::Amg::RowSum CouplingMetric; +#endif + +#if SYMMETRIC + typedef Dune::Amg::SymmetricCriterion CriterionBase; +#else + typedef Dune::Amg::UnSymmetricCriterion CriterionBase; +#endif + +#if SMOOTHER_ILU + typedef Dune::SeqILU0 Smoother; +#else + typedef Dune::SeqSOR Smoother; +#endif + typedef Dune::Amg::CoarsenCriterion Criterion; + typedef Dune::Amg::KAMG, std::allocator > Precond; + + Operator opA(A); + + // Construct preconditioner. + double relax = 1; + Precond::SmootherArgs smootherArgs; + smootherArgs.relaxationFactor = relax; + Criterion criterion; + criterion.setDebugLevel(verbosity); +#if ANISOTROPIC_3D + criterion.setDefaultValuesAnisotropic(3, 2); +#endif + Precond precond(opA, criterion, smootherArgs, 1,1,1,2,.2); + + // Construct linear solver. + Dune::GeneralizedPCGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + + // Solve system. + Dune::InverseOperatorResult result; + linsolve.apply(x, b, result); + + // Output results. + LinearSolverInterface::LinearSolverReport res; + res.converged = result.converged; + res.iterations = result.iterations; + res.residual_reduction = result.reduction; + return res; + } + +#ifdef HAS_DUNE_FAST_AMG + LinearSolverInterface::LinearSolverReport + solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) + { + // Solve with AMG solver. + +#if FIRST_DIAGONAL + typedef Dune::Amg::FirstDiagonal CouplingMetric; +#else + typedef Dune::Amg::RowSum CouplingMetric; +#endif + +#if SYMMETRIC + typedef Dune::Amg::AggregationCriterion > CriterionBase; +#else + typedef Dune::Amg::AggregationCriterion > CriterionBase; + #warn "Only symmetric matrices are supported currently. Computing anyway..." +#endif + + typedef Dune::Amg::CoarsenCriterion Criterion; + typedef Dune::Amg::FastAMG Precond; + + Operator opA(A); + + // Construct preconditioner. + Criterion criterion; + criterion.setDebugLevel(verbosity); +#if ANISOTROPIC_3D + criterion.setDefaultValuesAnisotropic(3, 2); +#else + criterion.setDefaultValuesIsotropic(3,2); +#endif + criterion.setProlongationDampingFactor(1.6); + Dune::Amg::Parameters parms; + parms.setDebugLevel(verbosity); + parms.setNoPreSmoothSteps(1); + parms.setNoPostSmoothSteps(1); + Precond precond(opA, criterion, parms); + + // Construct linear solver. + Dune::GeneralizedPCGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + + // Solve system. + Dune::InverseOperatorResult result; + linsolve.apply(x, b, result); + + // Output results. + LinearSolverInterface::LinearSolverReport res; + res.converged = result.converged; + res.iterations = result.iterations; + res.residual_reduction = result.reduction; + return res; + } +#endif + LinearSolverInterface::LinearSolverReport solveBiCGStab_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) @@ -279,13 +413,13 @@ namespace Opm Operator opA(A); // Construct preconditioner. - SeqILU0 precond(A, 1.0); + Dune::SeqILU0 precond(A, 1.0); // Construct linear solver. - BiCGSTABSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::BiCGSTABSolver linsolve(opA, precond, tolerance, maxit, verbosity); // Solve system. - InverseOperatorResult result; + Dune::InverseOperatorResult result; linsolve.apply(x, b, result); // Output results. diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp index 6de42e725..a12fb7eeb 100644 --- a/opm/core/linalg/LinearSolverIstl.hpp +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -40,6 +40,7 @@ namespace Opm /// linsolver_verbosity 0 /// linsolver_type 1 ( = CG_AMG), alternatives are: /// CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2 + /// FastAMG=3, KAMG=4 }; /// linsolver_save_system false /// linsolver_save_filename /// linsolver_max_iterations 0 (unlimited) @@ -83,11 +84,16 @@ namespace Opm private: double linsolver_residual_tolerance_; int linsolver_verbosity_; - enum LinsolverType { CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2 }; + enum LinsolverType { CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2, FastAMG=3, KAMG=4 }; LinsolverType linsolver_type_; bool linsolver_save_system_; std::string linsolver_save_filename_; int linsolver_max_iterations_; + /** \brief The number smoothing steps to apply in AMG. */ + int linsolver_smooth_steps_; + /** \brief The factor to scale the coarse grid correction with. */ + double linsolver_prolongate_factor_; + }; From a53d4f69c90cea6651a49352452d6f65346f9cd4 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 18 Mar 2013 10:10:38 +0100 Subject: [PATCH 2/2] Adapted the documentation. --- opm/core/linalg/LinearSolverIstl.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp index a12fb7eeb..b1f523484 100644 --- a/opm/core/linalg/LinearSolverIstl.hpp +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -43,7 +43,11 @@ namespace Opm /// FastAMG=3, KAMG=4 }; /// linsolver_save_system false /// linsolver_save_filename - /// linsolver_max_iterations 0 (unlimited) + /// linsolver_max_iterations 0 (unlimited=5000) + /// linsolver_residual_tolerance 1e-8 + /// linsolver_smooth_steps 2 + /// linsolver_prolongate_factor 1.6 + /// linsolver_verbosity 0 LinearSolverIstl(); /// Construct from parameters