diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index c32ad02ee..1871b8f8d 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 @@ -144,17 +150,33 @@ namespace Opm A[ri][ja[i]] = sa[i]; } } - // System RHS - Vector b(size); - std::copy(rhs, rhs + size, b.begin()); + + int maxit = linsolver_max_iterations_; + if (maxit == 0) { + maxit = 5000; + } + Dune::SeqScalarProduct sp; + Dune::Amg::SequentialInformation comm; + Operator opA(A); + return solveSystem(opA, solution, rhs, sp, comm, maxit); + } + + template + LinearSolverInterface::LinearSolverReport + LinearSolverIstl::solveSystem (O& opA, double* solution, const double* rhs, + S& sp, const C& comm, int maxit) const + { + // System RHS + Vector b(opA.getmat().N()); + std::copy(rhs, rhs+b.size(), b.begin()); // System solution - Vector x(size); + Vector x(opA.getmat().M()); x = 0.0; if (linsolver_save_system_) { // Save system to files. - writeMatrixToMatlab(A, linsolver_save_filename_ + "-mat"); + writeMatrixToMatlab(opA.getmat(), linsolver_save_filename_ + "-mat"); std::string rhsfile(linsolver_save_filename_ + "-rhs"); std::ofstream rhsf(rhsfile.c_str()); rhsf.precision(15); @@ -163,23 +185,18 @@ namespace Opm std::ostream_iterator(rhsf, "\n")); } - int maxit = linsolver_max_iterations_; - if (maxit == 0) { - maxit = 5000; - } - 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 +204,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 +303,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 +333,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 +354,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 +382,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 +389,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 +403,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 +426,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 +437,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 +452,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..4ab853cc7 100644 --- a/opm/core/linalg/LinearSolverIstl.hpp +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -86,6 +86,17 @@ 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[out] solution C array for storing the solution vector. + /// \param[in] rhs C array containing the right hand side. + /// \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, double* solution, const double *rhs, + 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 };