Poisson equation

Table of Contents

• Introduction
• Setup and initialization
• Problem data
• Grid creation
• Preparing thread-parallelism
• Create finite element basis
• Create matrix, vector, and constraints data structures
• Definition of the variational problem
• Global assembly of the problem
• Algebraic solution
• Write solution to VTK
• Error handling

Full code of this example This example on gitlab

Introduction

The following example shows how to solve a Poisson problem with Dirichlet boundary conditions using second order finite elements.

Setup and initialization

After including the required headers the main() function first creates a simple logger that allows to write log messages to std::cout. The logger will prefix log messages with the time passed since the program was started and since the last log message using the given format string. Afterwards MPI is initialized, passing the command line arguments.

int main (int argc, char *argv[]) try
{
 
  auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
 
  Dune::MPIHelper::instance(argc, argv);
 
  log("MPI initialized");

Problem data

Next we parse the command line arguments. To this end we make use of Dune::ParameterTree which is a tree-of-strings data structure. Here parameters.get<T>("key", defaultValue) reads the value with the respective key from the ParameterTree and converts it to the type T if it exists, otherwise the defaultValue is used. By parsing the command line into the ParameterTree we can pass parameters using command line arguments of the form -refinements 4. The supported arguments allow to specify the number of grid refinement steps, the number of in parallel algorithms, and the tolerance, maximal iteration count and verbosity of the iterative solver.

  auto parameters = Dune::ParameterTree();
  Dune::ParameterTreeParser::readOptions(argc, argv, parameters);
 
  auto refinements  = parameters.get<std::size_t>("refinements", 4);
  auto threadCount  = parameters.get<std::size_t>("threads", std::thread::hardware_concurrency());
  auto tol          = parameters.get<double>("tol", 1e-4);
  auto maxIter      = parameters.get<std::size_t>("maxIter", 1000);
  auto verbosity    = parameters.get<std::size_t>("verbosity", 2);
 
  log("Command line processed");

The example solves the Poisson equation with Dirichlet boundary conditions. The problem data consisting of the right-hand-side function and the Dirichlet values are defined using lambda expressions.

  auto rightHandSide = [] (const auto& x) { return 10;};
  auto dirichletValues = [](const auto& x){ return std::sin(2*std::numbers::pi*x[0]); };

Grid creation

Instead of creating a uniform grid using a factory or reading a grid from a file, the example manually creates a mixed grid of type UGGrid<2> using a GridFactory. The grid implementation UGGrid<2> supports unstructured grid in two dimensions. After inserting nine uniformly placed vertices on the unit square, the method inserts two quadrilateral elements in the lower half and four triangular elements in the upper half of the domain. Finally, the call to factory.createGrid() returns a std::unique_ptr referring to the newly creates grid. For simplicity the example also stores a reference to the grid. The coarse grid is then refined globally and copy of the leaf grid view is stored.

Note

A GridView is a view to a subset of the grid's elements, vertices, ... that should be stored by value and can be copied cheaply. Grids in dune are in general hierarchical and composed by elements on several levels. The discretization usually lives on the set of most refined elements that is denote the leaf GridView in dune.

  using Grid = Dune::UGGrid<2>;
 
  auto factory = Dune::GridFactory<Grid>();
  for(unsigned int k : Dune::range(9))
    factory.insertVertex({0.5*(k%3), 0.5*(k/3)});
  factory.insertElement(Dune::GeometryTypes::cube(2), {0, 1, 3, 4});
  factory.insertElement(Dune::GeometryTypes::cube(2), {1, 2, 4, 5});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {3, 4, 6});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {4, 7, 6});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {4, 5, 7});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {5, 8, 7});
 
  auto gridPtr = factory.createGrid();
  auto& grid = *gridPtr;
 
  log("Grid created");
 
  grid.globalRefine(refinements);
 
  log(Dune::formatString("Grid refined %d times", refinements));
 
  auto gridView = grid.leafGridView();

