System of PDEs

Introduction

In this example we will show how to implement a system of coupled PDEs. We define

\begin{eqnarray*} -\Delta u_1 &=& f &\text{ in } \Omega \\ u_0 - u_1 &=& 0 &\text{ in } \Omega \\ u_1 &=& 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 &=& -1.0 \\ g &=& 0. \end{eqnarray*}

The system can be written in matrix-vector form as

\[ \begin{pmatrix} I & -I \\ 0 & -\Delta \end{pmatrix} \begin{pmatrix} u_0\\ u_1 \end{pmatrix} = \begin{pmatrix} 0\\ f \end{pmatrix}. \]

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.

#include <config.h>
 
#include <iostream>
 
#include <amdis/AMDiS.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>

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. As every AMDiS program requires us to set up an Environment at the start we do that here.

int main(int argc, char** argv)
{
  Environment env(argc, argv);

We define the grid manager used in this example - a structured axis-aligned cube grid. Note that GRIDDIM is a constant passed by the buildsystem. By default we build the example for \( GRIDDIM = 2 \). Then we proceed to call the MeshCreator to build our grid.

  using Grid = Dune::YaspGrid<GRIDDIM>;
  auto grid = MeshCreator<Grid>{"elliptMesh"}.create();

Now we create and initialize the problem. We use the Dune::Functions::BasisFactory functions power and lagrange to define nodal (lagrange) basis functions of order 2 for both \( u_0 \text{ and } u_1 \). We also set up an AdaptInfo object that controls the adaptation procedure later on.

  using namespace Dune::Functions::BasisFactory;
  ProblemStat prob("ellipt", *grid, power<2>(lagrange<2>()));
  prob.initialize(INIT_ALL);
 
  AdaptInfo adaptInfo("adapt");

Next we add the operators that define our problem. Since we are using a basis with more than one basis node we need to specify the indices of the basis node an operator belongs to. This corresponds to the indices within the matrix-vector form above as follows

\begin{eqnarray*} (0,0) &\rightarrow& I \\ (0,1) &\rightarrow& -I \\ (1,0) &\rightarrow& 0 \\ (1,1) &\rightarrow& -\Delta. \end{eqnarray*}

  // ⟨∇u₁,∇v₁⟩
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
  prob.addMatrixOperator(opL, 1, 1);
 
  // ⟨u₀,v₀⟩ and ⟨u₁,v₀⟩
  auto opM = makeOperator(tag::test_trial{}, 1.0);
  prob.addMatrixOperator(opM, 0, 0);
  auto opM2 = makeOperator(tag::test_trial{}, -1.0);
  prob.addMatrixOperator(opM2, 0, 1);
 
  // (f,v₁) with f ≡ -1
  auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
  prob.addVectorOperator(opForce, 1);

Lastly we set the constraints: We have dirichlet boundary conditions on the left and lower part of the boundary, as indicated by the predicate function and prescribe a value of zero. Since we only want to apply the constraint to the first component \( u_1 \) we pass the index 1 as argument to the function.

  // set boundary condition
  auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
  auto dbcValues = [](auto const& x){ return 0.0; }; // set value
  prob.addDirichletBC(predicate, 1, 1, dbcValues);

Now we can assemble the system, solve it and write the results to a file.

  prob.buildAfterAdapt(adaptInfo, Flag(0));
  prob.solve(adaptInfo);
  prob.writeFiles(adaptInfo);

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 file vecellipt.dat.2d. For a full description of all options see the reference page.

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: 3
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". Note that the index in brackets relates to the index of the basis node, so [0] refers to the component \( u_0 \) while [1] refers to the component \( u_1 \).

• 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[0]->format: vtk
ellipt->output[0]->filename: vecellipt0.2d
ellipt->output[0]->name: u0
ellipt->output[0]->output directory: ./output
ellipt->output[1]->format: vtk
ellipt->output[1]->filename: vecellipt1.2d
ellipt->output[1]->name: u1
ellipt->output[1]->output directory: ./output

Running the example

Using a console we navigate to the examples directory where the vecellipt.2d files resides and run the file while also passing the necessary arguments.

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

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

Command line arguments

• Initfile, <amdis-root>/examples/init/vecellipt.dat.2d

Output

In this section we show some sample output of vecellipt.2d using the default arguments.

162 local DOFs
fill-in of assembled matrix: 4356
assemble needed 0.000800004 seconds
=== GeneralizedPCGSolver
=== rate=0.489578, T=0.000923321, TIT=2.71565e-05, IT=35
solution of discrete system needed 0.0066155 seconds
writeFiles needed 0.000301686 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:

• vecellipt.cc
• vecellipt.dat.2d