diff --git a/opm/core/linalg/LinearSolverInterface.cpp b/opm/core/linalg/LinearSolverInterface.cpp new file mode 100644 index 000000000..9e2d0bbfd --- /dev/null +++ b/opm/core/linalg/LinearSolverInterface.cpp @@ -0,0 +1,44 @@ +/* + 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 . +*/ + + +#include +#include +#include + +namespace Opm +{ + + LinearSolverInterface::~LinearSolverInterface() + { + } + + + + + LinearSolverInterface::LinearSolverReport + LinearSolverInterface::solve(const CSRMatrix* A, + const double* rhs, + double* solution) + { + return solve(A->m, A->nnz, A->ia, A->ja, A->sa, rhs, solution); + } + +} // namespace Opm + diff --git a/opm/core/linalg/LinearSolverInterface.hpp b/opm/core/linalg/LinearSolverInterface.hpp new file mode 100644 index 000000000..5835f11e3 --- /dev/null +++ b/opm/core/linalg/LinearSolverInterface.hpp @@ -0,0 +1,81 @@ +/* + 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 . +*/ + +#ifndef OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED +#define OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED + + +struct CSRMatrix; + +namespace Opm +{ + + + /// Abstract interface for linear solvers. + class LinearSolverInterface + { + public: + /// Virtual destructor. + virtual ~LinearSolverInterface(); + + /// Struct for reporting data about the solution process back + /// to the caller. The only field that is mandatory to set is + /// 'converged' (even for direct solvers) to indicate success. + struct LinearSolverReport + { + bool converged; + int iterations; + double residual_reduction; + }; + + /// Solve a linear system, with a matrix given in compressed sparse row format. + /// \param[in] A matrix in CSR format + /// \param[in] rhs array of length A->m containing the right hand side + /// \param[inout] solution array of length A->m to which the solution will be written, may also be used + /// as initial guess by iterative solvers. + /// Note: this method is a convenience method that calls the virtual solve() method. + LinearSolverReport solve(const CSRMatrix* A, + const double* rhs, + double* solution); + + /// Solve a linear system, with a matrix given in compressed sparse row format. + /// \param[in] size # of rows in matrix + /// \param[in] nonzeros # of nonzeros elements in matrix + /// \param[in] ia array of length (size + 1) containing start and end indices for each row + /// \param[in] ja array of length nonzeros containing column numbers for the nonzero elements + /// \param[in] sa array of length nonzeros containing the values of the nonzero elements + /// \param[in] rhs array of length size containing the right hand side + /// \param[inout] solution array of length size to which the solution will be written, may also be used + /// as initial guess by iterative solvers. + virtual LinearSolverReport solve(const int size, + const int nonzeros, + const int* ia, + const int* ja, + const double* sa, + const double* rhs, + double* solution) = 0; + + }; + + +} // namespace Opm + + + +#endif // OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp new file mode 100644 index 000000000..d8ac9c229 --- /dev/null +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -0,0 +1,290 @@ +/* + 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 + + +// TODO: clean up includes. +#define DUNE_DEPRECATED +#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) + { + // 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; + } + + + 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 + diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp new file mode 100644 index 000000000..1e883c285 --- /dev/null +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -0,0 +1,89 @@ +/* + 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 . +*/ + +#ifndef OPM_LINEARSOLVERISTL_HEADER_INCLUDED +#define OPM_LINEARSOLVERISTL_HEADER_INCLUDED + + +#include +#include +#include + + +namespace Opm +{ + + + /// Abstract interface for linear solvers. + class LinearSolverIstl : public LinearSolverInterface + { + public: + /// Default constructor. + /// All parameters controlling the solver are defaulted: + /// linsolver_residual_tolerance 1e-8 + /// linsolver_verbosity 0 + /// linsolver_type 1 ( = CG_AMG), alternatives are: + /// CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2 + /// linsolver_save_system false + /// linsolver_save_filename + /// linsolver_max_iterations 0 (unlimited) + LinearSolverIstl(); + + /// Construct from parameters + /// Accepted parameters are, with defaults, listed in the + /// default constructor. + LinearSolverIstl(const parameter::ParameterGroup& param); + + /// Destructor. + virtual ~LinearSolverIstl(); + + using LinearSolverInterface::solve; + + /// Solve a linear system, with a matrix given in compressed sparse row format. + /// \param[in] size # of rows in matrix + /// \param[in] nonzeros # of nonzeros elements in matrix + /// \param[in] ia array of length (size + 1) containing start and end indices for each row + /// \param[in] ja array of length nonzeros containing column numbers for the nonzero elements + /// \param[in] sa array of length nonzeros containing the values of the nonzero elements + /// \param[in] rhs array of length size containing the right hand side + /// \param[inout] solution array of length size to which the solution will be written, may also be used + /// as initial guess by iterative solvers. + virtual LinearSolverReport solve(const int size, + const int nonzeros, + const int* ia, + const int* ja, + const double* sa, + const double* rhs, + double* solution); + private: + double linsolver_residual_tolerance_; + int linsolver_verbosity_; + enum LinsolverType { CG_ILU0 = 0, CG_AMG = 1, BiCGStab_ILU0 = 2 }; + LinsolverType linsolver_type_; + bool linsolver_save_system_; + std::string linsolver_save_filename_; + int linsolver_max_iterations_; + }; + + +} // namespace Opm + + + +#endif // OPM_LINEARSOLVERISTL_HEADER_INCLUDED diff --git a/opm/core/linalg/LinearSolverUmfpack.cpp b/opm/core/linalg/LinearSolverUmfpack.cpp new file mode 100644 index 000000000..fd3d63bea --- /dev/null +++ b/opm/core/linalg/LinearSolverUmfpack.cpp @@ -0,0 +1,61 @@ +/* + 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 . +*/ + +#include +#include +#include + +namespace Opm +{ + + LinearSolverUmfpack::LinearSolverUmfpack() + { + } + + + + + LinearSolverUmfpack::~LinearSolverUmfpack() + { + } + + + + + LinearSolverInterface::LinearSolverReport + LinearSolverUmfpack::solve(const int size, + const int nonzeros, + const int* ia, + const int* ja, + const double* sa, + const double* rhs, + double* solution) + { + CSRMatrix A = { + size, + nonzeros, + const_cast(ia), + const_cast(ja), + const_cast(sa) + }; + call_UMFPACK(&A, rhs, solution); + } + +} // namespace Opm + diff --git a/opm/core/linalg/LinearSolverUmfpack.hpp b/opm/core/linalg/LinearSolverUmfpack.hpp new file mode 100644 index 000000000..3de9efbc9 --- /dev/null +++ b/opm/core/linalg/LinearSolverUmfpack.hpp @@ -0,0 +1,89 @@ +/* + 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 . +*/ + +#ifndef OPM_LINEARSOLVERUMFPACK_HEADER_INCLUDED +#define OPM_LINEARSOLVERUMFPACK_HEADER_INCLUDED +/* + 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 . +*/ + +#ifndef OPM_LINEARSOLVERISTL_HEADER_INCLUDED +#define OPM_LINEARSOLVERISTL_HEADER_INCLUDED + + +#include + + +namespace Opm +{ + + + /// Abstract interface for linear solvers. + class LinearSolverUmfpack : public LinearSolverInterface + { + public: + /// Default constructor. + LinearSolverUmfpack(); + + /// Destructor. + virtual ~LinearSolverUmfpack(); + + using LinearSolverInterface::solve; + + /// Solve a linear system, with a matrix given in compressed sparse row format. + /// \param[in] size # of rows in matrix + /// \param[in] nonzeros # of nonzeros elements in matrix + /// \param[in] ia array of length (size + 1) containing start and end indices for each row + /// \param[in] ja array of length nonzeros containing column numbers for the nonzero elements + /// \param[in] sa array of length nonzeros containing the values of the nonzero elements + /// \param[in] rhs array of length size containing the right hand side + /// \param[inout] solution array of length size to which the solution will be written, may also be used + /// as initial guess by iterative solvers. + virtual LinearSolverReport solve(const int size, + const int nonzeros, + const int* ia, + const int* ja, + const double* sa, + const double* rhs, + double* solution); + }; + + +} // namespace Opm + + + +#endif // OPM_LINEARSOLVERISTL_HEADER_INCLUDED + +#endif // OPM_LINEARSOLVERUMFPACK_HEADER_INCLUDED