Working interleaved solver implemented.

opm/autodiff/NewtonIterationBlackoilInterleaved.cpp

@@ -70,6 +70,17 @@
 #include <Eigen/SparseLU>
 #endif
 
+// namespace Dune
+// {
+//     typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
+//     typedef Dune::BCRSMatrix <MatrixBlockType>        Mat;
+
+//     template <>
+//     struct MatrixDimension<Opm::DuneMatrix> : public MatrixDimension<Mat>
+//     {
+//     };
+// }
+
 namespace Opm
 {
 
@@ -77,8 +88,8 @@ namespace Opm
     typedef AutoDiffBlock<double> ADB;
     typedef ADB::V V;
     typedef ADB::M M;
-    typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
+    // typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
-    typedef Dune::BCRSMatrix <MatrixBlockType>        Mat;
+    // typedef Dune::BCRSMatrix <MatrixBlockType>        Mat;
 
 
     namespace {
@@ -171,11 +182,67 @@ namespace Opm
             eqs[phase] = eqs[phase] * matbalscale[phase];
         }
 
-        // Form interleaved system.
+        // // Write matrices to disk for analysis.
+        // const std::string eqnames[] = { "wat", "oil", "gas" };
+        // const std::string varnames[] = { "p", "sw", "x" };
+        // for (int p1 = 0; p1 < np; ++p1) {
+        //     for (int p2 = 0; p2 < np; ++p2) {
+        //         DuneMatrix m(eqs[p1].derivative()[p2]);
+        //         std::string filename = "mat-";
+        //         filename += eqnames[p1];
+        //         filename += "-";
+        //         filename += varnames[p2];
+        //         writeMatrixToMatlab(m, filename.c_str());
+        //     }
+        // }
+
+        // Find sparsity structure as union of basic block sparsity structures.
+        // Use addition to get to the union structure.
+        Eigen::SparseMatrix<double> structure = eqs[0].derivative()[0];
+        for (int p1 = 0; p1 < np; ++p1) {
+            for (int p2 = 0; p2 < np; ++p2) {
+                structure += eqs[p1].derivative()[p2];
+            }
+        }
+        // DuneMatrix ms(structure);
+        // writeMatrixToMatlab(ms, "structurematrix");
+        // Get row major form.
+        Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure;
+
+        // Form modified system.
         Eigen::SparseMatrix<double, Eigen::RowMajor> A;
         V b;
         formInterleavedSystem(np, eqs, A, b);
 
+        // Create ISTL matrix.
+        assert(np == 3);
+        Mat istlA(s.rows(), s.cols(), s.nonZeros(), Mat::row_wise);
+        const int* ia = s.outerIndexPtr();
+        const int* ja = s.innerIndexPtr();
+        for (Mat::CreateIterator row = istlA.createbegin(); row != istlA.createend(); ++row) {
+            int ri = row.index();
+            for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
+                row.insert(ja[i]);
+            }
+        }
+        const int size = s.rows();
+        Span span[3] = { Span(size, 1, 0),
+                         Span(size, 1, size),
+                         Span(size, 1, 2*size) };
+        for (int row = 0; row < size; ++row) {
+            for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) {
+                const int col = ja[col_ix];
+                MatrixBlockType block;
+                for (int p1 = 0; p1 < np; ++p1) {
+                    for (int p2 = 0; p2 < np; ++p2) {
+                        block[p1][p2] = A.coeff(span[p1][row], span[p2][col]);
+                    }
+                }
+                istlA[row][col] = block;
+            }
+        }
+
+
         // // Scale pressure equation.
         // const double pscale = 200*unit::barsa;
         // const int nc = residual.material_balance_eq[0].size();
@@ -186,14 +253,19 @@ namespace Opm
         SolutionVector dx(SolutionVector::Zero(b.size()));
 
         // Create ISTL matrix.
-        DuneMatrix istlA( A );
+        // DuneMatrix istlA( A );
 
         // Create ISTL matrix for elliptic part.
         // DuneMatrix istlAe( A.topLeftCorner(nc, nc) );
 
         // Right hand side.
         Vector istlb(istlA.N());
-        std::copy_n(b.data(), istlb.size(), istlb.begin());
+        for (int i = 0; i < size; ++i) {
+            istlb[i][0] = b(i);
+            istlb[i][1] = b(size + i);
+            istlb[i][2] = b(2*size + i);
+        }
+
         // System solution
         Vector x(istlA.M());
         x = 0.0;
@@ -214,7 +286,11 @@ namespace Opm
         }
 
         // Copy solver output to dx.
-        std::copy(x.begin(), x.end(), dx.data());
+        for (int i = 0; i < size; ++i) {
+            dx(i)          = x[i][0];
+            dx(size + i)   = x[i][1];
+            dx(2*size + i) = x[i][2];
+        }
 
         if( hasWells )
         {

opm/autodiff/NewtonIterationBlackoilInterleaved.hpp

@@ -61,8 +61,8 @@ namespace Opm
     /// as a block-structured matrix (one block for all cell variables).
     class NewtonIterationBlackoilInterleaved : public NewtonIterationBlackoilInterface
     {
-        typedef Dune::FieldVector<double, 1   > VectorBlockType;
+        typedef Dune::FieldVector<double, 3   > VectorBlockType;
-        typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
+        typedef Dune::FieldMatrix<double, 3, 3> MatrixBlockType;
         typedef Dune::BCRSMatrix <MatrixBlockType>        Mat;
         typedef Dune::BlockVector<VectorBlockType>        Vector;