From e1d0995104c95d34a4f2d956eb870d51bdbccc59 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 19 Mar 2014 17:37:12 +0100 Subject: [PATCH] Explicitly construct the ScalarProduct, and SequentialInformation. 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. --- opm/core/linalg/LinearSolverIstl.cpp | 98 +++++++++++++++++++--------- opm/core/linalg/LinearSolverIstl.hpp | 13 ++++ 2 files changed, 80 insertions(+), 31 deletions(-) diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index c32ad02ee..181b9ed2f 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -39,6 +39,7 @@ #include #include #include +#include #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) #include @@ -59,25 +60,30 @@ namespace Opm typedef Dune::BlockVector Vector; typedef Dune::MatrixAdapter Operator; + template 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 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); #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) + template 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); + template 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); #endif + template 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 @@ -167,19 +173,29 @@ namespace Opm if (maxit == 0) { maxit = 5000; } + Dune::SeqScalarProduct sp; + Dune::Amg::SequentialInformation comm; + Operator opA(A); + return solveSystem(opA, x, solution, b, sp, comm, maxit); + } +template + LinearSolverInterface::LinearSolverReport + LinearSolverIstl::solveSystem (O& opA, V& x, double* solution, V& b, + S& sp, const C& comm, int maxit) const + { LinearSolverReport res; switch (linsolver_type_) { 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; 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_); break; case KAMG: #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_); #else throw std::runtime_error("KAMG not supported with this version of DUNE"); @@ -187,17 +203,17 @@ namespace Opm break; case FastAMG: #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_); #else if(linsolver_verbosity_) std::cerr<<"Fast AMG is not available; falling back to CG preconditioned with the normal one"< + struct PreconditionerTraits + { + typedef std::shared_ptr

PointerType; + }; + + template + typename PreconditionerTraits::PointerType + makePreconditioner(O& opA, double relax, const C& comm, int iterations=1) + { + typename Dune::Amg::SmootherTraits

::Arguments args; + typename Dune::Amg::ConstructionTraits

::Arguments cargs; + cargs.setMatrix(opA.getmat()); + args.iterations=iterations; + args.relaxationFactor=relax; + cargs.setArgs(args); + cargs.setComm(comm); + return std::shared_ptr

(Dune::Amg::ConstructionTraits

::construct(cargs)); + } + + template + 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. - Dune::SeqILU0 precond(A, 1.0); + typedef Dune::SeqILU0 Preconditioner; + auto precond = makePreconditioner(opA, 1.0, comm); // Construct linear solver. - Dune::CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::CGSolver linsolve(opA, sp, *precond, tolerance, maxit, verbosity); // Solve system. Dune::InverseOperatorResult result; @@ -266,8 +302,9 @@ namespace Opm criterion.setGamma(1); // V-cycle; this is the default } + template 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) { // Solve with AMG solver. @@ -295,13 +332,12 @@ namespace Opm // Construct preconditioner. Criterion criterion; Precond::SmootherArgs smootherArgs; - Operator opA(A); setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, linsolver_smooth_steps); Precond precond(opA, criterion, smootherArgs); // Construct linear solver. - Dune::CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::CGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity); // Solve system. Dune::InverseOperatorResult result; @@ -317,8 +353,9 @@ namespace Opm #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) + template 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) { // Solve with AMG solver. @@ -344,7 +381,6 @@ namespace Opm typedef Dune::Amg::KAMG Precond; // Construct preconditioner. - Operator opA(A); Precond::SmootherArgs smootherArgs; Criterion criterion; setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, @@ -352,7 +388,7 @@ namespace Opm Precond precond(opA, criterion, smootherArgs); // Construct linear solver. - Dune::GeneralizedPCGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::GeneralizedPCGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity); // Solve system. Dune::InverseOperatorResult result; @@ -366,8 +402,9 @@ namespace Opm return res; } + template 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) { // Solve with AMG solver. @@ -388,7 +425,6 @@ namespace Opm typedef Dune::Amg::FastAMG Precond; // Construct preconditioner. - Operator opA(A); Criterion criterion; const int smooth_steps = 1; setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, smooth_steps); @@ -400,7 +436,7 @@ namespace Opm Precond precond(opA, criterion, parms); // Construct linear solver. - Dune::GeneralizedPCGSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::GeneralizedPCGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity); // Solve system. Dune::InverseOperatorResult result; @@ -415,17 +451,17 @@ namespace Opm } #endif - + template 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. - Dune::SeqILU0 precond(A, 1.0); + typedef Dune::SeqILU0 Preconditioner; + auto precond = makePreconditioner(opA, 1.0, comm); // Construct linear solver. - Dune::BiCGSTABSolver linsolve(opA, precond, tolerance, maxit, verbosity); + Dune::BiCGSTABSolver linsolve(opA, sp, *precond, tolerance, maxit, verbosity); // Solve system. Dune::InverseOperatorResult result; diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp index da8e9ab39..df7b681e9 100644 --- a/opm/core/linalg/LinearSolverIstl.hpp +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -86,6 +86,19 @@ namespace Opm virtual double getTolerance() const; 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 + LinearSolverReport solveSystem(O& opA, V& x, double* solution, + V& b, S& sp, const C& comm, int maxit) const; + double linsolver_residual_tolerance_; int linsolver_verbosity_; enum LinsolverType { CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2, FastAMG=3, KAMG=4 };