Merge pull request #686 from totto82/sparsity_pattern

Add the option of computing the full sparity pattern based on all phases

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();
+        }
+
     };
 
 

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;

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;
         }
     };