Initial implementation of CPR preconditioner done.

With some caveats:
  - scaling factors for material balance equations and pressure are hardcoded.
  - pressure system is formed from sum of material balance equations, with
    no check for diagonal dominance.

opm/autodiff/CPRPreconditioner.hpp

@@ -76,7 +76,8 @@ namespace Opm
           \param w  The ILU0 relaxation factor.
         */
         CPRPreconditioner (const M& A, const M& Ae, const field_type relax)
-            : ILU_(A), // copy A (will be overwritten by ILU decomp)
+            : A_(A),
+              ILU_(A), // copy A (will be overwritten by ILU decomp)
               Ae_(Ae),
               relax_(relax)
         {
@@ -103,18 +104,26 @@ namespace Opm
         {
             // Extract part of d corresponding to elliptic part.
             Y de(Ae_.N());
+            // Note: Assumes that the elliptic part comes first.
             std::copy_n(d.begin(), Ae_.N(), de.begin());
 
             // Solve elliptic part, extend solution to full.
-            // "ve = Ae_ \ de;"
-            // "vfull = ve Extend full"
+            Y ve = solveElliptic(de);
+            Y vfull(ILU_.N());
+            vfull = 0.0;
+            // Again assuming that the elliptic part comes first.
+            std::copy(ve.begin(), ve.end(), vfull.begin());
 
             // Subtract elliptic residual from initial residual.
-            // "dmodified = d - A*vfull"
+            // dmodified = d - A * vfull
+            Y dmodified = d;
+            A_.mmv(vfull, dmodified);
 
             // Apply ILU0.
+            Y vilu(ILU_.N());
             Dune::bilu_backsolve(ILU_, vilu, dmodified);
-            // "v = vfull + vilu;"
+            v = vfull;
+            v += vilu;
             v *= relax_;
         }
 
@@ -129,6 +138,40 @@ namespace Opm
         }
 
     private:
+        Y solveElliptic(Y& de)
+        {
+            // std::cout << "solveElliptic()" << std::endl;
+            // Construct operator, scalar product and vectors needed.
+            typedef Dune::MatrixAdapter<M,X,X> Operator;
+            Operator opAe(Ae_);
+            Dune::SeqScalarProduct<X> sp;
+            // Right hand side.
+            // System solution
+            X x(opAe.getmat().M());
+            x = 0.0;
+
+            // Construct preconditioner.
+            typedef typename Dune::SeqILU0<M,X,X> Preconditioner;
+            const double relax = 1.0;
+            Preconditioner precond(Ae_, relax);
+
+            // Construct linear solver.
+            const double tolerance = 1e-4;
+            const int maxit = 5000;
+            const int verbosity = 0;
+            Dune::BiCGSTABSolver<X> linsolve(opAe, sp, precond, tolerance, maxit, verbosity);
+
+            // Solve system.
+            Dune::InverseOperatorResult result;
+            linsolve.apply(x, de, result);
+            if (result.converged) {
+                // std::cout << "solveElliptic() successful!" << std::endl;
+            }
+            return x;
+        }
+
+        //! \brief The matrix for the full linear problem.
+        const matrix_type& A_;
         //! \brief The ILU0 decomposition of the matrix.
         matrix_type ILU_;
         //! \brief The elliptic part of the matrix.

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -23,6 +23,7 @@
 #include <opm/autodiff/CPRPreconditioner.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/Units.hpp>
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 
 #include "disable_warning_pragmas.h"
@@ -133,6 +134,22 @@ namespace Opm
         eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns.
         assert(int(eqs.size()) == np);
 
+        // Scale material balance equations.
+        const double matbalscale[3] = { 1.1169, 1.0031, 0.0031 }; // HACK hardcoded instead of computed.
+        for (int phase = 0; phase < np; ++phase) {
+            eqs[phase] = eqs[phase] * matbalscale[phase];
+        }
+
+        // Add material balance equations to form pressure equation
+        // in place of first balance equation.
+        for (int phase = 1; phase < np; ++phase) {
+            eqs[0] += eqs[phase];
+        }
+
+        // Scale pressure equation.
+        const double pscale = 200*unit::barsa;
+        eqs[0] = eqs[0] * pscale;
+
         // Combine in single block.
         ADB total_residual = eqs[0];
         for (int phase = 1; phase < np; ++phase) {
@@ -143,6 +160,7 @@ namespace Opm
         // Solve reduced system.
         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(),
@@ -152,7 +170,7 @@ namespace Opm
                       "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
                       "Linear solver convergence failure.");
         }
-
+        */
 
 
         // Create ISTL matrix.
@@ -180,10 +198,12 @@ namespace Opm
         Preconditioner precond(A, Ae, relax);
 
         // Construct linear solver.
-        const double tolerance = 1e-8;
+        const double tolerance = 1e-3;
         const int maxit = 5000;
         const int verbosity = 1;
-        Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
+        // Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
+        const int restart = 40;
+        Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
 
         // Solve system.
         Dune::InverseOperatorResult result;
@@ -195,7 +215,7 @@ namespace Opm
         res.iterations = result.iterations;
         res.residual_reduction = result.reduction;
 
-        // std::copy(x.begin(), x.end(), dx.data());
+        std::copy(x.begin(), x.end(), dx.data());