Stationary problem with Dirichlet boundary condition

Introduction

In this example we will solve Poisson's equation with Dirichlet boundary condition

\begin{eqnarray*} - \Delta u &=& f &\text{ in } \Omega \\ u &=& g &\text{ on } \Gamma \end{eqnarray*}

whereby \( \Omega \subset \mathbb{R}^d, \quad 1 \le d \le 3 \) denotes the domain and \( \Gamma = \partial \Omega \) its boundary.

We set the right hand side functions to

\begin{eqnarray*} f &=& -(400 \| x \|^2 - 20 d) e^{-10 \|x\|^2} \\ g &=& e^{-10.0 \|x\|^2}. \end{eqnarray*}

The example code will solve the problem on a series of coarse to fine grids and compute the errors in \( L_2 \) and \( H_1 \) norm as well as experimental order of convergence (EOC).

Discretization

For using a Finite Element Method we derive a weak formulation

\[ \int_\Omega \nabla u \cdot \nabla v = \int_\Omega f v \quad \forall v \in V. \]

with solutions \( u \in \{ w \in H^1(\Omega) : tr_{\Gamma} w = g \} \) and trial space \( V = H^1_0(\Omega) = \{ w \in H^1(\Omega) : tr_{\Gamma} w = 0 \} \).

Implementation

We will now look at the source code of the example in detail. You can also find the complete source code for all files below.

Source code description

First we include the required AMDiS headers. AMDiS.hpp is a collection of headers that are commonly used. In this example we also use the ProblemStat class, which we explicitly include. We will also use the predefined operators in LocalOperators.hpp as well as the integrate function from Integrate.hpp.

#include <amdis/AMDiS.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>

On the following lines we include the grid manager used in this example - a structured axis-aligned cube grid. We set a name alias Grid for the type. Note that GRIDDIM is a constant passed by the buildsystem. By default we build the example for the values 2 and 3.

#include <dune/grid/yaspgrid.hh>
using Grid = Dune::YaspGrid<GRIDDIM>;

Most AMDiS classes and functions are found in the AMDiS namespace. To make life easier for us we therefore tell the compiler to look for our classes in this namespace.

using namespace AMDiS;

Our example begins with a call to the main function. Any command line argument we pass can be used with the argc and argv arguments.

int main(int argc, char** argv)
{

Every AMDiS program requires us to set up an Environment, so we do that here. We also start a timer to measure the execution time. Using the second command line argument we set the number of grid refinement steps. The next step is setting up our right hand side and boundary condition functors. We do that using a lambda function that takes a world coordinate as an argument.

  Environment env(argc, argv);
  Dune::Timer t;
 
  int numLevels = GRIDDIM == 2 ? 6 : 4;
  if (argc > 2)
    numLevels = std::atoi(argv[2]);
 
  auto f = [](auto const& x) {
    double r2 = dot(x,x);
    double ux = std::exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * x.size()) * ux;
  };
  auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
  auto grad_g = [](auto const& x) -> FieldMatrix<double,1,GRIDDIM> {
    return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
  };

Now we define the Finite Element space we want to work with - in this case using second order lagrange elements - and create and initialize the ProblemStat instance with the name ellipt. This class will contain all information about our problem as well as the means to solve it.

  using Param = LagrangeBasis<Grid, 2>;
  ProblemStat<Param> prob("ellipt");
  prob.initialize(INIT_ALL);

Next we tell the problem class what our PDE looks like. We add an operator of the form \( \int c \nabla u \cdot \nabla v\) with \( c = 1.0 \) as well as one of the form \( \int f v \) with a specified quadrature order of 6.

  // ⟨∇u,∇v⟩
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
  prob.addMatrixOperator(opL);
 
  // (f,v)
  auto opForce = makeOperator(tag::test{}, f, 6);
  prob.addVectorOperator(opForce);

Lastly we set the constraints: We have dirichlet boundary conditions on the complete boundary (since we did not specify any boundary segments the identifier 1 includes the complete boundary).

  // set boundary condition
  prob.addDirichletBC(1, g);

Before we start with the adapt-solve-loop we set up a helper object AdaptInfo that controls how the adaptation process is handled. We also initialize containers to store the error norms and element sizes in.

  AdaptInfo adaptInfo("adapt");
 
  std::vector<double> errL2; errL2.reserve(numLevels);
  std::vector<double> errH1; errH1.reserve(numLevels);
  std::vector<double> widths; widths.reserve(numLevels);

We now start with the adaptation loop by globally refining the domain.

  for (int l = 0; l < numLevels; ++l) {
    prob.globalRefine(1);

In this block we loop over all elements and calculate the maximum of the element widths.

    double h = 0;
    for (auto const& e : elements(prob.gridView(), Dune::Partitions::interior)) {
      auto refElem = Dune::referenceElement<double,GRIDDIM>(e.type());
      auto geo = e.geometry();
      for (int i = 0; i < refElem.size(GRIDDIM-1); ++i) { // edges
        auto v0 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,0,GRIDDIM), GRIDDIM));
        auto v1 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,1,GRIDDIM), GRIDDIM));
        h = std::max(h, (v0 - v1).two_norm());
      }
    }
    widths.push_back(h);

Here we tell the problem to assemble the operators we defined above and solve the resulting linear system of equations.

    prob.assemble(adaptInfo);
    prob.solve(adaptInfo);

