Interface and simple implementation of fully implicit system solver.

This is done in preparation for adding a cpr-preconditioning solver
for the fully implicit black-oil system. The existing implementation
that concatenates the whole system and passes it to some linear solver
has been moved from a private function of FullyImplicitBlackoilSolver
to the class FullyImplicitSolverSimple.

To enable this decoupling, the residual struct has been copied out
of the FullyImplicitBlackoilSolver class and is now an independent
struct: FullyImplicitBlackoilResidual. The opportunity has been used
to replace the field mass_balance with material_balance_eq, which is
more precise.

CMakeLists_files.cmake

@@ -29,6 +29,7 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/BlackoilPropsAd.cpp
 	opm/autodiff/BlackoilPropsAdInterface.cpp
 	opm/autodiff/FullyImplicitBlackoilSolver.cpp
+	opm/autodiff/FullyImplicitSystemSolverSimple.cpp
 	opm/autodiff/ImpesTPFAAD.cpp
 	opm/autodiff/SimulatorCompressibleAd.cpp
 	opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
@@ -83,7 +84,10 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/GeoProps.hpp
 	opm/autodiff/ImpesTPFAAD.hpp
+	opm/autodiff/FullyImplicitBlackoilResidual.hpp
 	opm/autodiff/FullyImplicitBlackoilSolver.hpp
+	opm/autodiff/FullyImplicitBlackoilSystemSolverInterface.hpp
+	opm/autodiff/FullyImplicitBlackoilSystemSolverSimple.hpp
 	opm/autodiff/SimulatorCompressibleAd.hpp
 	opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
 	opm/autodiff/SimulatorIncompTwophaseAd.hpp

opm/autodiff/FullyImplicitBlackoilResidual.hpp

@@ -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/>.
+*/
+
+#ifndef OPM_FULLYIMPLICITBLACKOILRESIDUAL_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITBLACKOILRESIDUAL_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+
+namespace Opm
+{
+
+    /// Residual structure of the fully implicit solver.
+    /// All equations are given as AD types, with multiple
+    /// jacobian blocks corresponding to the primary unknowns. The
+    /// primary unknowns are for a three-phase simulation, in order:
+    ///    p    (pressure)
+    ///    sw   (water saturation)
+    ///    xvar (gas saturation, gas-oil ratio or oil-gas ratio)
+    ///    qs   (well outflows by well and phase)
+    ///    bhp  (bottom hole pressures)
+    /// In the above, the xvar variable will have a different
+    /// meaning from cell to cell, corresponding to the state in
+    /// that cell (saturated, undersaturated oil or undersaturated
+    /// gas). In a two-phase simulation, either sw or xvar is not
+    /// used, depending on which face is missing.
+    ///
+    /// Note: this class is strongly coupled to the class
+    /// FullyImplicitBlackoilSolver, and is separated from that
+    /// class to facilitate the development of linear solver
+    /// strategies outside that class.
+    struct FullyImplicitBlackoilResidual {
+        /// A type alias for the automatic differentiation type.
+        typedef AutoDiffBlock<double> ADB;
+        /// The material_balance_eq vector has one element for each
+        /// active phase, each of which has size equal to the number
+        /// of cells. Each material balance equation is given in terms
+        /// of surface volumes.
+        std::vector<ADB> material_balance_eq;
+        /// The well_flux_eq has size equal to the number of wells
+        /// times the number of phases. It contains the well flow
+        /// equations, relating the total well flows to
+        /// bottom-hole pressures and reservoir conditions.
+        ADB well_flux_eq;
+        /// The well_eq has size equal to the number of wells. It
+        /// contains the well control equations, that is for each
+        /// well either a rate specification or bottom hole
+        /// pressure specification.
+        ADB well_eq;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_FULLYIMPLICITBLACKOILRESIDUAL_HEADER_INCLUDED

opm/autodiff/FullyImplicitSystemSolverInterface.hpp

@@ -0,0 +1,46 @@
+/*
+  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_FULLYIMPLICITSYSTEMSOLVERINTERFACE_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITSYSTEMSOLVERINTERFACE_HEADER_INCLUDED
+
+#include <opm/autodiff/FullyImplicitBlackoilResidual.hpp>
+
+namespace Opm
+{
+
+    /// Interface class for (linear) solvers for the fully implicit black-oil system.
+    class FullyImplicitSystemSolverInterface
+    {
+    public:
+        /// Return type for linearSolve(). A simple, non-ad vector type.
+        typedef FullyImplicitBlackoilResidual::ADB::V SolutionVector;
+
+        /// 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
+        virtual SolutionVector linearSolve(const FullyImplicitBlackoilResidual& residual) const = 0;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_FULLYIMPLICITSYSTEMSOLVERINTERFACE_HEADER_INCLUDED

opm/autodiff/FullyImplicitSystemSolverSimple.cpp

@@ -0,0 +1,68 @@
+/*
+  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 <opm/autodiff/FullyImplicitSystemSolverSimple.hpp>
+#include <opm/autodiff/AutodiffHelpers.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+namespace Opm
+{
+
+    /// Construct a system solver.
+    /// \param[in] linsolver   linear solver to use
+    FullyImplicitSystemSolverSimple::FullyImplicitSystemSolverSimple(const LinearSolverInterface& linsolver)
+        : linsolver_(linsolver)
+    {
+    }
+
+    /// 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
+    FullyImplicitSystemSolverSimple::SolutionVector
+    FullyImplicitSystemSolverSimple::linearSolve(const FullyImplicitBlackoilResidual& residual) const
+    {
+        typedef FullyImplicitBlackoilResidual::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_.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/FullyImplicitSystemSolverSimple.hpp

@@ -0,0 +1,55 @@
+/*
+  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_FULLYIMPLICITSYSTEMSOLVERSIMPLE_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITSYSTEMSOLVERSIMPLE_HEADER_INCLUDED
+
+
+#include <opm/autodiff/FullyImplicitSystemSolverInterface.hpp>
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+
+namespace Opm
+{
+
+    /// This class solves the fully implicit black-oil system by
+    /// simply concatenating the Jacobian matrices and passing the
+    /// resulting system to a linear solver. The linear solver used
+    /// can be passed in as a constructor argument.
+    class FullyImplicitSystemSolverSimple : public FullyImplicitSystemSolverInterface
+    {
+    public:
+        /// Construct a system solver.
+        /// \param[in] linsolver   linear solver to use
+        FullyImplicitSystemSolverSimple(const LinearSolverInterface& linsolver);
+
+        /// 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
+        virtual SolutionVector linearSolve(const FullyImplicitBlackoilResidual& residual) const;
+
+    private:
+        const LinearSolverInterface& linsolver_;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_FULLYIMPLICITSYSTEMSOLVERSIMPLE_HEADER_INCLUDED