From e1d0995104c95d34a4f2d956eb870d51bdbccc59 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 19 Mar 2014 17:37:12 +0100 Subject: [PATCH 1/3] 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 }; From 909d3c9019068741651bef9dec829759be79149a Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Mar 2014 09:37:41 +0100 Subject: [PATCH 2/3] Moved ISTL right and left hand side construction for ISTL to solveSystem --- opm/core/linalg/LinearSolverIstl.cpp | 43 ++++++++++++++-------------- opm/core/linalg/LinearSolverIstl.hpp | 10 +++---- 2 files changed, 26 insertions(+), 27 deletions(-) diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index 181b9ed2f..e679dad0a 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -150,24 +150,6 @@ namespace Opm A[ri][ja[i]] = sa[i]; } } - // System RHS - Vector b(size); - std::copy(rhs, rhs + size, b.begin()); - // System solution - Vector x(size); - x = 0.0; - - if (linsolver_save_system_) - { - // Save system to files. - writeMatrixToMatlab(A, linsolver_save_filename_ + "-mat"); - std::string rhsfile(linsolver_save_filename_ + "-rhs"); - std::ofstream rhsf(rhsfile.c_str()); - rhsf.precision(15); - rhsf.setf(std::ios::scientific | std::ios::showpos); - std::copy(b.begin(), b.end(), - std::ostream_iterator(rhsf, "\n")); - } int maxit = linsolver_max_iterations_; if (maxit == 0) { @@ -176,14 +158,33 @@ namespace Opm Dune::SeqScalarProduct sp; Dune::Amg::SequentialInformation comm; Operator opA(A); - return solveSystem(opA, x, solution, b, sp, comm, maxit); + return solveSystem(opA, solution, rhs, sp, comm, maxit); } -template + template LinearSolverInterface::LinearSolverReport - LinearSolverIstl::solveSystem (O& opA, V& x, double* solution, V& b, + 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 + size, b.begin()); + // System solution + Vector x(opA.getmat().M()); + x = 0.0; + + if (linsolver_save_system_) + { + // Save system to files. + writeMatrixToMatlab(opA.getmat(), linsolver_save_filename_ + "-mat"); + std::string rhsfile(linsolver_save_filename_ + "-rhs"); + std::ofstream rhsf(rhsfile.c_str()); + rhsf.precision(15); + rhsf.setf(std::ios::scientific | std::ios::showpos); + std::copy(b.begin(), b.end(), + std::ostream_iterator(rhsf, "\n")); + } + LinearSolverReport res; switch (linsolver_type_) { case CG_ILU0: diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp index df7b681e9..4ab853cc7 100644 --- a/opm/core/linalg/LinearSolverIstl.hpp +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -88,16 +88,14 @@ namespace Opm 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] 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, V& x, double* solution, - V& b, S& sp, const C& comm, int maxit) const; + 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_; From 7495bfa6c4d97c90bc59a9ae9798a82435622013 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 26 Mar 2014 10:23:11 +0100 Subject: [PATCH 3/3] [bugfix] Use the size of the vector for the copying. The last patch did not compile as there was no size member in scope. Therefore this patch resorts to using the size of the vector. --- opm/core/linalg/LinearSolverIstl.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index e679dad0a..1871b8f8d 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -168,7 +168,7 @@ namespace Opm { // System RHS Vector b(opA.getmat().N()); - std::copy(rhs, rhs + size, b.begin()); + std::copy(rhs, rhs+b.size(), b.begin()); // System solution Vector x(opA.getmat().M()); x = 0.0;