OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added a fast amg version of AMG (with one step of Gaus-Seidel smoothing) and AMG with Krylov-cycle.

Browse Source 
The former is only available when using the inofficial 2.2.1 cmake release.
The latter is currently not optimized.
This commit is contained in:
Markus Blatt 2013-03-18 10:04:03 +01:00
parent bb944b7e52
commit 2c4c4cbb2f
2 changed files with 173 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

198
opm/core/linalg/LinearSolverIstl.cpp
View File

@ -39,6 +39,7 @@
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
#include <stdexcept>
@ -46,14 +47,13 @@
namespace Opm
{
using namespace Dune; // While not great, it's okay in a cpp file like this.
namespace {
typedef FieldVector<double, 1 > VectorBlockType;
typedef FieldMatrix<double, 1, 1> MatrixBlockType;
typedef BCRSMatrix <MatrixBlockType> Mat;
typedef BlockVector<VectorBlockType> Vector;
typedef MatrixAdapter<Mat,Vector,Vector> Operator;
typedef Dune::FieldVector<double, 1 > VectorBlockType;
typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> Vector;
typedef Dune::MatrixAdapter<Mat,Vector,Vector> 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<Mat,Vector,Vector> precond(A, 1.0);
Dune::SeqILU0<Mat,Vector,Vector> precond(A, 1.0);
// Construct linear solver.
CGSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity);
Dune::CGSolver<Vector> 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<Mat,CouplingMetric> CriterionBase;
typedef Dune::Amg::SymmetricCriterion<Mat,CouplingMetric> CriterionBase;
#else
typedef Amg::UnSymmetricCriterion<Mat,CouplingMetric> CriterionBase;
typedef Dune::Amg::UnSymmetricCriterion<Mat,CouplingMetric> CriterionBase;
#endif
#if SMOOTHER_ILU
typedef SeqILU0<Mat,Vector,Vector> Smoother;
typedef Dune::SeqILU0<Mat,Vector,Vector> Smoother;
#else
typedef SeqSSOR<Mat,Vector,Vector> Smoother;
typedef Dune::SeqSOR<Mat,Vector,Vector> Smoother;
#endif
typedef Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Amg::AMG<Operator,Vector,Smoother> Precond;
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Dune::Amg::AMG<Operator,Vector,Smoother> 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<Vector> linsolve(opA, precond, tolerance, maxit, verbosity);
Dune::CGSolver<Vector> 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<Mat,CouplingMetric> CriterionBase;
#else
typedef Dune::Amg::UnSymmetricCriterion<Mat,CouplingMetric> CriterionBase;
#endif
#if SMOOTHER_ILU
typedef Dune::SeqILU0<Mat,Vector,Vector> Smoother;
#else
typedef Dune::SeqSOR<Mat,Vector,Vector> Smoother;
#endif
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation,Dune::GeneralizedPCGSolver<Vector>, std::allocator<Vector> > 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<Vector> 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<Dune::Amg::SymmetricMatrixDependency<Mat,CouplingMetric> > CriterionBase;
#else
typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Mat,CouplingMetric> > CriterionBase;
#warn "Only symmetric matrices are supported currently. Computing anyway..."
#endif
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Dune::Amg::FastAMG<Operator,Vector> 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<Vector> 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<Mat,Vector,Vector> precond(A, 1.0);
Dune::SeqILU0<Mat,Vector,Vector> precond(A, 1.0);
// Construct linear solver.
BiCGSTABSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity);
Dune::BiCGSTABSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity);
// Solve system.
InverseOperatorResult result;
Dune::InverseOperatorResult result;
linsolve.apply(x, b, result);
// Output results.

8
opm/core/linalg/LinearSolverIstl.hpp
View File

@ -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 <empty string>
/// 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_;
};