Fixed commented out functions and some warnings

opm/autodiff/AutoDiffMatrix.hpp

@@ -47,10 +47,10 @@ namespace Opm
 
 
 
-        AutoDiffMatrix(const int rows, const int cols)
+        AutoDiffMatrix(const int num_rows, const int num_cols)
             : type_(Z),
-              rows_(rows),
+              rows_(num_rows),
-              cols_(cols)
+              cols_(num_cols)
         {
         }
 
@@ -59,10 +59,10 @@ namespace Opm
         enum CreationType { ZeroMatrix, IdentityMatrix };
 
 
-        AutoDiffMatrix(const CreationType t, const int rows)
+        AutoDiffMatrix(const CreationType t, const int num_rows)
             : type_(t == ZeroMatrix ? Z : I),
-              rows_(rows),
+              rows_(num_rows),
-              cols_(rows)
+              cols_(num_rows)
         {
         }
 

opm/autodiff/NewtonIterationBlackoilSimple.cpp

@@ -47,7 +47,6 @@ namespace Opm
     NewtonIterationBlackoilSimple::SolutionVector
     NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
     {
-        /*
         typedef LinearisedBlackoilResidual::ADB ADB;
         const int np = residual.material_balance_eq.size();
         ADB mass_res = residual.material_balance_eq[0];
@@ -57,7 +56,8 @@ namespace Opm
         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];
+        Eigen::SparseMatrix<double, Eigen::RowMajor> matr;
+        total_residual.derivative()[0].toSparse(matr);
 
         SolutionVector dx(SolutionVector::Zero(total_residual.size()));
         Opm::LinearSolverInterface::LinearSolverReport rep
@@ -74,7 +74,6 @@ namespace Opm
                       "Linear solver convergence failure.");
         }
         return dx;
-        */
     }
 
     const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const

opm/autodiff/NewtonIterationUtilities.cpp

@@ -39,6 +39,7 @@ namespace Opm
     typedef AutoDiffBlock<double> ADB;
     typedef ADB::V V;
     typedef ADB::M M;
+    typedef Eigen::SparseMatrix<double> S;
 
 
     std::vector<ADB> eliminateVariable(const std::vector<ADB>& eqs, const int n)
@@ -75,7 +76,8 @@ namespace Opm
         const Sp Di = solver.solve(id);
 
         // compute inv(D)*bn for the update of the right hand side
-        const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix());
+        ADB::V eqs_n_v = eqs[n].value();
+        const Eigen::VectorXd& Dibn = solver.solve(eqs_n_v.matrix());
 
         std::vector<V> vals(num_eq);              // Number n will remain empty.
         std::vector<std::vector<M>> jacs(num_eq); // Number n will remain empty.
@@ -198,7 +200,6 @@ namespace Opm
                             Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
                             V& b)
     {
-        /*
         if (num_phases != 3) {
             OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases.");
         }
@@ -215,14 +216,17 @@ namespace Opm
         // The l1 block indicates if the equation for a given cell and phase is
         // sufficiently strong on the diagonal.
         Block l1 = Block::Zero(n, num_phases);
-        for (int phase = 0; phase < num_phases; ++phase) {
+        {
-            const M& J = eqs[phase].derivative()[0];
+            S J;
-            V dj = J.diagonal().cwiseAbs();
+            for (int phase = 0; phase < num_phases; ++phase) {
-            V sod = V::Zero(n);
+                eqs[phase].derivative()[0].toSparse(J);
-            for (int elem = 0; elem < n; ++elem) {
+                V dj = J.diagonal().cwiseAbs();
-                sod(elem) = J.col(elem).cwiseAbs().sum() - dj(elem);
+                V sod = V::Zero(n);
+                for (int elem = 0; elem < n; ++elem) {
+                    sod(elem) = J.col(elem).cwiseAbs().sum() - dj(elem);
+                }
+                l1.col(phase) = (dj/sod > ratio_limit).cast<double>();
             }
-            l1.col(phase) = (dj/sod > ratio_limit).cast<double>();
         }
 
         // By default, replace first equation with sum of all phase equations.
@@ -269,16 +273,18 @@ namespace Opm
             t.emplace_back(i3[ii], i1[ii], l31(ii));
             t.emplace_back(i3[ii], i3[ii], l33(ii));
         }
-        M L(3*n, 3*n);
+        S L(3*n, 3*n);
         L.setFromTriplets(t.begin(), t.end());
 
         // Combine in single block.
         ADB total_residual = vertcatCollapseJacs(eqs);
 
+        S derivative;
+        total_residual.derivative()[0].toSparse(derivative);
+
         // Create output as product of L with equations.
-        A = L * total_residual.derivative()[0];
+        A = L * derivative;
         b = L * total_residual.value().matrix();
-        */
     }
 
 

opm/autodiff/fastSparseProduct.hpp

@@ -151,6 +151,7 @@ inline void fastDiagSparseProduct(// const Eigen::DiagonalMatrix<double, Eigen::
 
     // Multiply rows by diagonal lhs.
     int n = res.cols();
+#pragma omp parallel for
     for (int col = 0; col < n; ++col) {
         typedef Eigen::SparseMatrix<double>::InnerIterator It;
         for (It it(res, col); it; ++it) {
@@ -171,6 +172,7 @@ inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
 
     // Multiply columns by diagonal rhs.
     int n = res.cols();
+#pragma omp parallel for
     for (int col = 0; col < n; ++col) {
         typedef Eigen::SparseMatrix<double>::InnerIterator It;
         for (It it(res, col); it; ++it) {