diff --git a/opm/core/linalg/LinearSolverAGMG.cpp b/opm/core/linalg/LinearSolverAGMG.cpp index 086c9316..2c774cda 100644 --- a/opm/core/linalg/LinearSolverAGMG.cpp +++ b/opm/core/linalg/LinearSolverAGMG.cpp @@ -69,12 +69,15 @@ namespace Opm LinearSolverAGMG::LinearSolverAGMG(const parameter::ParameterGroup& param) : max_it_(100) , - rtol_ (1.0e-6), - is_spd_(false) + rtol_ (1.0e-8), + is_spd_(false), + linsolver_verbosity_(0) { max_it_ = param.getDefault("max_it", max_it_); rtol_ = param.getDefault("rtol" , rtol_ ); is_spd_ = param.getDefault("is_spd", is_spd_); + linsolver_verbosity_ = param.getDefault("linsolver_verbosity", + linsolver_verbosity_); } LinearSolverAGMG::~LinearSolverAGMG() {} @@ -121,7 +124,7 @@ namespace Opm nrest = 10; // Suggested default number of GCR restarts. } - int iprint = 0; // Suppress most output + int iprint = linsolver_verbosity_; // Suppress most output DAGMG_(& size, & a[0], & j[0], & i[0], rhs, solution, & ijob, & iprint, & nrest, & rpt.iterations, & rtol_); diff --git a/opm/core/linalg/LinearSolverAGMG.hpp b/opm/core/linalg/LinearSolverAGMG.hpp index 9038ba75..82160d15 100644 --- a/opm/core/linalg/LinearSolverAGMG.hpp +++ b/opm/core/linalg/LinearSolverAGMG.hpp @@ -60,7 +60,7 @@ namespace Opm * systems are solved using GCR. */ LinearSolverAGMG(const int max_it = 100 , - const double rtol = 1.0e-6, + const double rtol = 1.0e-8, const bool is_spd = false); /** @@ -104,6 +104,7 @@ namespace Opm int max_it_; double rtol_ ; bool is_spd_; + int linsolver_verbosity_; }; diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index a3f9583e..8f6fabaf 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); +#ifdef HAS_DUNE_FAST_AMG + LinearSolverInterface::LinearSolverReport + solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); + + 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,22 @@ namespace Opm case CG_AMG: res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); break; + case KAMG: +#ifdef HAS_DUNE_FAST_AMG + res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); +#else + throw std::runtime_error("KAMG not supported with this version of DUNE"); +#endif + break; + case FastAMG: +#ifdef HAS_DUNE_FAST_AMG + res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); +#else + if(linsolver_verbosity_) + std::cerr<<"Fast AMG is not available; falling back to CG preconditioned with the normal one"< 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 +244,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 +285,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 +303,112 @@ namespace Opm } +#ifdef HAS_DUNE_FAST_AMG + 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 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; + } + + 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; +#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 +416,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 6de42e72..b1f52348 100644 --- a/opm/core/linalg/LinearSolverIstl.hpp +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -40,9 +40,14 @@ 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) + /// 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 @@ -83,11 +88,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_; + };