OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-23 07:53:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Explicitly construct the ScalarProduct, and SequentialInformation.

Browse Source 
This patch refactors the calls to the dune-istl solvers.
The SeqScalarProduct, and SequentialInformation is now explicitly
constructed and later used to construct the smoothers
generically. Aditionally the linear operator (MatrixAdapter) is
constructed before calling the various solver dependent solve
routines.

While this does not change the behaviour of the code it is a
preparatory step to support parallel solvers. These parallel
solvers only differ in the type of the scalarproduct and
linear operator used from the sequential ones.
This commit is contained in:
Markus Blatt 2014-03-19 17:37:12 +01:00
parent 2f3c17bd9d
commit e1d0995104
2 changed files with 80 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -39,6 +39,7 @@
#include <dune/istl/solvers.hh> #include <dune/istl/solvers.hh>
#include <dune/istl/paamg/amg.hh> #include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh> #include <dune/istl/paamg/kamg.hh>
#include <dune/istl/paamg/pinfo.hh>
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
#include <dune/istl/paamg/fastamg.hh> #include <dune/istl/paamg/fastamg.hh>
@ -59,25 +60,30 @@ namespace Opm
typedef Dune::BlockVector<VectorBlockType> Vector; typedef Dune::BlockVector<VectorBlockType> Vector;
typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator; typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); solveCG_ILU0(O& A, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity);
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, solveCG_AMG(O& A, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps); double prolongateFactor, int smoothsteps);
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, solveKAMG(O& A, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps); double prolongateFactor, int smoothsteps);
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, solveFastAMG(O& A, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity,
double prolongateFactor); double prolongateFactor);
#endif #endif
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveBiCGStab_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); solveBiCGStab_ILU0(O& A, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity);
} // anonymous namespace } // anonymous namespace
@ -167,19 +173,29 @@ namespace Opm
if (maxit == 0) { if (maxit == 0) {
maxit = 5000; maxit = 5000;
} }
Dune::SeqScalarProduct<Vector> sp;
Dune::Amg::SequentialInformation comm;
Operator opA(A);
return solveSystem(opA, x, solution, b, sp, comm, maxit);
}
template<class O, class V, class S, class C>
LinearSolverInterface::LinearSolverReport
LinearSolverIstl::solveSystem (O& opA, V& x, double* solution, V& b,
S& sp, const C& comm, int maxit) const
{
LinearSolverReport res; LinearSolverReport res;
switch (linsolver_type_) { switch (linsolver_type_) {
case CG_ILU0: case CG_ILU0:
res = solveCG_ILU0(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); res = solveCG_ILU0(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_);
break; break;
case CG_AMG: case CG_AMG:
res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_, res = solveCG_AMG(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_); linsolver_prolongate_factor_, linsolver_smooth_steps_);
break; break;
case KAMG: case KAMG:
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_, res = solveKAMG(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_); linsolver_prolongate_factor_, linsolver_smooth_steps_);
#else #else
throw std::runtime_error("KAMG not supported with this version of DUNE"); throw std::runtime_error("KAMG not supported with this version of DUNE");
@ -187,17 +203,17 @@ namespace Opm
break; break;
case FastAMG: case FastAMG:
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_, res = solveFastAMG(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_); linsolver_prolongate_factor_);
#else #else
if(linsolver_verbosity_) if(linsolver_verbosity_)
std::cerr<<"Fast AMG is not available; falling back to CG preconditioned with the normal one"<<std::endl; 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_, res = solveCG_AMG(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_); linsolver_prolongate_factor_, linsolver_smooth_steps_);
#endif #endif
break; break;
case BiCGStab_ILU0: case BiCGStab_ILU0:
res = solveBiCGStab_ILU0(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); res = solveBiCGStab_ILU0(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_);
break; break;
default: default:
std::cerr << "Unknown linsolver_type: " << int(linsolver_type_) << '\n'; std::cerr << "Unknown linsolver_type: " << int(linsolver_type_) << '\n';
@ -221,17 +237,37 @@ namespace Opm
namespace namespace
{ {
template<class P, class O>
LinearSolverInterface::LinearSolverReport struct PreconditionerTraits
solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) {
typedef std::shared_ptr<P> PointerType;
};
template<class P, class O, class C>
typename PreconditionerTraits<P,O>::PointerType
makePreconditioner(O& opA, double relax, const C& comm, int iterations=1)
{
typename Dune::Amg::SmootherTraits<P>::Arguments args;
typename Dune::Amg::ConstructionTraits<P>::Arguments cargs;
cargs.setMatrix(opA.getmat());
args.iterations=iterations;
args.relaxationFactor=relax;
cargs.setArgs(args);
cargs.setComm(comm);
return std::shared_ptr<P>(Dune::Amg::ConstructionTraits<P>::construct(cargs));
}
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport
solveCG_ILU0(O& opA, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity)
{ {
Operator opA(A);
// Construct preconditioner. // Construct preconditioner.
Dune::SeqILU0<Mat,Vector,Vector> precond(A, 1.0); typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
auto precond = makePreconditioner<Preconditioner>(opA, 1.0, comm);
// Construct linear solver. // Construct linear solver.
Dune::CGSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity); Dune::CGSolver<Vector> linsolve(opA, sp, *precond, tolerance, maxit, verbosity);
// Solve system. // Solve system.
Dune::InverseOperatorResult result; Dune::InverseOperatorResult result;
@ -266,8 +302,9 @@ namespace Opm
criterion.setGamma(1); // V-cycle; this is the default criterion.setGamma(1); // V-cycle; this is the default
} }
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, solveCG_AMG(O& opA, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor, int linsolver_smooth_steps) double linsolver_prolongate_factor, int linsolver_smooth_steps)
{ {
// Solve with AMG solver. // Solve with AMG solver.
@ -295,13 +332,12 @@ namespace Opm
// Construct preconditioner. // Construct preconditioner.
Criterion criterion; Criterion criterion;
Precond::SmootherArgs smootherArgs; Precond::SmootherArgs smootherArgs;
Operator opA(A);
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
linsolver_smooth_steps); linsolver_smooth_steps);
Precond precond(opA, criterion, smootherArgs); Precond precond(opA, criterion, smootherArgs);
// Construct linear solver. // Construct linear solver.
Dune::CGSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity); Dune::CGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
// Solve system. // Solve system.
Dune::InverseOperatorResult result; Dune::InverseOperatorResult result;
@ -317,8 +353,9 @@ namespace Opm
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, solveKAMG(O& opA, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor, int linsolver_smooth_steps) double linsolver_prolongate_factor, int linsolver_smooth_steps)
{ {
// Solve with AMG solver. // Solve with AMG solver.
@ -344,7 +381,6 @@ namespace Opm
typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond; typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
// Construct preconditioner. // Construct preconditioner.
Operator opA(A);
Precond::SmootherArgs smootherArgs; Precond::SmootherArgs smootherArgs;
Criterion criterion; Criterion criterion;
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
@ -352,7 +388,7 @@ namespace Opm
Precond precond(opA, criterion, smootherArgs); Precond precond(opA, criterion, smootherArgs);
// Construct linear solver. // Construct linear solver.
Dune::GeneralizedPCGSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity); Dune::GeneralizedPCGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
// Solve system. // Solve system.
Dune::InverseOperatorResult result; Dune::InverseOperatorResult result;
@ -366,8 +402,9 @@ namespace Opm
return res; return res;
} }
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveFastAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity, solveFastAMG(O& opA, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor) double linsolver_prolongate_factor)
{ {
// Solve with AMG solver. // Solve with AMG solver.
@ -388,7 +425,6 @@ namespace Opm
typedef Dune::Amg::FastAMG<Operator,Vector> Precond; typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
// Construct preconditioner. // Construct preconditioner.
Operator opA(A);
Criterion criterion; Criterion criterion;
const int smooth_steps = 1; const int smooth_steps = 1;
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, smooth_steps); setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, smooth_steps);
@ -400,7 +436,7 @@ namespace Opm
Precond precond(opA, criterion, parms); Precond precond(opA, criterion, parms);
// Construct linear solver. // Construct linear solver.
Dune::GeneralizedPCGSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity); Dune::GeneralizedPCGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
// Solve system. // Solve system.
Dune::InverseOperatorResult result; Dune::InverseOperatorResult result;
@ -415,17 +451,17 @@ namespace Opm
} }
#endif #endif
template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
solveBiCGStab_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) solveBiCGStab_ILU0(O& opA, Vector& x, Vector& b, S& sp, const C& comm, double tolerance, int maxit, int verbosity)
{ {
Operator opA(A);
// Construct preconditioner. // Construct preconditioner.
Dune::SeqILU0<Mat,Vector,Vector> precond(A, 1.0); typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
auto precond = makePreconditioner<Preconditioner>(opA, 1.0, comm);
// Construct linear solver. // Construct linear solver.
Dune::BiCGSTABSolver<Vector> linsolve(opA, precond, tolerance, maxit, verbosity); Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, *precond, tolerance, maxit, verbosity);
// Solve system. // Solve system.
Dune::InverseOperatorResult result; Dune::InverseOperatorResult result;

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

@ -86,6 +86,19 @@ namespace Opm
virtual double getTolerance() const; virtual double getTolerance() const;
private: private:
/// \brief Solve the linear system using ISTL
/// \param[in] opA The linear operator of the system to solve.
/// \param[in,out] x The vector with the initial guess that will store
/// the computed solution
/// \param[out] solution C array for storing the solution vector.
/// \param[in] b The vector containing the right hand side of the system.
/// \param[in] sp The scalar product to use.
/// \param[in] comm The information about the parallel domain decomposition.
/// \param[in] maxit The maximum number of iterations allowed.
template<class O, class V, class S, class C>
LinearSolverReport solveSystem(O& opA, V& x, double* solution,
V& b, S& sp, const C& comm, int maxit) const;
double linsolver_residual_tolerance_; double linsolver_residual_tolerance_;
int linsolver_verbosity_; int linsolver_verbosity_;
enum LinsolverType { CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2, FastAMG=3, KAMG=4 }; enum LinsolverType { CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2, FastAMG=3, KAMG=4 };