Added LinearSolverInterface and two subclasses, using Umfpack and Istl.

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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/linalg/sparse_sys.h>
+#include <opm/core/linalg/call_umfpack.h>
+
+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
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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

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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <opm/core/linalg/LinearSolverIstl.hpp>
+
+
+// TODO: clean up includes.
+#define DUNE_DEPRECATED
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/io.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/paamg/amg.hh>
+
+#include <stdexcept>
+
+
+namespace Opm
+{
+
+    using namespace Dune; // While not great, it's okay in a cpp file like this.
+
+    namespace {
+        typedef FieldVector<double, 1   > VectorBlockType;
+        typedef FieldMatrix<double, 1, 1> MatrixBlockType;
+        typedef BCRSMatrix <MatrixBlockType>        Mat;
+        typedef BlockVector<VectorBlockType>        Vector;
+        typedef MatrixAdapter<Mat,Vector,Vector> 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<VectorBlockType>(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<Mat,Vector,Vector> precond(A, 1.0);
+
+        // Construct linear solver.
+        CGSolver<Vector> 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<Mat,CouplingMetric>   CriterionBase;
+#else
+        typedef Amg::UnSymmetricCriterion<Mat,CouplingMetric> CriterionBase;
+#endif
+
+#if SMOOTHER_ILU
+        typedef SeqILU0<Mat,Vector,Vector>        Smoother;
+#else
+        typedef SeqSSOR<Mat,Vector,Vector>        Smoother;
+#endif
+        typedef Amg::CoarsenCriterion<CriterionBase> Criterion;
+        typedef Amg::AMG<Operator,Vector,Smoother>   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<Vector> 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<Mat,Vector,Vector> precond(A, 1.0);
+
+        // Construct linear solver.
+        BiCGSTABSolver<Vector> 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
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_LINEARSOLVERISTL_HEADER_INCLUDED
+#define OPM_LINEARSOLVERISTL_HEADER_INCLUDED
+
+
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <string>
+
+
+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       <empty string>
+        ///   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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/core/linalg/LinearSolverUmfpack.hpp>
+#include <opm/core/linalg/sparse_sys.h>
+#include <opm/core/linalg/call_umfpack.h>
+
+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<int*>(ia),
+            const_cast<int*>(ja),
+            const_cast<double*>(sa)
+        };
+        call_UMFPACK(&A, rhs, solution);
+    }
+
+} // namespace Opm
+

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_LINEARSOLVERISTL_HEADER_INCLUDED
+#define OPM_LINEARSOLVERISTL_HEADER_INCLUDED
+
+
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+
+
+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