diff --git a/examples/poisson-pq2.cc b/examples/poisson-pq2.cc
index aa127d2fd9662d1af29beba7c16899dbd88c2012..9d42013b901894f456aea05b319f76504d7b6f1f 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,33 @@ int main (int argc, char *argv[]) try
// Technicality: turn the matrix into a linear operator
MatrixAdapter<Matrix,Vector,Vector> op(stiffnessMatrix);
+// 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
- 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 +252,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) {
diff --git a/examples/stokes-taylorhood.cc b/examples/stokes-taylorhood.cc
index e24c91bf77a1294c0f4647ecb7187e576e16051c..adf3c5cdf3d8c13904f5ab520d0b158fbe6f0328 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,26 +60,7 @@ 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);
-}
-
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 10)
template<class V, class... P>
class BlockDiagonalPreconditioner : public Preconditioner<V,V> {
@@ -329,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-
@@ -355,8 +330,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 523f453f659d3f49f6c237bdecf65080ed35c5b1..ceb1b407aea023dada0056ba3c7f8250219b4c6e 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_GTE(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,53 @@ 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.
+ * \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);
+}
+#endif
+
+
+
// 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)