At the end of the loop we calculate the error norms. Note that g is the exact solution of Poisson's equation.

    double errorL2 = integrate(sqr(g - prob.solution()), prob.gridView(), 6);
    errL2.push_back(std::sqrt(errorL2));
    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientOf(prob.solution())), prob.gridView(), 6);
    errH1.push_back(std::sqrt(errorH1));
  }

After we finished the adaptation loop for the last time we write the final solution to an output file defined in the initfile. We also print a table of the error norms and the EOC as well as the total execution time.

  // write last solution to file
  prob.writeFiles(adaptInfo);
 
  msg("");
  msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
      "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
  msg_("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
 
  std::vector<double> eocL2(numLevels, 0.0), eocH1(numLevels, 0.0);
  for (int i = 1; i < numLevels; ++i) {
    eocL2[i] = std::log(errL2[i]/errL2[i-1]) / std::log(widths[i]/widths[i-1]);
    eocH1[i] = std::log(errH1[i]/errH1[i-1]) / std::log(widths[i]/widths[i-1]);
  }
 
  for (int i = 0; i < numLevels; ++i)
    msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
        i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
 
  msg("elapsed time: {} seconds", t.elapsed());

With that the example ends. We return 0 to indicate a successful run.

  return 0;
}

Parameter file

We describe the options used by the example's parameter files. For a full description of all options see the reference page.

ellipt.dat.2d

Parameters for the grid "elliptMesh":

• global refinements: number of global refinement steps done as part of initialization
• num cells: number of grid cells in x- and y-direction
• overlap: size of the interface between neighbouring subdomains in parallel computations

elliptMesh->global refinements: 0
elliptMesh->num cells: 8 8
elliptMesh->overlap: 1

ellipt->mesh specifies the name of the grid used by the problem named "ellipt".

ellipt->mesh: elliptMesh

This parameter specifies the symmetry structure of the assembled matrix.

ellipt->symmetry: spd

Parameters for the solver employed by the problem "ellipt":

• solver: type of the solver, pcg means preconditioned conjugate gradient method
• info: verbosity of the output generated by the solver
• max iteration: maximum number of iterations until the solver aborts
• relative tolerance: factor by which the initial defect is reduced before the solver ends the iteration
• precon: preconditioner used in the computation, ilu means incomplete LU-deconposition

ellipt->solver: pcg
ellipt->solver->info: 1
ellipt->solver->max iteration: 10000
ellipt->solver->relative tolerance: 1.e-10
ellipt->solver->precon: ilu

Parameters for the output generated by the problem "ellipt":

• format: output format
• filename: name of the output file
• name: name of the solution within the output file
• output directory: directory where the output files are generated

ellipt->output->format: vtk
ellipt->output->filename: ellipt.2d
ellipt->output->name: u
ellipt->output->output directory: ./output

ellipt.dat.3d

See ellipt.dat.2d, except num cells has three entries since the grid is 3-dimensional.

Running the example

Using a console we navigate to the examples directory where the ellipt.#d files reside (# can be either 2 or 3) and run the file while also passing the necessary arguments.

Running ellipt.2d in bash could look like this:

<amdis-root>/build-cmake/examples $ ./ellipt.2d "<amdis-root>/examples/init/ellipt.dat.2d" 4

Command line arguments

• Initfile, <amdis-root>/examples/init/ellipt.dat.#d
• (optional) Number grid levels used for calculation of errors and EOC

Output

In this section we show some sample output of ellipt.2d using the default arguments. The text output of ellipt.3d is similar.

Every time the adapt loop is run the time for a global refinement step, the number of DOFs, the number of nonzero entries inserted into the matrix and the assembly time is printed.

globalRefine needed 0.000197335 seconds
289 local DOFs
fill-in of assembled matrix: 4225
assemble needed 0.000797755 seconds

After each assemble step the solver runs and prints convergence rate, elapsed time (total and per iteration) and total iterations until the tolerance is reached. Lastly the total time for the solve step is printed.

=== GeneralizedPCGSolver
=== rate=0.0718531, T=0.000249728, TIT=2.77476e-05, IT=10
solution of discrete system needed 0.0034118 seconds

After the last loop iteration finishes the output file is generated an a table with error norms and EOC as well as the total execution time is printed.

writeFiles needed 0.195444 seconds

level | h            | |u - u*|_L2  | |u - u*|_H1  | eoc_L2       | eoc_H1
------+--------------+--------------+--------------+--------------+-------------
1     | 0.125        | 0.000378951  | 0.0197187    | 0            | 0
2     | 0.0625       | 4.79736e-05  | 0.00497982   | 2.9817       | 1.9854
3     | 0.03125      | 6.01667e-06  | 0.00124812   | 2.9952       | 1.99634
4     | 0.015625     | 7.52723e-07  | 0.000312228  | 2.99877      | 1.99908
5     | 0.0078125    | 9.41105e-08  | 7.80695e-05  | 2.99969      | 1.99977
6     | 0.00390625   | 1.17645e-08  | 1.95181e-05  | 2.99992      | 1.99994
elapsed time: 12.9328 seconds

Using ParaView we can inspect the output file and obtain a visualization that looks like this.

See also

The complete source code for the example and the parameter files can be found here:

• ellipt.cc
• ellipt.dat.2d
• ellipt.dat.3d