Add initial skeleton of CPR preconditioning class.

Implementation is identical to NewtonIterationBlackoilSimple for now.

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -0,0 +1,70 @@
+/*
+  Copyright 2014 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 <config.h>
+
+#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/linalg/LinearSolverUmfpack.hpp>
+
+namespace Opm
+{
+
+    /// Construct a system solver.
+    /// \param[in] linsolver   linear solver to use
+    NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
+    {
+        linsolver_full_.reset(new LinearSolverUmfpack);
+    }
+
+    /// Solve the linear system Ax = b, with A being the
+    /// combined derivative matrix of the residual and b
+    /// being the residual itself.
+    /// \param[in] residual   residual object containing A and b.
+    /// \return               the solution x
+    NewtonIterationBlackoilCPR::SolutionVector
+    NewtonIterationBlackoilCPR::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
+    {
+        typedef LinearisedBlackoilResidual::ADB ADB;
+        const int np = residual.material_balance_eq.size();
+        ADB mass_res = residual.material_balance_eq[0];
+        for (int phase = 1; phase < np; ++phase) {
+            mass_res = vertcat(mass_res, residual.material_balance_eq[phase]);
+        }
+        const ADB well_res = vertcat(residual.well_flux_eq, residual.well_eq);
+        const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));
+
+        const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual.derivative()[0];
+
+        SolutionVector dx(SolutionVector::Zero(total_residual.size()));
+        Opm::LinearSolverInterface::LinearSolverReport rep
+            = linsolver_full_->solve(matr.rows(), matr.nonZeros(),
+                                     matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
+                                     total_residual.value().data(), dx.data());
+        if (!rep.converged) {
+            OPM_THROW(std::runtime_error,
+                      "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
+                      "Linear solver convergence failure.");
+        }
+        return dx;
+    }
+
+} // namespace Opm
+

opm/autodiff/NewtonIterationBlackoilCPR.hpp

@@ -0,0 +1,58 @@
+/*
+  Copyright 2014 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_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
+#define OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
+
+
+#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <memory>
+
+namespace Opm
+{
+
+    /// This class solves the fully implicit black-oil system by
+    /// applying a Constrained Pressure Residual preconditioning
+    /// strategy.
+    class NewtonIterationBlackoilCPR : public NewtonIterationBlackoilInterface
+    {
+    public:
+        /// Construct a system solver.
+        /// \param[in] param   parameters controlling the behaviour of
+        ///                    the preconditioning and choice of
+        ///                    linear solvers.
+        NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param);
+
+        /// Solve the system of linear equations Ax = b, with A being the
+        /// combined derivative matrix of the residual and b
+        /// being the residual itself.
+        /// \param[in] residual   residual object containing A and b.
+        /// \return               the solution x
+        virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
+
+    private:
+        std::unique_ptr<LinearSolverInterface> linsolver_elliptic_;
+        std::unique_ptr<LinearSolverInterface> linsolver_full_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED