Merge branch 'feature/umfpack-multitype' into 'master'

UMFPackSolver for MultiTypeBlockMatrix

See merge request !74

CHANGELOG.md

 # Master (will become release 2.10)
 
+- `UMFPackSolver` accepts each correctly formed blocked matrix. The user has to make sure that the vector types of `x` and `rhs` are compatible to the matrix.
+  The main advantage is that it is now possible to use `MultiTypeBlockMatrix`.
+
 - A new solver `ProximalNewtonSolver` is added which solves non-smooth minimization problems.
 
 # Release 2.9

dune/solvers/solvers/umfpacksolver.hh

@@ -61,8 +61,6 @@ public:
 
   void solve() override
   {
-    // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
-    VectorType mutableRhs = *rhs_;
 
     if (not this->hasIgnore())
     {
@@ -72,83 +70,86 @@ public:
       Dune::InverseOperatorResult statistics;
       Dune::UMFPack<MatrixType> solver(*matrix_);
       solver.setOption(UMFPACK_PRL, 0);   // no output at all
+
+      // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
+      VectorType mutableRhs = *rhs_;
       solver.apply(*x_, mutableRhs, statistics);
 
     }
     else
     {
-      ///////////////////////////////////////////////////////////////////////////////////////////
-      //  Extract the set of matrix rows that do not correspond to ignored degrees of freedom.
-      //  Unfortunately, not all cases are handled by the ISTL UMFPack solver.  Currently,
-      //  you can only remove complete block rows.  If a block is only partially ignored,
-      //  the dune-istl UMFPack solver cannot do it, and the code here will throw an exception.
-      //  All this can be fixed, but it needs going into the istl UMFPack code.
-      ///////////////////////////////////////////////////////////////////////////////////////////
-      std::set<std::size_t> nonIgnoreRows;
-      for (size_t i=0; i<matrix_->N(); i++)
-      {
-        auto const &ignore = this->ignore()[i];
-        if constexpr (std::is_convertible<decltype(ignore), bool>::value)
-        {
-          if (!ignore)
-            nonIgnoreRows.insert(i);
-        } else {
-          if (ignore.none())
-            nonIgnoreRows.insert(i);
-          else if (not ignore.all())
-            DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
-        }
-      }
-
       // Construct the solver
       Dune::InverseOperatorResult statistics;
       Dune::UMFPack<MatrixType> solver;
       solver.setOption(UMFPACK_PRL, 0);   // no output at all
       // We eliminate all rows and columns(!) from the matrix that correspond to ignored degrees of freedom.
       // Here is where the sparse LU decomposition is happenening.
-      solver.setSubMatrix(*matrix_,nonIgnoreRows);
+      solver.setMatrix(*matrix_,this->ignore());
+
+
+      // total number of dofs
+      auto N = flatVectorForEach(*rhs_, [](auto&&, auto&&){});
 
       // We need to modify the rhs vector by static condensation.
-      VectorType shortRhs(nonIgnoreRows.size());
-      int shortRowCount=0;
-      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
-        shortRhs[shortRowCount] = mutableRhs[*it];
+      std::vector<bool> flatIgnore(N);
+      size_t numberOfIgnoredDofs = 0;
+
+      flatVectorForEach(this->ignore(), [&](auto&& entry, auto&& offset)
+      {
+        flatIgnore[offset] = entry;
+        if ( entry )
+        {
+          numberOfIgnoredDofs++;
+        }
+      });
+
 
-      shortRowCount = 0;
-      for (size_t i=0; i<matrix_->N(); i++)
+      using field_type = typename MatrixType::field_type;
+      std::vector<field_type> shortRhs(N-numberOfIgnoredDofs);
+
+      // mapping of long indices to short indices
+      std::vector<size_t> subIndices(N,std::numeric_limits<size_t>::max());
+
+      int shortRhsCount=0;
+      flatVectorForEach(*rhs_, [&](auto&& entry, auto&& offset)
       {
-        if constexpr (std::is_convertible<decltype(this->ignore()[i]), const bool>::value) {
-          if (this->ignore()[i])
-            continue;
-        } else {
-          if (this->ignore()[i].all())
-            continue;
+        if ( not flatIgnore[offset] )
+        {
+          shortRhs[shortRhsCount] = entry;
+          subIndices[offset] = shortRhsCount;
+          shortRhsCount++;
         }
+      });
 
-        auto cIt    = (*matrix_)[i].begin();
-        auto cEndIt = (*matrix_)[i].end();
+      std::vector<field_type> flatX(N);
+      flatVectorForEach(*x_, [&](auto&& entry, auto&& offset)
+      {
+        flatX[offset] = entry;
+      });
 
-        for (; cIt!=cEndIt; ++cIt)
-          if constexpr (std::is_convertible<decltype(this->ignore()[cIt.index()]), const bool>::value) {
-            if (this->ignore()[cIt.index()])
-              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
-          } else {
-            if (this->ignore()[cIt.index()].all())
-              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
-          }
 
-        shortRowCount++;
-      }
+      flatMatrixForEach(*matrix_, [&](auto&& entry, auto&& row, auto&& col)
+      {
+        if ( flatIgnore[col] and not flatIgnore[row] )
+        {
+          shortRhs[ subIndices[ row ] ] -= entry * flatX[ col ];
+        }
+      });
 
       // Solve the reduced system
-      VectorType shortX(nonIgnoreRows.size());
-      solver.apply(shortX, shortRhs, statistics);
+      std::vector<field_type> shortX(N-numberOfIgnoredDofs);
 
-      // Blow up the solution vector
-      shortRowCount=0;
-      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
-        (*x_)[*it] = shortX[shortRowCount];
+      // call the raw-pointer variant of the ISTL UMPACK solver
+      solver.apply(shortX.data(), shortRhs.data());
 
+      // Blow up the solution vector
+      flatVectorForEach(*x_, [&](auto&& entry, auto&& offset)
+      {
+        if ( not flatIgnore[ offset ] )
+        {
+          entry = shortX[ subIndices[ offset ] ];
+        }
+      });
     }
   }
 

