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 };