Time dependant scalar problem

Introduction

In this example we will solve the heat equation with Dirichlet boundary condition

\begin{eqnarray*} - \partial_t u - \Delta u &=& f &\text{ in } \Omega \times (t_{begin}, t_{end}) \\ u &=& g &\text{ on } \Gamma \times (t_{begin}, t_{end}) \\ u &=& u_0 &\text{ on } \Omega \times t_{begin} \end{eqnarray*}

whereby \( \Omega \subset \mathbb{R}^d, \quad 1 \le d \le 3 \) denotes the domain and \( \Gamma = \{ x \in \partial \Omega | \exists i \in \{1, ..., d\}: x_i = 0 \} \) a subset of its boundary. We solve the problem in the time interval \( (t_{begin}, t_{end}) \) with Dirichlet boundary conditions on \( \Gamma \) and set

\begin{eqnarray*} f(x,t) &=& -1 \\ g(x,t) &=& 0 \\ u_0(x) &=& 0. \end{eqnarray*}

Discretization

In this example we use the implicit Euler scheme

\[ \frac{u^{new} - u^{old}}{\tau} - \Delta u^{new} = f(\cdot, t^{old} + \tau) \]

where \( \tau = t^{new} - t^{old} \) is the timestep size between the old and new problem time. \( u^{new} \) is the searched solution at \( t = t^{new} \), \( u^{old} \) is the solution at \( t = t^{old} \), which is already known from the last timestep. If we bring all terms that depend on \( u^{old} \) to the right hand side, the equation reads

\[ \frac{u^{new}}{\tau} - \Delta u^{new} = - \frac{u^{old}}{\tau} + f(\cdot, t^{old} + \tau). \]

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 and ProblemInstat classes, 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/AdaptInstationary.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemInstat.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/GridFunctions.hpp>

On the following lines we define the grid manager and finite element space used in this example - a structured axis-aligned cube grid with nodal (lagrange) basis functions of order 2. We set name aliases for ease of notation for the basis and problem types. Note that GRIDDIM is a constant passed by the buildsystem. By default we build the example for \( GRIDDIM = 2 \).

// 1 component with polynomial degree 1
//using Grid = Dune::AlbertaGrid<GRIDDIM, WORLDDIM>;
using HeatParam   = YaspGridBasis<GRIDDIM, 2>;
using HeatProblem = ProblemStat<HeatParam>;
using HeatProblemInstat = ProblemInstat<HeatParam>;

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);

Now we create the stationary base problem using the ProblemStat class, as well as the instationary problem using the ProblemInstat class and initialize both. Note that we name both problem instances heat and initialize everything. We also set up an AdaptInfo object that controls the adaptation procedure later on.

  HeatProblem prob("heat");
  prob.initialize(INIT_ALL);
 
  HeatProblemInstat probInstat("heat", prob);
  probInstat.initialize(INIT_UH_OLD);
 
  AdaptInfo adaptInfo("adapt");

Next we add the operators that define our problem. Note that we store references to the inverse of the time step size invTau and the last solution prob.solution() so we use the correct values every timestep. See the operator reference page for descriptions of the operators used here.

  auto invTau = std::ref(probInstat.invTau());
 
  // 1/τ*⟨u,v⟩
  auto opTimeLhs = makeOperator(tag::test_trial{}, invTau);
  prob.addMatrixOperator(opTimeLhs);
 
  // ⟨∇u,∇v⟩
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
  prob.addMatrixOperator(opL);
 
  // 1/τ*⟨u_old,v⟩
  auto opTimeRhs = makeOperator(tag::test{},
    invokeAtQP([invTau](double u) { return u * invTau.get(); }, prob.solution()), 2);
  prob.addVectorOperator(opTimeRhs);
 
  // (-1.0,v)
  auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
  prob.addVectorOperator(opForce);

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.

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

With everything set up we can now create a helper object AdaptInstationary that performs the time-stepping scheme automatically using the parameters given by the AdaptInfo object we created earlier. We only need to call the adapt() function.

  AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
  adapt.adapt();

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

Parameters for the grid "heatMesh":

• global refinements: number of global refinement steps done as part of initialization
• min corner: lower left corner of the domain
• max corner: upper right corner of the domain
• num cells: number of grid cells in x- and y-direction

heatMesh->global refinements: 2
heatMesh->min corner: 0 0
heatMesh->max corner: 2 2
heatMesh->num cells: 8 8

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

heat->mesh: heatMesh

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

• solver: type of the solver, pcg means preconditioned conjugate gradient method
• info: verbosity of the output generated by the solver
• absolute tolerance: error tolerance that needs to be reached before the solver ends the iteration
• max iteration: maximum number of iterations until the solver aborts
• precon: preconditioner used in the computation, ilu means incomplete LU-deconposition

heat->solver: pcg
heat->solver->max iteration: 1000
heat->solver->absolute tolerance: 1e-6
heat->solver->info: 1
heat->solver->precon: ilu

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

• format: output format, here we generate a backup file
• filename: name of the output file
• name: name of the solution within the output file
• output directory: directory where the output files are generated
• animation: write an animated file opposed to just a single timestep

heat->output->format: vtk
heat->output->name: u
heat->output->filename: heat.2d
heat->output->output directory: output
heat->output->animation: 1

Parameters for the time stepping scheme:

• timestep: timestep size
• start time: \( t_{begin}\)
• end time: \( t_{end}\)

adapt->timestep:    0.1
adapt->start time:  0.0
adapt->end time:    1.0

Running the example

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

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

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

Command line arguments

• Initfile, <amdis-root>/examples/init/heat.dat.#d

Output

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

The example produces a series of output blocks like the one below for each time step. In each of them the current simulation time and time step size is printed before the inner iteration begins. Some performance information is printed while the iteration is in progress. After one step the inner iteration finishes and an output file is generated.

time = 0.1, timestep = 0.1

begin of iteration number: 1
=============================
markElements needed 7.7e-08 seconds
4225 local DOFs
fill-in of assembled matrix: 66049
assemble needed 0.0135029 seconds
=== GeneralizedPCGSolver
=== rate=0.395649, T=0.00769608, TIT=0.000481005, IT=17
solution of discrete system needed 0.0489367 seconds

end of iteration number: 1
=============================
writeFiles needed 0.00137441 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 file can be found here:

• heat.cc
• heat.dat.2d