Skip to content
Snippets Groups Projects
Commit f6af07e2 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

Merge branch 'feature/extend-example-utilities' into 'master' 

[examples] Add convenience utility for FastAMG preconditioner and use it in Poisson example

See merge request !249
parents c4921b6d 0558c96c
Branches
No related tags found
1 merge request!249[examples] Add convenience utility for FastAMG preconditioner and use it in Poisson example
Pipeline #76352 passed
Stage: .pre
Stage: test
Hide whitespace changes
Inline Side-by-side
examples/poisson-pq2.cc
+ 42
− 13
View file @ f6af07e2
......@@ -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) {
......
examples/stokes-taylorhood.cc
+ 4
− 29
View file @ f6af07e2
......@@ -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(
......
examples/utilities.hh
+ 51
− 0
View file @ f6af07e2
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment