Reference documentation for deal.II version Git ab1cc5b 2017-03-24 06:23:49 -0600
The step-26 tutorial program

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

Note
The material presented here is also discussed in video lecture 29, video lecture 30. (All video lectures are also available here.) (See also video lecture 31.7.)

This program implements the heat equation

\begin{align*} \frac{\partial u(\mathbf x, t)}{\partial t} - \Delta u(\mathbf x, t) &= f(\mathbf x, t), \qquad\qquad && \forall \mathbf x \in \Omega, t\in (0,T), \\ u(\mathbf x, 0) &= u_0(\mathbf x) && \forall \mathbf x \in \Omega, \\ \\ u(\mathbf x, t) &= g(\mathbf x,t) && \forall \mathbf x \in \partial\Omega, t \in (0,T). \end{align*}

In some sense, this equation is simpler than the ones we have discussed in the preceding programs step-23, step-24, step-25, namely the wave equation. This is due to the fact that the heat equation smoothes out the solution over time, and is consequently more forgiving in many regards. For example, when using implicit time stepping methods, we can actually take large time steps, we have less trouble with the small disturbances we introduce through adapting the mesh every few time steps, etc.

Our goal here will be to solve the equations above using the theta-scheme that discretizes the equation in time using the following approach, where we would like \(u^n(\mathbf x)\) to approximate \(u(\mathbf x, t_n)\) at some time \(t_n\):

\begin{align*} \frac{u^n(\mathbf x)-u^{n-1}(\mathbf x)}{k_n} - \left[ (1-\theta)\Delta u^{n-1}(\mathbf x) + \theta\Delta u^n(\mathbf x) \right] &= \left[ (1-\theta)f(\mathbf x, t_{n-1}) + \theta f(\mathbf x, t_n) \right]. \end{align*}

Here, \(k_n=t_n-t_{n-1}\) is the time step size. The theta-scheme generalizes the explicit Euler ( \(\theta=0\)), implicit Euler ( \(\theta=1\)) and Crank-Nicolson ( \(\theta=\frac 12\)) time discretizations. Since the latter has the highest convergence order, we will choose \(\theta=\frac 12\) in the program below, but make it so that playing with this parameter remains simple. (If you are interested in playing with higher order methods, take a look at step-52.)

Given this time discretization, space discretization happens as it always does, by multiplying with test functions, integrating by parts, and then restricting everything to a finite dimensional subspace. This yields the following set of fully discrete equations after multiplying through with \(k_n\):

\begin{align*} M U^n-MU^{n-1} + k_n \left[ (1-\theta)A U^{n-1} + \theta A U^n \right] &= k_n \left[ (1-\theta)F^{n-1} + \theta F^n \right], \end{align*}

where \(M\) is the mass matrix and \(A\) is the stiffness matrix that results from discretizing the Laplacian. Bringing all known quantities to the right hand side yields the linear system we have to solve in every step:

\begin{align*} (M + k_n \theta A) U^n &= MU^{n-1} - k_n (1-\theta)A U^{n-1} + k_n \left[ (1-\theta)F^{n-1} + \theta F^n \right]. \end{align*}

The linear system on the left hand side is symmetric and positive definite, so we should have no trouble solving it with the Conjugate Gradient method.

We can start the iteration above if we have the set of nodal coefficients \(U^0\) at the initial time. Here, we take the ones we get by interpolating the initial values \(u_0(\mathbf x)\) onto the mesh used for the first time step. We will also need to choose a time step; we will here just choose it as fixed, but clearly advanced simulators will want to choose it adaptively. We will briefly come back to this in the results section below.

Adapting meshes for time dependent problems

When solving the wave equation and its variants in the previous few programs, we kept the mesh fixed. Just as for stationary equations, one can make a good case that this is not the smartest approach and that significant savings can be had by adapting the mesh. There are, however, significant difficulties compared to the stationary case. Let us go through them in turn:

What could possibly go wrong? Verifying whether the code is correct

There are a number of things one can typically get wrong when implementing a finite element code. In particular, for time dependent problems, the following are common sources of bugs:

A less common problem is getting the initial conditions wrong because one can typically see that it is wrong by just outputting the first time step. In any case, in order to verify the correctness of the code, it is helpful to have a testing protocol that allows us to verify each of these components separately. This means:

