From 09ab53067458c07a35ccffe7566114debb956667 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Wed, 11 May 2016 15:13:52 +0200
Subject: [PATCH] Add the option of computing the full sparity pattern based on
 all phases

- For some cases (for instance involving solvent flow) the reasoning for
only adding the pressure derivatives seems to fail. As getting the
sparsity pattern is non-trivial, in terms of work, the full sparsity
pattern is only added when specified by the parameter
"require_full_sparsity_pattern"
- For solvent runs "require_full_sparsity_pattern" defaults to true for
all other runs the default is to only extract the sparsity pattern from
the pressure derivatives.
---
 opm/autodiff/FlowMainSolvent.hpp                    |  9 +++++++++
 opm/autodiff/NewtonIterationBlackoilInterleaved.cpp | 13 ++++++++++++-
 opm/autodiff/NewtonIterationBlackoilInterleaved.hpp |  3 +++
 3 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/opm/autodiff/FlowMainSolvent.hpp b/opm/autodiff/FlowMainSolvent.hpp
index eaac004a8..4715f7c77 100644
--- a/opm/autodiff/FlowMainSolvent.hpp
+++ b/opm/autodiff/FlowMainSolvent.hpp
@@ -107,6 +107,15 @@ namespace Opm
                                                  Base::deck_->hasKeyword("SOLVENT")));
         }
 
+        void setupLinearSolver()
+        {
+            // require_full_sparsity_pattern as default for solvent runs
+            if (Base::deck_->hasKeyword("SOLVENT") && !Base::param_.has("require_full_sparsity_pattern") )  {
+                Base::param_.insertParameter("require_full_sparsity_pattern","true");
+            }
+            Base::setupLinearSolver();
+        }
+
     };
 
 
diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
index 88e16d224..f99cf36bd 100644
--- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
@@ -296,13 +296,24 @@ namespace Opm
             // Find sparsity structure as union of basic block sparsity structures,
             // corresponding to the jacobians with respect to pressure.
             // Use our custom PointOneOp to get to the union structure.
-            // Note that we only iterate over the pressure derivatives on purpose.
+            // As default we only iterate over the pressure derivatives.
             Eigen::SparseMatrix<double, Eigen::ColMajor> col_major = eqs[0].derivative()[0].getSparse();
             detail::PointOneOp<double> point_one;
             for (int phase = 1; phase < np; ++phase) {
                 const AutoDiffMatrix::SparseRep& mat = eqs[phase].derivative()[0].getSparse();
                 col_major = col_major.binaryExpr(mat, point_one);
             }
+            // For some cases (for instance involving Solvent flow) the reasoning for only adding
+            // the pressure derivatives fails. As getting the sparsity pattern is non-trivial, in terms
+            // of work, the full sparsity pattern is only added when required.
+            if (parameters_.require_full_sparsity_pattern_) {
+                for (int p1 = 0; p1 < np; ++p1) {
+                    for (int p2 = 1; p2 < np; ++p2) { // pressure is already added
+                        const AutoDiffMatrix::SparseRep& mat = eqs[p1].derivative()[p2].getSparse();
+                        col_major = col_major.binaryExpr(mat, point_one);
+                    }
+                }
+            }
 
             // Automatically convert the column major structure to a row-major structure
             Eigen::SparseMatrix<double, Eigen::RowMajor> row_major = col_major;
diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
index 0ce740fa8..f407cbdf6 100644
--- a/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
+++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
@@ -41,6 +41,7 @@ namespace Opm
         int    linear_solver_restart_;
         int    linear_solver_verbosity_;
         bool   newton_use_gmres_;
+        bool   require_full_sparsity_pattern_;
 
         NewtonIterationBlackoilInterleavedParameters() { reset(); }
         // read values from parameter class
@@ -55,6 +56,7 @@ namespace Opm
             linear_solver_maxiter_   = param.getDefault("linear_solver_maxiter", linear_solver_maxiter_);
             linear_solver_restart_   = param.getDefault("linear_solver_restart", linear_solver_restart_);
             linear_solver_verbosity_ = param.getDefault("linear_solver_verbosity", linear_solver_verbosity_);
+            require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_);
         }
 
         // set default values
@@ -65,6 +67,7 @@ namespace Opm
             linear_solver_maxiter_   = 75;
             linear_solver_restart_   = 40;
             linear_solver_verbosity_ = 0;
+            require_full_sparsity_pattern_ = false;
         }
     };