Merge branch 'cpr-preconditioning' into cpr-preconditioning-again

Conflicts:
	examples/sim_fibo_ad.cpp
	examples/test_implicit_ad.cpp
	opm/autodiff/FullyImplicitBlackoilSolver.cpp
	opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
	opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

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
@@ -100,7 +101,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

examples/sim_fibo_ad.cpp

@@ -37,6 +37,7 @@
 #include <opm/core/props/rock/RockCompressibility.hpp>
 
 #include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/autodiff/FullyImplicitSystemSolverSimple.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@@ -141,6 +142,7 @@ try
 
     // Linear solver.
     LinearSolverFactory linsolver(param);
+    FullyImplicitSystemSolverSimple fis_solver(linsolver);
 
     // Write parameters used for later reference.
     bool output = param.getDefault("output", true);
@@ -210,7 +212,7 @@ try
                                                  *new_props,
                                                  rock_comp->isActive() ? rock_comp.get() : 0,
                                                  wells,
-                                                 linsolver,
+                                                 fis_solver,
                                                  grav);
         SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
 

examples/test_implicit_ad.cpp

@@ -24,6 +24,7 @@
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/BlackoilPropsAd.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/autodiff/FullyImplicitSystemSolverSimple.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -103,8 +104,9 @@ try
     Opm::DerivedGeology geo(*g, props, grav);
 
     Opm::LinearSolverFactory linsolver(param);
+    Opm::FullyImplicitSystemSolverSimple fis_solver(linsolver);
 
-    Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, linsolver);
+    Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, fis_solver);
 
     Opm::BlackoilState state;
     initStateBasic(*g, props0, param, 0.0, state);

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/FullyImplicitBlackoilSolver.cpp

@@ -195,7 +195,7 @@ namespace {
                                 const DerivedGeology&           geo  ,
                                 const RockCompressibility*      rock_comp_props,
                                 const Wells&                    wells,
-                                const LinearSolverInterface&    linsolver)
+                                const FullyImplicitSystemSolverInterface&    linsolver)
         : grid_  (grid)
         , fluid_ (fluid)
         , geo_   (geo)
@@ -643,13 +643,13 @@ namespace {
             // std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl;
             // std::cout << rq_[phase].mflux;
 
-            residual_.mass_balance[ phaseIdx ] =
+            residual_.material_balance_eq[ phaseIdx ] =
                 pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
                 + ops_.div*rq_[phaseIdx].mflux;
 
 
             // DUMP(ops_.div*rq_[phase].mflux);
-            // DUMP(residual_.mass_balance[phase]);
+            // DUMP(residual_.material_balance_eq[phase]);
         }
 
         // -------- Extra (optional) rs and rv contributions to the mass balance equations --------
@@ -663,16 +663,16 @@ namespace {
                                                 rq_[po].head.value());
             const ADB rs_face = upwindOil.select(state.rs);
 
-            residual_.mass_balance[ Gas ] += ops_.div * (rs_face * rq_[po].mflux);
+            residual_.material_balance_eq[ Gas ] += ops_.div * (rs_face * rq_[po].mflux);
 
             const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
             const UpwindSelector<double> upwindGas(grid_, ops_,
                                                 rq_[pg].head.value());
             const ADB rv_face = upwindGas.select(state.rv);
 
-            residual_.mass_balance[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux);
+            residual_.material_balance_eq[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux);
 
-            // DUMP(residual_.mass_balance[ Gas ]);
+            // DUMP(residual_.material_balance_eq[ Gas ]);
 
         }
 
@@ -850,7 +850,7 @@ namespace {
 
         // Add well contributions to mass balance equations
         for (int phase = 0; phase < np; ++phase) {
-            residual_.mass_balance[phase] -= superset(cq_s[phase],well_cells,nc);
+            residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc);
         }
 
         // Add WELL EQUATIONS
@@ -1133,7 +1133,7 @@ namespace {
             // const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
             well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc);
             // DUMP(well_contribs[phase]);
-            residual_.mass_balance[phase] += well_contribs[phase];
+            residual_.material_balance_eq[phase] += well_contribs[phase];
         }
         if (active_[Gas] && active_[Oil]) {
             const int oilpos = pu.phase_pos[Oil];
@@ -1143,8 +1143,8 @@ namespace {
             well_rates_all += superset(wops_.p2w * (well_perf_rates[oilpos]*rs_perf), Span(nw, 1, gaspos*nw), nw*np);
             well_rates_all += superset(wops_.p2w * (well_perf_rates[gaspos]*rv_perf), Span(nw, 1, oilpos*nw), nw*np);
             // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs);
-            residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs;
+            residual_.material_balance_eq[gaspos] += well_contribs[oilpos]*state.rs;
-            residual_.mass_balance[oilpos] += well_contribs[gaspos]*state.rv;
+            residual_.material_balance_eq[oilpos] += well_contribs[gaspos]*state.rv;
         }
 
         // Set the well flux equation
@@ -1158,40 +1158,7 @@ namespace {
 
     V FullyImplicitBlackoilSolver::solveJacobianSystem() const
     {
-        const int np = fluid_.numPhases();
+        return linsolver_.linearSolve(residual_);
-        ADB mass_res = residual_.mass_balance[0];
-        for (int phase = 1; phase < np; ++phase) {
-            mass_res = vertcat(mass_res, residual_.mass_balance[phase]);
-        }
-        const ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq);
-        const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));
-        // DUMP(total_residual);
-
-        const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual.derivative()[0];
-
-        V dx(V::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());
-        /*
-        std::ofstream filestream("matrix.out");
-        filestream << matr;
-        filestream.close();
-        std::ofstream filestream2("sol.out");
-        filestream2 << dx;
-        filestream2.close();
-        std::ofstream filestream3("r.out");
-        filestream3 << total_residual.value();
-        filestream3.close(); */
-
-
-        if (!rep.converged) {
-            OPM_THROW(std::runtime_error,
-                      "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
-                      "Linear solver convergence failure.");
-        }
-        return dx;
     }
 
 
