From 0007a32fd1d959612abb11e2052e14c3ed0fa67a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@math.fau.de>
Date: Mon, 20 Jan 2025 15:13:11 +0100
Subject: [PATCH 1/4] [examples] Add convenience functions for creating a
Dune::FastAMG preconditioner
---
examples/stokes-taylorhood.cc | 23 ++--------------
examples/utilities.hh | 49 +++++++++++++++++++++++++++++++++++
2 files changed, 51 insertions(+), 21 deletions(-)
diff --git a/examples/stokes-taylorhood.cc b/examples/stokes-taylorhood.cc
index e24c91bf..e4591689 100644
--- a/examples/stokes-taylorhood.cc
+++ b/examples/stokes-taylorhood.cc
@@ -67,25 +67,6 @@ using namespace Dune;
#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
-template<class V, class M>
-auto makeAMGPreconditioner(const V& x, const M& m)
-{
- using Operator = Dune::MatrixAdapter<M,V,V>;
- using Criterion = Dune::Amg::SymmetricCriterion<M, Dune::Amg::RowSum>;
- using AMG = Dune::Amg::FastAMG<Operator, V>;
-
- auto parameters = Dune::Amg::Parameters();
- parameters.setDebugLevel(0);
- parameters.setNoPreSmoothSteps(1);
- parameters.setNoPostSmoothSteps(1);
-
- auto criterion = Criterion();
- criterion.setDebugLevel(0);
-
- auto opPtr = std::make_shared<Operator>(m);
- return AMG(opPtr, criterion, parameters);
-}
-
template<class V, class... P>
class BlockDiagonalPreconditioner : public Preconditioner<V,V> {
@@ -355,8 +336,8 @@ int main (int argc, char *argv[]) try
using namespace Dune::Indices;
auto preconditioner = BlockDiagonalPreconditioner(x,
- makeAMGPreconditioner(x[_0], stiffnessMatrix[_0][_0]),
- makeAMGPreconditioner(x[_1], pressureMassMatrix));
+ makeAMGPreconditioner(stiffnessMatrix[_0][_0], x[_0]),
+ makeAMGPreconditioner(pressureMassMatrix, x[_1]));
// Construct the actual iterative solver
MINRESSolver<Vector> solver(
diff --git a/examples/utilities.hh b/examples/utilities.hh
index 523f453f..cd48ae52 100644
--- a/examples/utilities.hh
+++ b/examples/utilities.hh
@@ -10,6 +10,10 @@
#include <dune/common/typetraits.hh>
#include <dune/common/concept.hh>
+#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
+#include <dune/istl/paamg/fastamg.hh>
+#endif
+
#include <dune/grid/common/rangegenerators.hh>
#include <dune/functions/backends/concepts.hh>
@@ -141,6 +145,51 @@ void incorporateEssentialConstraints(const Basis& basis, Matrix& matrix, Vector&
+/**
+ * \brief Create a FastAMG preconditioner
+ * \tparam Matrix Type of the matrix the preconditioner should act on.
+ * \tparam Vector Type of vectors the preconditioner should act on.
+ * \param matrix The matrix used in the preconditioner
+ *
+ * This constructs a Dune::FastAMG preconditioner. With some
+ * standard arguments that are not defaulted in dune-istl.
+ */
+template<class Vector, class Matrix>
+auto makeAMGPreconditioner(const Matrix& matrix)
+{
+ using Operator = Dune::MatrixAdapter<Matrix,Vector,Vector>;
+ using Criterion = Dune::Amg::SymmetricCriterion<Matrix, Dune::Amg::RowSum>;
+ using AMG = Dune::Amg::FastAMG<Operator, Vector>;
+
+ auto parameters = Dune::Amg::Parameters();
+ parameters.setDebugLevel(0);
+ parameters.setNoPreSmoothSteps(1);
+ parameters.setNoPostSmoothSteps(1);
+
+ auto criterion = Criterion();
+ criterion.setDebugLevel(0);
+
+ auto opPtr = std::make_shared<Operator>(matrix);
+ return AMG(opPtr, criterion, parameters);
+}
+
+/**
+ * \brief Create a FastAMG preconditioner
+ * \tparam Matrix Type of the matrix the preconditioner should act on.
+ * \tparam Vector Type of vectors the preconditioner should act on.
+ * \param matrix The matrix used in the preconditioner
+ * \param dummy A dummy vector argument to enable type deduction for the Vector type.
+ *
+ * This constructs a Dune::FastAMG preconditioner.
+ */
+template<class Vector, class Matrix>
+auto makeAMGPreconditioner(const Matrix& matrix, const Vector& dummy)
+{
+ return makeAMGPreconditioner<Vector>(matrix);
+}
+
+
+
// Compute the (i,j)-th cofactor of F
template<class K, int d>
constexpr K cofactor(const Dune::FieldMatrix<K,d,d>& F, std::size_t i, std::size_t j)
--
GitLab
From e208f1d007f1255077578bb711fb75df5271b2ea Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@math.fau.de>
Date: Mon, 20 Jan 2025 17:33:20 +0100
Subject: [PATCH 2/4] [examples] Improve poisson-pq2.cc
* Use logger to unify messages.
* Fix matrix and vector types (no nesting since we have flat indices).
* Use AMG instead of ILDL preconditioner if available.
---
examples/poisson-pq2.cc | 54 +++++++++++++++++++++++++++++++----------
1 file changed, 41 insertions(+), 13 deletions(-)
diff --git a/examples/poisson-pq2.cc b/examples/poisson-pq2.cc
index aa127d2f..24e4e844 100644
--- a/examples/poisson-pq2.cc
+++ b/examples/poisson-pq2.cc
@@ -7,6 +7,7 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/version.hh>
+#include <dune/common/stringutility.hh>
#include <dune/geometry/quadraturerules.hh>
@@ -32,6 +33,7 @@
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#endif
+#include <dune/fufem/logger.hh>
#include <dune/fufem/parallel/elementcoloring.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
@@ -71,6 +73,9 @@ auto createMixedGrid()
int main (int argc, char *argv[]) try
{
+
+ auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
+
std::size_t threadCount = 4;
if (argc>1)
threadCount = std::stoul(std::string(argv[1]));
@@ -79,10 +84,12 @@ int main (int argc, char *argv[]) try
if (argc>2)
refinements = std::stoul(std::string(argv[2]));
+ log("Command line processed");
- // Set up MPI, if available
MPIHelper::instance(argc, argv);
+ log("MPI initialized");
+
///////////////////////////////////
// Generate the grid
///////////////////////////////////
@@ -93,8 +100,12 @@ int main (int argc, char *argv[]) try
auto grid = createUniformCubeGrid<2>();
#endif
+ log("Grid created");
+
grid->globalRefine(refinements);
+ log(Dune::formatString("Grid refined %d times", refinements));
+
auto gridView = grid->leafGridView();
using GridView = decltype(gridView);
@@ -105,13 +116,15 @@ int main (int argc, char *argv[]) try
using Basis = Dune::Functions::LagrangeBasis<GridView,2>;
Basis basis(gridView);
+ log(Dune::formatString("Basis created with dimension %d", basis.dimension()));
+
/////////////////////////////////////////////////////////
// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
- using Vector = Dune::BlockVector<Dune::FieldVector<double,1>>;
- using BitVector = Dune::BlockVector<Dune::FieldVector<char,1>>;
- using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
+ using Vector = Dune::BlockVector<double>;
+ using BitVector = std::vector<bool>;
+ using Matrix = Dune::BCRSMatrix<double>;
Vector rhs;
Vector sol;
@@ -122,8 +135,6 @@ int main (int argc, char *argv[]) try
// Assemble the system
/////////////////////////////////////////////////////////
- std::cout << "Number of DOFs is " << basis.dimension() << std::endl;
-
auto dirichletIndicatorFunction = [](auto x) { return true; };
auto dirichletPatch = BoundaryPatch(basis.gridView());
@@ -137,9 +148,11 @@ int main (int argc, char *argv[]) try
// Disable parallel assembly for UGGrid before 2.10
#if DUNE_VERSION_GT(DUNE_FUNCTIONS, 2, 9)
auto gridViewPartition = Dune::Fufem::coloredGridViewPartition(gridView);
- auto globalAssembler = GlobalAssembler(basis, basis, 1, std::cref(gridViewPartition), threadCount);
+ log("Grid coloring computed");
+
+ auto globalAssembler = GlobalAssembler(basis, basis, 0, std::cref(gridViewPartition), threadCount);
#else
- auto globalAssembler = GlobalAssembler(basis, basis, 1);
+ auto globalAssembler = GlobalAssembler(basis, basis, 0);
#endif
{
@@ -153,8 +166,11 @@ int main (int argc, char *argv[]) try
auto a = integrate(dot(grad(u), grad(v)));
auto b = integrate(f*v);
+ log("Assembler set up");
globalAssembler.assembleOperator(stiffnessMatrix, a);
+ log("Matrix assembled");
globalAssembler.assembleFunctional(rhs, b);
+ log("RHS assembled");
}
// *********************************************
@@ -171,11 +187,13 @@ int main (int argc, char *argv[]) try
Dune::Functions::istlVectorBackend(isConstrained).resize(basis);
Dune::Functions::istlVectorBackend(isConstrained) = 0;
Dune::Fufem::markBoundaryPatchDofs(dirichletPatch, basis, isConstrained);
+ log("Boundary DOFs marked");
interpolate(basis, sol, dirichletValues, isConstrained);
+ log("Boundary DOFs interpolated");
incorporateEssentialConstraints(basis, stiffnessMatrix, rhs, isConstrained, sol);
-
+ log("Boundary condition incorporated into system");
////////////////////////////
// Compute solution
@@ -184,22 +202,32 @@ int main (int argc, char *argv[]) try
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> op(stiffnessMatrix);
+#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
+ auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
+#else
// Sequential incomplete LU decomposition as the preconditioner
- SeqILDL<Matrix,Vector,Vector> ildl(stiffnessMatrix,1.0);
+ auto preconditioner = SeqILDL<Matrix,Vector,Vector>(stiffnessMatrix,1.0);
+#endif
+
+ log("Preconditioner created");
// Preconditioned conjugate-gradient solver
CGSolver<Vector> cg(op,
- ildl, // preconditioner
+ preconditioner, // preconditioner
1e-4, // desired residual reduction factor
- 100, // maximum number of iterations
+ 1000, // maximum number of iterations
2); // verbosity of the solver
+ log("Solver created");
+
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
// Solve!
cg.apply(sol, rhs, statistics);
+ log("Linear system solved");
+
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
////////////////////////////////////////////////////////////////////////////
@@ -223,7 +251,7 @@ int main (int argc, char *argv[]) try
vtkWriter.addVertexData(solFunction, VTK::FieldInfo("sol", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("poisson-pq2");
#endif
-
+ log("Solution written to vtk file");
}
// Error handling
catch (Exception& e) {
--
GitLab
From b92f96cc89cf3f4e2f68685f498244e6253302ba Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@math.fau.de>
Date: Tue, 21 Jan 2025 11:19:05 +0100
Subject: [PATCH 3/4] [examples] Only enable makeFastAMG in version >=2.10
---
examples/utilities.hh | 2 ++
1 file changed, 2 insertions(+)
diff --git a/examples/utilities.hh b/examples/utilities.hh
index cd48ae52..ef9083b0 100644
--- a/examples/utilities.hh
+++ b/examples/utilities.hh
@@ -145,6 +145,7 @@ void incorporateEssentialConstraints(const Basis& basis, Matrix& matrix, Vector&
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 10)
/**
* \brief Create a FastAMG preconditioner
* \tparam Matrix Type of the matrix the preconditioner should act on.
@@ -187,6 +188,7 @@ auto makeAMGPreconditioner(const Matrix& matrix, const Vector& dummy)
{
return makeAMGPreconditioner<Vector>(matrix);
}
+#endif
--
GitLab
From 0558c96c8dd07aab72765b1f5431073927389485 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@math.fau.de>
Date: Tue, 21 Jan 2025 13:31:16 +0100
Subject: [PATCH 4/4] [examples] Use more fine-grained version testing for
FastAMG
* The way the `FastAMG` preconditioner is set up only works in 2.10.
* Using non-blocked matrices only works in master (>2.10) which is
now tested as >=2.11, because in the 2.10 CI job, checking fo >2.10
succeeds.
---
examples/poisson-pq2.cc | 3 ++-
examples/stokes-taylorhood.cc | 10 ++--------
examples/utilities.hh | 2 +-
3 files changed, 5 insertions(+), 10 deletions(-)
diff --git a/examples/poisson-pq2.cc b/examples/poisson-pq2.cc
index 24e4e844..9d42013b 100644
--- a/examples/poisson-pq2.cc
+++ b/examples/poisson-pq2.cc
@@ -202,7 +202,8 @@ int main (int argc, char *argv[]) try
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> op(stiffnessMatrix);
-#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
+// FastAMG is not working for non-blocked matrices in 2.10
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 11)
auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
#else
// Sequential incomplete LU decomposition as the preconditioner
diff --git a/examples/stokes-taylorhood.cc b/examples/stokes-taylorhood.cc
index e4591689..adf3c5cd 100644
--- a/examples/stokes-taylorhood.cc
+++ b/examples/stokes-taylorhood.cc
@@ -22,12 +22,6 @@
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
-
-#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
-#include <dune/istl/paamg/fastamg.hh>
-#endif
-
-
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/taylorhoodbasis.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
@@ -66,7 +60,7 @@ using namespace Dune;
-#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 10)
template<class V, class... P>
class BlockDiagonalPreconditioner : public Preconditioner<V,V> {
@@ -310,7 +304,7 @@ int main (int argc, char *argv[]) try
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> stiffnessOperator(stiffnessMatrix);
-#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 10)
// As iterative solver we use MinRes with a block-diagonal preconditioner
// for the saddle point problem. An optimal preconditioner would be obtained,
// if the (0,0)-velocity-block is exact, while the (1,1)-block is the Schur-
diff --git a/examples/utilities.hh b/examples/utilities.hh
index ef9083b0..ceb1b407 100644
--- a/examples/utilities.hh
+++ b/examples/utilities.hh
@@ -10,7 +10,7 @@
#include <dune/common/typetraits.hh>
#include <dune/common/concept.hh>
-#if DUNE_VERSION_GT(DUNE_ISTL, 2, 10)
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 10)
#include <dune/istl/paamg/fastamg.hh>
#endif
--
GitLab