This sounds complicated, but fortunately, for linear partial differential equations without coefficients (or constant coefficients) like the one here, there is a fairly standard protocol that rests on the following observation: if you choose as your domain a square \([0,1]^2\) (or, with slight modifications, a rectangle), then the exact solution can be written as

\begin{align*} u(x,y,t) = a(t) \sin(n_x \pi x) \sin(n_y \pi y) \end{align*}

(with integer constants \(n_x,n_y\)) if only the initial condition, right hand side and boundary values are all of the form \(\sin(n_x \pi x) \sin(n_y \pi y)\) as well. This is due to the fact that the function \(\sin(n_x \pi x) \sin(n_y \pi y)\) is an eigenfunction of the Laplace operator and allows us to compute things like the time factor \(a(t)\) analytically and, consequently, compare with what we get numerically.

As an example, let us consider the situation where we have \(u_0(x,y)=\sin(n_x \pi x) \sin(n_x \pi y)\) and \(f(x,y,t)=0\). With the claim (ansatz) of the form for \(u(x,y,t)\) above, we get that

\begin{align*} \left(\frac{\partial}{\partial t} -\Delta\right) u(x,y,t) &= \left(\frac{\partial}{\partial t} -\Delta\right) a(t) \sin(n_x \pi x) \sin(n_y \pi y) \\ &= \left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y). \end{align*}

For this to be equal to \(f(x,y,t)=0\), we need that

\begin{align*} a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 0 \end{align*}

and due to the initial conditions, \(a(0)=1\). This differential equation can be integrated to yield

\begin{align*} a(t) = - e^{-(n_x^2+n_y^2)\pi^2 t}. \end{align*}

In other words, if the initial condition is a product of sines, then the solution has exactly the same shape of a product of sines that decays to zero with a known time dependence. This is something that is easy to test if you have a sufficiently fine mesh and sufficiently small time step.

What is typically going to happen if you get the time integration scheme wrong (e.g., by having the wrong factors of \(\theta\) or \(k\) in front of the various terms) is that you don't get the right temporal behavior of the solution. Double check the various factors until you get the right behavior. You may also want to verify that the temporal decay rate (as determined, for example, by plotting the value of the solution at a fixed point) does not double or halve each time you double or halve the time step or mesh size. You know that it's not the handling of the boundary conditions or right hand side because these were both zero.

If you have so verified that the time integrator is correct, take the situation where the right hand side is nonzero but the initial conditions are zero: \(u_0(x,y)=0\) and \(f(x,y,t)=\sin(n_x \pi x) \sin(n_x \pi y)\). Again,