Preparing thread-parallelism

Dune-fufem supports thread-parallel assembly based on a coloring of the grid view. To use this feature we need to compute such a coloring. The coloredElementRange() function from the dune-assembler module computes such a coloring and returns a partition of the grid view's elements according to the coloring. The ColoredRangeExecutor class then encapsulates the information for executing parallel algorithms based on this partition.

  auto coloredRange = Dune::Assembler::Experimental::coloredElementRange(gridView);
 
  log("Grid coloring computed");
 
  auto executor = Dune::Assembler::Experimental::ColoredRangeExecutor(coloredRange, threadCount);
 
  log(Dune::formatString("Using parallelism with %d threads", threadCount));

Create finite element basis

As a next step the program creates a global finite element basis on the GridView. The helper function lagrange<2>() declares that we want to use Lagrange finite elements of order two, while the makeBasis() function creates the basis of this type on the given GridView. Both functions are from the namespace Dune::Functions::BasisFactory which is included for simplicity.

  using namespace Dune::Functions::BasisFactory;
  auto basis = makeBasis(gridView, lagrange<2>());
 
  log(Dune::formatString("Basis created with dimension %d", basis.dimension()));

Create matrix, vector, and constraints data structures

Since the basis' indices are plain integers we next define types for flat vectors and matrices of reals numbers and create a solution and a right-hand-side vector and a stiffness matrix.

  using Vector = Dune::BlockVector<double>;
  using Matrix = Dune::BCRSMatrix<double>;
 
  auto sol = Vector();
  auto rhs = Vector();
  auto stiffnessMatrix = Matrix();

To handle the Dirichlet boundary conditions we then create an object for storing affine constraints with respect to the basis and then fill the object by computing the actual constraints values on the boundary from the Dirichlet value function.

  auto constraints = Dune::Fufem::makeAffineConstraints(basis);
 
  Dune::Fufem::computeBoundaryConstraints(constraints, basis, dirichletValues);
 
  log("Boundary constraints computed");

Definition of the variational problem

Inspired by the UFL language of the Fenics project project, dune-fufem allows to define local assemblers for linear and bilinear forms using simple expressions from the Dune::Fufem::Forms namespace. After including the namespace, we first obtain expressions for trial and test functions associated to the basis. The only difference between trialFunction() and testFunction() is that the resulting objects later refer to column and row indices, respectively. In order to use the right-hand-side lambda functions within a form expression we first wrap it into a GridView-functions and then declare this as coefficient function, i.e. a function not depending on test and trial functions. Now we can define the bilinear and linear form in a straight forward way. The defined bilinear form and linear form objects satisfy the interfaces for local matrix and vector assemblers to be passed to a global assembler.

  using namespace Dune::Fufem::Forms;
 
  auto u = trialFunction(basis);
  auto v = testFunction(basis);
  auto f = Coefficient(Dune::Functions::makeGridViewFunction(rightHandSide, gridView));
 
  auto a = integrate(dot(grad(u), grad(v)));
  auto l = integrate(f*v);

Global assembly of the problem

For global assembly we use the dune-assembler module. This module allows to use different linear algebra frameworks through a unified matrix- and vector-backend interface. To use the matrix- and vector-classes from dune-istl, we wrap them into corresponding ISTL-backends.

  auto matrixBackend = Dune::Assembler::ISTLMatrixBackend(stiffnessMatrix);
  auto rhsBackend = Dune::Assembler::ISTLVectorBackend(rhs);

Next we create a global assembler object using the same basis for the test and trial space. In order to use thread-parallel assembly, we additionally pass the parallel executor which encapsulates the grid coloring and requested thread count to the assembler.

  auto assembler = Dune::Assembler::Assembler(basis, basis, executor);
 
  log("Assembler set up");

