OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #262 from blattms/unify-amg-calls

Browse Source 
Unify amg calls
This commit is contained in:
Bård Skaflestad 2013-06-21 04:01:03 -07:00
commit 6afd914756
1 changed files with 56 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -59,16 +59,19 @@ namespace Opm
solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity);
LinearSolverInterface::LinearSolverReport
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity);
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps);
#ifdef HAS_DUNE_FAST_AMG
LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity);
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps);
LinearSolverInterface::LinearSolverReport
solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity);
solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double prolongateFactor);
#endif
LinearSolverInterface::LinearSolverReport
solveBiCGStab_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity);
} // anonymous namespace
@ -109,18 +112,10 @@ namespace Opm
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_);
}
LinearSolverIstl::~LinearSolverIstl()
{
}
{}
LinearSolverInterface::LinearSolverReport
LinearSolverIstl::solve(const int size,
@ -163,34 +158,38 @@ namespace Opm
std::copy(b.begin(), b.end(),
std::ostream_iterator<VectorBlockType>(rhsf, "\n"));
}
int maxit = linsolver_max_iterations_;
if (maxit == 0) {
maxit = 5000;
}
LinearSolverReport res;
LinearSolverReport res;
switch (linsolver_type_) {
case CG_ILU0:
res = solveCG_ILU0(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_);
break;
case CG_AMG:
res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_);
res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_);
break;
case KAMG:
#ifdef HAS_DUNE_FAST_AMG
res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_);
res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_);
#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_);
res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_);
#else
if(linsolver_verbosity_)
std::cerr<<"Fast AMG is not available; falling back to CG preconditioned with the normal one"<<std::endl;
res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_);
if(linsolver_verbosity_)
std::cerr<<"Fast AMG is not available; falling back to CG preconditioned with the normal one"<<std::endl;
res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_);
#endif
break;
case BiCGStab_ILU0:
@ -249,8 +248,20 @@ namespace Opm
#define SMOOTHER_ILU 0
#define ANISOTROPIC_3D 0
template<typename C>
void setUpCriterion(C& criterion, double linsolver_prolongate_factor,
int verbosity)
{
criterion.setDebugLevel(verbosity);
#if ANISOTROPIC_3D
criterion.setDefaultValuesAnisotropic(3, 2);
#endif
criterion.setProlongationDampingFactor(linsolver_prolongate_factor);
}
LinearSolverInterface::LinearSolverReport
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity)
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor, int linsolver_smooth_steps)
{
// Solve with AMG solver.
@ -274,19 +285,14 @@ namespace Opm
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Dune::Amg::AMG<Operator,Vector,Smoother> 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);
Precond::SmootherArgs smootherArgs;
Operator opA(A);
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity);
Precond precond(opA, criterion, smootherArgs, 1, linsolver_smooth_steps,
linsolver_smooth_steps);
// Construct linear solver.
Dune::CGSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity);
@ -305,7 +311,8 @@ namespace Opm
#ifdef HAS_DUNE_FAST_AMG
LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity)
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor, int linsolver_smooth_steps)
{
// Solve with AMG solver.
@ -328,18 +335,14 @@ namespace Opm
#endif
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
Operator opA(A);
// Construct preconditioner.
double relax = 1;
Operator opA(A);
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);
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity);
Precond precond(opA, criterion, smootherArgs, 1, linsolver_smooth_steps,
linsolver_smooth_steps);
// Construct linear solver.
Dune::GeneralizedPCGSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity);
@ -357,7 +360,8 @@ namespace Opm
}
LinearSolverInterface::LinearSolverReport
solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity)
solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor)
{
// Solve with AMG solver.
@ -370,27 +374,21 @@ namespace Opm
#if SYMMETRIC
typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Mat,CouplingMetric> > CriterionBase;
#else
typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Mat,CouplingMetric> > CriterionBase;
typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Mat,CouplingMetric> > CriterionBase;
#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);
Operator opA(A);
Criterion criterion;
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity);
Dune::Amg::Parameters parms;
parms.setDebugLevel(verbosity);
parms.setNoPreSmoothSteps(1);
parms.setNoPostSmoothSteps(1);
parms.setProlongationDampingFactor(linsolver_prolongate_factor);
Precond precond(opA, criterion, parms);
// Construct linear solver.
@ -440,4 +438,3 @@ namespace Opm
} // namespace Opm