[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

@@ -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) {