\begin{align*} \left(\frac{\partial}{\partial t} -\Delta\right) u(x,y,t) &= \left(\frac{\partial}{\partial t} -\Delta\right) a(t) \sin(n_x \pi x) \sin(n_y \pi y) \\ &= \left(a'(t) + (n_x^2+n_y^2)\pi^2 a(t) \right) \sin(n_x \pi x) \sin(n_y \pi y), \end{align*}

and for this to be equal to \(f(x,y,t)\), we need that

\begin{align*} a'(t) + (n_x^2+n_y^2)\pi^2 a(t) = 1 \end{align*}

and due to the initial conditions, \(a(0)=0\). Integrating this equation in time yields

\begin{align*} a(t) = \frac{1}{(n_x^2+n_y^2)\pi^2} \left[ 1 - e^{-(n_x^2+n_y^2)\pi^2 t} \right]. \end{align*}

Again, if you have the wrong factors of \(\theta\) or \(k\) in front of the right hand side terms you will either not get the right temporal behavior of the solution, or it will converge to a maximum value other than \(\frac{1}{(n_x^2+n_y^2)\pi^2}\).

Once we have verified that the time integration and right hand side handling are correct using this scheme, we can go on to verifying that we have the boundary values correct, using a very similar approach.

The testcase

Solving the heat equation on a simple domain with a simple right hand side almost always leads to solutions that are exceedingly boring, since they become very smooth very quickly and then do not move very much any more. Rather, we here solve the equation on the L-shaped domain with zero Dirichlet boundary values and zero initial conditions, but as right hand side we choose

\begin{align*} f(\mathbf x, t) = \left\{ \begin{array}{ll} \chi_1(\mathbf x) & \text{if \(0\le t \le 0.2\tau\) or \(\tau\le t \le 1.2\tau\) or \(2\tau\le t \le 2.2\tau\), etc} \\ \chi_2(\mathbf x) & \text{if \(0.5\le t \le 0.7\tau\) or \(1.5\tau\le t \le 1.7\tau\) or \(2.5\tau\le t \le 2.7\tau\), etc} \\ 0 & \text{otherwise} \end{array} \right. \end{align*}

Here,

\begin{align*} \chi_1(\mathbf x) &= \left\{ \begin{array}{ll} 1 & \text{if \(x>0.5\) and \(y>-0.5\)} \\ 0 & \text{otherwise} \end{array} \right. \\ \chi_2(\mathbf x) &= \left\{ \begin{array}{ll} 1 & \text{if \(x>-0.5\) and \(y>0.5\)} \\ 0 & \text{otherwise} \end{array} \right. \end{align*}

In other words, in every period of lenght \(\tau\), the right hand side first flashes on in domain 1, then off completely, then on in domain 2, then off completely again. This pattern is probably best observed via the little animation of the solution shown in the results section.

If you interpret the heat equation as finding the spatially and temporally variable temperature distribution of a conducting solid, then the test case above corresponds to an L-shaped body where we keep the boundary at zero temperature, and heat alternatingly in two parts of the domain. While heating is in effect, the temperature rises in these places, after which it diffuses and diminishes again. The point of these initial conditions is that they provide us with a solution that has singularities both in time (when sources switch on and off) as well as time (at the reentrant corner as well as at the edges and corners of the regions where the source acts).

The commented program

The program starts with the usual include files, all of which you should have seen before by now:

#include <deal.II/base/utilities.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/matrix_tools.h>
#include <fstream>
#include <iostream>

Then the usual placing of all content of this program into a namespace and the importation of the deal.II namespace into the one we will work in:

namespace Step26
{
using namespace dealii;

The HeatEquation class

The next piece is the declaration of the main class of this program. It follows the well trodden path of previous examples. If you have looked at step-6, for example, the only thing worth noting here is that we need to build two matrices (the mass and Laplace matrix) and keep the current and previous time step's solution. We then also need to store the current time, the size of the time step, and the number of the current time step. The last of the member variables denotes the theta parameter discussed in the introduction that allows us to treat the explicit and implicit Euler methods as well as the Crank-Nicolson method and other generalizations all in one program.

As far as member functions are concerned, the only possible surprise is that the refine_mesh function takes arguments for the minimal and maximal mesh refinement level. The purpose of this is discussed in the introduction.

template<int dim>
class HeatEquation
{
public:
HeatEquation();
void run();
private:
void setup_system();
void solve_time_step();
void output_results() const;
void refine_mesh (const unsigned int min_grid_level,
const unsigned int max_grid_level);
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> old_solution;
Vector<double> system_rhs;
double time;
double time_step;
unsigned int timestep_number;
const double theta;
};

Equation data

In the following classes and functions, we implement the various pieces of data that define this problem (right hand side and boundary values) that are used in this program and for which we need function objects. The right hand side is chosen as discussed at the end of the introduction. For boundary values, we choose zero values, but this is easily changed below.

template<int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide ()
:
Function<dim>(),
period (0.2)
{}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
private:
const double period;
};
template<int dim>
double RightHandSide<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
(void) component;
Assert (component == 0, ExcIndexRange(component, 0, 1));
Assert (dim == 2, ExcNotImplemented());
const double time = this->get_time();
const double point_within_period = (time/period - std::floor(time/period));
if ((point_within_period >= 0.0) && (point_within_period <= 0.2))
{
if ((p[0] > 0.5) && (p[1] > -0.5))
return 1;
else
return 0;
}
else if ((point_within_period >= 0.5) && (point_within_period <= 0.7))
{
if ((p[0] > -0.5) && (p[1] > 0.5))
return 1;
else
return 0;
}
else
return 0;
}
template<int dim>
class BoundaryValues : public Function<dim>
{
public:
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template<int dim>
double BoundaryValues<dim>::value (const Point<dim> &/ *p* /,
const unsigned int component) const
{
(void) component;
Assert (component == 0, ExcIndexRange(component, 0, 1));
return 0;
}

The HeatEquation implementation

It is time now for the implementation of the main class. Let's start with the constructor which selects a linear element, a time step constant at 1/500 (remember that one period of the source on the right hand side was set to 0.2 above, so we resolve each period with 100 time steps) and chooses the Crank Nicolson method by setting \(\theta=1/2\).

template<int dim>
HeatEquation<dim>::HeatEquation ()
:
fe(1),
dof_handler(triangulation),
time (0.0),
time_step(1. / 500),
timestep_number (0),
theta(0.5)
{}

HeatEquation::setup_system

The next function is the one that sets up the DoFHandler object, computes the constraints, and sets the linear algebra objects to their correct sizes. We also compute the mass and Laplace matrix here by simply calling two functions in the library.

Note that we do not take the hanging node constraints into account when assembling the matrices (both functions have a ConstraintMatrix argument that defaults to an empty object). This is because we are going to condense the constraints in run() after combining the matrices for the current time-step.

template<int dim>
void HeatEquation<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << std::endl
<< "==========================================="
<< std::endl
<< "Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl
<< std::endl;
constraints.clear ();
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
dsp,
constraints,
/ *keep_constrained_dofs = * / true);
sparsity_pattern.copy_from(dsp);
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
system_matrix.reinit(sparsity_pattern);
mass_matrix);
laplace_matrix);
solution.reinit(dof_handler.n_dofs());
old_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}

