/* Copyright 2012 SINTEF ICT, Applied Mathematics. This file is part of the Open Porous Media project (OPM). OPM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OPM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OPM. If not, see . */ #if HAVE_CONFIG_H #include "config.h" #endif #include // Work around the fact that istl headers expect // HAVE_BOOST to be 1, and not just defined. #undef HAVE_BOOST #define HAVE_BOOST 1 // TODO: clean up includes. #include #include #include #include #include #include #include #include #include namespace Opm { using namespace Dune; // While not great, it's okay in a cpp file like this. namespace { typedef FieldVector VectorBlockType; typedef FieldMatrix MatrixBlockType; typedef BCRSMatrix Mat; typedef BlockVector Vector; typedef MatrixAdapter Operator; LinearSolverInterface::LinearSolverReport solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); LinearSolverInterface::LinearSolverReport solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); LinearSolverInterface::LinearSolverReport solveBiCGStab_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity); } // anonymous namespace LinearSolverIstl::LinearSolverIstl() : linsolver_residual_tolerance_(1e-8), linsolver_verbosity_(0), linsolver_type_(CG_AMG), linsolver_save_system_(false), linsolver_max_iterations_(0) { } LinearSolverIstl::LinearSolverIstl(const parameter::ParameterGroup& param) : linsolver_residual_tolerance_(1e-8), linsolver_verbosity_(0), linsolver_type_(CG_AMG), linsolver_save_system_(false), linsolver_max_iterations_(0) { linsolver_residual_tolerance_ = param.getDefault("linsolver_residual_tolerance", linsolver_residual_tolerance_); linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_); linsolver_type_ = LinsolverType(param.getDefault("linsolver_type", int(linsolver_type_))); linsolver_save_system_ = param.getDefault("linsolver_save_system", linsolver_save_system_); if (linsolver_save_system_) { linsolver_save_filename_ = param.getDefault("linsolver_save_filename", std::string("linsys")); } linsolver_max_iterations_ = param.getDefault("linsolver_max_iterations", linsolver_max_iterations_); } LinearSolverIstl::~LinearSolverIstl() { } LinearSolverInterface::LinearSolverReport LinearSolverIstl::solve(const int size, const int nonzeros, const int* ia, const int* ja, const double* sa, const double* rhs, double* solution) const { // Build Istl structures from input. // System matrix Mat A(size, size, nonzeros, Mat::row_wise); for (Mat::CreateIterator row = A.createbegin(); row != A.createend(); ++row) { int ri = row.index(); for (int i = ia[ri]; i < ia[ri + 1]; ++i) { row.insert(ja[i]); } } for (int ri = 0; ri < size; ++ri) { for (int i = ia[ri]; i < ia[ri + 1]; ++i) { 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) { maxit = A.N(); } LinearSolverReport res; switch (linsolver_type_) { case CG_ILU0: res = solveCG_ILU0(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); break; case CG_AMG: res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); break; case BiCGStab_ILU0: res = solveBiCGStab_ILU0(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_); break; default: std::cerr << "Unknown linsolver_type: " << int(linsolver_type_) << '\n'; throw std::runtime_error("Unknown linsolver_type"); } std::copy(x.begin(), x.end(), solution); return res; } void LinearSolverIstl::setTolerance(const double tol) { linsolver_residual_tolerance_ = tol; } double LinearSolverIstl::getTolerance() const { return linsolver_residual_tolerance_; } namespace { LinearSolverInterface::LinearSolverReport solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) { Operator opA(A); // Construct preconditioner. SeqILU0 precond(A, 1.0); // Construct linear solver. CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); // Solve system. InverseOperatorResult result; linsolve.apply(x, b, result); // Output results. LinearSolverInterface::LinearSolverReport res; res.converged = result.converged; res.iterations = result.iterations; res.residual_reduction = result.reduction; return res; } LinearSolverInterface::LinearSolverReport solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) { // Solve with AMG solver. #define FIRST_DIAGONAL 1 #define SYMMETRIC 1 #define SMOOTHER_ILU 1 #define ANISOTROPIC_3D 0 #if FIRST_DIAGONAL typedef Amg::FirstDiagonal CouplingMetric; #else typedef Amg::RowSum CouplingMetric; #endif #if SYMMETRIC typedef Amg::SymmetricCriterion CriterionBase; #else typedef Amg::UnSymmetricCriterion CriterionBase; #endif #if SMOOTHER_ILU typedef SeqILU0 Smoother; #else typedef SeqSSOR Smoother; #endif typedef Amg::CoarsenCriterion Criterion; typedef Amg::AMG Precond; Operator opA(A); // Construct preconditioner. double relax = 1; Precond::SmootherArgs smootherArgs; smootherArgs.relaxationFactor = relax; Criterion criterion; criterion.setDebugLevel(verbosity); #if ANISOTROPIC_3D criterion.setDefaultValuesAnisotropic(3, 2); #endif Precond precond(opA, criterion, smootherArgs); // Construct linear solver. CGSolver linsolve(opA, precond, tolerance, maxit, verbosity); // Solve system. InverseOperatorResult result; linsolve.apply(x, b, result); // Output results. LinearSolverInterface::LinearSolverReport res; res.converged = result.converged; res.iterations = result.iterations; res.residual_reduction = result.reduction; return res; } LinearSolverInterface::LinearSolverReport solveBiCGStab_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity) { Operator opA(A); // Construct preconditioner. SeqILU0 precond(A, 1.0); // Construct linear solver. BiCGSTABSolver linsolve(opA, precond, tolerance, maxit, verbosity); // Solve system. InverseOperatorResult result; linsolve.apply(x, b, result); // Output results. LinearSolverInterface::LinearSolverReport res; res.converged = result.converged; res.iterations = result.iterations; res.residual_reduction = result.reduction; return res; } } // anonymous namespace } // namespace Opm