]> https://gitweb.dealii.org/ - dealii.git/commitdiff
step-86: A tutorial for PETScWrappers::TimeStepper.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 28 Jun 2024 22:14:59 +0000 (16:14 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 17 Jul 2024 17:58:35 +0000 (11:58 -0600)
examples/step-86/CMakeLists.txt [new file with mode: 0644]
examples/step-86/doc/builds-on [new file with mode: 0644]
examples/step-86/doc/intro.dox [new file with mode: 0644]
examples/step-86/doc/kind [new file with mode: 0644]
examples/step-86/doc/results.dox [new file with mode: 0644]
examples/step-86/doc/tooltip [new file with mode: 0644]
examples/step-86/heat_equation.prm [new file with mode: 0644]
examples/step-86/step-86.cc [new file with mode: 0644]

diff --git a/examples/step-86/CMakeLists.txt b/examples/step-86/CMakeLists.txt
new file mode 100644 (file)
index 0000000..acf8b26
--- /dev/null
@@ -0,0 +1,54 @@
+##
+#  CMake script for the step-86 tutorial program:
+##
+
+# Set the name of the project and target:
+set(TARGET "step-86")
+
+# 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
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+cmake_minimum_required(VERSION 3.13.4)
+
+find_package(deal.II 9.6.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_PETSC OR DEAL_II_PETSC_WITH_COMPLEX) # keep in one line
+  message(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following options:
+    DEAL_II_WITH_PETSC = ON
+    DEAL_II_PETSC_WITH_COMPLEX = OFF
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options:
+    DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC}
+    DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX}
+This conflicts with the requirements."
+    )
+endif()
+
+deal_ii_initialize_cached_variables()
+project(${TARGET})
+deal_ii_invoke_autopilot()
diff --git a/examples/step-86/doc/builds-on b/examples/step-86/doc/builds-on
new file mode 100644 (file)
index 0000000..fa766ce
--- /dev/null
@@ -0,0 +1 @@
+step-26 step-40
diff --git a/examples/step-86/doc/intro.dox b/examples/step-86/doc/intro.dox
new file mode 100644 (file)
index 0000000..7fe39d2
--- /dev/null
@@ -0,0 +1,442 @@
+<br>
+
+<i>
+This program was contributed by Wolfgang Bangerth (Colorado State University),
+Stefano Zampini (King Abdullah University of Science and Technology), and
+Luca Heltai (University of Pisa).
+
+This material is based upon work partially supported by National Science
+Foundation grants OAC-1835673, DMS-1821210, and EAR-1925595.
+</i>
+<br>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+step-26 solved the simple heat equation, one of the prototypical examples
+of time dependent problems:
+@f{align*}{
+  \frac{\partial u(\mathbf x, t)}{\partial t}
+  -
+  \Delta u(\mathbf x, t)
+  &=
+  f(\mathbf x, t),
+  \qquad\qquad &&
+  \forall \mathbf x \in \Omega, t\in (0,T),
+  \\
+  u(\mathbf x, 0) &= u_0(\mathbf x) &&
+  \forall \mathbf x \in \Omega,
+  \\
+  u(\mathbf x, t) &= g(\mathbf x,t) &&
+  \forall \mathbf x \in \partial\Omega, t \in (0,T).
+@f}
+While that program showed a number of advanced techniques such as
+using adaptive mesh refinement for time-dependent problems, it did not address one big issue:
+It hand-rolls its own time stepping scheme, which in that program
+is the simple
+<a href="https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method">Crank-Nicolson</a>
+method with a fixed time step. This is neither accurate nor efficient: We
+should be using a higher-order time stepping algorithm, and we should
+use one of the many ways to efficiently and automatically choose the
+length of the time step in response to the accuracy obtained.
+
+This would of course require quite a lot of development effort -- unless,
+of course, you do what we always advise: You build on what others have
+already done and have likely done in a way far superior to what one can
+do by oneself. In the current case, deal.II has interfaces to two
+such libraries: SUNDIALS, the *SUite of Nonlinear and DIfferential/ALgebraic
+equation Solvers* (and here specifically the Runge-Kutta-type solvers
+wrapped in the SUNDIALS::ARKode class), and PETSc's TS sub-package
+(wrapped in the PETScWrappers::TimeStepper class). In this program, we will
+use the PETSc TS interfaces.
+
+While we're at it with updating step-26, we will also make the program run
+in parallel -- a minor change if you've read step-40, for example.
+
+
+<h3> %Mapping the heat equation onto an ordinary differential equation formulation </h3>
+
+Both the PETSc TS and the SUNDIALS interfaces require that we first write the partial differential equation
+in the form of a system of ordinary differential equations. To this end, let us turn
+around the approach we used in step-26. There, we first discretized in time,
+obtaining a PDE to be solved at each time step that we could then discretize
+using the finite element method. This approach is called the "Rothe method".
+Instead, here, we use what's called the "Method of Lines" where we first
+discretize in space, obtaining a system of ordinary differential equations
+to which we can apply traditional time steppers. (There are some trade-offs
+between these two strategies, principally around using dynamically changing
+meshes; we will get back to this issue later on.)
+
+To get this started, we take the equation above and multiply it by a test
+function $\varphi(\mathbf x)$ and integrate by parts to get a weak form:
+We seek a function $u(\mathbf x, t)$ that for all test functions
+$\varphi \in H^1_0(\Omega)$
+satisfies
+@f{align*}{
+\left(\varphi(\mathbf x),
+  \frac{\partial u(\mathbf x, t)}{\partial t} \right)_\Omega
+  +
+\left(\nabla \varphi(\mathbf x),
+  \nabla u(\mathbf x, t) \right)_\Omega
+  &=
+\left(\varphi(\mathbf x),
+  f(\mathbf x, t) \right)_\Omega,
+  \\
+\left(\varphi(\mathbf x),  u(\mathbf x, 0)\right)_\Omega &=
+\left(\varphi(\mathbf x), u_0(\mathbf x)\right)_\Omega, &&
+  \\
+  u(\mathbf x, t) &= g(\mathbf x,t) &&
+  \forall \mathbf x \in \partial\Omega, t \in (0,T).
+@f}
+(Integration by parts would ordinarily result in boundary terms unless
+one has Dirichlet boundary conditions -- possibly nonzero -- all
+around the boundary. We will assume that this is indeed so herein.)
+
+We then discretize by restricting ourself to finite element functions
+of the form
+@f{align*}{
+u_h(\mathbf x,t) = \sum_j U_j(t) \varphi_j(\mathbf x),
+@f}
+which leads to the problem of finding a function $u_h(\mathbf x, t)$ that for all
+discrete test functions $\varphi \in V_h(\Omega)$ satisfies
+@f{align*}{
+\left(\varphi_i(\mathbf x),
+  \frac{\partial u_h(\mathbf x, t)}{\partial t} \right)_\Omega
+  +
+\left(\nabla \varphi_j(\mathbf x),
+  \nabla u_h(\mathbf x, t) \right)_\Omega
+  &=
+\left(\varphi_i(\mathbf x),
+  f(\mathbf x, t) \right)_\Omega,
+  \\
+\left(\varphi_i(\mathbf x),  u_h(\mathbf x, 0)\right)_\Omega &=
+\left(\varphi_i(\mathbf x), u_0(\mathbf x)\right)_\Omega, &&
+  \\
+  u_h(\mathbf x, t) &= g_h(\mathbf x,t) &&
+  \forall \mathbf x \in \partial\Omega, t \in (0,T),
+@f}
+where $g_h$ is an interpolant of the function $g$ on the boundary.
+
+This equation can be rewritten in matrix form in the usual way, by
+expanding $u_h$ into its coefficients times shape function form,
+pulling the sum over $j$ out of the integrals, and then considering
+that choosing test function $\varphi_i$ leads to the $i$th row
+of the linear system. This then gives us
+@f{align*}{
+  M
+  \frac{\partial U(t)}{\partial t}
+  +
+  AU(t)
+  &=
+  F(t),
+  \\
+  U(0) = U_0,
+@f}
+plus appropriate boundary conditions.
+
+There are now two perspectives on how one should proceed. If we
+were to use the SUNDIALS::ARKode wrappers to solve this linear system,
+we would bring the $AU$ term to the right hand side and consider
+the ODE
+@f{align*}{
+  M
+  \frac{\partial U(t)}{\partial t}
+  &=
+  -
+  AU(t)
+  +
+  F(t),
+@f}
+which matches the form stated in the documentation of SUNDIALS::ARKode.
+In particular, ARKode is able to deal with the fact that the time
+derivative is multiplied by the mass matrix $M$, which is always
+there when using finite elements.
+
+On the other hand, when using the PETScWrappers::TimeStepper class,
+we can solve ODEs that are stated in a general "implicit" form, and in that
+case we simply bring everything to the left hand side and obtain
+@f{align*}{
+  \underbrace{
+    M
+    \frac{\partial U(t)}{\partial t}
+    +
+    AU(t)
+    -
+    F(t)
+  }_{=:R(t,U,\dot U)}
+  =
+  0.
+@f}
+This matches the form $R(t,U,\dot U) = 0$ you can find in the
+documentation of PETScWrappers::TimeStepper if you identify the time
+dependent function $y=y(t)$ used there with our solution vector
+$U(t)$, and our notation $R(t,U,\dot U)$ instead of the $F(t,y,\dot y)$
+used there and which we rename because we want to use $F$ as the right
+hand side vector of the ODE indicating forcing terms.
+
+This program uses the PETScWrappers::TimeStepper class, and so we will
+take the latter viewpoint. (It is worth noting that SUNDIALS also has
+a package that can solve ODEs in implicit form, wrapped by the
+SUNDIALS::IDA class.) In what follows, we will continue to use $U(t)$
+as the function we seek, even though the documentation of the class
+uses $y(t)$.
+
+
+<h3> %Mapping the differential equation formulation to the time stepper</h3>
+
+Having identified how we want to see the problem (namely, as an "implicit"
+ODE), the question is how we describe the problem to the time stepper.
+Conceptually, all of the wrappers for time stepping packages we support
+in deal.II only requires us to provide them with a very limited set of
+operations. Specifically, for the implicit formulation used by
+PETScWrappers::TimeStepper, all we need to implement are functions
+that provide the following:
+- A way, for a given $t,U,\dot U$, to evaluate the residual vector
+  $R(t,U,\dot U)$.
+- A way, for a given $t,U,\dot U, \alpha$, to set up a matrix
+  $J := \dfrac{\partial R}{\partial y} +
+  \alpha \dfrac{\partial R}{\partial \dot y}$. This is often
+  called the "Jacobian" of the implicit function $R$, perhaps with
+  a small abuse of terminology. In the current case, this matrix
+  is $J=A + \alpha M$. If you have read through step-26, it is probably
+  not lost on you that this matrix appears prominently there as well --
+  with $\alpha$ being a multiple of the inverse of the time step (which
+  there we had denoted by $k_n$).
+  Importantly, for the linear problem we consider here, $J$ is a linear
+  combination of matrices that do not depend on $U$.
+- A way to solve a linear system with this matrix $J$.
+
+That's really it. If we can provide these three functions, PETSc will do
+the rest (as would, for example, SUNDIALS::ARKode or, if you prefer
+the implicit form, SUNDIALS::IDA). It will not be
+very difficult to set these things up. In practice, the way this will
+work is that inside the `run()` function, we will set up lambda functions
+that can access the information of the surrounding scopes and that
+return the requested information.
+
+In practice, we often want to provide a fourth function:
+- A callback that is called at the end of each time step
+  and that is provided with the current solution
+  and other information that can be used to "monitor" the progress
+  of computations. One of the ways in which this can be used is to
+  output visualization data every few time steps.
+
+
+<h3> Complication 1: Dirichlet boundary values </h3>
+
+While we like to say that all nodes in a finite element mesh are
+"degrees of freedom", this is not actually true if we have Dirichlet
+boundary conditions: Degrees of "freedom" located along Dirichlet
+boundaries have specific values that are prescribed by the boundary
+conditions and are, consequently, not "free". Moreover, while the form
+@f{align*}{
+  M
+  \frac{\partial U(t)}{\partial t}
+  &=
+  -
+  AU(t)
+  +
+  F(t)
+@f}
+suggests that *all* elements of the vector $U(t)$ satisfy a
+differential equation, this is not actually true for those components
+of $U$ that correspond to boundary nodes; rather, their values are
+prescribed and will in general not satisfy the equation. On second
+thought, you will also find that the same sort of issue happens as
+well with handing nodes: These, too, are not really "free" but
+constrained to values that are derived from the values of neighboring
+nodes.
+
+For the same reason as we will also discuss in the next section, all
+of this is easier using the Rothe method we have used in all previous
+tutorials. There, we end up with a PDE at every time step for which we
+can independently prescribe boundary conditions and hanging node
+constraints, and we then deal with those by modifying the matrix and
+right hand side appropriately. Here, with the method of lines, things
+are slightly more complicated.
+
+Not too complicated, however: Like with the mesh refinement issues of
+the next section, Dirichlet boundary conditions (and constrained
+degrees of freedom in general) are something every PDE solver has to
+deal with, and because the people who write ODE solvers also have PDEs
+in mind, they needed to address these cases too and so the interfaces
+we use are prepared for it. Specifically, what we need to do is mark
+which entries of the solution vector (i.e., which degrees of freedom)
+are "algebraic" -- that is, satisfy an algebraic, rather than a
+differential, equation. The way we will do this is that the ODE
+integrator interface requires us to provide a "callback" function that
+it can call and that needs to return an IndexSet object in which all
+the algebraic degrees of freedom are listed. At the end of each
+solution stage, a second callback then needs to
+provide the ability to take a solution vector and set its constrained
+(algebraic) entries to their correct values; for us, that will mean
+setting boundary values and hanging node constraints correctly.
+
+
+<h3> Complication 2: Mesh refinement </h3>
+
+When stating an ODE in the form
+@f{align*}{
+  M
+  \frac{\partial U(t)}{\partial t}
+  &=
+  -
+  AU(t)
+  +
+  F(t),
+@f}
+or one of the reformulations discussed above, there is an implicit
+assumption that the number of entries in the vector $U$ stays constant
+and that each entry continues to correspond to the same quantity. But
+if you use mesh refinement, this is not the case: The number of unknowns
+will go up or down whenever you refine or coarsen the mesh, and the
+42nd (or any other) degree of freedom may be located at an entirely different
+physical location after mesh refinement than where it was located
+below. In other words, the size of vectors and what individual vector
+entries mean changes when we do mesh refinement. The ODE form we derived above
+after spatial discretization simply ceases to be meaningful at these times.
+
+This was precisely why, in all previous time-dependent tutorial
+programs, we have adopted the Rothe approach. There, one first
+discretizes in time, and obtains a PDE at each time step. This PDE can
+then be discretized anew -- if one really wanted to, with an entirely
+different mesh in each time step, though we typically don't go that
+far. On the other hand, being able to use external ODE integrators
+*is* undoubtedly very useful, and so let us see if we can shoehorn the
+mesh refinement complication into what external ODE integrators
+do. This is, in practice, not as difficult as it may at first sound
+because, perhaps not surprisingly, ODE integrators are written by
+people who want to solve problems like ours, and so they have had to
+deal with the same kinds of complications we are discussing here.
+
+The way we approach the situation from a *conceptual* perspective is
+that we break things into "time slabs". Let's say we want to solve on
+the time interval $[0,T]$, then we break things into slabs
+$[0=\tau_0,\tau_1], [\tau_1,\tau_2], \ldots [\tau_{n-1},\tau_n=T]$
+where the break points satisfy $\tau_{k-1}<\tau_k$.  On each time
+slab, we keep the mesh the same, and so we can call into our time
+integrator. At the end of a time slab, we then save the solution,
+refine the mesh, set up other data structures, and restore the
+solution on the new mesh; then we start the time integrator again at
+the start of the new time slab.  This approach guarantees that for the
+purposes of ODE solvers, we really only ever give them something that
+can rightfully be considered an ODE system. A disadvantage is that we
+typically want to refine or coarsen the mesh relatively frequently (in
+large-scale codes one often chooses to refine and coarsen the mesh
+every 10-20 time steps), and that limits the efficiency of time
+integrators: They gain much of their advantage from being able to
+choose the time step length automatically, but there is often a cost
+associated with starting up; if the slabs are too short, then neither
+the start-up cost nor the benefit of potentially long time steps are
+realized.
+
+In *practice*, good integrators such as those in PETSc TS can deal
+with this transparently. We just have to give them a way to call back
+into our code at the end of each time step
+to ask whether we want to refine the mesh and do some
+prep work; and a second function that the integrator can then call to
+do the actual refinement and interpolate solution vectors from old to
+new mesh. You will see in the `run()` function that none of this is
+difficult to explain to the ODE integrator.
+
+
+<h3> Structure of the code </h3>
+
+Compared to the code structure of step-26, the program the current tutorial
+is principally based on, there are essentially two sets of changes:
+
+- The program runs in parallel via MPI. This may, at first, seem like
+  a major change, but the reality is that it does not show up in a very
+  large number of places. If you compare step-6 (a sequential Laplace
+  solver) with step-40 (its parallel version), you will see that it takes
+  maybe 20 or 30 extra lines of code to make a simple code run in parallel.
+  These are principally related to keeping track which cells and degrees
+  of freedom are locally owned, or are on ghost cells. We will also have
+  to use matrix and vector data types that are MPI-aware, and we will
+  get those from the PETScWrappers namespace given that we are already
+  using PETSc for the time stepping.
+
+- In step-26 (and most other tutorials), the logic that drives the program's
+  execution lives in the `run()` function: You can see the loop over all
+  time steps, and which functions are called where within the loop. Here,
+  however, this is no longer the case. In essence, in `run()`, we create
+  an object of type PETScWrappers::TimeStepper, and after some set-up,
+  we turn over control to that object's PETScWrappers::TimeStepper::solve()
+  function that contains the loop over time steps and the logic that
+  decides how large the time step needs to be, what needs to happen when,
+  etc. In other words, the *details* of the program's logic are no longer
+  visible. Instead, what you have to provide to the PETScWrappers::TimeStepper
+  object is a series of "callbacks": Functions that the time stepper
+  can call whenever appropriate. These callbacks are typically small
+  [lambda functions](https://en.cppreference.com/w/cpp/language/lambda) that,
+  if the functionality required only takes a few lines of code do exactly
+  that or, otherwise, call larger member functions of the main class.
+
+
+<h3> The test case </h3>
+
+The program solves the heat equation, which with all right hand sides,
+initial, and boundary values reads as
+@f{align*}{
+  \frac{\partial u(\mathbf x, t)}{\partial t}
+  -
+  \Delta u(\mathbf x, t)
+  &=
+  f(\mathbf x, t),
+  \qquad\qquad &&
+  \forall \mathbf x \in \Omega, t\in (0,T),
+  \\
+  u(\mathbf x, 0) &= u_0(\mathbf x) &&
+  \forall \mathbf x \in \Omega,
+  \\
+  u(\mathbf x, t) &= g(\mathbf x,t) &&
+  \forall \mathbf x \in \partial\Omega, t \in (0,T).
+@f}
+The right hand side $f$, initial conditions $u_0$, and Dirichlet boundary
+values $g$ are all specified in an input file `heat_equation.prm` in which
+these functions are provided as expressions that are parsed and evaluated
+at run time using the Functions::ParsedFunction<dim> class. The version of
+this file that is distributed with the library uses
+@f{align*}{
+  f(\mathbf x,t) &= 0, \\
+  u_0(\mathbf x) &= 0, \\
+  g(\mathbf x,t) &= \begin{cases}
+    \cos(4\pi t) & \text{if $x=-1$}, \\
+    -\cos(4\pi t) & \text{if $x=1$}, \\
+    0 & \text{otherwise}
+@f}
+but this is easily changed.
+
+The program's input file also contains two sections that control the
+time stepper:
+@code
+  subsection Time stepper
+    subsection Running parameters
+      set final time              = 5
+      set initial step size       = 0.025
+      set initial time            = 0
+      set match final time        = false
+      set maximum number of steps = -1
+      set options prefix          =
+      set solver type             = beuler
+    end
+
+    subsection Error control
+      set absolute error tolerance = -1
+      set relative error tolerance = -1
+      set adaptor type             = none
+      set ignore algebraic lte     = true
+      set maximum step size        = -1
+      set minimum step size        = -1
+    end
+  end
+@endcode
+The first of these two sections describes things such as the end time up to
+which we want to run the program, the initial time step size, and the type
+of the time stepper (where `beuler` indicates "backward Euler"; other
+choices are
+[listed here](https://petsc.org/release/overview/integrator_table/#integrator-table).
+We will play with some of these parameters in the results section.
+As usual when using PETSc solvers, these runtime configuration
+options can always be complemented (or overridden) via command
+line options.
diff --git a/examples/step-86/doc/kind b/examples/step-86/doc/kind
new file mode 100644 (file)
index 0000000..86a44aa
--- /dev/null
@@ -0,0 +1 @@
+time dependent
diff --git a/examples/step-86/doc/results.dox b/examples/step-86/doc/results.dox
new file mode 100644 (file)
index 0000000..5e3ef5d
--- /dev/null
@@ -0,0 +1,254 @@
+<h1>Results</h1>
+
+When you run this program with the input file as is, you get output as follows:
+@code
+
+Number of active cells: 768
+Number of degrees of freedom: 833
+
+Time step 0 at t=0
+     5 linear iterations.
+     8 linear iterations.
+Time step 1 at t=0.025
+     6 linear iterations.
+Time step 2 at t=0.05
+     5 linear iterations.
+     8 linear iterations.
+Time step 3 at t=0.075
+     6 linear iterations.
+Time step 4 at t=0.1
+     6 linear iterations.
+Time step 5 at t=0.125
+     6 linear iterations.
+Time step 6 at t=0.15
+     6 linear iterations.
+Time step 7 at t=0.175
+     6 linear iterations.
+Time step 8 at t=0.2
+     6 linear iterations.
+Time step 9 at t=0.225
+     6 linear iterations.
+
+Adapting the mesh...
+
+Number of active cells: 1050
+Number of degrees of freedom: 1155
+
+Time step 10 at t=0.25
+     5 linear iterations.
+     8 linear iterations.
+Time step 11 at t=0.275
+     5 linear iterations.
+     7 linear iterations.
+
+[...]
+
+Time step 195 at t=4.875
+     6 linear iterations.
+Time step 196 at t=4.9
+     6 linear iterations.
+Time step 197 at t=4.925
+     6 linear iterations.
+Time step 198 at t=4.95
+     6 linear iterations.
+Time step 199 at t=4.975
+     5 linear iterations.
+
+Adapting the mesh...
+
+Number of active cells: 1380
+Number of degrees of freedom: 1547
+
+Time step 200 at t=5
+
+
++---------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start    |      43.2s |            |
+|                                             |            |            |
+| Section                         | no. calls |  wall time | % of total |
++---------------------------------+-----------+------------+------------+
+| assemble implicit Jacobian      |       226 |      9.93s |        23% |
+| implicit function               |       426 |      16.2s |        37% |
+| output results                  |       201 |      9.74s |        23% |
+| set algebraic components        |       200 |    0.0496s |      0.11% |
+| setup system                    |        21 |     0.799s |       1.8% |
+| solve with Jacobian             |       226 |      0.56s |       1.3% |
+| update current constraints      |       201 |      1.53s |       3.5% |
++---------------------------------+-----------+------------+------------+
+
+@endcode
+
+We can generate output for this in the form of a video:
+@htmlonly
+<p align="center">
+  <iframe width="560" height="315" src="https://www.youtube.com/embed/fhJVkcdRksM"
+   frameborder="0"
+   allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
+   allowfullscreen></iframe>
+ </p>
+@endhtmlonly
+The solution here is driven by boundary values (the initial conditions are zero,
+and so is the right hand side of the equation). It takes a little bit of time
+for the boundary values to diffuse into the domain, and so the temperature
+(the solution of the heat equation) in the interior of the domain has a slight
+lag compared to the temperature at the boundary.
+
+
+The more interesting component of this program is how easy it is to play with
+the details of the time stepping algorithm. Recall that the solution above
+is controlled by the following parameters:
+@code
+  subsection Time stepper
+    subsection Running parameters
+      set final time              = 5
+      set initial step size       = 0.025
+      set initial time            = 0
+      set match final time        = false
+      set maximum number of steps = -1
+      set options prefix          =
+      set solver type             = beuler
+    end
+
+    subsection Error control
+      set absolute error tolerance = -1
+      set adaptor type             = none
+      set ignore algebraic lte     = true
+      set maximum step size        = -1
+      set minimum step size        = -1
+      set relative error tolerance = -1
+    end
+  end
+@endcode
+Of particular interest for us here is to set the time stepping algorithm
+and the adaptive time step control method. The latter is set to "none"
+above, but there are
+[several alternative choices](https://petsc.org/release/manualpages/TS/TSAdaptType/)
+for this parameter. For example, we can set parameters as follows,
+@code
+  subsection Time stepper
+    subsection Running parameters
+      set final time              = 5
+      set initial step size       = 0.025
+      set initial time            = 0
+      set match final time        = false
+      set maximum number of steps = -1
+      set options prefix          =
+      set solver type             = bdf
+    end
+
+    subsection Error control
+      set absolute error tolerance = 1e-2
+      set relative error tolerance = 1e-2
+      set adaptor type             = basic
+      set ignore algebraic lte     = true
+      set maximum step size        = 1
+      set minimum step size        = 0.01
+    end
+  end
+@endcode
+What we do here is set the initial time step size to 0.025, and choose relatively
+large absolute and relative error tolerances of 0.01 for the time step size
+adaptation algorithm (for which we choose "basic"). We ask PETSc TS to use a
+[Backward Differentiation Formula (BDF)](https://en.wikipedia.org/wiki/Backward_differentiation_formula)
+method, and we get the following as output:
+@code
+
+===========================================
+Number of active cells: 768
+Number of degrees of freedom: 833
+
+Time step 0 at t=0
+     5 linear iterations.
+     5 linear iterations.
+     5 linear iterations.
+     4 linear iterations.
+     6 linear iterations.
+Time step 1 at t=0.01
+     5 linear iterations.
+     5 linear iterations.
+Time step 2 at t=0.02
+     5 linear iterations.
+     5 linear iterations.
+Time step 3 at t=0.03
+     5 linear iterations.
+     5 linear iterations.
+Time step 4 at t=0.042574
+     5 linear iterations.
+     5 linear iterations.
+Time step 5 at t=0.0588392
+     5 linear iterations.
+     5 linear iterations.
+Time step 6 at t=0.0783573
+     5 linear iterations.
+     7 linear iterations.
+     5 linear iterations.
+     7 linear iterations.
+Time step 7 at t=0.100456
+     5 linear iterations.
+     7 linear iterations.
+     5 linear iterations.
+     7 linear iterations.
+Time step 8 at t=0.124982
+     5 linear iterations.
+     8 linear iterations.
+     5 linear iterations.
+     7 linear iterations.
+Time step 9 at t=0.153156
+     6 linear iterations.
+     5 linear iterations.
+
+[...]
+
+Time step 206 at t=4.96911
+     5 linear iterations.
+     5 linear iterations.
+Time step 207 at t=4.99398
+     5 linear iterations.
+     5 linear iterations.
+Time step 208 at t=5.01723
+
+
++---------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start    |       117s |            |
+|                                             |            |            |
+| Section                         | no. calls |  wall time | % of total |
++---------------------------------+-----------+------------+------------+
+| assemble implicit Jacobian      |       593 |      35.6s |        31% |
+| implicit function               |      1101 |      56.3s |        48% |
+| output results                  |       209 |      12.5s |        11% |
+| set algebraic components        |       508 |     0.156s |      0.13% |
+| setup system                    |        21 |      1.11s |      0.95% |
+| solve with Jacobian             |       593 |      1.97s |       1.7% |
+| update current constraints      |       509 |      4.53s |       3.9% |
++---------------------------------+-----------+------------+------------+
+@endcode
+What is happening here is that apparently PETSc TS is not happy with our
+choice of initial time step size, and after several linear solves has
+reduced it to the minimum we allow it to, 0.01. The following time steps
+then run at a time step size of 0.01 until it decides to make it slightly
+larger again and (apparently) switches to a higher order method that
+requires more linear solves per time step but allows for a larger time
+step closer to our initial choice 0.025 again. It does not quite hit the
+final time of $T=5$ with its time step choices, but we've got only
+ourselves to blame for that by setting
+@code
+      set match final time        = false
+@endcode
+in the input file.
+
+Not all combinations of methods, time step adaptation algorithms, and
+other parameters are valid, but the main messages from the experiment
+above that you should take away are:
+- It would undoubtedly be quite time consuming to implement many of the
+  methods that PETSc TS offers for time stepping -- but with a program
+  such as the one here, we don't need to: We can just select from the
+  many methods PETSc TS already has.
+- Adaptive time step control is difficult; adaptive choice of which
+  method or order to choose is perhaps even more difficult. None of the
+  time dependent programs that came before the current one (say, step-23,
+  step-26, step-31, step-58, and a good number of others) have either.
+  Moreover, while deal.II is good at spatial adaptation of meshes, it
+  is not a library written by experts in time step adaptation, and so will
+  likely not gain this ability either. But, again, it doesn't
+  have to: We can rely on a library written by experts in that area.
+
diff --git a/examples/step-86/doc/tooltip b/examples/step-86/doc/tooltip
new file mode 100644 (file)
index 0000000..aca910f
--- /dev/null
@@ -0,0 +1 @@
+A parallel solver for the heat equation, using better time steppers.
diff --git a/examples/step-86/heat_equation.prm b/examples/step-86/heat_equation.prm
new file mode 100644 (file)
index 0000000..da176a9
--- /dev/null
@@ -0,0 +1,44 @@
+subsection Heat Equation
+  set Initial global refinement      = 4
+  set Maximum delta refinement level = 2
+  set Mesh adaptation frequency      = 10
+
+  subsection Right hand side
+    set Function constants  = 
+    set Function expression = 0
+    set Variable names      = x,y,t
+  end
+  
+  subsection Initial value
+    set Function constants  = 
+    set Function expression = 0
+    set Variable names      = x,y,t
+  end
+  
+  subsection Boundary values
+    set Function constants  = 
+    set Function expression = (if(x<-0.999, 1, 0) + if(x>0.999, -1, 0)) * cos(4*pi*t)
+    set Variable names      = x,y,t
+  end
+
+  subsection Time stepper
+    subsection Running parameters
+      set final time              = 5
+      set initial step size       = 0.025
+      set initial time            = 0
+      set match final time        = false
+      set maximum number of steps = -1
+      set options prefix          = 
+      set solver type             = beuler
+    end
+
+    subsection Error control
+      set absolute error tolerance = -1
+      set relative error tolerance = -1
+      set adaptor type             = none
+      set ignore algebraic lte     = true
+      set maximum step size        = -1
+      set minimum step size        = -1
+    end    
+  end
+end
diff --git a/examples/step-86/step-86.cc b/examples/step-86/step-86.cc
new file mode 100644 (file)
index 0000000..d2149d4
--- /dev/null
@@ -0,0 +1,1145 @@
+/* ------------------------------------------------------------------------
+ *
+ * SPDX-License-Identifier: LGPL-2.1-or-later
+ * Copyright (C) 2000 - 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:
+ *   Wolfgang Bangerth, Colorado State University, 2024
+ *   Stefano Zampini, King Abdullah University of Science and Technology, 2024
+ */
+
+
+// The program starts with the usual prerequisites, all of which you will know
+// by now from either step-26 (the heat equation solver) or step-40 (the
+// parallel Laplace solver):
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/parsed_function.h>
+#include <deal.II/base/parameter_acceptor.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/solution_transfer.h>
+#include <deal.II/grid/grid_generator.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/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+// The only new include file relevant here is the one that provides us with the
+// PETScWrappers::TimerStepper class:
+#include <deal.II/lac/petsc_ts.h>
+
+#include <fstream>
+#include <iostream>
+
+
+namespace Step86
+{
+  using namespace dealii;
+
+  // @sect3{The HeatEquation class}
+  //
+  // At its core, this program's principal structure can be understood quite
+  // easily if you know step-26 (for the heat equation) and step-40 (for how
+  // a parallel solver looks like). It has many of the usual member functions
+  // and member variables that for convenience of documentation we will list
+  // towards the top of the class so that we can document the remainder
+  // separately below.
+  //
+  // We derive the main class from ParameterAcceptor to make dealing with run
+  // time parameters and reading them a parameter file easier. step-60 and
+  // step-70 have already explained how this works.
+  template <int dim>
+  class HeatEquation : public ParameterAcceptor
+  {
+  public:
+    HeatEquation(const MPI_Comm mpi_communicator);
+    void run();
+
+  private:
+    const MPI_Comm mpi_communicator;
+
+    ConditionalOStream pcout;
+    TimerOutput        computing_timer;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+    FE_Q<dim>                                 fe;
+    DoFHandler<dim>                           dof_handler;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+
+    void setup_system(const double time);
+
+    void output_results(const double                      time,
+                        const unsigned int                timestep_number,
+                        const PETScWrappers::MPI::Vector &solution);
+
+    // At this point, we start to deviate from the "manual" way of solving time
+    // dependent problems. High level packages for the solution of ODEs
+    // usually expect the user to provide a function that computes the residual
+    // of the equation and the time derivative of the solution.
+    //
+    // This allows those packages to abstract away the details of how the time
+    // derivative is defined (e.g., backward Euler, Crank-Nicolson, etc.) and
+    // allow the user to provide a uniform interface, irrespective of the time
+    // stepper used. PETSc TS and SUNDIALS are two examples of such packages,
+    // and they require the user to provide C style call-backs to compute
+    // residuals, Jacobians, etc. In deal.II, we wrap these C style libraries
+    // with our own wrappers that use a style closer to c++, and that allows us
+    // to simply define the callbacks via lambda functions. Several of these
+    // lambda functions will simply call member functions to do the actual work.
+    //
+    // To make it clear what we are doing here, we start by defining functions
+    // with the same name of the interface that we will use later on (and
+    // documented both in the introduction of this program as well in the
+    // documentation of the PETScWrappers::TimeStepper class). These are
+    // the functions that compute the right hand side, the "Jacobian matrix"
+    // (see the introduction for how that is defined), and that can solve
+    // a linear system with the Jacobian matrix. At the bottom of the following
+    // block, we also declare the matrix object that will store this Jacobian.
+    //
+    // Note that all of these functions receive the current solution (and,
+    // where necessary, its time derivative) as inputs. This is because *the*
+    // solution vector -- i.e., the variable that stores the current state
+    // of the solution -- is kept inside the time integrator object. This
+    // is useful: If we kept a copy of the solution vector as a member variable
+    // of the current class, one would continuously have to wonder whether it
+    // is still in sync with the version of the solution vector the time
+    // integrator object internally believes is the currently correct
+    // version of this vector. As a consequence, we do not store such a
+    // copy here: Whenever a function requires access to the current value
+    // of the solution, it receives it from the time integrator as a const
+    // argument. The same observation can be made about the variable that stores
+    // the current time, or the current length of the time step, or the number
+    // of time steps performed so far: They are all kept inside the time
+    // stepping object.
+    void implicit_function(const double                      time,
+                           const PETScWrappers::MPI::Vector &solution,
+                           const PETScWrappers::MPI::Vector &solution_dot,
+                           PETScWrappers::MPI::Vector       &residual);
+
+    void
+    assemble_implicit_jacobian(const double                      time,
+                               const PETScWrappers::MPI::Vector &solution,
+                               const PETScWrappers::MPI::Vector &solution_dot,
+                               const double                      shift);
+
+    void solve_with_jacobian(const PETScWrappers::MPI::Vector &src,
+                             PETScWrappers::MPI::Vector       &residual);
+
+    PETScWrappers::MPI::SparseMatrix jacobian_matrix;
+
+
+    // In this tutorial program, similar to what we did in step-26, we
+    // want to adapt the mesh at regular time intervals. However, if
+    // we are using an external time stepper, we need to make sure
+    // that the time stepper is aware of the mesh changes -- see the
+    // discussion on this topic in the introduction. In particular,
+    // the time stepper must support the fact that the number of
+    // degrees of freedom may change, and it must support the fact
+    // that all internal vectors (e.g., the solution vector and all
+    // intermediate stages used to compute the current time
+    // derivative) will need to be transferred to the new mesh. For
+    // reasons that will be discussed in the implementation below, we
+    // split the mesh refinement operation into two functions: The
+    // first will mark which cells to refine, based on a given
+    // solution vector from which we can compute error indicators; the
+    // second one will do the actual mesh refinement and transfer a
+    // set of vectors from the old to the new mesh.
+    void prepare_for_coarsening_and_refinement(
+      const PETScWrappers::MPI::Vector &solution);
+
+    void interpolate(const std::vector<PETScWrappers::MPI::Vector> &all_in,
+                     std::vector<PETScWrappers::MPI::Vector>       &all_out);
+
+    // As also discussed in the introduction, we also have to deal with
+    // "algebraic" solution components, i.e., degrees of freedom that are not
+    // really free but instead have values determined either by boundary
+    // conditions or by hanging node constraints. While the values of the
+    // boundary conditions can change from time step to time step (because the
+    // function $g(\mathbf x,t)$ may indeed depend on time), hanging node
+    // constraints remain the same as long as the mesh remains the same. As a
+    // consequence, we will keep an AffineConstraints object that stores the
+    // hanging node constraints and that is only updated when the mesh changes,
+    // and then an AffineConstraints object `current_constraints` that we will
+    // initialize with the hanging node constraints and then add the constraints
+    // due to Dirichlet boundary values at a specified time. This evaluation of
+    // boundary values and combining of constraints happens in the
+    // `update_current_constraints()` function.
+    //
+    // At one place in the program, we will also need an object that constrains
+    // the same degrees of freedom, but with zero values even if the boundary
+    // values for the solution are non-zero; we will keep this modified set of
+    // constraints in `homogeneous_constraints`.
+    AffineConstraints<double> hanging_nodes_constraints;
+    AffineConstraints<double> current_constraints;
+    AffineConstraints<double> homogeneous_constraints;
+
+    void update_current_constraints(const double time);
+
+    // The remainder of the class is simply consumed with objects that
+    // describe either the functioning of the time stepper, when and how to
+    // do mesh refinement, and the objects that describe right hand side,
+    // initial conditions, and boundary conditions for the PDE.
+    // Since we already derive our class from ParameterAcceptor, we exploit its
+    // facilities to parse also the parameters of the functions that define the
+    // initial value, the right hand side, and the boundary values. We therefore
+    // create three ParameterAcceptorProxy objects, which wrap the actual
+    // ParsedFunction class into objects that ParameterAcceptor can handle.
+    PETScWrappers::TimeStepperData time_stepper_data;
+
+    unsigned int initial_global_refinement;
+    unsigned int max_delta_refinement_level;
+    unsigned int mesh_adaptation_frequency;
+
+    ParameterAcceptorProxy<Functions::ParsedFunction<dim>>
+      right_hand_side_function;
+    ParameterAcceptorProxy<Functions::ParsedFunction<dim>>
+      initial_value_function;
+    ParameterAcceptorProxy<Functions::ParsedFunction<dim>>
+      boundary_values_function;
+  };
+
+
+  // @sect4{The HeatEquation constructor}
+  //
+  // The constructor is responsible for initializing all member variables of
+  // the class. This is relatively straightforward, and includes setting up
+  // parameters to be set upon reading from an input file via the
+  // ParameterAcceptor mechanism previously detailed in step-60 and step-70.
+  template <int dim>
+  HeatEquation<dim>::HeatEquation(const MPI_Comm mpi_communicator)
+    : ParameterAcceptor("/Heat Equation/")
+    , mpi_communicator(mpi_communicator)
+    , pcout(std::cout,
+            (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+    , computing_timer(mpi_communicator,
+                      pcout,
+                      TimerOutput::summary,
+                      TimerOutput::wall_times)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , fe(1)
+    , dof_handler(triangulation)
+    , time_stepper_data("",
+                        "beuler",
+                        /* start time */ 0.0,
+                        /* end time */ 1.0,
+                        /* initial time step */ 0.025)
+    , initial_global_refinement(5)
+    , max_delta_refinement_level(2)
+    , mesh_adaptation_frequency(0)
+    , right_hand_side_function("/Heat Equation/Right hand side")
+    , initial_value_function("/Heat Equation/Initial value")
+    , boundary_values_function("/Heat Equation/Boundary values")
+  {
+    enter_subsection("Time stepper");
+    {
+      enter_my_subsection(this->prm);
+      {
+        time_stepper_data.add_parameters(this->prm);
+      }
+      leave_my_subsection(this->prm);
+    }
+    leave_subsection();
+
+    add_parameter("Initial global refinement",
+                  initial_global_refinement,
+                  "Number of times the mesh is refined globally before "
+                  "starting the time stepping.");
+    add_parameter("Maximum delta refinement level",
+                  max_delta_refinement_level,
+                  "Maximum number of local refinement levels.");
+    add_parameter("Mesh adaptation frequency",
+                  mesh_adaptation_frequency,
+                  "When to adapt the mesh.");
+  }
+
+
+
+  // @sect4{The HeatEquation::setup_system() function}
+  //
+  // This function is not very different from what we do in many other programs
+  // (including step-26). We enumerate degrees of freedom, output some
+  // information about then, build constraint objects (recalling that we
+  // put hanging node constraints into their separate object), and then
+  // also build an AffineConstraint object that contains both the hanging
+  // node constraints as well as constraints corresponding to zero Dirichlet
+  // boundary conditions. This last object is needed since we impose the
+  // constraints through algebraic equations. While technically it would be
+  // possible to use the time derivative of the boundary function as a boundary
+  // conditions for the time derivative of the solution, this is not done here.
+  // Instead, we impose the boundary conditions through algebraic equations, and
+  // therefore the time derivative of the boundary conditions is not part of
+  // the algebraic system, and we need zero boundary conditions on the time
+  // derivative of the solution when computing the residual. We use the
+  // `homogeneous_constraints` object for this purpose.
+  //
+  // Finally, we create the actual non-homogeneous `current_constraints` by
+  // calling `update_current_constraints). These are also used during the
+  // assembly and during the residual evaluation.
+  template <int dim>
+  void HeatEquation<dim>::setup_system(const double time)
+  {
+    TimerOutput::Scope t(computing_timer, "setup system");
+
+    dof_handler.distribute_dofs(fe);
+    pcout << std::endl
+          << "Number of active cells: " << triangulation.n_active_cells()
+          << std::endl
+          << "Number of degrees of freedom: " << dof_handler.n_dofs()
+          << std::endl
+          << std::endl;
+
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    locally_relevant_dofs =
+      DoFTools::extract_locally_relevant_dofs(dof_handler);
+
+
+    hanging_nodes_constraints.clear();
+    hanging_nodes_constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints(dof_handler,
+                                            hanging_nodes_constraints);
+    hanging_nodes_constraints.make_consistent_in_parallel(locally_owned_dofs,
+                                                          locally_relevant_dofs,
+                                                          mpi_communicator);
+    hanging_nodes_constraints.close();
+
+
+    homogeneous_constraints.clear();
+    homogeneous_constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
+    homogeneous_constraints.merge(hanging_nodes_constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(),
+                                             homogeneous_constraints);
+    homogeneous_constraints.make_consistent_in_parallel(locally_owned_dofs,
+                                                        locally_relevant_dofs,
+                                                        mpi_communicator);
+    homogeneous_constraints.close();
+
+
+    update_current_constraints(time);
+
+
+    // The final block of code resets and initializes the matrix object with
+    // the appropriate sparsity pattern. Recall that we do not store solution
+    // vectors in this class (the time integrator object does that internally)
+    // and so do not have to resize and initialize them either.
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    dsp,
+                                    homogeneous_constraints,
+                                    false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               locally_owned_dofs,
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    jacobian_matrix.reinit(locally_owned_dofs,
+                           locally_owned_dofs,
+                           dsp,
+                           mpi_communicator);
+  }
+
+
+  // @sect4{The HeatEquation::output_results() function}
+  //
+  // This function is called from "monitor" function that is called in turns
+  // by the time stepper in each time step. We use it to write the solution to a
+  // file, and provide graphical output through paraview or visit. We also write
+  // a pvd file, which groups all metadata about the `.vtu` files into a single
+  // file that can be used to load the full time dependent solution in paraview.
+  template <int dim>
+  void HeatEquation<dim>::output_results(const double       time,
+                                         const unsigned int timestep_number,
+                                         const PETScWrappers::MPI::Vector &y)
+  {
+    TimerOutput::Scope t(computing_timer, "output results");
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(y, "U");
+    data_out.build_patches();
+
+    data_out.set_flags(DataOutBase::VtkFlags(time, timestep_number));
+
+    const std::string filename =
+      "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
+    data_out.write_vtu_in_parallel(filename, mpi_communicator);
+
+    if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+      {
+        static std::vector<std::pair<double, std::string>> times_and_names;
+        times_and_names.emplace_back(time, filename);
+
+        std::ofstream pvd_output("solution.pvd");
+        DataOutBase::write_pvd_record(pvd_output, times_and_names);
+      }
+  }
+
+
+
+  // @sect4{The HeatEquation::implicit_function() function}
+  //
+  // As discussed in the introduction, we describe the ODE system to the time
+  // stepper via its residual,
+  // @f[
+  //   R(t,U,\dot U) = M \frac{\partial U(t)}{\partial t} + AU(t) - F(t).
+  // @f]
+  // The following function computes it, given vectors for $U,\dot U$ that
+  // we will denote by `y` and `y_dot` because that's how they are called
+  // in the documentation of the PETScWrappers::TimeStepper class.
+  //
+  // At the top of the function, we do the usual set up when computing
+  // integrals. We face two minor difficulties here: the first is that
+  // the `y` and `y_dot` vectors we get as input are read only, but we
+  // need to make sure they satisfy the correct boundary conditions and so
+  // have to set elements in these vectors. The second is that we need to
+  // compute the residual, and therefore in general we need to evaluate solution
+  // values and gradients inside locally owned cells, and for this need access
+  // to degrees of freedom which may be owned by neighboring processors. To
+  // address these issues, we create (non-ghosted) writable copies of the input
+  // vectors, apply boundary conditions and hanging node current_constraints;
+  // and then copy these vectors to ghosted vectors before we can do anything
+  // sensible with them.
+  template <int dim>
+  void
+  HeatEquation<dim>::implicit_function(const double                      time,
+                                       const PETScWrappers::MPI::Vector &y,
+                                       const PETScWrappers::MPI::Vector &y_dot,
+                                       PETScWrappers::MPI::Vector &residual)
+  {
+    TimerOutput::Scope t(computing_timer, "implicit function");
+
+    PETScWrappers::MPI::Vector tmp_solution(locally_owned_dofs,
+                                            mpi_communicator);
+    PETScWrappers::MPI::Vector tmp_solution_dot(locally_owned_dofs,
+                                                mpi_communicator);
+    tmp_solution     = y;
+    tmp_solution_dot = y_dot;
+
+    update_current_constraints(time);
+    current_constraints.distribute(tmp_solution);
+    homogeneous_constraints.distribute(tmp_solution_dot);
+
+    PETScWrappers::MPI::Vector locally_relevant_solution(locally_owned_dofs,
+                                                         locally_relevant_dofs,
+                                                         mpi_communicator);
+    PETScWrappers::MPI::Vector locally_relevant_solution_dot(
+      locally_owned_dofs, locally_relevant_dofs, mpi_communicator);
+    locally_relevant_solution     = tmp_solution;
+    locally_relevant_solution_dot = tmp_solution_dot;
+
+
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+    FEValues<dim>     fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              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();
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    std::vector<Tensor<1, dim>> solution_gradients(n_q_points);
+    std::vector<double>         solution_dot_values(n_q_points);
+
+    Vector<double> cell_residual(dofs_per_cell);
+
+    right_hand_side_function.set_time(time);
+
+    // Now for computing the actual residual. Recall that we wan to compute the
+    // vector
+    // @f[
+    //   R(t,U,\dot U) = M \frac{\partial U(t)}{\partial t} + AU(t) - F(t).
+    // @f]
+    // We could do that by actually forming the matrices $M$ and $A$, but this
+    // is not efficient. Instead, recall (by writing out how the elements of
+    // $M$ and $A$ are defined, and exchanging integrals and sums) that the
+    // $i$th element of the residual vector is given by
+    // @f{align*}{
+    //   R(t,U,\dot U)_i
+    //     &= \sum_j \int_\Omega \varphi_i(\mathbf x, t) \varphi_j(\mathbf x, t)
+    //     {\partial U_j(t)}{\partial t}
+    //       + \sum_j \int_\Omega \nabla \varphi_i(\mathbf x, t) \cdot \nabla
+    //       \varphi_j(\mathbf x, t) U_j(t)
+    //       - \int_\Omega \varphi_i f(\mathbf x, t)
+    //     \\ &=
+    //     \int_\Omega \varphi_i(\mathbf x, t) u_h(\mathbf x, t)
+    //       + \int_\Omega \nabla \varphi_i(\mathbf x, t) \cdot \nabla
+    //         u_h(\mathbf x, t)
+    //       - \int_\Omega \varphi_i f(\mathbf x, t).
+    // @f}
+    // We can compute these integrals efficiently by breaking them up into
+    // a sum over all cells and then applying quadrature. For the integrand,
+    // we need to evaluate the solution and its gradient at the quadrature
+    // points within each locally owned cell, and for this, we need also
+    // access to degrees of freedom that may be owned by neighboring
+    // processors. We therefore use the locally_relevant_solution and and
+    // locally_relevant_solution_dot vectors.
+    residual = 0;
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          fe_values.reinit(cell);
+
+          fe_values.get_function_gradients(locally_relevant_solution,
+                                           solution_gradients);
+          fe_values.get_function_values(locally_relevant_solution_dot,
+                                        solution_dot_values);
+
+          cell->get_dof_indices(local_dof_indices);
+
+          cell_residual = 0;
+          for (const unsigned int q : fe_values.quadrature_point_indices())
+            for (const unsigned int i : fe_values.dof_indices())
+              {
+                cell_residual(i) +=
+                  (fe_values.shape_value(i, q) *       // [phi_i(x_q) *
+                     solution_dot_values[q]            //  u(x_q)
+                   +                                   //  +
+                   fe_values.shape_grad(i, q) *        //  grad phi_i(x_q) *
+                     solution_gradients[q]             //  grad u(x_q)
+                   -                                   //  -
+                   fe_values.shape_value(i, q) *       //  phi_i(x_q) *
+                     right_hand_side_function.value(   //
+                       fe_values.quadrature_point(q))) //  f(x_q)]
+                  * fe_values.JxW(q);                  // * dx
+              }
+          current_constraints.distribute_local_to_global(cell_residual,
+                                                         local_dof_indices,
+                                                         residual);
+        }
+    residual.compress(VectorOperation::add);
+
+    // The end result of the operations above is a vector that contains the
+    // residual vector, having taken into account the constraints due to
+    // hanging nodes and Dirichlet boundary conditions (by virtue of having
+    // used `current_constraints.distribute_local_to_global()` to add the
+    // local contributions to the global vector. At the end of the day, the
+    // residual vector $r$ will be used in the solution of linear systems
+    // of the form $J z = r$ with the "Jacobian" matrix that we define
+    // below. We want to achieve that for algebraic components, the algebraic
+    // components of $z$ have predictable values that achieve the purposes
+    // discussed in the following. We do this by ensuring that the entries
+    // corresponding to algebraic components in the residual $r$ have specific
+    // values, and then we will do the same in the next function for the
+    // matrix; for this, you will have to know that the rows and columns
+    // of the matrix corresponding to constrained entries are zero with the
+    // exception of the diagonal entries. We will manually set that diagonal
+    // entry to one, and so $z_i=r_i$ for algebraic components.
+    //
+    // From the point of view of the residual vector, if the input `y`
+    // vector does not contain the correct values on constrained degrees of
+    // freedom (hanging nodes or boundary conditions), we need to communicate
+    // this to the time stepper, and we do so by setting the residual to the
+    // actual difference between the input `y` vector and the our local copy of
+    // it, in which we have applied the constraints (see the top of the
+    // function where we called `current_constraints.distribute(tmp_solution)`
+    // and a similar operation on the time derivative). Since we have made a
+    // copy of the input vector for this purpose, we use it to compute the
+    // residual value. However, there is a difference between hanging nodes
+    // constraints and boundary conditions: we do not want to make hanging node
+    // constraints actually depend on their dependent degrees of freedom, since
+    // this would imply that we are actually solving for the dependent degrees
+    // of freedom. This is not what we are actually doing, however, since
+    // hanging nodes are not actually solved for. They are eliminated from the
+    // system by the call to AffineConstraints::distribute_local_to_global()
+    // above. From the point of view of the Jacobian matrix, we are effectively
+    // setting hanging nodes to an artificial value (usually zero), and we
+    // simply want to make sure that we solve for those degrees of freedom a
+    // posteriori, by calling the function AffineConstraints::distribute().
+    //
+    // Here we therefore check that the residual is equal to the input value on
+    // the constrained dofs corresponding to hanging nodes (i.e., those for
+    // which the lines of the `current_constraints` contain at least one other
+    // entry), and to the difference between the input vector and the actual
+    // solution on those constraints that correspond to boundary conditions. 
+    for (const auto &c : current_constraints.get_lines())
+      if (locally_owned_dofs.is_element(c.index))
+        {
+          if (c.entries.empty()) /* no dependencies -> a Dirichlet node */
+            residual[c.index] = y[c.index] - tmp_solution[c.index];
+          else /* has dependencies -> a hanging node */
+            residual[c.index] = y[c.index];
+        }
+    residual.compress(VectorOperation::insert);
+  }
+
+
+  // @sect4{The HeatEquation::assemble_implicit_jacobian() function}
+  //
+  // The next operation is to compute the "Jacobian", which PETSc TS defines
+  // as the matrix
+  // @f[
+  //   J_\alpha = \dfrac{\partial R}{\partial y} + \alpha \dfrac{\partial
+  //   R}{\partial \dot y}
+  // @f]
+  // which, for the current linear problem, is simply
+  // @f[
+  //   J_\alpha = A + \alpha M
+  // @f]
+  // and which is in particular independent of time and the current solution
+  // vectors $y$ and $\dot y$.
+  //
+  // Having seen the assembly of matrices before, there is little that should
+  // surprise you in the actual assembly here:
+  template <int dim>
+  void HeatEquation<dim>::assemble_implicit_jacobian(
+    const double /* time */,
+    const PETScWrappers::MPI::Vector & /* y */,
+    const PETScWrappers::MPI::Vector & /* y_dot */,
+    const double alpha)
+  {
+    TimerOutput::Scope t(computing_timer, "assemble implicit Jacobian");
+
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+    FEValues<dim>     fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+
+    jacobian_matrix = 0;
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          fe_values.reinit(cell);
+
+          cell->get_dof_indices(local_dof_indices);
+
+          cell_matrix = 0;
+          for (const unsigned int q : fe_values.quadrature_point_indices())
+            for (const unsigned int i : fe_values.dof_indices())
+              for (const unsigned int j : fe_values.dof_indices())
+                {
+                  cell_matrix(i, j) +=
+                    (fe_values.shape_grad(i, q) *      // grad phi_i(x_q) *
+                       fe_values.shape_grad(j, q)      // grad phi_j(x_q)
+                     + alpha *                         //
+                         fe_values.shape_value(i, q) * // phi_i(x_q) *
+                         fe_values.shape_value(j, q)   // phi_j(x_q)
+                     ) *
+                    fe_values.JxW(q); // * dx
+                }
+          current_constraints.distribute_local_to_global(cell_matrix,
+                                                         local_dof_indices,
+                                                         jacobian_matrix);
+        }
+    jacobian_matrix.compress(VectorOperation::add);
+
+    // The only interesting part is the following. Recall that we modified the
+    // residual vector's entries corresponding to the algebraic components of
+    // the solution in the previous function. The outcome of calling
+    // `current_constraints.distribute_local_to_global()` a few lines
+    // above is that the global matrix has zero rows and columns for the
+    // algebraic (constrained) components of the solution; the function
+    // puts a value on the diagonal that is nonzero and has about the
+    // same size as the remaining diagonal entries of the matrix. What
+    // this diagonal value is is unknown to us -- in other cases where
+    // we call `current_constraints.distribute_local_to_global()` on
+    // both the left side matrix and the right side vector, as in most
+    // other tutorial programs, the matrix diagonal entries and right
+    // hand side values are chosen in such a way that the result of
+    // solving a linear system is what we want it to be, but the scaling
+    // is done automatically.
+    //
+    // This is not good enough for us here, because we are building the
+    // right hand side independently from the matrix in different functions.
+    // Thus, for any constrained degree of freedom, we set the diagonal of the
+    // Jacobian to be one. This leaves the Jacobian matrix invertible,
+    // consistent with what the time stepper expects, and it also makes sure
+    // that if we did not make a mistake in the residual and/or in the Jacbian
+    // matrix, then asking the time stepper to check the Jacobian with a finite
+    // difference method will produce the correct result. This can be activated
+    // at run time via passing the `-snes_test_jacobian` option on the command
+    // line.
+    for (const auto &c : current_constraints.get_lines())
+      jacobian_matrix.set(c.index, c.index, 1.0);
+    jacobian_matrix.compress(VectorOperation::insert);
+  }
+
+
+  // @sect4{The HeatEquation::solve_with_jacobian() function}
+  //
+  // This is the function that actually solves the linear system with the
+  // Jacobian matrix we have previously built (in a call to the previous
+  // function during the current time step or another earlier one -- time
+  // steppers are quite sophisticated in determining internally whether it is
+  // necessary to update the Jacobian matrix, or whether one can reuse it for
+  // another time step without rebuilding it; this is similar to how one can
+  // re-use the Newton matrix for several Newton steps, see for example the
+  // discussion in step-77). We could in principle not provide this function to
+  // the time stepper, and instead select a specific solver on the command line
+  // by using the `-ksp_*` options of PETSc. However, by providing this
+  // function, we can use a specific solver and preconditioner for the linear
+  // system, and still have the possibility to change them on the command line.
+  //
+  // Providing a specific solver is more in line with the way we usually do
+  // things in other deal.II examples, while letting PETSc choose a generic
+  // solver, and changing it on the command line via `-ksp_type` is more in line
+  // with the way PETSc is usually used, and it can be a convenient approach
+  // when we are experimenting to find an optimal solver for our problem. Both
+  // options are available here, since we can still change both the solver and
+  // the preconditioner on the command line via `-user_ksp_type` and
+  // `-user_pc_type` options.
+  //
+  // In any case, recall that the Jacobian we built in the previous function is
+  // always of the form
+  // @f[
+  //   J_\alpha = \alpha M + A
+  // @f]
+  // where $M$ is a mass matrix and $A$ a Laplace matrix. $M$ is symmetric and
+  // positive definite; $A$ is symmetric and at least positive semidefinite;
+  // $\alpha> 0$. As a consequence, the Jacobian matrix is a symmetric and
+  // positive definite matrix, which we can efficiently solve with the Conjugate
+  // Gradient method, along with either SSOR or (if available) the algebraic
+  // multigrid implementation provided by PETSc (via the Hypre package) as
+  // preconditioner. In practice, if you wanted to solve "real" problems, one
+  // would spend some time finding which preconditioner is optimal, perhaps
+  // using PETSc's ability to read solver and preconditioner choices from the
+  // command line. But this is not the focus of this tutorial program, and so
+  // we just go with the following:
+  template <int dim>
+  void
+  HeatEquation<dim>::solve_with_jacobian(const PETScWrappers::MPI::Vector &src,
+                                         PETScWrappers::MPI::Vector       &dst)
+  {
+    TimerOutput::Scope t(computing_timer, "solve with Jacobian");
+
+#if defined(PETSC_HAVE_HYPRE)
+    PETScWrappers::PreconditionBoomerAMG preconditioner;
+    preconditioner.initialize(jacobian_matrix);
+#else
+    PETScWrappers::PreconditionSSOR preconditioner;
+    preconditioner.initialize(
+      jacobian_matrix, PETScWrappers::PreconditionSSOR::AdditionalData(1.0));
+#endif
+
+    SolverControl           solver_control(1000, 1e-8 * src.l2_norm());
+    PETScWrappers::SolverCG cg(solver_control);
+    cg.set_prefix("user_");
+
+    cg.solve(jacobian_matrix, dst, src, preconditioner);
+
+    pcout << "     " << solver_control.last_step() << " linear iterations."
+          << std::endl;
+  }
+
+
+  // @sect4{The HeatEquation::prepare_for_coarsening_and_refinement() function}
+  //
+  // The next block of functions deals with mesh refinement. We split this
+  // process up into a "decide whether and what you want to refine" and a
+  // "please transfer these vectors from old to new mesh" phase, where the first
+  // also deals with marking cells for refinement. (The decision whether or
+  // not to refine is done in the lambda function that calls the current
+  // function.)
+  //
+  // Breaking things into a "mark cells" function and into a "execute mesh
+  // adaptation and transfer solution vectors" function is awkward, though
+  // conceptually not difficult to understand. These two pieces of code
+  // should really be part of the same function, as they are in step-26.
+  // The issue is with what PETScWrappers::TimeStepper provides us with in these
+  // callbacks. Specifically, the "decide whether and what you want to refine"
+  // callback has access to the current solution, and so can evaluate
+  // (spatial) error estimators to decide which cells to refine. The
+  // second callback that transfers vectors from old to new mesh gets a bunch
+  // of vectors, but without the semantic information on which of these is
+  // the current solution vector. As a consequence, it cannot do the marking
+  // of cells for refinement or coarsening, and we have to do that from the
+  // first callback.
+  //
+  // In practice, however, the problem is minor. The first of these
+  // two functions is essentially identical to the
+  // first half of the corresponding function in step-26, with the only
+  // difference that it uses the parallel::distributed::GridRefinement namespace
+  // function instead of the serial one. Once again, we make sure that we never
+  // fall below the minimum refinement level, and above the maximum one, that we
+  // can select from the parameter file.
+  template <int dim>
+  void HeatEquation<dim>::prepare_for_coarsening_and_refinement(
+    const PETScWrappers::MPI::Vector &y)
+  {
+    PETScWrappers::MPI::Vector locally_relevant_solution(locally_owned_dofs,
+                                                         locally_relevant_dofs,
+                                                         mpi_communicator);
+    locally_relevant_solution = y;
+
+    Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+    KellyErrorEstimator<dim>::estimate(dof_handler,
+                                       QGauss<dim - 1>(fe.degree + 1),
+                                       {},
+                                       locally_relevant_solution,
+                                       estimated_error_per_cell);
+
+    parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
+      triangulation, estimated_error_per_cell, 0.6, 0.4);
+
+    const unsigned int max_grid_level =
+      initial_global_refinement + max_delta_refinement_level;
+    const unsigned int min_grid_level = initial_global_refinement;
+
+    if (triangulation.n_levels() > max_grid_level)
+      for (const auto &cell :
+           triangulation.active_cell_iterators_on_level(max_grid_level))
+        cell->clear_refine_flag();
+    for (const auto &cell :
+         triangulation.active_cell_iterators_on_level(min_grid_level))
+      cell->clear_coarsen_flag();
+  }
+
+
+  // @sect4{The HeatEquation::interpolate() function}
+  //
+  // The following function then is the second half of the correspond function
+  // in step-26. It is called by the time stepper whenever it requires to
+  // transfer the solution and any intermediate stage vectors to a new mesh. We
+  // must make sure that all input vectors are transformed into ghosted vectors
+  // before the actual transfer is executed, and that we distribute the hanging
+  // node constraints on the output vectors as soon as we have interpolated the
+  // vectors to the new mesh -- i.e., that all constraints are satisfied on
+  // the vectors we transfer.
+  //
+  // We have no way to enforce boundary conditions at this stage, however. This
+  // is because the various vectors may correspond to solutions at previous time
+  // steps if the method used here is a multistep time integrator, and so may
+  // correspond to different time points that we are not privy to.
+  //
+  // While this could be a problem if we used the values of the solution and of
+  // the intermediate stages on the constrained degrees of freedom to compute
+  // the errors, we do not do so. Instead, we compute the errors on the
+  // differential equation, and ignore algebraic constraints; therefore we do no
+  // need to guarantee that the boundary conditions are satisfied also in the
+  // intermediate stages.
+  //
+  // We have at our disposal the hanging node current_constraints alone, though,
+  // and therefore we enforce them on the output vectors, even if this is not
+  // really needed.
+  template <int dim>
+  void HeatEquation<dim>::interpolate(
+    const std::vector<PETScWrappers::MPI::Vector> &all_in,
+    std::vector<PETScWrappers::MPI::Vector>       &all_out)
+  {
+    parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector>
+      solution_trans(dof_handler);
+
+    std::vector<PETScWrappers::MPI::Vector> all_in_ghosted(all_in.size());
+    std::vector<const PETScWrappers::MPI::Vector *> all_in_ghosted_ptr(
+      all_in.size());
+    std::vector<PETScWrappers::MPI::Vector *> all_out_ptr(all_in.size());
+    for (unsigned int i = 0; i < all_in.size(); ++i)
+      {
+        all_in_ghosted[i].reinit(locally_owned_dofs,
+                                 locally_relevant_dofs,
+                                 mpi_communicator);
+        all_in_ghosted[i]     = all_in[i];
+        all_in_ghosted_ptr[i] = &all_in_ghosted[i];
+      }
+
+    triangulation.prepare_coarsening_and_refinement();
+    solution_trans.prepare_for_coarsening_and_refinement(all_in_ghosted_ptr);
+    triangulation.execute_coarsening_and_refinement();
+
+    setup_system(boundary_values_function.get_time());
+
+    all_out.resize(all_in.size());
+    for (unsigned int i = 0; i < all_in.size(); ++i)
+      {
+        all_out[i].reinit(locally_owned_dofs, mpi_communicator);
+        all_out_ptr[i] = &all_out[i];
+      }
+    solution_trans.interpolate(all_out_ptr);
+
+    for (PETScWrappers::MPI::Vector &v : all_out)
+      hanging_nodes_constraints.distribute(v);
+  }
+
+
+  // @sect4{The HeatEquation::update_current_constraints() function}
+  //
+  // Since regenerating the constraints at each time step may be expensive, we
+  // make sure that we only do so when the time changes. We track time change by
+  // checking if the time of the boundary_values_function has changed, with
+  // respect to the time of the last call to this function. This will work most
+  // of the times, but not the very first time we call this function, since the
+  // time then may be zero and the time of the `boundary_values_function` is
+  // zero at construction time. We therefore also check if the number
+  // constraints in `current_constraints`, and if these are empty, we regenerate
+  // the constraints regardless of the time variable.
+  template <int dim>
+  void HeatEquation<dim>::update_current_constraints(const double time)
+  {
+    if (current_constraints.n_constraints() == 0 ||
+        time != boundary_values_function.get_time())
+      {
+        TimerOutput::Scope t(computing_timer, "update current constraints");
+
+        boundary_values_function.set_time(time);
+        current_constraints.clear();
+        current_constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
+        current_constraints.merge(hanging_nodes_constraints);
+        VectorTools::interpolate_boundary_values(dof_handler,
+                                                 0,
+                                                 boundary_values_function,
+                                                 current_constraints);
+        current_constraints.make_consistent_in_parallel(locally_owned_dofs,
+                                                        locally_relevant_dofs,
+                                                        mpi_communicator);
+        current_constraints.close();
+      }
+  }
+
+
+  // @sect4{The HeatEquation::run() function}
+  //
+  // We have finally arrived at the main function of the class. At the top, it
+  // creates the mesh and sets up the variables that make up the linear system:
+  template <int dim>
+  void HeatEquation<dim>::run()
+  {
+    GridGenerator::hyper_L(triangulation);
+    triangulation.refine_global(initial_global_refinement);
+
+    setup_system(/* time */ 0);
+
+    // We then set up the time stepping object and associate the matrix we will
+    // build whenever requested for both the Jacobian matrix (see the definition
+    // above of what the "Jacobian" actually refers to) and for the matrix
+    // that will be used as a preconditioner for the Jacobian.
+    PETScWrappers::TimeStepper<PETScWrappers::MPI::Vector,
+                               PETScWrappers::MPI::SparseMatrix>
+      petsc_ts(time_stepper_data);
+
+    petsc_ts.set_matrices(jacobian_matrix, jacobian_matrix);
+
+
+    // The real work setting up the time stepping object starts here. As
+    // discussed in the introduction, the way the PETScWrappers::TimeStepper
+    // class is used is by inverting control: At the end of this function, we
+    // will call PETScWrappers::TimeStepper::solve() which internally will
+    // run the loop over time steps, and at the appropriate places *call back*
+    // into user code for whatever functionality is required. What we need to
+    // do is to hook up these callback functions by assigning appropriate
+    // lambda functions to member variables of the `petsc_ts` object.
+    //
+    // We start by creating lambda functions that provide information about
+    // the "implicit function" (i.e., that part of the right hand side of the
+    // ODE system that we want to treat implicitly -- which in our case is
+    // the entire right hand side), a function that assembles the Jacobian
+    // matrix, and a function that solves a linear system with the Jacobian.
+    petsc_ts.implicit_function = [&](const double                      time,
+                                     const PETScWrappers::MPI::Vector &y,
+                                     const PETScWrappers::MPI::Vector &y_dot,
+                                     PETScWrappers::MPI::Vector       &res) {
+      this->implicit_function(time, y, y_dot, res);
+    };
+
+    petsc_ts.setup_jacobian = [&](const double                      time,
+                                  const PETScWrappers::MPI::Vector &y,
+                                  const PETScWrappers::MPI::Vector &y_dot,
+                                  const double                      alpha) {
+      this->assemble_implicit_jacobian(time, y, y_dot, alpha);
+    };
+
+    petsc_ts.solve_with_jacobian = [&](const PETScWrappers::MPI::Vector &src,
+                                       PETScWrappers::MPI::Vector       &dst) {
+      this->solve_with_jacobian(src, dst);
+    };
+
+    // The next two callbacks deal with identifying and setting variables
+    // that are considered "algebraic" (rather than "differential"), i.e., for
+    // which we know what values they are supposed to have rather than letting
+    // their values be determined by the differential equation. We need to
+    // instruct the time stepper to ignore these components when computing the
+    // error in the residuals, and we do so by first providing a function that
+    // returns an IndexSet with the indices of these algebraic components of the
+    // solution vector (or rather, that subset of the locally-owned part of the
+    // vector that is algebraic, in case we are running in parallel). This first
+    // of the following two functions does that. Specifically, both nodes at
+    // which Dirichlet boundary conditions are applied, and hanging nodes are
+    // algebraically constrained. This function then returns a set of
+    // indices that is initially empty (but knows about the size of the index
+    // space) and which we then construct as the union of boundary and hanging
+    // node indices.
+    //
+    // Following this, we then also need a function that, given a solution
+    // vector `y` and the current time, sets the algebraic components of that
+    // vector to their correct value. This comes down to ensuring that we have
+    // up to date constraints in the `constraints` variable, and then applying
+    // these constraints to the solution vector via
+    // AffineConstraints::distribute(). (It is perhaps worth noting that we
+    // *could* have achieved the same in `solve_with_jacobian()`. Whenever the
+    // time stepper solves a linear system, it follows up the call to the solver
+    // by calling the callback to set algebraic components correct. We could
+    // also have put the calls to `update_current_constraints()` and
+    // `distribute()` into the `solve_with_jacobian` function, but by not doing
+    // so, we can also replace the `solve_with_jacobian` function with a call to
+    // a PETSc solver, and we would still have the current_constraints correctly
+    // applied to the solution vector.)
+    petsc_ts.algebraic_components = [&]() {
+      IndexSet algebraic_set(dof_handler.n_dofs());
+      algebraic_set.add_indices(DoFTools::extract_boundary_dofs(dof_handler));
+      algebraic_set.add_indices(
+        DoFTools::extract_hanging_node_dofs(dof_handler));
+      return algebraic_set;
+    };
+
+    petsc_ts.update_constrained_components =
+      [&](const double time, PETScWrappers::MPI::Vector &y) {
+        TimerOutput::Scope t(computing_timer, "set algebraic components");
+        update_current_constraints(time);
+        current_constraints.distribute(y);
+      };
+
+
+    // The next two callbacks relate to mesh refinement. As discussed in the
+    // introduction, PETScWrappers::TimeStepper knows how to deal with the
+    // situation where we want to change the mesh. All we have to provide
+    // is a callback that returns `true` if we are at a point where we want
+    // to refine the mesh (and `false` otherwise) and that if we want to
+    // do mesh refinement does some prep work for that in the form of
+    // calling the `prepare_for_coarsening_and_refinement` function.
+    //
+    // If the first callback below returns `true`, then PETSc TS will
+    // do some clean-up operations, and call the second of the callback
+    // functions (`petsc_ts.interpolate`) with a collection of vectors that
+    // need to be interpolated from the old to the new mesh. This may include
+    // the current solution, perhaps the current time derivative of the
+    // solution, and in the case of
+    // [multistep time
+    // integrators](https://en.wikipedia.org/wiki/Linear_multistep_method) also
+    // the solutions of some previous time steps. We hand all of these over to
+    // the `interpolate()` member function of this class.
+    petsc_ts.decide_and_prepare_for_remeshing =
+      [&](const double /* time */,
+          const unsigned int                step_number,
+          const PETScWrappers::MPI::Vector &y) -> bool {
+      if (step_number > 0 && this->mesh_adaptation_frequency > 0 &&
+          step_number % this->mesh_adaptation_frequency == 0)
+        {
+          pcout << std::endl << "Adapting the mesh..." << std::endl;
+          this->prepare_for_coarsening_and_refinement(y);
+          return true;
+        }
+      else
+        return false;
+    };
+
+    petsc_ts.interpolate =
+      [&](const std::vector<PETScWrappers::MPI::Vector> &all_in,
+          std::vector<PETScWrappers::MPI::Vector>       &all_out) {
+        this->interpolate(all_in, all_out);
+      };
+
+    // The final callback is a "monitor" that is called in each
+    // time step. Here we use it to create a graphical output. Perhaps a better
+    // scheme would output the solution at fixed time intervals, rather
+    // than in every time step, but this is not the main point of
+    // this program and so we go with the easy approach:
+    petsc_ts.monitor = [&](const double                      time,
+                           const PETScWrappers::MPI::Vector &y,
+                           const unsigned int                step_number) {
+      pcout << "Time step " << step_number << " at t=" << time << std::endl;
+      this->output_results(time, step_number, y);
+    };
+
+
+    // With all of this out of the way, the rest of the function is
+    // anticlimactic: We just have to initialize the solution vector with
+    // the initial conditions and call the function that does the time
+    // stepping, and everything else will happen automatically:
+    PETScWrappers::MPI::Vector solution(locally_owned_dofs, mpi_communicator);
+    VectorTools::interpolate(dof_handler, initial_value_function, solution);
+
+    petsc_ts.solve(solution);
+  }
+} // namespace Step86
+
+
+// @sect3{The main() function}
+//
+// The rest of the program is as it always looks. We read run-time parameters
+// from an input file via the ParameterAcceptor class in the same way as
+// we showed in step-60 and step-70.
+int main(int argc, char **argv)
+{
+  try
+    {
+      using namespace Step86;
+
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      HeatEquation<2>                  heat_equation_solver(MPI_COMM_WORLD);
+      ParameterAcceptor::initialize("heat_equation.prm",
+                                    "heat_equation_used.prm");
+      heat_equation_solver.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}

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.