]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove step-52
authorDaniel Arndt <arndtd@ornl.gov>
Fri, 26 Jul 2024 15:33:09 +0000 (11:33 -0400)
committerDaniel Arndt <arndtd@ornl.gov>
Fri, 23 Aug 2024 13:30:32 +0000 (09:30 -0400)
14 files changed:
doc/doxygen/tutorial/tutorial.h.in
doc/news/8.1.0-vs-8.2.0.h
doc/news/9.1.1-vs-9.2.0.h
examples/step-26/doc/intro.dox
examples/step-26/doc/results.dox
examples/step-52/CMakeLists.txt [deleted file]
examples/step-52/doc/builds-on [deleted file]
examples/step-52/doc/intro.dox [deleted file]
examples/step-52/doc/kind [deleted file]
examples/step-52/doc/results.dox [deleted file]
examples/step-52/doc/tooltip [deleted file]
examples/step-52/step-52.cc [deleted file]
tests/examples/step-52.diff [deleted file]
tests/examples/step-52.with_umfpack=true.output [deleted file]

index 9b8f63bfc7c6f85282226d580017b42da58f9a57..a8a00514541fd83a356862b8e7fcadb6fb26d039 100644 (file)
  *       </td></tr>
  *
  *   <tr valign="top">
- *       <td>step-52</td>
- *       <td> Solving the time dependent neutron diffusion equation using
- *       Runge-Kutta methods. Explicit and implicit time stepping.
- *       </td></tr>
- *
- *   <tr valign="top">
  *       <td>step-53</td>
  *       <td> Describing the geometry of complex domains and curved boundaries.
  *       </td></tr>
  *   </tr>
  *
  *   <tr valign="top">
- *     <td> Time dependent neutron diffusion equation
- *     </td>
- *     <td>step-52</td>
- *   </tr>
- *
- *   <tr valign="top">
  *     <td> Time dependent fluid structure interaction problems
  *     </td>
  *     <td>step-70</td>
index 708f8e840d8d9abab1456edb2806ec79e4be5969..5a6dddc369c6bb36c913f73936d962f8e1bbdd17 100644 (file)
@@ -163,7 +163,7 @@ inconvenience this causes.
   (Luca Heltai, 2014/09/13)
   </li>
 
-  <li> New: The new tutorial program step-52 explains how to use the
+  <li> New: The new tutorial program step 52 explains how to use the
   new time stepping methods.
   <br>
   (Bruno Turcksin, Damien Lebrun-Grandie, 2014/09/12)
index cbcc535151cb55a9010fae1906f0f2d454e98d25..aa1c1e96512cc06f772d42c9848fdf7f570733ee 100644 (file)
@@ -989,7 +989,7 @@ inconvenience this causes.
  </li>
 
  <li>
-  Fixed: step-26 and step-52 now also output the simulation time
+  Fixed: step-26 and step 52 now also output the simulation time
   in the VTK/VTU output files.
   <br>
   (Wolfgang Bangerth, 2020/03/27)
index 44115e69987f9446d3089b76a80e28be5d798efb..ae3564e24684c09734c472d9e69388e45ab7f3db 100644 (file)
@@ -52,7 +52,7 @@ 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.)
+are interested in playing with higher order methods, take a look at step-86.)
 
 Given this time discretization, space discretization happens as it always
 does, by multiplying with test functions, integrating by parts, and then
index 36cf16c28c0687fbad243190505708ddaa8fdc61..d8b4f2172b829ef6b6f776fe43d0a899165bc827 100644 (file)
@@ -90,7 +90,7 @@ more widely used (and even less accurate) implicit (=backward) Euler method,
 but many other, more accurate methods such as BDF or
 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, if one wanted; a more systematic treatment is also given in step-52.
+program, if one wanted; a more systematic treatment is also given in step-86.
 
 As a general rule, however, one should not be implementing time stepping
 methods by hand, as we do here, for problems that do not require