HeatEquation::solve_time_step

The next function is the one that solves the actual linear system for a single time step. There is nothing surprising here:

template<int dim>
void HeatEquation<dim>::solve_time_step()
{
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
SolverCG<> cg(solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.0);
cg.solve(system_matrix, solution, system_rhs,
preconditioner);
constraints.distribute(solution);
std::cout << " " << solver_control.last_step()
<< " CG iterations." << std::endl;
}

HeatEquation::output_results

Neither is there anything new in generating graphical output:

template<int dim>
void HeatEquation<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "U");
data_out.build_patches();
const std::string filename = "solution-"
+ Utilities::int_to_string(timestep_number, 3) +
".vtk";
std::ofstream output(filename.c_str());
data_out.write_vtk(output);
}

HeatEquation::refine_mesh

This function is the interesting part of the program. It takes care of the adaptive mesh refinement. The three tasks this function performs is to first find out which cells to refine/coarsen, then to actually do the refinement and eventually transfer the solution vectors between the two different grids. The first task is simply achieved by using the well-established Kelly error estimator on the solution. The second task is to actually do the remeshing. That involves only basic functions as well, such as the refine_and_coarsen_fixed_fraction that refines those cells with the largest estimated error that together make up 60 per cent of the error, and coarsens those cells with the smallest error that make up for a combined 40 per cent of the error. Note that for problems such as the current one where the areas where something is going on are shifting around, we want to aggressively coarsen so that we can move cells around to where it is necessary.

As already discussed in the introduction, too small a mesh leads to too small a time step, whereas too large a mesh leads to too little resolution. Consequently, after the first two steps, we have two loops that limit refinement and coarsening to an allowable range of cells:

template <int dim>
void HeatEquation<dim>::refine_mesh (const unsigned int min_grid_level,
const unsigned int max_grid_level)
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
solution,
estimated_error_per_cell);
estimated_error_per_cell,
0.6, 0.4);
if (triangulation.n_levels() > max_grid_level)
cell = triangulation.begin_active(max_grid_level);
cell != triangulation.end(); ++cell)
cell->clear_refine_flag ();
cell = triangulation.begin_active(min_grid_level);
cell != triangulation.end_active(min_grid_level); ++cell)
cell->clear_coarsen_flag ();

These two loops above are slightly different but this is easily explained. In the first loop, instead of calling triangulation.end() we may as well have called triangulation.end_active(max_grid_level). The two calls should yield the same iterator since iterators are sorted by level and there should not be any cells on levels higher than on level max_grid_level. In fact, this very piece of code makes sure that this is the case.

As part of mesh refinement we need to transfer the solution vectors from the old mesh to the new one. To this end we use the SolutionTransfer class and we have to prepare the solution vectors that should be transferred to the new grid (we will lose the old grid once we have done the refinement so the transfer has to happen concurrently with refinement). At the point where we call this function, we will have just computed the solution, so we no longer need the old_solution variable (it will be overwritten by the solution just after the mesh may have been refined, i.e., at the end of the time step; see below). In other words, we only need the one solution vector, and we copy it to a temporary object where it is safe from being reset when we further down below call setup_system().

Consequently, we initialize a SolutionTransfer object by attaching it to the old DoF handler. We then prepare the triangulation and the data vector for refinement (in this order).

