This tutorial depends on step-26.
This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie.
- Note
- In order to run this program, deal.II must be configured to use the UMFPACK sparse direct solver. Refer to the ReadMe for instructions how to do this.
Introduction
This program shows how to use Runge-Kutta methods to solve a time-dependent problem. It solves a small variation of the heat equation discussed first in step-26 but, since the purpose of this program is only to demonstrate using more advanced ways to interface with deal.II's time stepping algorithms, only solves a simple problem on a uniformly refined mesh.
- Note
- At the end of the day, time stepping is a problem that is only solved efficiently if you use adaptive time step selection and error control. This is implemented in many external libraries, but not in deal.II itself. You may want to consider looking at step-86 for a worked out example. In the meantime, this program is deprecated and will be removed after the deal.II 9.6 release.
Problem statement
In this example, we solve the one-group time-dependent diffusion approximation of the neutron transport equation (see step-28 for the time-independent multigroup diffusion). This is a model for how neutrons move around highly scattering media, and consequently it is a variant of the time-dependent diffusion equation – which is just a different name for the heat equation discussed in step-26, plus some extra terms. We assume that the medium is not fissible and therefore, the neutron flux satisfies the following equation:
\begin{eqnarray*}
\frac{1}{v}\frac{\partial \phi(x,t)}{\partial t} = \nabla \cdot D(x) \nabla \phi(x,t)
- \Sigma_a(x) \phi(x,t) + S(x,t)
\end{eqnarray*}
augmented by appropriate boundary conditions. Here, \(v\) is the velocity of neutrons (for simplicity we assume it is equal to 1 which can be achieved by simply scaling the time variable), \(D\) is the diffusion coefficient, \(\Sigma_a\) is the absorption cross section, and \(S\) is a source. Because we are only interested in the time dependence, we assume that \(D\) and \(\Sigma_a\) are constant.
Since this program only intends to demonstrate how to use advanced time stepping algorithms, we will only look for the solutions of relatively simple problems. Specifically, we are looking for a solution on a square domain \([0,b]\times[0,b]\) of the form
\begin{eqnarray*}
\phi(x,t) = A\sin(\omega t)(bx-x^2).
\end{eqnarray*}
By using quadratic finite elements, we can represent this function exactly at any particular time, and all the error will be due to the time discretization. We do this because it is then easy to observe the order of convergence of the various time stepping schemes we will consider, without having to separate spatial and temporal errors.
We impose the following boundary conditions: homogeneous Dirichlet for \(x=0\) and \(x=b\) and homogeneous Neumann conditions for \(y=0\) and \(y=b\). We choose the source term so that the corresponding solution is in fact of the form stated above:
\begin{eqnarray*}
S=A\left(\frac{1}{v}\omega \cos(\omega t)(bx -x^2) + \sin(\omega t)
\left(\Sigma_a (bx-x^2)+2D\right) \right).
\end{eqnarray*}
Because the solution is a sine in time, we know that the exact solution satisfies \(\phi\left(x,\frac{\pi}{\omega}\right) = 0\). Therefore, the error at time \(t=\frac{\pi}{\omega}\) is simply the norm of the numerical solution, i.e., \(\|e(\cdot,t=\frac{\pi}{\omega})\|_{L_2} = \|\phi_h(\cdot,t=\frac{\pi}{\omega})\|_{L_2}\), and is particularly easily evaluated. In the code, we evaluate the \(l_2\) norm of the vector of nodal values of \(\phi_h\) instead of the \(L_2\) norm of the associated spatial function, since the former is simpler to compute; however, on uniform meshes, the two are just related by a constant and we can consequently observe the temporal convergence order with either.
Runge-Kutta methods
The Runge-Kutta methods implemented in deal.II assume that the equation to be solved can be written as:
\begin{eqnarray*}
\frac{dy}{dt} = g(t,y).
\end{eqnarray*}
On the other hand, when using finite elements, discretized time derivatives always result in the presence of a mass matrix on the left hand side. This can easily be seen by considering that if the solution vector \(y(t)\) in the equation above is in fact the vector of nodal coefficients \(U(t)\) for a variable of the form
\begin{eqnarray*}
u_h(x,t) = \sum_j U_j(t) \varphi_j(x)
\end{eqnarray*}
with spatial shape functions \(\varphi_j(x)\), then multiplying an equation of the form
\begin{eqnarray*}
\frac{\partial u(x,t)}{\partial t} = q(t,u(x,t))
\end{eqnarray*}
by test functions, integrating over \(\Omega\), substituting \(u\rightarrow u_h\) and restricting the test functions to the \(\varphi_i(x)\) from above, then this spatially discretized equation has the form
\begin{eqnarray*}
M\frac{dU}{dt} = f(t,U),
\end{eqnarray*}
where \(M\) is the mass matrix and \(f(t,U)\) is the spatially discretized version of \(q(t,u(x,t))\) (where \(q\) is typically the place where spatial derivatives appear, but this is not of much concern for the moment given that we only consider time derivatives). In other words, this form fits the general scheme above if we write
\begin{eqnarray*}
\frac{dy}{dt} = g(t,y) = M^{-1}f(t,y).
\end{eqnarray*}
Runke-Kutta methods are time stepping schemes that approximate \(y(t_n)\approx
y_{n}\) through a particular one-step approach. They are typically written in the form
\begin{eqnarray*}
y_{n+1} = y_n + \sum_{i=1}^s b_i k_i
\end{eqnarray*}
where for the form of the right hand side above
\begin{eqnarray*}
k_i = \Delta t \, M^{-1} f\left(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j\right).
\end{eqnarray*}
Here \(a_{ij}\), \(b_i\), and \(c_i\) are known coefficients that identify which particular Runge-Kutta scheme you want to use, and \(\Delta t=t_{n+1}-t_n\) is the time step used. Different time stepping methods of the Runge-Kutta class differ in the number of stages \(s\) and the values they use for the coefficients \(a_{ij}\), \(b_i\), and \(c_i\) but are otherwise easy to implement since one can look up tabulated values for these coefficients. (These tables are often called Butcher tableaus.)
At the time of the writing of this tutorial, the methods implemented in deal.II can be divided in three categories:
-
Explicit Runge-Kutta; in order for a method to be explicit, it is necessary that in the formula above defining \(k_i\), \(k_i\) does not appear on the right hand side. In other words, these methods have to satisfy \(a_{ii}=0, i=1,\ldots,s\).
-
Embedded (or adaptive) Runge-Kutta; we will discuss their properties below.
-
Implicit Runge-Kutta; this class of methods require the solution of a possibly nonlinear system the stages \(k_i\) above, i.e., they have \(a_{ii}\neq 0\) for at least one of the stages \(i=1,\ldots,s\).
Many well known time stepping schemes that one does not typically associate with the names Runge or Kutta can in fact be written in a way so that they, too, can be expressed in these categories. They oftentimes represent the lowest-order members of these families; one example is the simple explicit Euler method.
Explicit Runge-Kutta methods
These methods, only require a function to evaluate \(M^{-1}f(t,y)\) but not (as implicit methods) to solve an equation that involves \(f(t,y)\) for \(y\). As all explicit time stepping methods, they become unstable when the time step chosen is too large.
Well known methods in this class include forward Euler, third order Runge-Kutta, and fourth order Runge-Kutta (often abbreviated as RK4).
Embedded Runge-Kutta methods
These methods use both a lower and a higher order method to estimate the error and decide if the time step needs to be shortened or can be increased. The term "embedded" refers to the fact that the lower-order method does not require additional evaluates of the function \(M^{-1}f(\cdot,\cdot)\) but reuses data that has to be computed for the high order method anyway. It is, in other words, essentially free, and we get the error estimate as a side product of using the higher order method.
This class of methods include Heun-Euler, Bogacki-Shampine, Dormand-Prince (ode45 in Matlab and often abbreviated as RK45 to indicate that the lower and higher order methods used here are 4th and 5th order Runge-Kutta methods, respectively), Fehlberg, and Cash-Karp.
At the time of the writing, only embedded explicit methods have been implemented.
Implicit Runge-Kutta methods
Implicit methods require the solution of (possibly nonlinear) systems of the form \(\alpha y = f(t,y)\) for \(y\) in each (sub-)timestep. Internally, this is done using a Newton-type method and, consequently, they require that the user provide functions that can evaluate \(M^{-1}f(t,y)\) and \(\left(I-\tau M^{-1} \frac{\partial f}{\partial y}\right)^{-1}\) or equivalently \(\left(M - \tau \frac{\partial f}{\partial y}\right)^{-1} M\).
The particular form of this operator results from the fact that each Newton step requires the solution of an equation of the form
\begin{align*}
\left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y
= -M h(t,y)
\end{align*}
for some (given) \(h(t,y)\). Implicit methods are always stable, regardless of the time step size, but too large time steps of course affect the accuracy of the solution, even if the numerical solution remains stable and bounded.
Methods in this class include backward Euler, implicit midpoint, Crank-Nicolson, and the two stage SDIRK method (short for "singly diagonally
implicit Runge-Kutta", a term coined to indicate that the diagonal elements \(a_{ii}\) defining the time stepping method are all equal; this property allows for the Newton matrix \(I-\tau M^{-1}\frac{\partial f}{\partial y}\) to be re-used between stages because \(\tau\) is the same every time).
Spatially discrete formulation
By expanding the solution of our model problem as always using shape functions \(\psi_j\) and writing
\begin{eqnarray*}
\phi_h(x,t) = \sum_j U_j(t) \psi_j(x),
\end{eqnarray*}
we immediately get the spatially discretized version of the diffusion equation as
\begin{eqnarray*}
M \frac{dU(t)}{dt}
= -{\cal D} U(t) - {\cal A} U(t) + {\cal S}(t)
\end{eqnarray*}
where
\begin{eqnarray*}
M_{ij} &=& (\psi_i,\psi_j), \\
{\cal D}_{ij} &=& (D\nabla\psi_i,\nabla\psi_j)_\Omega, \\
{\cal A}_{ij} &=& (\Sigma_a\psi_i,\psi_j)_\Omega, \\
{\cal S}_{i}(t) &=& (\psi_i,S(x,t))_\Omega.
\end{eqnarray*}
See also step-24 and step-26 to understand how we arrive here. Boundary terms are not necessary due to the chosen boundary conditions for the current problem. To use the Runge-Kutta methods, we recast this as follows:
\begin{eqnarray*}
f(y) = -{\cal D}y - {\cal A}y + {\cal S}.
\end{eqnarray*}
In the code, we will need to be able to evaluate this function \(f(U)\) along with its derivative,
\begin{eqnarray*}
\frac{\partial f}{\partial y} = -{\cal D} - {\cal A}.
\end{eqnarray*}
Notes on the testcase
To simplify the problem, the domain is two dimensional and the mesh is uniformly refined (there is no need to adapt the mesh since we use quadratic finite elements and the exact solution is quadratic). Going from a two dimensional domain to a three dimensional domain is not very challenging. However if you intend to solve more complex problems where the mesh must be adapted (as is done, for example, in step-26), then it is important to remember the following issues:
-
You will need to project the solution to the new mesh when the mesh is changed. Of course, the mesh used should be the same from the beginning to the end of each time step, a question that arises because Runge-Kutta methods use multiple evaluations of the equations within each time step.
-
You will need to update the mass matrix and its inverse every time the mesh is changed.
The techniques for these steps are readily available by looking at step-26.
The commented program
Include files
The first task as usual is to include the functionality of these well-known deal.II library files and some C++ header files.
#include <deal.II/base/discrete_time.h>
#include <deal.II/base/function.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_handler.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/lac/affine_constraints.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <cmath>
#include <map>
This is the only include file that is new: It includes all the Runge-Kutta methods.
#include <deal.II/base/time_stepping.h>
The next step is like in all previous tutorial programs: We put everything into a namespace of its own and then import the deal.II classes and functions into it.
The Diffusion class
The next piece is the declaration of the main class. Most of the functions in this class are not new and have been explained in previous tutorials. The only interesting functions are evaluate_diffusion() and id_minus_tau_J_inverse(). evaluate_diffusion() evaluates the diffusion equation, \(M^{-1}(f(t,y))\), at a given time and a given \(y\). id_minus_tau_J_inverse() evaluates \(\left(I-\tau
M^{-1} \frac{\partial f(t,y)}{\partial y}\right)^{-1}\) or equivalently \(\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M\) at a given time, for a given \(\tau\) and \(y\). This function is needed when an implicit method is used.
class Diffusion
{
public:
Diffusion();
void run();
private:
void setup_system();
void assemble_system();
double get_source(
const double time,
const Point<2> &point)
const;
const double tau,
void output_results(const double time,
const unsigned int time_step,
The next three functions are the drivers for the explicit methods, the implicit methods, and the embedded explicit methods respectively. The driver function for embedded explicit methods returns the number of steps executed given that it only takes the number of time steps passed as an argument as a hint, but internally computed the optimal time step itself.
const unsigned int n_time_steps,
const double initial_time,
const double final_time);
const unsigned int n_time_steps,
const double initial_time,
const double final_time);
unsigned int
const unsigned int n_time_steps,
const double initial_time,
const double final_time);
const unsigned int fe_degree;
const double diffusion_coefficient;
const double absorption_cross_section;
};
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
We choose quadratic finite elements and we initialize the parameters.
Diffusion::Diffusion()
: fe_degree(2)
, diffusion_coefficient(1. / 30.)
, absorption_cross_section(1.)
, fe(fe_degree)
{}
Diffusion::setup_system
Now, we create the constraint matrix and the sparsity pattern. Then, we initialize the matrices and the solution vector.
void Diffusion::setup_system()
{
dof_handler.distribute_dofs(fe);
1,
constraint_matrix);
constraint_matrix.close();
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
mass_matrix.reinit(sparsity_pattern);
mass_minus_tau_Jacobian.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
}
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Diffusion::assemble_system
In this function, we compute \(-\int D \nabla b_i \cdot \nabla b_j
d\boldsymbol{r} - \int \Sigma_a b_i b_j d\boldsymbol{r}\) and the mass matrix \(\int b_i b_j d\boldsymbol{r}\). The mass matrix is then inverted using a direct solver; the inverse_mass_matrix variable will then store the inverse of the mass matrix so that \(M^{-1}\) can be applied to a vector using the vmult() function of that object. (Internally, UMFPACK does not really store the inverse of the matrix, but its LU factors; applying the inverse matrix is then equivalent to doing one forward and one backward solves with these two factors, which has the same complexity as applying an explicit inverse of the matrix).
void Diffusion::assemble_system()
{
system_matrix = 0.;
mass_matrix = 0.;
const QGauss<2> quadrature_formula(fe_degree + 1);
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0.;
cell_mass_matrix = 0.;
fe_values.reinit(cell);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
cell_matrix(i, j) +=
((-diffusion_coefficient *
fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point)
- absorption_cross_section *
fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point))
* fe_values.JxW(q_point));
cell_mass_matrix(i, j) += fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point) *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraint_matrix.distribute_local_to_global(cell_matrix,
local_dof_indices,
system_matrix);
constraint_matrix.distribute_local_to_global(cell_mass_matrix,
local_dof_indices,
mass_matrix);
}
inverse_mass_matrix.initialize(mass_matrix);
}
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
Diffusion::get_source
In this function, the source term of the equation for a given time and a given point is computed.
double Diffusion::get_source(
const double time,
const Point<2> &point)
const
{
const double intensity = 10.;
const double b = 5.;
const double x = point[0];
return intensity *
(frequency *
std::cos(frequency * time) * (b * x - x * x) +
(absorption_cross_section * (b * x - x * x) +
2. * diffusion_coefficient));
}
static constexpr double PI
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Diffusion::evaluate_diffusion
Next, we evaluate the weak form of the diffusion equation at a given time \(t\) and for a given vector \(y\). In other words, as outlined in the introduction, we evaluate \(M^{-1}(-{\cal D}y - {\cal A}y + {\cal
S})\). For this, we have to apply the matrix \(-{\cal D} - {\cal A}\) (previously computed and stored in the variable system_matrix) to \(y\) and then add the source term which we integrate as we usually do. (Integrating up the solution could be done using VectorTools::create_right_hand_side() if you wanted to save a few lines of code, or wanted to take advantage of doing the integration in parallel.) The result is then multiplied by \(M^{-1}\).
{
tmp = 0.;
system_matrix.vmult(tmp, y);
const QGauss<2> quadrature_formula(fe_degree + 1);
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_source = 0.;
fe_values.reinit(cell);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const double source =
get_source(time, fe_values.quadrature_point(q_point));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_source(i) += fe_values.shape_value(i, q_point) *
source *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraint_matrix.distribute_local_to_global(cell_source,
local_dof_indices,
tmp);
}
inverse_mass_matrix.vmult(value, tmp);
}
@ update_quadrature_points
Transformed quadrature points.
Diffusion::id_minus_tau_J_inverse
We compute \(\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} M\). This is done in several steps:
- compute \(M-\tau \frac{\partial f}{\partial y}\)
- invert the matrix to get \(\left(M-\tau \frac{\partial f}
{\partial y}\right)^{-1}\)
- compute \(tmp=My\)
- compute \(z=\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} tmp =
\left(M-\tau \frac{\partial f}{\partial y}\right)^{-1} My\)
- return z.
const double tau,
{
mass_minus_tau_Jacobian.copy_from(mass_matrix);
mass_minus_tau_Jacobian.add(-tau, system_matrix);
inverse_mass_minus_tau_Jacobian.
initialize(mass_minus_tau_Jacobian);
mass_matrix.vmult(tmp, y);
inverse_mass_minus_tau_Jacobian.vmult(result, tmp);
return result;
}
void initialize(const SparsityPattern &sparsity_pattern)
Diffusion::output_results
The following function then outputs the solution in vtu files indexed by the number of the time step and the name of the time stepping method. Of course, the (exact) result should really be the same for all time stepping method, but the output here at least allows us to compare them.
void Diffusion::output_results(const double time,
const unsigned int time_step,
{
std::string method_name;
switch (method)
{
{
method_name = "forward_euler";
break;
}
{
method_name = "rk3";
break;
}
{
method_name = "rk4";
break;
}
{
method_name = "backward_euler";
break;
}
{
method_name = "implicit_midpoint";
break;
}
{
method_name = "sdirk";
break;
}
{
method_name = "heun_euler";
break;
}
{
method_name = "bogacki_shampine";
break;
}
{
method_name = "dopri";
break;
}
{
method_name = "fehlberg";
break;
}
{
method_name = "cash_karp";
break;
}
default:
{
break;
}
}
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
const std::string filename = "solution_" + method_name + "-" +
".vtu";
std::ofstream output(filename);
data_out.write_vtu(output);
static std::vector<std::pair<double, std::string>> times_and_names;
static std::string method_name_prev = "";
static std::string pvd_filename;
if (method_name_prev != method_name)
{
times_and_names.clear();
method_name_prev = method_name;
pvd_filename = "solution_" + method_name + ".pvd";
}
times_and_names.emplace_back(time, filename);
std::ofstream pvd_output(pvd_filename);
}
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > ×_and_names)
@ RK_CLASSIC_FOURTH_ORDER
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Diffusion::explicit_method
This function is the driver for all the explicit methods. At the top it initializes the time stepping and the solution (by setting it to zero and then ensuring that boundary value and hanging node constraints are respected; of course, with the mesh we use here, hanging node constraints are not in fact an issue). It then calls evolve_one_time_step which performs one time step. Time is stored and incremented through a DiscreteTime object.
For explicit methods, evolve_one_time_step needs to evaluate \(M^{-1}(f(t,y))\), i.e, it needs evaluate_diffusion. Because evaluate_diffusion is a member function, it needs to be bound to this. After each evolution step, we again apply the correct boundary values and hanging node constraints.
Finally, the solution is output every 10 time steps.
const unsigned int n_time_steps,
const double initial_time,
const double final_time)
{
const double time_step =
(final_time - initial_time) / static_cast<double>(n_time_steps);
solution = 0.;
constraint_matrix.distribute(solution);
method);
output_results(initial_time, 0, method);
while (time.is_at_end() == false)
{
explicit_runge_kutta.evolve_one_time_step(
return this->evaluate_diffusion(time, y);
},
time.get_current_time(),
time.get_next_step_size(),
solution);
time.advance_time();
constraint_matrix.distribute(solution);
if (time.get_step_number() % 10 == 0)
output_results(time.get_current_time(),
time.get_step_number(),
method);
}
}
Diffusion::implicit_method
This function is equivalent to explicit_method but for implicit methods. When using implicit methods, we need to evaluate \(M^{-1}(f(t,y))\) and \(\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial
y}\right)^{-1}\) for which we use the two member functions previously introduced.
const unsigned int n_time_steps,
const double initial_time,
const double final_time)
{
const double time_step =
(final_time - initial_time) / static_cast<double>(n_time_steps);
solution = 0.;
constraint_matrix.distribute(solution);
method);
output_results(initial_time, 0, method);
while (time.is_at_end() == false)
{
implicit_runge_kutta.evolve_one_time_step(
return this->evaluate_diffusion(time, y);
},
[
this](
const double time,
const double tau,
const Vector<double> &y) {
return this->id_minus_tau_J_inverse(time, tau, y);
},
time.get_current_time(),
time.get_next_step_size(),
solution);
time.advance_time();
constraint_matrix.distribute(solution);
if (time.get_step_number() % 10 == 0)
output_results(time.get_current_time(),
time.get_step_number(),
method);
}
}
Diffusion::embedded_explicit_method
This function is the driver for the embedded explicit methods. It requires more parameters:
- coarsen_param: factor multiplying the current time step when the error is below the threshold.
- refine_param: factor multiplying the current time step when the error is above the threshold.
- min_delta: smallest time step acceptable.
- max_delta: largest time step acceptable.
- refine_tol: threshold above which the time step is refined.
- coarsen_tol: threshold below which the time step is coarsen.
Embedded methods use a guessed time step. If the error using this time step is too large, the time step will be reduced. If the error is below the threshold, a larger time step will be tried for the next time step. delta_t_guess is the guessed time step produced by the embedded method. In summary, time step size is potentially modified in three ways:
unsigned int Diffusion::embedded_explicit_method(
const unsigned int n_time_steps,
const double initial_time,
const double final_time)
{
const double time_step =
(final_time - initial_time) / static_cast<double>(n_time_steps);
const double coarsen_param = 1.2;
const double refine_param = 0.8;
const double min_delta = 1e-8;
const double max_delta = 10 * time_step;
const double refine_tol = 1e-1;
const double coarsen_tol = 1e-5;
solution = 0.;
constraint_matrix.distribute(solution);
embedded_explicit_runge_kutta(method,
coarsen_param,
refine_param,
min_delta,
max_delta,
refine_tol,
coarsen_tol);
output_results(initial_time, 0, method);
while (time.is_at_end() == false)
{
const double new_time =
embedded_explicit_runge_kutta.evolve_one_time_step(
return this->evaluate_diffusion(time, y);
},
time.get_current_time(),
time.get_next_step_size(),
solution);
time.set_next_step_size(new_time - time.get_current_time());
time.advance_time();
constraint_matrix.distribute(solution);
if (time.get_step_number() % 10 == 0)
output_results(time.get_current_time(),
time.get_step_number(),
method);
time.set_desired_next_step_size(
embedded_explicit_runge_kutta.get_status().delta_t_guess);
}
return time.get_step_number();
}
Diffusion::run
The following is the main function of the program. At the top, we create the grid (a \([0,5]\times [0,5]\) square) and refine it four times to get a mesh that has 16 by 16 cells, for a total of 256. We then set the boundary indicator to 1 for those parts of the boundary where \(x=0\) and \(x=5\).
void Diffusion::run()
{
for (const auto &face : cell->face_iterators())
if (face->at_boundary())
{
face->set_boundary_id(1);
else
face->set_boundary_id(0);
}
void refine_global(const unsigned int times=1)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
Next, we set up the linear systems and fill them with content so that they can be used throughout the time stepping process:
setup_system();
assemble_system();
Finally, we solve the diffusion problem using several of the Runge-Kutta methods implemented in namespace TimeStepping, each time outputting the error at the end time. (As explained in the introduction, since the exact solution is zero at the final time, the error equals the numerical solution and can be computed by just taking the \(l_2\) norm of the solution vector.)
unsigned int n_steps = 0;
const unsigned int n_time_steps = 200;
const double initial_time = 0.;
const double final_time = 10.;
std::cout << "Explicit methods:" << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Forward Euler: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Third order Runge-Kutta: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Fourth order Runge-Kutta: error=" << solution.l2_norm()
<< std::endl;
std::cout << std::endl;
std::cout << "Implicit methods:" << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Backward Euler: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Implicit Midpoint: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Crank-Nicolson: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " SDIRK: error=" << solution.l2_norm()
<< std::endl;
std::cout << std::endl;
std::cout << "Embedded explicit methods:" << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Heun-Euler: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Bogacki-Shampine: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Dopri: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Fehlberg: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Cash-Karp: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
}
}
The main() function
The following main function is similar to previous examples and need not be commented on.
int main()
{
try
{
Step52::Diffusion diffusion;
diffusion.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
The point of this program is less to show particular results, but instead to show how it is done. This we have already demonstrated simply by discussing the code above. Consequently, the output the program yields is relatively sparse and consists only of the console output and the solutions given in VTU format for visualization.
The console output contains both errors and, for some of the methods, the number of steps they performed:
Explicit methods:
Forward Euler: error=1.00883
Third order Runge-Kutta: error=0.000227982
Fourth order Runge-Kutta: error=1.90541e-06
Implicit methods:
Backward Euler: error=1.03428
Implicit Midpoint: error=0.00862702
Crank-Nicolson: error=0.00862675
SDIRK: error=0.0042349
Embedded explicit methods:
Heun-Euler: error=0.0073012
steps performed=284
Bogacki-Shampine: error=0.000408407
steps performed=181
Dopri: error=0.000836695
steps performed=120
Fehlberg: error=0.00248922
steps performed=106
Cash-Karp: error=0.0787735
steps performed=106
As expected the higher order methods give (much) more accurate solutions. We also see that the (rather inaccurate) Heun-Euler method increased the number of time steps in order to satisfy the tolerance. On the other hand, the other embedded methods used a lot less time steps than what was prescribed.
The plain program
#include <fstream>
#include <iostream>
#include <cmath>
#include <map>
namespace Step52
{
class Diffusion
{
public:
Diffusion();
void run();
private:
void setup_system();
void assemble_system();
double get_source(
const double time,
const Point<2> &point)
const;
const double tau,
void output_results(const double time,
const unsigned int time_step,
const unsigned int n_time_steps,
const double initial_time,
const double final_time);
const unsigned int n_time_steps,
const double initial_time,
const double final_time);
unsigned int
const unsigned int n_time_steps,
const double initial_time,
const double final_time);
const unsigned int fe_degree;
const double diffusion_coefficient;
const double absorption_cross_section;
};
Diffusion::Diffusion()
: fe_degree(2)
, diffusion_coefficient(1. / 30.)
, absorption_cross_section(1.)
, fe(fe_degree)
{}
void Diffusion::setup_system()
{
dof_handler.distribute_dofs(fe);
1,
constraint_matrix);
constraint_matrix.close();
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
mass_minus_tau_Jacobian.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
}
void Diffusion::assemble_system()
{
system_matrix = 0.;
const QGauss<2> quadrature_formula(fe_degree + 1);
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_mass_matrix = 0.;
fe_values.reinit(cell);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
((-diffusion_coefficient *
fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point)
- absorption_cross_section *
fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point))
* fe_values.JxW(q_point));
cell_mass_matrix(i, j) += fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point) *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraint_matrix.distribute_local_to_global(cell_matrix,
local_dof_indices,
system_matrix);
constraint_matrix.distribute_local_to_global(cell_mass_matrix,
local_dof_indices,
mass_matrix);
}
inverse_mass_matrix.initialize(mass_matrix);
}
double Diffusion::get_source(
const double time,
const Point<2> &point)
const
{
const double intensity = 10.;
const double x =
point[0];
return intensity *
(frequency *
std::cos(frequency * time) * (
b * x - x * x) +
(absorption_cross_section * (
b * x - x * x) +
2. * diffusion_coefficient));
}
{
tmp = 0.;
system_matrix.vmult(tmp, y);
const QGauss<2> quadrature_formula(fe_degree + 1);
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_source = 0.;
fe_values.reinit(cell);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const double source =
get_source(time, fe_values.quadrature_point(q_point));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_source(i) += fe_values.shape_value(i, q_point) *
source *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraint_matrix.distribute_local_to_global(cell_source,
local_dof_indices,
tmp);
}
inverse_mass_matrix.vmult(value, tmp);
}
const double tau,
{
mass_minus_tau_Jacobian.copy_from(mass_matrix);
mass_minus_tau_Jacobian.add(-tau, system_matrix);
inverse_mass_minus_tau_Jacobian.
initialize(mass_minus_tau_Jacobian);
inverse_mass_minus_tau_Jacobian.
vmult(result, tmp);
return result;
}
void Diffusion::output_results(const double time,
const unsigned int time_step,
{
std::string method_name;
switch (method)
{
{
method_name = "forward_euler";
break;
}
{
method_name = "rk3";
break;
}
{
method_name = "rk4";
break;
}
{
method_name = "backward_euler";
break;
}
{
method_name = "implicit_midpoint";
break;
}
{
method_name = "sdirk";
break;
}
{
method_name = "heun_euler";
break;
}
{
method_name = "bogacki_shampine";
break;
}
{
method_name = "dopri";
break;
}
{
method_name = "fehlberg";
break;
}
{
method_name = "cash_karp";
break;
}
default:
{
break;
}
}
const std::string filename = "solution_" + method_name + "-" +
".vtu";
std::ofstream output(filename);
static std::vector<std::pair<double, std::string>> times_and_names;
static std::string method_name_prev = "";
static std::string pvd_filename;
if (method_name_prev != method_name)
{
times_and_names.clear();
method_name_prev = method_name;
pvd_filename = "solution_" + method_name + ".pvd";
}
times_and_names.emplace_back(time, filename);
std::ofstream pvd_output(pvd_filename);
}
const unsigned int n_time_steps,
const double initial_time,
const double final_time)
{
const double time_step =
(final_time - initial_time) / static_cast<double>(n_time_steps);
solution = 0.;
constraint_matrix.distribute(solution);
method);
output_results(initial_time, 0, method);
while (time.is_at_end() == false)
{
explicit_runge_kutta.evolve_one_time_step(
return this->evaluate_diffusion(time, y);
},
time.get_current_time(),
time.get_next_step_size(),
solution);
time.advance_time();
constraint_matrix.distribute(solution);
if (time.get_step_number() % 10 == 0)
output_results(time.get_current_time(),
time.get_step_number(),
method);
}
}
const unsigned int n_time_steps,
const double initial_time,
const double final_time)
{
const double time_step =
(final_time - initial_time) / static_cast<double>(n_time_steps);
solution = 0.;
constraint_matrix.distribute(solution);
method);
output_results(initial_time, 0, method);
while (time.is_at_end() == false)
{
implicit_runge_kutta.evolve_one_time_step(
return this->evaluate_diffusion(time, y);
},
[
this](
const double time,
const double tau,
const Vector<double> &y) {
return this->id_minus_tau_J_inverse(time, tau, y);
},
time.get_current_time(),
time.get_next_step_size(),
solution);
time.advance_time();
constraint_matrix.distribute(solution);
if (time.get_step_number() % 10 == 0)
output_results(time.get_current_time(),
time.get_step_number(),
method);
}
}
unsigned int Diffusion::embedded_explicit_method(
const unsigned int n_time_steps,
const double initial_time,
const double final_time)
{
const double time_step =
(final_time - initial_time) / static_cast<double>(n_time_steps);
const double coarsen_param = 1.2;
const double refine_param = 0.8;
const double min_delta = 1
e-8;
const double max_delta = 10 * time_step;
const double refine_tol = 1
e-1;
const double coarsen_tol = 1
e-5;
solution = 0.;
constraint_matrix.distribute(solution);
embedded_explicit_runge_kutta(method,
coarsen_param,
refine_param,
min_delta,
max_delta,
refine_tol,
coarsen_tol);
output_results(initial_time, 0, method);
while (time.is_at_end() == false)
{
embedded_explicit_runge_kutta.evolve_one_time_step(
return this->evaluate_diffusion(time, y);
},
time.get_current_time(),
time.get_next_step_size(),
solution);
time.set_next_step_size(new_time - time.get_current_time());
time.advance_time();
constraint_matrix.distribute(solution);
if (time.get_step_number() % 10 == 0)
output_results(time.get_current_time(),
time.get_step_number(),
method);
time.set_desired_next_step_size(
embedded_explicit_runge_kutta.get_status().delta_t_guess);
}
return time.get_step_number();
}
void Diffusion::run()
{
for (const auto &face : cell->face_iterators())
if (face->at_boundary())
{
if ((face->center()[0] == 0.) || (face->center()[0] == 5.))
face->set_boundary_id(1);
else
face->set_boundary_id(0);
}
setup_system();
assemble_system();
unsigned int n_steps = 0;
const unsigned int n_time_steps = 200;
const double initial_time = 0.;
const double final_time = 10.;
std::cout << "Explicit methods:" << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Forward Euler: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Third order Runge-Kutta: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Fourth order Runge-Kutta: error=" << solution.l2_norm()
<< std::endl;
std::cout << std::endl;
std::cout << "Implicit methods:" << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Backward Euler: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Implicit Midpoint: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Crank-Nicolson: error=" << solution.l2_norm()
<< std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " SDIRK: error=" << solution.l2_norm()
<< std::endl;
std::cout << std::endl;
std::cout << "Embedded explicit methods:" << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Heun-Euler: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Bogacki-Shampine: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Dopri: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Fehlberg: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
n_time_steps,
initial_time,
final_time);
std::cout << " Cash-Karp: error=" << solution.l2_norm()
<< std::endl;
std::cout << " steps performed=" << n_steps << std::endl;
}
}
int main()
{
try
{
Step52::Diffusion diffusion;
diffusion.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;
}
void write_vtu(std::ostream &out) const
void set_flags(const FlagType &flags)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
void vmult(Vector< double > &dst, const Vector< double > &src) const
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
1.11.0