CHANGELOG.md

@@ -5,6 +5,8 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 # Master (will become release 2.10)
 
+- `setupLaplacian` has now an optional diagonal regularization parameter.
+
 - `UMFPACK` is extended to arbitrary blocked matrix structures. This includes `MultiTypeBlockMatrix`. The external interface is unchanged.
   However, internally the `BCCSMatrixInitializer` is replaced by direct calls of `flatMatrixForEach` similar to `Cholmod`. This requires a
   compatible vector field of "ignored" degrees of freedom. The method `setSubMatrix` with a top-level index set is preserved for compatibility.

dune/istl/test/cholmodtest.cc

@@ -40,7 +40,8 @@ bool test()
 
     // fill matrix with external method
     BCRSMatrix<FieldMatrix<double,bs,bs>> A;
-    setupLaplacian(A, N);
+    // no regularization, A may be singular!
+    setupLaplacian(A, N, 0.0);
 
     BlockVector<FieldVector<double,bs>> b,x;
     b.resize(A.N());

dune/istl/test/laplacian.hh

@@ -4,6 +4,9 @@
 // vi: set et ts=4 sw=2 sts=2:
 #ifndef LAPLACIAN_HH
 #define LAPLACIAN_HH
+
+#include <optional>
+
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/scalarmatrixview.hh>
@@ -38,15 +41,22 @@ void setupSparsityPattern(Dune::BCRSMatrix<B,Alloc>& A, int N)
   }
 }
 
-
+/* \brief Setup a FD Laplace matrix on a square 2D grid
+ *
+ * \param A The resulting matrix
+ * \param N Number of grid nodes per direction
+ * \param reg (optional) Regularization added to the diagonal to ensure a good condition number
+ */
 template<class B, class Alloc>
-void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A, int N)
+void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A,
+                    int N,
+                    std::optional<typename Dune::BCRSMatrix<B,Alloc>::field_type> reg = {} )
 {
-  typedef typename Dune::BCRSMatrix<B,Alloc>::field_type FieldType;
+  using field_type = typename Dune::BCRSMatrix<B,Alloc>::field_type;
 
   setupSparsityPattern(A,N);
 
-  B diagonal(static_cast<FieldType>(0)), bone(static_cast<FieldType>(0));
+  B diagonal(static_cast<field_type>(0)), bone(static_cast<field_type>(0));
 
   auto setDiagonal = [](auto&& scalarOrMatrix, const auto& value) {
     auto&& matrix = Dune::Impl::asMatrix(scalarOrMatrix);
@@ -57,42 +67,42 @@ void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A, int N)
   setDiagonal(diagonal, 4.0);
   setDiagonal(bone, -1.0);
 
-  for (typename Dune::BCRSMatrix<B,Alloc>::RowIterator i = A.begin(); i != A.end(); ++i) {
+  for (auto i = A.begin(); i != A.end(); ++i)
+  {
     int x = i.index()%N; // x coordinate in the 2d field
     int y = i.index()/N;  // y coordinate in the 2d field
 
-    /*    if(x==0 || x==N-1 || y==0||y==N-1){
-
-       i->operator[](i.index())=1.0;
+    i->operator[](i.index())=diagonal;
 
-       if(y>0)
-       i->operator[](i.index()-N)=0;
+    if(y>0)
+      i->operator[](i.index()-N)=bone;
 
-       if(y<N-1)
-       i->operator[](i.index()+N)=0.0;
+    if(y<N-1)
+      i->operator[](i.index()+N)=bone;
 
-       if(x>0)
-       i->operator[](i.index()-1)=0.0;
+    if(x>0)
+      i->operator[](i.index()-1)=bone;
 
-       if(x < N-1)
-       i->operator[](i.index()+1)=0.0;
+    if(x < N-1)
+      i->operator[](i.index()+1)=bone;
+  }
 
-       }else*/
+  // add regularization to the diagonal
+  if ( reg.has_value() )
+  {
+    for ( auto rowIt=A.begin(); rowIt!=A.end(); rowIt++ )
     {
-
-      i->operator[](i.index())=diagonal;
-
-      if(y>0)
-        i->operator[](i.index()-N)=bone;
-
-      if(y<N-1)
-        i->operator[](i.index()+N)=bone;
-
-      if(x>0)
-        i->operator[](i.index()-1)=bone;
-
-      if(x < N-1)
-        i->operator[](i.index()+1)=bone;
+      for ( auto entryIt=rowIt->begin(); entryIt!=rowIt->end(); entryIt++ )
+      {
+        if ( entryIt.index() == rowIt.index() )
+        {
+          auto&& block = Dune::Impl::asMatrix(*entryIt);
+          for ( auto it=block.begin(); it!=block.end(); ++it )
+          {
+            (*it)[it.index()] += reg.value();
+          }
+        }
+      }
     }
   }
 }

dune/istl/test/umfpacktest.cc

@@ -18,24 +18,6 @@
 
 using namespace Dune;
 
-// helper to add something to the diagonal of a BCRSMatrix with blocks
-auto addToDiagonal = [](auto&& matrix, const auto& value) {
-  for( auto rowIt=matrix.begin(); rowIt!=matrix.end(); rowIt++ )
-  {
-    for( auto entryIt=rowIt->begin(); entryIt!=rowIt->end(); entryIt++ )
-    {
-      if ( entryIt.index() == rowIt.index() )
-      {
-        auto&& block = Dune::Impl::asMatrix(*entryIt);
-        for (auto it=block.begin(); it!=block.end(); ++it)
-        {
-          (*it)[it.index()] += value;
-        }
-      }
-    }
-  }
-};
-
 template<typename Matrix, typename Vector, typename BitVector>
 TestSuite runUMFPack(std::size_t N)
 {
@@ -44,10 +26,8 @@ TestSuite runUMFPack(std::size_t N)
   Matrix mat;
   Vector b(N*N), x(N*N), b1(N/2), x1(N/2);
 
-  setupLaplacian(mat,N);
-
-  // make this regular
-  addToDiagonal(mat,100.0);
+  // with additional diagonal regularization
+  setupLaplacian(mat,N,100.0);
 
   b=1;
   b1=1;
@@ -150,14 +130,10 @@ TestSuite runUMFPackMultitype2x2(std::size_t N)
 
 
   // initialize all 4 matrix blocks
-  setupLaplacian(mat[_0][_0],N);
+  setupLaplacian(mat[_0][_0],N,100.0); // regularize the diagonal block
   setupLaplacian(mat[_0][_1],N);
   setupLaplacian(mat[_1][_0],N);
-  setupLaplacian(mat[_1][_1],N);
-
-  // make this regular
-  addToDiagonal(mat[_0][_0],200.0);
-  addToDiagonal(mat[_1][_1],100.0);
+  setupLaplacian(mat[_1][_1],N,200.0); // regularize the diagonal block
 
   b=1.0;
   b1=1.0;