dune/solvers/test/umfpacksolvertest.cc

@@ -62,6 +62,99 @@ struct UMFPackSolverTestSuite
 
 
 
+/**  \brief test for the UMFPackSolver class with multitype matrix */
+template <size_t blocksize>
+struct UMFPackMultiTypeSolverTestSuite
+{
+  template <class GridType>
+  bool check(const GridType& grid)
+  {
+    double tol = 1e-11;
+
+    bool passed = true;
+
+    using Problem =
+    SymmetricSampleProblem<blocksize, typename GridType::LeafGridView>;
+    Problem p(grid.leafGridView());
+
+    // create a multitype problem
+    using BlockRow = MultiTypeBlockVector<typename Problem::Matrix,typename Problem::Matrix>;
+    using Matrix = MultiTypeBlockMatrix<BlockRow,BlockRow>;
+
+    using Vector = MultiTypeBlockVector<typename Problem::Vector,typename Problem::Vector>;
+
+    Matrix A;
+    Vector u, rhs;
+
+    using namespace Dune::Indices;
+    // initialize each block
+    A[_0][_0] = p.A;
+    A[_0][_1] = p.A;
+    A[_1][_0] = p.A;
+    A[_1][_1] = p.A;
+
+    u[_0] = p.u;
+    u[_1] = p.u;
+    rhs[_0] = p.rhs;
+    rhs[_1] = p.rhs;
+
+    using IgnoreBlock = decltype(p.ignore);
+    using Ignore = Solvers::TupleVector<IgnoreBlock,IgnoreBlock>;
+
+    Ignore ignore;
+    ignore[_0] = p.ignore;
+    ignore[_1] = p.ignore;
+
+    typedef Solvers::UMFPackSolver<Matrix,Vector> Solver;
+    Solver solver(A,u,rhs);
+
+    solver.setIgnore(ignore);
+
+    // solve blockdiagonal
+    A[_0][_1] *= 0.0;
+    A[_1][_0] *= 0.0;
+    solver.solve();
+
+    // the solution is ( u_ex, u_ex )
+    if (p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) > tol)
+    {
+      std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) << std::endl;
+      passed = false;
+    }
+
+    // solve a problem more generic problem without ignore field
+
+    A[_0][_1] = p.A;
+    A[_1][_0] = p.A;
+
+    // make A regular
+    A[_0][_0] *= 2.0;
+
+    u[_0] = p.u;
+    u[_1] = p.u;
+
+    // solve problem
+    Solver solver2(A,u,rhs);
+    solver2.preprocess();
+    solver2.solve();
+
+    // check residual
+    auto resisual = rhs;
+    A.mmv(u,resisual);
+
+    if (p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) > tol)
+    {
+      std::cout << "The UMFPackSolver did not produce a satisfactory result. ||residual||=" << p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) << std::endl;
+      passed = false;
+    }
+
+
+    return passed;
+  }
+};
+
+
+
 int main(int argc, char** argv)
 {
     Dune::MPIHelper::instance(argc, argv);
@@ -71,9 +164,13 @@ int main(int argc, char** argv)
     UMFPackSolverTestSuite<0> testsuite0;
     UMFPackSolverTestSuite<1> testsuite1;
     UMFPackSolverTestSuite<2> testsuite2;
+
+    UMFPackMultiTypeSolverTestSuite<3> testsuite3;
+
     passed = checkWithStandardGrids(testsuite0);
     passed = checkWithStandardGrids(testsuite1);
     passed = checkWithStandardGrids(testsuite2);
+    passed = checkWithStandardGrids(testsuite3);
 
     return passed ? 0 : 1;
 }