As a first step we need to assemble the matrix pattern. To this end we obtain a patternBuilder object from the matrix backend and initialize its size according to the test and trial basis. Then we use the global assemblers assembleMatrixPattern() method to assemble the pattern. While local matrices are often fully populated, this is not necessarily the case. Hence the local assembler also defines a local pattern and thus has to be passed when assembling the pattern, too.

  auto patternBuilder = matrixBackend.patternBuilder();
 
  patternBuilder.resize(basis, basis);
 
  assembler.assembleMatrixPattern(a, patternBuilder);
 
  log("Matrix pattern assembled");

Note

There are cases where the constraints object requires to enlarge the computed pattern. This is especially the case for constraints enforcing \(H^1\)-conforming finite element spaces on grids with hanging nodes. In the present example this is not the case. Hence we can skip this step here.

Now we are ready to assemble the problem. To this end we let the patternBuilder initialize the matrix size and pattern and then call the actual assembly method for the matrix.

  patternBuilder.setupMatrix();
 
  log("Matrix set up");
 
  assembler.assembleMatrixEntries(a, matrixBackend);
 
  log("Matrix assembled");

Similarly we assemble the right hand side by first resizing the vector according to the test basis and initializing it with zero and then calling the corresponding assembly method.

  rhsBackend.resize(basis);
  rhsBackend.setZero();
 
  assembler.assembleVectorEntries(l, rhsBackend);
 
  log("RHS assembled");

Once the system is assembled, we can constrain it to the affine subspace satisfying the constraints that we defined earlier.

  constraints.constrainLinearSystem(stiffnessMatrix, rhs);
 
  log("Boundary condition incorporated into system");

Algebraic solution

Before solving the linear system, we have to initialize the solution vector. While we could directly use the corresponding resize method and assignment operator, we here demonstrate how the ISTLVectorBackend from dune-assembler can be used, since this approach would also work for basis using multi-indices and corresponding nested vector containers.

Note

There is also a utility Dune::Functions::istlVectorBackend() provided by the dune-functions module. The Dune::Assemble::ISTLVectorBackend extends the latter by a few methods required for assembly. Here we could have used the version from the dune-functions module as well with a slightly different syntax.

  Dune::Assembler::ISTLVectorBackend(sol).resize(basis);
  Dune::Assembler::ISTLVectorBackend(sol).setZero();

To use the matrix in the solvers from the dune-istl module we first wrap it into a suitable operator interface. Then we create an AMG preconditioner and pass the operator and preconditioner to a CG solver. The apply() method of the latter solves the linear system up to the given tolerance and writes some solver statistics to the respective object that needs to be passed additionally.

  auto matrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>(stiffnessMatrix);
 
  auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
 
  log("Preconditioner created");
 
  auto solver = Dune::CGSolver<Vector>(
                  matrixOperator,
                  preconditioner,
                  tol,
                  maxIter,
                  verbosity
                );
 
  log("Solver created");
 
  auto solverStatistics = Dune::InverseOperatorResult();
 
  solver.apply(sol, rhs, solverStatistics);
 
  log("Linear system solved");

Write solution to VTK

Finally we can write the solution to a VTK file understood by the ParaView visualization program. To this end we need to create a finite element function by combining the basis with the coefficient vector. Here this is obtained using the expression u|sol which binds the linear operator u (i.e. the identity on the finite element space) to the coefficient vector sol. We create a writer object capable of writing unstructured grids and tell it to write polynomial data of order 2 and pass the function to the writer, specifying the name and data type of the resulting field in the output.

  auto vtkWriter = Dune::Vtk::UnstructuredGridWriter(Dune::Vtk::LagrangeDataCollector(basis.gridView(), 2));
  vtkWriter.addPointData(u|sol, Dune::VTK::FieldInfo("sol", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.write("poisson-pq2");
 
  log("Solution written to vtk file");

Error handling

In case an exception is thrown within main() it is followed by a catch block that will print the error message and exit.

}
catch (Dune::Exception& e)
{
  std::cout << e.what() << std::endl;
  return 1;
}