SolutionTransfer<dim> solution_trans(dof_handler);
Vector<double> previous_solution;
previous_solution = solution;
triangulation.prepare_coarsening_and_refinement();
solution_trans.prepare_for_coarsening_and_refinement(previous_solution);

Now everything is ready, so do the refinement and recreate the DoF structure on the new grid, and finally initialize the matrix structures and the new vectors in the setup_system function. Next, we actually perform the interpolation of the solution from old to new grid. The final step is to apply the hanging node constraints to the solution vector, i.e., to make sure that the values of degrees of freedom located on hanging nodes are so that the solution is continuous. This is necessary since SolutionTransfer only operates on cells locally, without regard to the neighborhoof.

triangulation.execute_coarsening_and_refinement ();
setup_system ();
solution_trans.interpolate(previous_solution, solution);
constraints.distribute (solution);
}

HeatEquation::run

This is the main driver of the program, where we loop over all time steps. At the top of the function, we set the number of initial global mesh refinements and the number of initial cycles of adaptive mesh refinement by repeating the first time step a few times. Then we create a mesh, initialize the various objects we will work with, set a label for where we should start when re-running the first time step, and interpolate the initial solution onto out mesh (we choose the zero function here, which of course we could do in a simpler way by just setting the solution vector to zero). We also output the initial time step once.

Note
If you're an experienced programmer, you may be surprised that we use a goto statement in this piece of code! goto statements are not particularly well liked any more since Edsgar Dijkstra, one of the greats of computer science, wrote a letter in 1968 called "Go To Statement considered harmful" (see here). The author of this code subscribes to this notion whole-heartedly: goto is hard to understand. In fact, deal.II contains virtually no occurrences: excluding code that was essentially transcribed from books and not counting duplicated code pieces, there are 3 locations in about 600,000 lines of code; we also use it in 4 tutorial programs, in exactly the same context as here. Instead of trying to justify the occurrence here, let's first look at the code and we'll come back to the issue at the end of function.
template<int dim>
void HeatEquation<dim>::run()
{
const unsigned int initial_global_refinement = 2;
const unsigned int n_adaptive_pre_refinement_steps = 4;
GridGenerator::hyper_L (triangulation);
triangulation.refine_global (initial_global_refinement);
setup_system();
unsigned int pre_refinement_step = 0;
Vector<double> forcing_terms;
start_time_iteration:
tmp.reinit (solution.size());
forcing_terms.reinit (solution.size());
old_solution);
solution = old_solution;
output_results();

Then we start the main loop until the computed time exceeds our end time of 0.5. The first task is to build the right hand side of the linear system we need to solve in each time step. Recall that it contains the term \(MU^{n-1}-(1-\theta)k_n AU^{n-1}\). We put these terms into the variable system_rhs, with the help of a temporary vector:

while (time <= 0.5)
{
time += time_step;
++timestep_number;
std::cout << "Time step " << timestep_number << " at t=" << time
<< std::endl;
mass_matrix.vmult(system_rhs, old_solution);
laplace_matrix.vmult(tmp, old_solution);
system_rhs.add(-(1 - theta) * time_step, tmp);

The second piece is to compute the contributions of the source terms. This corresponds to the term \(k_n \left[ (1-\theta)F^{n-1} + \theta F^n \right]\). The following code calls VectorTools::create_right_hand_side to compute the vectors \(F\), where we set the time of the right hand side (source) function before we evaluate it. The result of this all ends up in the forcing_terms variable:

RightHandSide<dim> rhs_function;
rhs_function.set_time(time);
rhs_function,
tmp);
forcing_terms = tmp;
forcing_terms *= time_step * theta;
rhs_function.set_time(time - time_step);
rhs_function,
tmp);
forcing_terms.add(time_step * (1 - theta), tmp);

Next, we add the forcing terms to the ones that come from the time stepping, and also build the matrix \(M+k_n\theta A\) that we have to invert in each time step. The final piece of these operations is to eliminate hanging node constrained degrees of freedom from the linear system:

system_rhs += forcing_terms;
system_matrix.copy_from(mass_matrix);
system_matrix.add(theta * time_step, laplace_matrix);
constraints.condense (system_matrix, system_rhs);

There is one more operation we need to do before we can solve it: boundary values. To this end, we create a boundary value object, set the proper time to the one of the current time step, and evaluate it as we have done many times before. The result is used to also set the correct boundary values in the linear system:

{
BoundaryValues<dim> boundary_values_function;
boundary_values_function.set_time(time);
std::map<types::global_dof_index, double> boundary_values;
0,
boundary_values_function,
boundary_values);
system_matrix,
solution,
system_rhs);
}

With this out of the way, all we have to do is solve the system, generate graphical data, and...

solve_time_step();
output_results();

...take care of mesh refinement. Here, what we want to do is (i) refine the requested number of times at the very beginning of the solution procedure, after which we jump to the top to restart the time iteration, (ii) refine every fifth time step after that.

The time loop and, indeed, the main part of the program ends with starting into the next time step by setting old_solution to the solution we have just computed.

if ((timestep_number == 1) &&
(pre_refinement_step < n_adaptive_pre_refinement_steps))
{
refine_mesh (initial_global_refinement,
initial_global_refinement + n_adaptive_pre_refinement_steps);
++pre_refinement_step;
tmp.reinit (solution.size());
forcing_terms.reinit (solution.size());
std::cout << std::endl;
goto start_time_iteration;
}
else if ((timestep_number > 0) && (timestep_number % 5 == 0))
{
refine_mesh (initial_global_refinement,
initial_global_refinement + n_adaptive_pre_refinement_steps);
tmp.reinit (solution.size());
forcing_terms.reinit (solution.size());
}
old_solution = solution;
}
}
}

Now that you have seen what the function does, let us come back to the issue of the goto. In essence, what the code does is something like this:

void run ()
{
initialize;
start_time_iteration:
for (timestep=1...)
{
solve timestep;
if (timestep==1 && not happy with the result)
{
adjust some data structures;
goto start_time_iteration; // simply try again
}
postprocess;
}
}

Here, the condition "happy with the result" is whether we'd like to keep the current mesh or would rather refine the mesh and start over on the new mesh. We could of course replace the use of the goto by the following:

void run ()
{
initialize;
while (true)
{
solve timestep;
if (not happy with the result)
adjust some data structures;
else
break;
}
postprocess;
for (timestep=2...)
{
solve timestep;
postprocess;
}
}

This has the advantage of getting rid of the goto but the disadvantage of having to duplicate the code that implements the "solve timestep" and "postprocess" operations in two different places. This could be countered by putting these parts of the code (sizable chunks in the actual implementation above) into their own functions, but a while(true) loop with a break statement is not really all that much easier to read or understand than a goto.

In the end, one might simply agree that in general goto statements are a bad idea but be pragmatic and state that there may be occasions where they can help avoid code duplication and awkward control flow. This may be one of these places.

The main function

Having made it this far, there is, again, nothing much to discuss for the main function of this program: it looks like all such functions since step-6.