diff --git a/examples/step-52/CMakeLists.txt b/examples/step-52/CMakeLists.txt
deleted file mode 100644 (file)
index 4d39b41..0000000
+++ /dev/null
@@ -1,55 +0,0 @@
-##
-#  CMake script for the step-52 tutorial program:
-##
-
-# Set the name of the project and target:
-set(TARGET "step-52")
-
-# Declare all source files the target consists of. Here, this is only
-# the one step-X.cc file, but as you expand your project you may wish
-# to add other source files as well. If your project becomes much larger,
-# you may want to either replace the following statement by something like
-#  file(GLOB_RECURSE TARGET_SRC  "source/*.cc")
-#  file(GLOB_RECURSE TARGET_INC  "include/*.h")
-#  set(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
-# or switch altogether to the large project CMakeLists.txt file discussed
-# in the "CMake in user projects" page accessible from the "User info"
-# page of the documentation.
-set(TARGET_SRC
-  ${TARGET}.cc
-  )
-
-# Define the output that should be cleaned:
-set(CLEAN_UP_FILES *.vtu)
-
-# Usually, you will not need to modify anything beyond this point...
-
-cmake_minimum_required(VERSION 3.13.4)
-
-find_package(deal.II 9.7.0
-  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
-  )
-if(NOT ${deal.II_FOUND})
-  message(FATAL_ERROR "\n"
-    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
-    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
-    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
-    )
-endif()
-
-#
-# Are all dependencies fulfilled?
-#
-if(NOT DEAL_II_WITH_UMFPACK)
-  message(FATAL_ERROR "
-Error! This tutorial requires a deal.II library that was configured with the following options:
-    DEAL_II_WITH_UMFPACK = ON
-However, the deal.II library found at ${DEAL_II_PATH} was configured with these options:
-    DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK}
-This conflicts with the requirements."
-    )
-endif()
-
-deal_ii_initialize_cached_variables()
-project(${TARGET})
-deal_ii_invoke_autopilot()
diff --git a/examples/step-52/doc/builds-on b/examples/step-52/doc/builds-on
deleted file mode 100644 (file)
index 1aabbdf..0000000
+++ /dev/null
@@ -1 +0,0 @@
-step-26
diff --git a/examples/step-52/doc/intro.dox b/examples/step-52/doc/intro.dox
deleted file mode 100644 (file)
index dc7b1d4..0000000
+++ /dev/null
@@ -1,260 +0,0 @@
-<br>
-
-<i>This program was contributed by Bruno Turcksin and Damien Lebrun-Grandie.</i>
-
-@note In order to run this program, deal.II must be configured to use
-the UMFPACK sparse direct solver. Refer to the <a
-href="../../readme.html#umfpack">ReadMe</a> for instructions how to do this.
-
-<a name="step-52-Intro"></a>
-<h1>Introduction</h1>
-
-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.
-
-
-
-<h3>Problem statement</h3>
-
-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:
-@f{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)
-@f}
-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
-@f{eqnarray*}{
-\phi(x,t) = A\sin(\omega t)(bx-x^2).
-@f}
-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:
-@f{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).
-@f}
-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.
-
-
-<h3>Runge-Kutta methods</h3>
-
-The Runge-Kutta methods implemented in deal.II assume that the equation to be
-solved can be written as:
-@f{eqnarray*}{
-\frac{dy}{dt} = g(t,y).
-@f}
-On the other hand, when using finite elements, discretized time derivatives always result in the
-presence of a @ref GlossMassMatrix "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
-@f{eqnarray*}{
-  u_h(x,t) = \sum_j U_j(t) \varphi_j(x)
-@f}
-with spatial shape functions $\varphi_j(x)$, then multiplying an equation of
-the form
-@f{eqnarray*}{
-  \frac{\partial u(x,t)}{\partial t} = q(t,u(x,t))
-@f}
-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
-@f{eqnarray*}{
-M\frac{dU}{dt} = f(t,U),
-@f}
-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
-@f{eqnarray*}{
-\frac{dy}{dt} = g(t,y) = M^{-1}f(t,y).
-@f}
-
-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
-@f{eqnarray*}{
-y_{n+1} = y_n + \sum_{i=1}^s b_i k_i
-@f}
-where for the form of the right hand side above
-@f{eqnarray*}{
-k_i = \Delta t \, M^{-1} f\left(t_n+c_ih,y_n+\sum_{j=1}^sa_{ij}k_j\right).
-@f}
-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:
-<ol>
-<li> 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$.
-<li> Embedded (or adaptive) Runge-Kutta; we will discuss their properties below.
-<li> 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$.
-</ol>
-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.
-
-
-<h4>Explicit Runge-Kutta methods</h4>
-
-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).
-
-
-<h4>Embedded Runge-Kutta methods</h4>
-
-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.
-
-
-<h4>Implicit Runge-Kutta methods</h4>
-
-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
-@f{align*}{
-  \left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y
-  = -M h(t,y)
-@f}
-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 <i>accuracy</i> 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).
-
-
-<h3>Spatially discrete formulation</h3>
-
-By expanding the solution of our model problem
-as always using shape functions $\psi_j$ and writing
-@f{eqnarray*}{
-\phi_h(x,t) = \sum_j U_j(t) \psi_j(x),
-@f}
-we immediately get the spatially discretized version of the diffusion equation as
-@f{eqnarray*}{
-  M \frac{dU(t)}{dt}
-  = -{\cal D} U(t) - {\cal A} U(t) + {\cal S}(t)
-@f}
-where
-@f{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.
-@f}
-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:
-@f{eqnarray*}{
-f(y) = -{\cal D}y - {\cal A}y + {\cal S}.
-@f}
-In the code, we will need to be able to evaluate this function $f(U)$ along
-with its derivative,
-@f{eqnarray*}{
-\frac{\partial f}{\partial y} = -{\cal D} - {\cal A}.
-@f}
-
-
-<h3>Notes on the testcase</h3>
-
-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:
-
-<ol>
-<li> 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.
-<li> You will need to update the mass matrix and its inverse every time the
-     mesh is changed.
-</ol>
-The techniques for these steps are readily available by looking at step-26.
diff --git a/examples/step-52/doc/kind b/examples/step-52/doc/kind
deleted file mode 100644 (file)
index 86a44aa..0000000
+++ /dev/null
@@ -1 +0,0 @@
-time dependent
diff --git a/examples/step-52/doc/results.dox b/examples/step-52/doc/results.dox
deleted file mode 100644 (file)
index 03ea68c..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-<h1>Results</h1>
-
-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:
-@code
-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
-@endcode
-
-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.
diff --git a/examples/step-52/doc/tooltip b/examples/step-52/doc/tooltip
deleted file mode 100644 (file)
index 7754b44..0000000
+++ /dev/null
@@ -1 +0,0 @@
-Time-dependent diffusion equation. Time stepping methods.
diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc
deleted file mode 100644 (file)
index 6c61089..0000000
+++ /dev/null
@@ -1,836 +0,0 @@
-/* ------------------------------------------------------------------------
- *
- * SPDX-License-Identifier: LGPL-2.1-or-later
- * Copyright (C) 2014 - 2024 by the deal.II authors
- *
- * This file is part of the deal.II library.
- *
- * Part of the source code is dual licensed under Apache-2.0 WITH
- * LLVM-exception OR LGPL-2.1-or-later. Detailed license information
- * governing the source code and code contributions can be found in
- * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
- *
- * ------------------------------------------------------------------------
- *
- * Authors: Damien Lebrun-Grandie, Bruno Turcksin, 2014
- */
-
-// @sect3{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.
-namespace Step52
-{
-  using namespace dealii;
-
-  // @sect3{The <code>Diffusion</code> 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
-  // <code>evaluate_diffusion()</code> and
-  // <code>id_minus_tau_J_inverse()</code>. <code>evaluate_diffusion()</code>
-  // evaluates the diffusion equation, $M^{-1}(f(t,y))$, at a given time and a
-  // given $y$. <code>id_minus_tau_J_inverse()</code> 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;
-
-    Vector<double> evaluate_diffusion(const double          time,
-                                      const Vector<double> &y) const;
-
-    Vector<double> id_minus_tau_J_inverse(const double          time,
-                                          const double          tau,
-                                          const Vector<double> &y);
-
-    void output_results(const double                     time,
-                        const unsigned int               time_step,
-                        TimeStepping::runge_kutta_method method) const;
-
-    // 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.
-    void explicit_method(const TimeStepping::runge_kutta_method method,
-                         const unsigned int                     n_time_steps,
-                         const double                           initial_time,
-                         const double                           final_time);
-
-    void implicit_method(const TimeStepping::runge_kutta_method method,
-                         const unsigned int                     n_time_steps,
-                         const double                           initial_time,
-                         const double                           final_time);
-
-    unsigned int
-    embedded_explicit_method(const TimeStepping::runge_kutta_method method,
-                             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;
-
-    Triangulation<2> triangulation;
-
-    const FE_Q<2> fe;
-
-    DoFHandler<2> dof_handler;
-
-    AffineConstraints<double> constraint_matrix;
-
-    SparsityPattern sparsity_pattern;
-
-    SparseMatrix<double> system_matrix;
-    SparseMatrix<double> mass_matrix;
-    SparseMatrix<double> mass_minus_tau_Jacobian;
-
-    SparseDirectUMFPACK inverse_mass_matrix;
-
-    Vector<double> solution;
-  };
-
-
-
-  // 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)
-    , dof_handler(triangulation)
-  {}
-
-
-
-  // @sect4{<code>Diffusion::setup_system</code>}
-  // 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);
-
-    VectorTools::interpolate_boundary_values(dof_handler,
-                                             1,
-                                             Functions::ZeroFunction<2>(),
-                                             constraint_matrix);
-    constraint_matrix.close();
-
-    DynamicSparsityPattern dsp(dof_handler.n_dofs());
-    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraint_matrix);
-    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());
-  }
-
-
-
-  // @sect4{<code>Diffusion::assemble_system</code>}
-  // 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 @ref GlossMassMatrix "mass matrix" is then
-  // inverted using a direct solver; the <code>inverse_mass_matrix</code>
-  // variable will then store the inverse of the mass matrix so that
-  // $M^{-1}$ can be applied to a vector using the <code>vmult()</code>
-  // 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);
-
-    FEValues<2> fe_values(fe,
-                          quadrature_formula,
-                          update_values | update_gradients | update_JxW_values);
-
-
-    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
-    const unsigned int n_q_points    = quadrature_formula.size();
-
-    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
-    FullMatrix<double> cell_mass_matrix(dofs_per_cell, dofs_per_cell);
-
-    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 *                // (-D
-                      fe_values.shape_grad(i, q_point) *    //  * grad phi_i
-                      fe_values.shape_grad(j, q_point)      //  * grad phi_j
-                    - absorption_cross_section *            //  -Sigma
-                        fe_values.shape_value(i, q_point) * //  * phi_i
-                        fe_values.shape_value(j, q_point))  //  * phi_j)
-                   * fe_values.JxW(q_point));               // * dx
-                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);
-  }
-
-
-
-  // @sect4{<code>Diffusion::get_source</code>}
-  //
-  // 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 frequency = numbers::PI / 10.;
-    const double b         = 5.;
-    const double x         = point[0];
-
-    return intensity *
-           (frequency * std::cos(frequency * time) * (b * x - x * x) +
-            std::sin(frequency * time) *
-              (absorption_cross_section * (b * x - x * x) +
-               2. * diffusion_coefficient));
-  }
-
-
-
-  // @sect4{<code>Diffusion::evaluate_diffusion</code>}
-  //
-  // 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
-  // <code>system_matrix</code>) 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}$.
-  Vector<double> Diffusion::evaluate_diffusion(const double          time,
-                                               const Vector<double> &y) const
-  {
-    Vector<double> tmp(dof_handler.n_dofs());
-    tmp = 0.;
-    system_matrix.vmult(tmp, y);
-
-    const QGauss<2> quadrature_formula(fe_degree + 1);
-
-    FEValues<2> fe_values(fe,
-                          quadrature_formula,
-                          update_values | update_quadrature_points |
-                            update_JxW_values);
-
-
-    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
-    const unsigned int n_q_points    = quadrature_formula.size();
-
-    Vector<double> cell_source(dofs_per_cell);
-
-    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) * // phi_i(x)
-                                source *                            // * S(x)
-                                fe_values.JxW(q_point);             // * dx
-          }
-
-        cell->get_dof_indices(local_dof_indices);
-
-        constraint_matrix.distribute_local_to_global(cell_source,
-                                                     local_dof_indices,
-                                                     tmp);
-      }
-
-    Vector<double> value(dof_handler.n_dofs());
-    inverse_mass_matrix.vmult(value, tmp);
-
-    return value;
-  }
-
-
-  // @sect4{<code>Diffusion::id_minus_tau_J_inverse</code>}
-  //
-  // 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.
-  Vector<double> Diffusion::id_minus_tau_J_inverse(const double /*time*/,
-                                                   const double          tau,
-                                                   const Vector<double> &y)
-  {
-    SparseDirectUMFPACK inverse_mass_minus_tau_Jacobian;
-
-    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);
-
-    Vector<double> tmp(dof_handler.n_dofs());
-    mass_matrix.vmult(tmp, y);
-
-    Vector<double> result(y);
-    inverse_mass_minus_tau_Jacobian.vmult(result, tmp);
-
-    return result;
-  }
-
-
-
-  // @sect4{<code>Diffusion::output_results</code>}
-  //
-  // 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,
-                                 TimeStepping::runge_kutta_method method) const
-  {
-    std::string method_name;
-
-    switch (method)
-      {
-        case TimeStepping::FORWARD_EULER:
-          {
-            method_name = "forward_euler";
-            break;
-          }
-        case TimeStepping::RK_THIRD_ORDER:
-          {
-            method_name = "rk3";
-            break;
-          }
-        case TimeStepping::RK_CLASSIC_FOURTH_ORDER:
-          {
-            method_name = "rk4";
-            break;
-          }
-        case TimeStepping::BACKWARD_EULER:
-          {
-            method_name = "backward_euler";
-            break;
-          }
-        case TimeStepping::IMPLICIT_MIDPOINT:
-          {
-            method_name = "implicit_midpoint";
-            break;
-          }
-        case TimeStepping::SDIRK_TWO_STAGES:
-          {
-            method_name = "sdirk";
-            break;
-          }
-        case TimeStepping::HEUN_EULER:
-          {
-            method_name = "heun_euler";
-            break;
-          }
-        case TimeStepping::BOGACKI_SHAMPINE:
-          {
-            method_name = "bogacki_shampine";
-            break;
-          }
-        case TimeStepping::DOPRI:
-          {
-            method_name = "dopri";
-            break;
-          }
-        case TimeStepping::FEHLBERG:
-          {
-            method_name = "fehlberg";
-            break;
-          }
-        case TimeStepping::CASH_KARP:
-          {
-            method_name = "cash_karp";
-            break;
-          }
-        default:
-          {
-            break;
-          }
-      }
-
-    DataOut<2> data_out;
-
-    data_out.attach_dof_handler(dof_handler);
-    data_out.add_data_vector(solution, "solution");
-
-    data_out.build_patches();
-
-    data_out.set_flags(DataOutBase::VtkFlags(time, time_step));
-
-    const std::string filename = "solution_" + method_name + "-" +
-                                 Utilities::int_to_string(time_step, 3) +
-                                 ".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);
-    DataOutBase::write_pvd_record(pvd_output, times_and_names);
-  }
-
-
-  // @sect4{<code>Diffusion::explicit_method</code>}
-  //
-  // 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
-  // <code>evolve_one_time_step</code> which performs one time step.
-  // Time is stored and incremented through a DiscreteTime object.
-  //
-  // For explicit methods, <code>evolve_one_time_step</code> needs to
-  // evaluate $M^{-1}(f(t,y))$, i.e, it needs
-  // <code>evaluate_diffusion</code>. Because
-  // <code>evaluate_diffusion</code> is a member function, it needs to
-  // be bound to <code>this</code>. After each evolution step, we
-  // again apply the correct boundary values and hanging node
-  // constraints.
-  //
-  // Finally, the solution is output
-  // every 10 time steps.
-  void Diffusion::explicit_method(const TimeStepping::runge_kutta_method 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);
-
-    TimeStepping::ExplicitRungeKutta<Vector<double>> explicit_runge_kutta(
-      method);
-    output_results(initial_time, 0, method);
-    DiscreteTime time(initial_time, final_time, time_step);
-    while (time.is_at_end() == false)
-      {
-        explicit_runge_kutta.evolve_one_time_step(
-          [this](const double time, const Vector<double> &y) {
-            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);
-      }
-  }
-
-
-
-  // @sect4{<code>Diffusion::implicit_method</code>}
-  // This function is equivalent to <code>explicit_method</code> 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.
-  void Diffusion::implicit_method(const TimeStepping::runge_kutta_method 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);
-
-    TimeStepping::ImplicitRungeKutta<Vector<double>> implicit_runge_kutta(
-      method);
-    output_results(initial_time, 0, method);
-    DiscreteTime time(initial_time, final_time, time_step);
-    while (time.is_at_end() == false)
-      {
-        implicit_runge_kutta.evolve_one_time_step(
-          [this](const double time, const Vector<double> &y) {
-            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);
-      }
-  }
-
-
-
-  // @sect4{<code>Diffusion::embedded_explicit_method</code>}
-  // 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.
-  // <code>delta_t_guess</code> is the guessed time step produced by the
-  // embedded method. In summary, time step size is potentially modified in
-  // three ways:
-  //   - Reducing or increasing time step size within
-  //     TimeStepping::EmbeddedExplicitRungeKutta::evolve_one_time_step().
-  //   - Using the calculated <code>delta_t_guess</code>.
-  //   - Automatically adjusting the step size of the last time step to ensure
-  //     simulation ends precisely at <code>final_time</code>. This adjustment
-  //     is handled inside the DiscreteTime instance.
-  unsigned int Diffusion::embedded_explicit_method(
-    const TimeStepping::runge_kutta_method 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);
-
-    TimeStepping::EmbeddedExplicitRungeKutta<Vector<double>>
-      embedded_explicit_runge_kutta(method,
-                                    coarsen_param,
-                                    refine_param,
-                                    min_delta,
-                                    max_delta,
-                                    refine_tol,
-                                    coarsen_tol);
-    output_results(initial_time, 0, method);
-    DiscreteTime time(initial_time, final_time, time_step);
-    while (time.is_at_end() == false)
-      {
-        const double new_time =
-          embedded_explicit_runge_kutta.evolve_one_time_step(
-            [this](const double time, const Vector<double> &y) {
-              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();
-  }
-
-
-
-  // @sect4{<code>Diffusion::run</code>}
-  //
-  // 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()
-  {
-    GridGenerator::hyper_cube(triangulation, 0., 5.);
-    triangulation.refine_global(4);
-
-    for (const auto &cell : triangulation.active_cell_iterators())
-      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);
-          }
-
-    // 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;
-    explicit_method(TimeStepping::FORWARD_EULER,
-                    n_time_steps,
-                    initial_time,
-                    final_time);
-    std::cout << "   Forward Euler:            error=" << solution.l2_norm()
-              << std::endl;
-
-    explicit_method(TimeStepping::RK_THIRD_ORDER,
-                    n_time_steps,
-                    initial_time,
-                    final_time);
-    std::cout << "   Third order Runge-Kutta:  error=" << solution.l2_norm()
-              << std::endl;
-
-    explicit_method(TimeStepping::RK_CLASSIC_FOURTH_ORDER,
-                    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;
-    implicit_method(TimeStepping::BACKWARD_EULER,
-                    n_time_steps,
-                    initial_time,
-                    final_time);
-    std::cout << "   Backward Euler:           error=" << solution.l2_norm()
-              << std::endl;
-
-    implicit_method(TimeStepping::IMPLICIT_MIDPOINT,
-                    n_time_steps,
-                    initial_time,
-                    final_time);
-    std::cout << "   Implicit Midpoint:        error=" << solution.l2_norm()
-              << std::endl;
-
-    implicit_method(TimeStepping::CRANK_NICOLSON,
-                    n_time_steps,
-                    initial_time,
-                    final_time);
-    std::cout << "   Crank-Nicolson:           error=" << solution.l2_norm()
-              << std::endl;
-
-    implicit_method(TimeStepping::SDIRK_TWO_STAGES,
-                    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_steps = embedded_explicit_method(TimeStepping::HEUN_EULER,
-                                       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_steps = embedded_explicit_method(TimeStepping::BOGACKI_SHAMPINE,
-                                       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_steps = embedded_explicit_method(TimeStepping::DOPRI,
-                                       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_steps = embedded_explicit_method(TimeStepping::FEHLBERG,
-                                       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_steps = embedded_explicit_method(TimeStepping::CASH_KARP,
-                                       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;
-  }
-} // namespace Step52
-
-
-
-// @sect3{The <code>main()</code> function}
-//
-// The following <code>main</code> 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;
-}
diff --git a/tests/examples/step-52.diff b/tests/examples/step-52.diff
deleted file mode 100644 (file)
index cd5758a..0000000
+++ /dev/null
@@ -1,74 +0,0 @@
-669c669
-<     triangulation.refine_global(4);
----
->     triangulation.refine_global(3);
-693c693
-<     unsigned int       n_steps      = 0;
----
->     unsigned int       n_steps [[maybe_unused]] = 0;
-703,704c703
-<     std::cout << "   Forward Euler:            error=" << solution.l2_norm()
-<               << std::endl;
----
->     std::cout << "   Forward Euler: [no error reported]" << std::endl;
-710,711c709
-<     std::cout << "   Third order Runge-Kutta:  error=" << solution.l2_norm()
-<               << std::endl;
----
->     std::cout << "   Third order Runge-Kutta: [no error reported]" << std::endl;
-717,719c715
-<     std::cout << "   Fourth order Runge-Kutta: error=" << solution.l2_norm()
-<               << std::endl;
-<     std::cout << std::endl;
----
->     std::cout << "   Fourth order Runge-Kutta: [no error reported]" << std::endl;
-727,728c723
-<     std::cout << "   Backward Euler:           error=" << solution.l2_norm()
-<               << std::endl;
----
->     std::cout << "   Backward Euler: [no error reported]" << std::endl;
-734,735c729
-<     std::cout << "   Implicit Midpoint:        error=" << solution.l2_norm()
-<               << std::endl;
----
->     std::cout << "   Implicit Midpoint: [no error reported]" << std::endl;
-741,742c735
-<     std::cout << "   Crank-Nicolson:           error=" << solution.l2_norm()
-<               << std::endl;
----
->     std::cout << "   Crank-Nicolson: [no error reported]" << std::endl;
-748,749c741
-<     std::cout << "   SDIRK:                    error=" << solution.l2_norm()
-<               << std::endl;
----
->     std::cout << "   SDIRK: [no error reported]" << std::endl;
-758,760c750
-<     std::cout << "   Heun-Euler:               error=" << solution.l2_norm()
-<               << std::endl;
-<     std::cout << "                   steps performed=" << n_steps << std::endl;
----
->     std::cout << "   Heun-Euler: [no error reported]" << std::endl;
-766,768c756
-<     std::cout << "   Bogacki-Shampine:         error=" << solution.l2_norm()
-<               << std::endl;
-<     std::cout << "                   steps performed=" << n_steps << std::endl;
----
->     std::cout << "   Bogacki-Shampine: [no error reported]" << std::endl;
-774,776c762
-<     std::cout << "   Dopri:                    error=" << solution.l2_norm()
-<               << std::endl;
-<     std::cout << "                   steps performed=" << n_steps << std::endl;
----
->     std::cout << "   Dopri: [no error reported]" << std::endl;
-782,784c768
-<     std::cout << "   Fehlberg:                 error=" << solution.l2_norm()
-<               << std::endl;
-<     std::cout << "                   steps performed=" << n_steps << std::endl;
----
->     std::cout << "   Fehlberg: [no error reported]" << std::endl;
-790,792c774
-<     std::cout << "   Cash-Karp:                error=" << solution.l2_norm()
-<               << std::endl;
-<     std::cout << "                   steps performed=" << n_steps << std::endl;
----
->     std::cout << "   Cash-Karp: [no error reported]" << std::endl;
diff --git a/tests/examples/step-52.with_umfpack=true.output b/tests/examples/step-52.with_umfpack=true.output
deleted file mode 100644 (file)
index 92c3651..0000000
+++ /dev/null
@@ -1,16 +0,0 @@
-Explicit methods:
-   Forward Euler: [no error reported]
-   Third order Runge-Kutta: [no error reported]
-   Fourth order Runge-Kutta: [no error reported]
-Implicit methods:
-   Backward Euler: [no error reported]
-   Implicit Midpoint: [no error reported]
-   Crank-Nicolson: [no error reported]
-   SDIRK: [no error reported]
-
-Embedded explicit methods:
-   Heun-Euler: [no error reported]
-   Bogacki-Shampine: [no error reported]
-   Dopri: [no error reported]
-   Fehlberg: [no error reported]
-   Cash-Karp: [no error reported]

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.