@@ -1599,10 +1566,9 @@ namespace {
     FullyImplicitBlackoilSolver::residualNorm() const
     {
         double globalNorm = 0;
-        std::vector<ADB>::const_iterator quantityIt = residual_.mass_balance.begin();
+        std::vector<ADB>::const_iterator quantityIt = residual_.material_balance_eq.begin();
-        const std::vector<ADB>::const_iterator endQuantityIt = residual_.mass_balance.end();
+        const std::vector<ADB>::const_iterator endQuantityIt = residual_.material_balance_eq.end();
-        for (; quantityIt != endQuantityIt; ++quantityIt)
+        for (; quantityIt != endQuantityIt; ++quantityIt) {
-        {
             const double quantityResid = (*quantityIt).value().matrix().norm();
             if (!std::isfinite(quantityResid)) {
                 OPM_THROW(Opm::NumericalProblem,

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -23,6 +23,8 @@
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/FullyImplicitBlackoilResidual.hpp>
+#include <opm/autodiff/FullyImplicitSystemSolverInterface.hpp>
 
 struct UnstructuredGrid;
 struct Wells;
@@ -31,7 +33,7 @@ namespace Opm {
 
     class DerivedGeology;
     class RockCompressibility;
-    class LinearSolverInterface;
+    class FullyImplicitSystemSolverInterface;
     class BlackoilState;
     class WellStateFullyImplicitBlackoil;
 
@@ -62,7 +64,7 @@ namespace Opm {
                                     const DerivedGeology&           geo  ,
                                     const RockCompressibility*      rock_comp_props,
                                     const Wells&                    wells,
-                                    const LinearSolverInterface&    linsolver);
+                                    const FullyImplicitSystemSolverInterface& linsolver);
 
         /// Take a single forward step, modifiying
         ///   state.pressure()
@@ -123,7 +125,7 @@ namespace Opm {
         const DerivedGeology&           geo_;
         const RockCompressibility*      rock_comp_props_;
         const Wells&                    wells_;
-        const LinearSolverInterface&    linsolver_;
+        const FullyImplicitSystemSolverInterface&    linsolver_;
         // For each canonical phase -> true if active
         const std::vector<bool>         active_;
         // Size = # active faces. Maps active -> canonical phase indices.
@@ -137,14 +139,7 @@ namespace Opm {
         std::vector<PhasePresence> phaseCondition_;
         V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
 
-        // The mass_balance vector has one element for each active phase,
+        FullyImplicitBlackoilResidual residual_;
-        // each of which has size equal to the number of cells.
-        // The well_eq has size equal to the number of wells.
-        struct {
-            std::vector<ADB> mass_balance;
-            ADB well_flux_eq;
-            ADB well_eq;
-        } residual_;
 
         // Private methods.
         SolutionState

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

opm/autodiff/SimulatorFullyImplicitBlackoil.cpp

@@ -70,7 +70,7 @@ namespace Opm
              BlackoilPropsAdInterface& props,
              const RockCompressibility* rock_comp_props,
              WellsManager& wells_manager,
-             LinearSolverInterface& linsolver,
+             FullyImplicitSystemSolverInterface& linsolver,
              const double* gravity);
 
         SimulatorReport run(SimulatorTimer& timer,
@@ -110,7 +110,7 @@ namespace Opm
                                                                    BlackoilPropsAdInterface& props,
                                                                    const RockCompressibility* rock_comp_props,
                                                                    WellsManager& wells_manager,
-                                                                   LinearSolverInterface& linsolver,
+                                                                   FullyImplicitSystemSolverInterface& linsolver,
                                                                    const double* gravity)
 
     {
@@ -257,7 +257,7 @@ namespace Opm
                                                BlackoilPropsAdInterface& props,
                                                const RockCompressibility* rock_comp_props,
                                                WellsManager& wells_manager,
-                                               LinearSolverInterface& linsolver,
+                                               FullyImplicitSystemSolverInterface& linsolver,
                                                const double* gravity)
         : grid_(grid),
           props_(props),

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -33,7 +33,7 @@ namespace Opm
     class BlackoilPropsAdInterface;
     class RockCompressibility;
     class WellsManager;
-    class LinearSolverInterface;
+    class FullyImplicitSystemSolverInterface;
     class SimulatorTimer;
     class BlackoilState;
     class WellStateFullyImplicitBlackoil;
@@ -70,7 +70,7 @@ namespace Opm
                                        BlackoilPropsAdInterface& props,
                                        const RockCompressibility* rock_comp_props,
                                        WellsManager& wells_manager,
-                                       LinearSolverInterface& linsolver,
+                                       FullyImplicitSystemSolverInterface& linsolver,
                                        const double* gravity);
 
         /// Run the simulation.