int main()
{
try
{
using namespace dealii;
using namespace Step26;
HeatEquation<2> heat_equation_solver;
heat_equation_solver.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl << exc.what()
<< std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl << "Aborting!"
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}

Results

As in many of the tutorials, the actual output of the program matters less than how we arrived there. Nonetheless, here it is:

===========================================
Number of active cells: 48
Number of degrees of freedom: 65
Time step 1 at t=0.002
7 CG iterations.
===========================================
Number of active cells: 60
Number of degrees of freedom: 81
Time step 1 at t=0.002
7 CG iterations.
===========================================
Number of active cells: 105
Number of degrees of freedom: 136
Time step 1 at t=0.002
7 CG iterations.
[...]
Time step 249 at t=0.498
13 CG iterations.
Time step 250 at t=0.5
14 CG iterations.
===========================================
Number of active cells: 1803
Number of degrees of freedom: 2109

Maybe of more interest is a visualization of the solution and the mesh on which it was computed:

Animation of the solution of step 26.

The movie shows how the two sources switch on and off and how the mesh reacts to this. It is quite obvious that the mesh as is is probably not the best we could come up with. We'll get back to this in the next section.

Possibilities for extensions

There are at least two areas where one can improve this program significantly: adaptive time stepping and a better choice of the mesh.

Adaptive time stepping

Having chosen an implicit time stepping scheme, we are not bound by any CFL-like condition on the time step. Furthermore, because the time scales on which change happens on a given cell in the heat equation are not bound to the cells diameter (unlike the case with the wave equation, where we had a fixed speed of information transport that couples the temporal and spatial scales), we can choose the time step as we please. Or, better, choose it as we deem necessary for accuracy.

Looking at the solution, it is clear that the action does not happen uniformly over time: a lot is changing around the time we switch on a source, things become less dramatic once a source is on for a little while, and we enter a long phase of decline when both sources are off. During these times, we could surely get away with a larger time step than before without sacrificing too much accuracy.

The literature has many suggestions on how to choose the time step size adaptively. Much can be learned, for example, from the way ODE solvers choose their time steps. One can also be inspired by a posteriori error estimators that can, ideally, be written in a way that the consist of a temporal and a spatial contribution to the overall error. If the temporal one is too large, we should choose a smaller time step. Ideas in this direction can be found, for example, in the PhD thesis of a former principal developer of deal.II, Ralf Hartmann, published by the University of Heidelberg, Germany, in 2002.

Better time stepping methods

We here use one of the simpler time stepping methods, namely the second order in time Crank-Nicolson method. However, more accurate methods such as Runge-Kutta methods are available and should be used as they do not represent much additional effort. It is not difficult to implement this for the current program, but a more systematic treatment is also given in step-52.

Better refinement criteria

If you look at the meshes in the movie above, it is clear that they are not particularly well suited to the task at hand. In fact, they look rather random.

There are two factors at play. First, there are some islands where cells have been refined but that are surrounded by non-refined cells (and there are probably also a few occasional coarsened islands). These are not terrible, as they most of the time do not affect the approximation quality of the mesh, but they also don't help because so many of their additional degrees of freedom are in fact constrained by hanging node constraints. That said, this is easy to fix: the Triangulation class takes an argument to its constructor indicating a level of "mesh smoothing". Passing one of many possible flags, this instructs the triangulation to refine some additional cells, or not to refine some cells, so that the resulting mesh does not have these artifacts.

The second problem is more severe: the mesh appears to lag the solution. The underlying reason is that we only adapt the mesh once every fifth time step, and only allow for a single refinement in these cases. Whenever a source switches on, the solution had been very smooth in this area before and the mesh was consequently rather coarse. This implies that the next time step when we refine the mesh, we will get one refinement level more in this area, and five time steps later another level, etc. But this is not enough: first, we should refine immediately when a source switches on (after all, in the current context we at least know what the right hand side is), and we should allow for more than one refinement level. Of course, all of this can be done using deal.II, it just requires a bit of algorithmic thinking in how to make this work!

The plain program

/* ---------------------------------------------------------------------
*
* Copyright (C) 2013 - 2017 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, Texas A&M University, 2013
*/
#include <deal.II/base/utilities.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/matrix_tools.h>
#include <fstream>
#include <iostream>
namespace Step26
{
using namespace dealii;
template<int dim>
class HeatEquation
{
public:
HeatEquation();
void run();
private:
void setup_system();
void solve_time_step();
void output_results() const;
void refine_mesh (const unsigned int min_grid_level,
const unsigned int max_grid_level);
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> old_solution;
Vector<double> system_rhs;
double time;
double time_step;
unsigned int timestep_number;
const double theta;
};
template<int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide ()
:
Function<dim>(),
period (0.2)
{}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
private:
const double period;
};
template<int dim>
double RightHandSide<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
(void) component;
Assert (component == 0, ExcIndexRange(component, 0, 1));
Assert (dim == 2, ExcNotImplemented());
const double time = this->get_time();
const double point_within_period = (time/period - std::floor(time/period));
if ((point_within_period >= 0.0) && (point_within_period <= 0.2))
{
if ((p[0] > 0.5) && (p[1] > -0.5))
return 1;
else
return 0;
}
else if ((point_within_period >= 0.5) && (point_within_period <= 0.7))
{
if ((p[0] > -0.5) && (p[1] > 0.5))
return 1;
else
return 0;
}
else
return 0;
}
template<int dim>
class BoundaryValues : public Function<dim>
{
public:
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template<int dim>
double BoundaryValues<dim>::value (const Point<dim> &/*p*/,
const unsigned int component) const
{
(void) component;
Assert (component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
template<int dim>
HeatEquation<dim>::HeatEquation ()
:
fe(1),
dof_handler(triangulation),
time (0.0),
time_step(1. / 500),
timestep_number (0),
theta(0.5)
{}
template<int dim>
void HeatEquation<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << std::endl
<< "==========================================="
<< std::endl
<< "Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl
<< std::endl;
constraints.clear ();
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
dsp,
constraints,
/*keep_constrained_dofs = */ true);
sparsity_pattern.copy_from(dsp);
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
system_matrix.reinit(sparsity_pattern);
mass_matrix);
laplace_matrix);
solution.reinit(dof_handler.n_dofs());
old_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template<int dim>
void HeatEquation<dim>::solve_time_step()
{
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
SolverCG<> cg(solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.0);
cg.solve(system_matrix, solution, system_rhs,
preconditioner);
constraints.distribute(solution);
std::cout << " " << solver_control.last_step()
<< " CG iterations." << std::endl;
}
template<int dim>
void HeatEquation<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "U");
data_out.build_patches();
const std::string filename = "solution-"
+ Utilities::int_to_string(timestep_number, 3) +
".vtk";
std::ofstream output(filename.c_str());
data_out.write_vtk(output);
}
template <int dim>
void HeatEquation<dim>::refine_mesh (const unsigned int min_grid_level,
const unsigned int max_grid_level)
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
solution,
estimated_error_per_cell);
estimated_error_per_cell,
0.6, 0.4);
if (triangulation.n_levels() > max_grid_level)
cell = triangulation.begin_active(max_grid_level);
cell != triangulation.end(); ++cell)
cell->clear_refine_flag ();
cell = triangulation.begin_active(min_grid_level);
cell != triangulation.end_active(min_grid_level); ++cell)
cell->clear_coarsen_flag ();
SolutionTransfer<dim> solution_trans(dof_handler);
Vector<double> previous_solution;
previous_solution = solution;
triangulation.prepare_coarsening_and_refinement();
solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
triangulation.execute_coarsening_and_refinement ();
setup_system ();
solution_trans.interpolate(previous_solution, solution);
constraints.distribute (solution);
}
template<int dim>
void HeatEquation<dim>::run()
{
const unsigned int initial_global_refinement = 2;
const unsigned int n_adaptive_pre_refinement_steps = 4;
GridGenerator::hyper_L (triangulation);
triangulation.refine_global (initial_global_refinement);
setup_system();
unsigned int pre_refinement_step = 0;
Vector<double> forcing_terms;
start_time_iteration:
tmp.reinit (solution.size());
forcing_terms.reinit (solution.size());
old_solution);
solution = old_solution;
output_results();
while (time <= 0.5)
{
time += time_step;
++timestep_number;
std::cout << "Time step " << timestep_number << " at t=" << time
<< std::endl;
mass_matrix.vmult(system_rhs, old_solution);
laplace_matrix.vmult(tmp, old_solution);
system_rhs.add(-(1 - theta) * time_step, tmp);
RightHandSide<dim> rhs_function;
rhs_function.set_time(time);
rhs_function,
tmp);
forcing_terms = tmp;
forcing_terms *= time_step * theta;
rhs_function.set_time(time - time_step);
rhs_function,
tmp);
forcing_terms.add(time_step * (1 - theta), tmp);
system_rhs += forcing_terms;
system_matrix.copy_from(mass_matrix);
system_matrix.add(theta * time_step, laplace_matrix);
constraints.condense (system_matrix, system_rhs);
{
BoundaryValues<dim> boundary_values_function;
boundary_values_function.set_time(time);
std::map<types::global_dof_index, double> boundary_values;
0,
boundary_values_function,
boundary_values);
system_matrix,
solution,
system_rhs);
}
solve_time_step();
output_results();
if ((timestep_number == 1) &&
(pre_refinement_step < n_adaptive_pre_refinement_steps))
{
refine_mesh (initial_global_refinement,
initial_global_refinement + n_adaptive_pre_refinement_steps);
++pre_refinement_step;
tmp.reinit (solution.size());
forcing_terms.reinit (solution.size());
std::cout << std::endl;
goto start_time_iteration;
}
else if ((timestep_number > 0) && (timestep_number % 5 == 0))
{
refine_mesh (initial_global_refinement,
initial_global_refinement + n_adaptive_pre_refinement_steps);
tmp.reinit (solution.size());
forcing_terms.reinit (solution.size());
}
old_solution = solution;
}
}
}
int main()
{
try
{
using namespace dealii;
using namespace Step26;
HeatEquation<2> heat_equation_solver;
heat_equation_solver.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl << exc.what()
<< std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl << "Aborting!"
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}