]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-83: A program demonstrating checkpoint/restart.
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 19 Jul 2024 00:03:02 +0000 (18:03 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Sun, 21 Jul 2024 18:03:16 +0000 (12:03 -0600)
examples/step-83/CMakeLists.txt [new file with mode: 0644]
examples/step-83/doc/builds-on [new file with mode: 0644]
examples/step-83/doc/intro.dox [new file with mode: 0644]
examples/step-83/doc/kind [new file with mode: 0644]
examples/step-83/doc/results.dox [new file with mode: 0644]
examples/step-83/doc/tooltip [new file with mode: 0644]
examples/step-83/step-83.cc [new file with mode: 0644]

diff --git a/examples/step-83/CMakeLists.txt b/examples/step-83/CMakeLists.txt
new file mode 100644 (file)
index 0000000..f78956b
--- /dev/null
@@ -0,0 +1,54 @@
+##
+#  CMake script for the step-83 tutorial program:
+##
+
+# Set the name of the project and target:
+set(TARGET "step-83")
+
+# 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_MPI) # keep in one line
+  message(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following option:
+    DEAL_II_WITH_MPI = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options:
+    DEAL_II_WITH_MPI = ${DEAL_II_WITH_MPI}
+This conflicts with the requirements."
+    )
+endif()
+
+
+deal_ii_initialize_cached_variables()
+set(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.ucd *.d2 *.vtu *.pvtu)
+project(${TARGET})
+deal_ii_invoke_autopilot()
diff --git a/examples/step-83/doc/builds-on b/examples/step-83/doc/builds-on
new file mode 100644 (file)
index 0000000..4060861
--- /dev/null
@@ -0,0 +1 @@
+step-19
diff --git a/examples/step-83/doc/intro.dox b/examples/step-83/doc/intro.dox
new file mode 100644 (file)
index 0000000..6b0a356
--- /dev/null
@@ -0,0 +1,645 @@
+<br>
+
+<i>
+This program was contributed by
+Pasquale Africa (SISSA),
+Wolfgang Bangerth (Colorado State University), and
+Bruno Blais (Polytechnique Montreal).
+
+This material is based upon work partially supported by National Science
+Foundation grants OAC-1835673, DMS-1821210, and EAR-1925595;
+and by the Computational Infrastructure in
+Geodynamics initiative (CIG), through the National Science Foundation under
+Award No. EAR-1550901 and The University of California-Davis.
+</i>
+<br>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+
+<h3>Checkpoint/restart</h3>
+
+A common problem in scientific computing is that a program runs such a
+long time that one would like to be able to periodically save the
+state of the program, and then have the possibility to restart at the
+same point where the program was when the state was saved. There are a
+number of reasons why one would want to do that:
+
+- On machines that are shared with many other users, say on
+  supercomputers, it is common to use queues that have a time limit on
+  how long programs can run. For example, one can only run programs
+  that run at most 48 hours. But what if a simulation hasn't done the
+  requisite number of time steps when those 48 hours are up, or for
+  some other reason isn't finished yet? If we could save the state of
+  the program every hour, that wouldn't matter so much: If the queuing
+  system terminates the program after 48 hours, we could simply let
+  the program run again and pick up where the last state was saved
+  somewhere between hour 47 and 48 -- at most an hour of run time
+  would have been lost.
+
+- For some of the very largest computations one can do today, with
+  tens or hundreds of thousands of processors (this being written in
+  2023), it is not entirely uncommon that one has to expect that a
+  node in a supercomputer dies over the course of long run. This could
+  be because of hardware failures, system software failures, power
+  failures, or any other reason that might affect the stability of the
+  system. It is not entirely unreasonable that one would have to
+  expect that: A computation with, say, 100,000 cores on machines with
+  32 cores each uses 3,000 nodes; a computation that takes 24 hours on
+  such a machine uses nearly 10 years of machine time, a time frame
+  within which we would certainly expect that a typical machine might
+  die. In such cases, losing a simulation bears a very substsntial
+  cost, given how expensive it is to get 100,000 cores for a whole
+  day, and being able to periodically save the state will guard
+  against the loss of this investment.
+
+- There are sometimes cases where one would like to compute an
+  expensive step once and then do many inexpensive computations on top
+  of it. For example, one might be interested in an expensive
+  simulation of the weather pattern today, and then consider where
+  rain falls by looking at many stochastic perturbations. In such a
+  case, one would save the state at the end of the expensive
+  computation, and re-start from it many times over with different
+  perturbations.
+
+The ability to save the state of a program and restart from it is
+often called *checkpoint/restart*. "Checkpointing" refers to the step
+of writing the current state of the program into one or several files
+on disk (or, in more general schemes, into some kind of temporary
+storage: disk files may be one option, but it could also be
+non-volatile memory on other nodes of a parallel
+machine). "Restarting" means starting a program not from scratch, but
+from a previously stored "checkpoint". This tutorial discusses both how
+one should *conceptually* think about checkpoint/restart, as well as
+how this is typically implemented.
+
+The program is informed by how this has been implemented in the
+[ASPECT](https://aspect.geodynamics.org) code, in which we have used
+this strategy for many years successfully. Many of the ideas that can
+be found in ASPECT are based on approaches that other codes have used
+long before that.
+
+
+<h3>Serialization/deserialization</h3>
+
+The second term that we should introduce is "serialization" and
+"deserialization". A key piece of checkpoint/restart is that we need
+to "dump" the *state of the program* (i.e., the data structures
+currently stored in memory, and where the program currently is in the
+sequence of what operations it is doing one after the other) into a
+file or some other temporary storage.
+
+The thing, however, is that we are often using very complicated data
+structures. Think, for example, of a `std::map<unsigned
+int,std::string>` object: This isn't just an array of `double`
+numbers, but is in general stored as a tree data structure where every
+node will have to consist of the key (an integer) and the value (a
+string, which is itself a non-trivial object). How would one store
+this on disk?
+
+The solution is that one has to convert *every* piece of data that
+corresponds to a *linear array of bytes* that can then be saved in a
+file -- because what files store are ultimately just a sequence of
+bytes. This process is called "serialization": It converts the data
+that represents the program into a series of bytes. The way this
+conceptually works is that one has to define bottom-up serialization
+functions: If we have a function that converts an `unsigned int` into
+a series of bytes (for example by just storing the four bytes that
+represent an integer on most current architectures) and a function
+that converts a `std::string` into a series of bytes (for example, by
+storing first the length of the string as an `unsigned int` in the way
+described before, and then the individual characters that make up that
+string), then we can define a function that stores a `std::map`: For
+example, we could first store the number of entries in the map, and
+then for each key/value pair, we first serialize the key and then the
+value. We can do this for the larger-and-large classes and eventually
+the whole program: We convert each member variable in turn into a
+sequence of bytes if we have already defined how each of the member
+variables individually would do that.
+
+From this kind of information, we can then also re-load the state of
+the program: Starting with the top-level object of our program, we
+read each member variable in turn from the file, where each member
+variable may again be a class of its own that read its member
+variables, etc. This process is then called "deserialization":
+Building data structures from a serial representation.
+
+It seems like a monumental task to write functions that serialize and
+again deserialize all possible data structures. One would have to do
+this for built-in types like `double`, `int`, etc., but then also for
+all of the C++ containers such as `std::vector`, `std::list`,
+`std::map` along with C-style arrays. From there, one might work one's
+way up to classes such as Tensor, Point, and eventually Triangulation,
+Vector, or ParticleHandler. This would be a lot of work.
+
+Fortunately, help exists: There are a number of external libraries
+that make this task relatively straightforward. The one that deal.II
+uses for serialization and deserialization is the
+[BOOST serialization
+library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html)
+that has everything related to (de)serialization of built-in and
+`std::` data types already available. One then only has to write
+functions for each class that needs to be (de)serialized, but this is
+also made relatively simple by BOOST serialization: One only has to
+add a single function in which one can use operator overloading. For
+example, here is a very short excerpt of how the Tensor class might
+do this:
+@code
+template <int rank, int dim, typename Number>
+class Tensor
+{
+public:
+  //
+  // Read or write the data of this object to or from a stream for the purpose
+  // of serialization using the BOOST serialization
+  // library.
+  //
+  template <class Archive>
+  void
+  serialize(Archive &ar, const unsigned int version);
+
+private:
+  //
+  // Array of tensors holding the subelements.
+  //
+  Tensor<rank-1, dim, Number> values[dim];
+};
+
+
+template <int rank, int dim, typename Number>
+template <class Archive>
+inline void
+Tensor<rank, dim, Number>::serialize(Archive &ar, const unsigned int /*version*/)
+{
+  ar & values;
+}
+@endcode
+
+In other words, the Tensor class stores an array of `dim` objects of a
+tensor of one rank lower, and in order to serialize this array, the
+function simply uses the overloaded `operator&` which recognizes that
+`values` is a C-style array of fixed length, and what the data type of
+the elements of this array are. It will then in turn call the
+`serialize()` member function of `Tensor<rank-1,dim>`. (The recursion
+has a stopping point which is not of importance for this discussion. The
+actual implementation of class Tensor looks different, but you get the
+general idea here.)
+
+Depending on whether the `Archive` type used for the first argument of
+this function is an input or output archive, `operator&` reads from or
+writes into the `values` member variable. If the class had more member
+variables, one would list them one after the other. The same is true
+if a class has base classes. For example, here is the corresponding
+function for the Quadrature class that stores quadrature points and
+weights, and is derived from the Subscriptor base class:
+@code
+template <int dim>
+template <class Archive>
+inline void
+Quadrature<dim>::serialize(Archive &ar, const unsigned int)
+{
+  // Forward to the (de)serialization function in the base class:
+  ar & static_cast<Subscriptor &>(*this);
+
+  // Then (de)serialize the member variables:
+  ar & quadrature_points & weights;
+}
+@endcode
+The beauty of having to write only one function for both serialization
+and deserialization is that obviously one has to read data from the
+archive in the same order as one writes into it, and that is easy
+to get wrong if one had to write two separate functions -- but that
+is automatically correct if one has the same function for both
+purposes!
+
+The BOOST serialization library also has mechanisms if one wants to
+write separate save and load functions. This is useful if a class
+stores a lot of data in one array and then has a cache of commonly
+accessed data on the side for fast access; one would then only store
+and load the large array, and in the load function re-build the cache
+from the large array. (See also the discussion below on what one
+actually wants to save.)
+
+Based on these principles, let us consider how one would  (de)serialize
+a program such as step-19. Recall that its principal class looked
+essentially as follows:
+@code
+  template <int dim>
+  class CathodeRaySimulator
+  {
+  public:
+    CathodeRaySimulator();
+
+    void run();
+
+  private:
+    void make_grid();
+    void setup_system();
+    void assemble_system();
+    void solve_field();
+    void refine_grid();
+
+    void create_particles();
+    void move_particles();
+    void track_lost_particle(const typename Particles::ParticleIterator<dim> &        particle,
+                             const typename Triangulation<dim>::active_cell_iterator &cell);
+
+
+    void update_timestep_size();
+    void output_results() const;
+
+    Triangulation<dim>        triangulation;
+    MappingQGeneric<dim>      mapping;
+    FE_Q<dim>                 fe;
+    DoFHandler<dim>           dof_handler;
+    AffineConstraints<double> constraints;
+
+    SparseMatrix<double>      system_matrix;
+    SparsityPattern           sparsity_pattern;
+
+    Vector<double>            solution;
+    Vector<double>            system_rhs;
+
+    Particles::ParticleHandler<dim> particle_handler;
+    types::particle_index           next_unused_particle_id;
+    types::particle_index           n_recently_lost_particles;
+    types::particle_index           n_total_lost_particles;
+    types::particle_index           n_particles_lost_through_anode;
+
+    DiscreteTime time;
+  };
+@endcode
+One would then first write a
+function that allows (de)serialization of the data of this class
+as follows:
+@code
+  template <int dim>
+  template <class Archive>
+  void CathodeRaySimulator<dim>::serialize(Archive &ar,
+                                           const unsigned int /* version */)
+  {
+    ar & triangulation;
+    ar & solution;
+    ar & particle_handler;
+    ar & next_unused_particle_id;
+    ar & n_recently_lost_particles;
+    ar & n_total_lost_particles;
+    ar & n_particles_lost_through_anode;
+    ar & time;
+  }
+@endcode
+As discussed below, this is not the entire list of member variables: We
+only want to serialize those that we cannot otherwise re-generate.
+
+Next, we need a function that writes a checkpoint by creating an "output
+archive" and using the usual `operator<<` to put an object into this archive. In
+the code of this program, this will then look as follows (with some
+minor modifications we will discuss later):
+@code
+  template <int dim>
+  void CathodeRaySimulator<dim>::checkpoint()
+  {
+    std::ofstream                 checkpoint_file("checkpoint");
+    boost::archive::text_oarchive archive(checkpoint_file,
+                                          boost::archive::no_header);
+
+    archive << *this;
+  }
+@endcode
+
+What `operator<<` does here is to call the serialization functions of the right hand
+operand (here, the `serialize()` function described above, with an output
+archive as template argument), and create a serialized representation of the data. In the
+end, serialization has put everything into the "archive" from which one
+can extract a (sometimes very long) string that one can save in bulk
+into a file, and that is exactly what happens when the destructor of the
+`archive` variable is called.
+
+BOOST serialization offers different archives, including ones that
+store the data in text format (as we do above), in binary format, or
+already compressed with something like gzip to minimize the amount of
+space necessary. The specific type of archive to be used is selected
+by the type of the `archive` variable above (and the corresponding
+variable in the `restart()` function of course). This program uses a
+text archive so that you can look at how a serialized representation
+would actually look at, though a "real" program would of course try to
+be more space efficient by using binary (and possibly compressed)
+representations of the data.
+
+In any case, the data we have thus created is read in very similarly
+in the following function (again with some minor modifications):
+@code
+  template <int dim>
+  void CathodeRaySimulator<dim>::restart()
+  {
+    std::ifstream checkpoint_file("checkpoint");
+
+    boost::archive::text_iarchive archive(checkpoint_file,
+                                          boost::archive::no_header);
+
+    archive >> *this;
+  }
+@endcode
+
+The magic of this approach is that one doesn't actually have to write very
+much code to checkpoint or restart programs: It is all hidden in the
+serialization functions of classes such as Triangulation,
+Particles::ParticleHandler, etc., which the deal.II library provides for you.
+
+
+<h4> Serialization of parallel programs </h4>
+
+The program as provided here is sequential, i.e., it runs on a single
+processor just like the original step-19 did. But what if you had a parallel
+program -- say, something like step-40 -- that runs in parallel with MPI?
+In that case, serialization and checkpoint/restart becomes more complicated.
+While parallel execution is not something that is of concern to us
+in this program, it is an issue that has influenced the design of how
+deal.II does serialization; as a consequence, we need to talk through
+what makes serialization of parallel programs more difficult in order
+to understand why this program does things the way it does.
+
+Intuitively, one might simply want to use the same idea as we used
+here except that we let every MPI process serialize its own data, and
+read its own data. This works, but there are some drawbacks:
+- There is a certain subset of data that is replicated across all
+  MPI processes and that would be then written by all processes. An example
+  is the `time` data structure that stores the current time, time step
+  size, and other information, and that better be the same on all processes.
+  Typically, the replicated data isn't a large fraction of a program's
+  overall memory usage, and perhaps writing it more than once isn't
+  going to be too much of a problem, but it is unsatisfactory anyway
+  to have this kind of data on disk more than once.
+- One would have to think in detail how exactly one wants to represent
+  the data on disk. One possibility would be for every MPI process to write
+  its own file. On the other hand, checkpointing is most useful for large
+  programs, using large numbers of processes -- it is not uncommon to use
+  checkpointing on programs that run on 10,000 or more processors in parallel.
+  This would then lead to 10,000 or more files on disk. That's unpleasant, and
+  probably inefficient as well. We could address this by first letting all
+  the processes serialize into a string in memory (using `std::ostringstream`)
+  and then collating all of these strings into one file. The MPI I/O sub-system
+  has facilities to make that happen, for example, but it will require a bit
+  of thought not the least because the serialized data from each process
+  will likely result in strings of different sizes.
+- Perhaps the most important reason to rethink how one does things in parallel
+  is because, with a little bit of thought, it is possible to checkpoint a
+  program running with $N$ MPI processes and restart it with $M\neq N$
+  processes. This may, at first, seem like a pointless exercise, but it is
+  useful if one had, for example, a program that repeatedly refines the mesh
+  and where it is inefficient to run the early refinement steps with a
+  coarse mesh on too many processes, whereas it is too slow to run the later
+  refinement steps with a fine mesh on too few processes.
+
+In order to address these issues, in particular the last one, the right approach
+is to deviate a bit from the simple scheme of having a `serialize()` function
+that simply serializes/deserializes *everything* into an archive, and then have two
+functions `checkpoint()` and `restart()` that for all practical purposes defer
+all the work to the `serialize()` function. Instead, one splits all data into
+two categories:
+- Data that is tied to the cells of a triangulation. This includes the mesh itself,
+  but also the particles in the Particles::ParticleHandler class, and most
+  importantly the solution vector(s). The way to serialize such data is to
+  attach the data to cells and then let the Triangulation class serialize the attached
+  data along with its own data. If this is done in a way so that we can re-load
+  a triangulation on a different number of processes than the data was written,
+  then this automatically also ensures that we can restore solution vectors
+  and Particles::ParticleHandler objects (and everything else we can attach
+  to the cells of a triangulation) on a different number of processes.
+- Other data. In finite element programs, this data is almost always replicated
+  across processes, and so it is enough if the "root" process (typically the
+  process with MPI rank zero) writes it to disk. Upon restart, the root
+  process reads the data from disk, sends it to all other processes (however many
+  of them there may be), and these then initialize their own copies of the
+  replicated data structures.
+
+These sorts of considerations have influenced the design of the Triangulation and
+Particles::ParticleHandler classes. In particular, Particles::ParticleHandler's
+`serialize()` function only serializes the "other data" category, but not the
+actual particles; these can instead be attached to the triangulation by calling
+Particles::ParticleHandler::prepare_for_serialization(), and then one can
+call Triangulation::save() to actually write this information into a set of
+files that become part of the checkpoint. Upon restart, we then first call
+Triangulation::load(), followed by Particles::ParticleHandler::deserialize()
+to retrieve the particles from the cells they are attached to.
+
+(We could, with relatively minimal effort, use the same scheme for the solution
+vector: The SolutionTransfer class can be used to attach the values of degrees
+of freedom to cells, and then Triangulation::save() also writes these into
+checkpoint files. SolutionTransfer would then be able to re-create the solution
+vector upon restart in a similar way. However, in contrast to
+Particles::ParticleHandler, the Vector class we use for the solution vector
+can actually serialize itself completely, and so we will go with this
+approach and save ourselves the dozen or so additional lines of code.)
+
+Finally, even though we wrote the `serialize()` function above in such
+a way that it also serializes the `triangulation` member variable, in practice
+the call to Triangulation::save() we needed to deal with the particles *also*
+saves the same kind of information, and Triangulation::load() reads it.
+In other words, we are saving redundant information; in the actual
+implementation of the program, we therefore skip the call to
+@code
+  ar & triangulation;
+@endcode
+We do still need to say
+@code
+  ar & particle_handler;
+@endcode
+because the information attached to the cells of the triangulation only contains
+information about the particles themselves, whereas the previous line is
+necessary to store information such as how many particles there are, what the
+next unused particle index is, and other internal information about the class.
+
+
+<h3>Checkpointing strategies</h3>
+
+Having discussed the general idea of checkpoint/restart, let us turn
+to some more specific questions one has to answer: (i) What do we
+actually want to save/restore? (ii) How often do we want to write
+checkpoints?
+
+
+<h4>What to save/restore</h4>
+
+We will base this tutorial on step-19, and so let us use it as an
+example in this section. Recall that that program simulates an
+electric field in which particles move from the electrode on one
+side to the other side of the domain, i.e., we have both field-based
+and particle-based information to store.
+
+Recall the main class of step-19, which
+had quite a lot of member variables one might want to
+(de)serialize:
+@code
+  template <int dim>
+  class CathodeRaySimulator
+  {
+  public:
+    CathodeRaySimulator();
+
+    void run();
+
+  private:
+    [... member functions ...]
+
+    Triangulation<dim>        triangulation;
+    MappingQGeneric<dim>      mapping;
+    FE_Q<dim>                 fe;
+    DoFHandler<dim>           dof_handler;
+    AffineConstraints<double> constraints;
+
+    SparseMatrix<double>      system_matrix;
+    SparsityPattern           sparsity_pattern;
+
+    Vector<double>            solution;
+    Vector<double>            system_rhs;
+
+    Particles::ParticleHandler<dim> particle_handler;
+    types::particle_index           next_unused_particle_id;
+    types::particle_index           n_recently_lost_particles;
+    types::particle_index           n_total_lost_particles;
+    types::particle_index           n_particles_lost_through_anode;
+
+    DiscreteTime time;
+  };
+@endcode
+Do we really need to save all of these to disk? That would presumably
+lead to quite a lot of data that needs to be stored and, if necessary,
+re-loaded.
+
+In practice, one does not save all of this information, but only what
+cannot be reasonably re-computed in different ways. What is saved
+should also depend on also *where* in the program's algorithm one
+currently is, and one generally finds a convenient point at which not
+so much data needs to be stored. For the
+current example of step-19, a time dependent problem, one could apply
+the following considerations:
+
+- The program runs with the same finite element every time, so there
+  is no need to actually save the element: We know what polynomial
+  degree we want, and can just re-generate the element upon
+  restart. If the polynomial degree was a run-time parameter, then
+  maybe we should serialize the polynomial degree rather than all of
+  the myriad data structures that characterize a FE_Q object, given
+  that we can always re-generate the object by just knowing its
+  polynomial degree. This is the classical trade-off of space vs time:
+  We can achieve the same outcome by saving far less data if we are
+  willing to offer a bit of CPU time to regenerate all of the internal
+  data structures of the FE_Q given the polynomial degree.
+
+- We rebuild the matrix and sparsity pattern in each time step from
+  the DoFHandler and the finite element. These are quite large data
+  structures, but they are conceptually easy to re-create again as
+  necessary. So they do not need to be saved to disk, and this is
+  going to save quite a lot of space. Furthermore, we really only need
+  the matrix for the linear solve; once we are done with the linear
+  solve in the `solve_field()` function, the contents of the matrix
+  are no longer used and are, indeed, overwritten in the next time
+  step. As a consequence, there would really only be a point in saving
+  the matrix if we did the checkpointing between the assembly and the
+  linear solve -- but maybe that is just not a convenient point for
+  this operation, and we should pick a better location. In practice,
+  one generally puts the checkpointing at either the very end or the
+  very beginning of the time stepping loop, given that this is the
+  point where the number of variables whose values are currently
+  active is minimal.
+
+- We also do not need to save the DoFHandler object: If we know the
+  triangulation, we can always just create a DoFHandler object during
+  restart to enumerate degrees of freedom in the same way as we did
+  the last time before a previous program run was checkpointed. In
+  fact, the example implementation of the `checkpoint()` function
+  shown above did not serialize the DoFHandler object for this very
+  reason. On the other hand, we probably do want to save the
+  Triangulation here given that the triangulation is not statically
+  generated once at the beginning of the program and then never
+  changed, but is dynamically adapted every few time steps. In order
+  to re-generate the triangulation, we would therefore have to save
+  which cells were refined/coarsened and when (i.e., the *history* of
+  the triangulation), and this would likely cost substantially more
+  disk space for long-running computations than just saving the
+  triangulation itself.
+
+Similar considerations can be applied to all member variables: Can we
+re-generate their values with relatively little effort (in which case
+they do not have to be saved) or is their state difficult or
+impossible to re-generate if it is not saved to disk (in which case
+the variable needs to be serialized)?
+
+@note If you have carefully read step-19, you might now realize that
+  strictly speaking, we do not *need* to checkpoint to solution vector.
+  This is because the solution vector represents the electric field,
+  which the program solves for at the beginning of each timestep and
+  that this solve does not make reference to the electric field at
+  previous time steps -- in other words, the electric field is not
+  a "history variable": If we know the domain, the mesh, the finite
+  element, and the positions of the particles, we can recompute the
+  solution vector, and consequently we would not have to save it
+  into the checkpoint file. However, this is perhaps more work than
+  we want to do for checkpointing (which you will see is otherwise
+  rather little code) and so, for pedagological purposes, we simply
+  save the solution vector along with the other variables that
+  actually do represent the history of the program.
+
+
+<h4>How precisely should we save the data of a checkpoint</h4>
+
+Recall that the goal of checkpointing is to end up with a safe copy of
+where the program currently is in its computations. As a consequence,
+we need to make sure that we do not end up in a situation where, for
+example, we start overwriting the previous checkpoint file and
+somewhere halfway through the serialization process, the machine
+crashes and we end up with an aborted program and no functional
+checkpoint file.
+
+Instead, the procedure one generally follows to guard against this
+kind of scenario is that checkpoints are written into a file that is
+*separate* from the previous checkpoint file; only once we are past
+the writing process and the file is safely on disk can we replace the
+previous checkpoint file by the one just written -- that is, we *move*
+the new file into place of the old one. You will see in the code how
+this two-step process is implemented in practice.
+
+The situation is made slightly more complicated by the fact that
+in the program below, a "checkpoint" actually consists of a number
+of files -- one file into which we write the program's member
+variables, and several into which the triangulation puts its
+information. We then would have to rename several files,
+preferrably as a single, "atomic" operation that cannot be
+interrupted. Implementing this is tricky and non-trivial (though
+possible), and so we will not show this part and instead just
+assume that nothing will happen between renaming the first
+and the last of the files -- maybe not a great strategy in
+general, but good enough for this tutorial program.
+
+
+<h4>How often to save/restore</h4>
+
+Now that we know *what* we want to save and how we want to restore it,
+we need to answer the question *how often* we want to checkpoint the
+program. At least theoretically, this question has been answered many
+decades ago already, see @cite Young1974 and @cite Daly2006. In
+practice (as actually also in these theoretical derivations), it comes
+down to (i) how long it takes to checkpoint data, and (ii) how
+frequently we expect that the stored data will have to be used, i.e.,
+how often the system crashes.
+
+For example, if it takes five minutes to save the state of the
+program, then we probably do not want to write a checkpoint every ten
+minutes. On the other hand, if it only takes five seconds, then maybe
+ten minutes is a reasonable frequency if we run on a modest 100 cores
+and the machine does not crash very often, given that in that case the
+overhead is only approximately 1%. Finally, if it takes five seconds
+to save the state, but we are running on 100,000 processes (i.e., a
+very expensive simulation) and the machine frequently crashes, then
+maybe we are willing to offer a 5% penalty in the overall run time and
+write a checkpoint every minute and a half given that we lose far less
+work this way on average if the machine crashes at an unpredictable
+moment in our computations. The papers cited above essentially just
+formalize this sort of consideration.
+
+In the program, we will not dwell on this and simply choose an ad-hoc
+value of saving the state every ten time steps: That's too often in
+practice, but is useful for experiencing how this works in
+practice without having to run the program too long.
diff --git a/examples/step-83/doc/kind b/examples/step-83/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-83/doc/results.dox b/examples/step-83/doc/results.dox
new file mode 100644 (file)
index 0000000..75c2676
--- /dev/null
@@ -0,0 +1,226 @@
+<h1>Results</h1>
+
+When you run this program, it produces the following output that is
+almost exactly identical to what you get from step-19:
+@code
+Timestep 1
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  20
+  Number of particles lost this time step:  0
+
+  Now at t=2.12647e-07, dt=2.12647e-07.
+
+Timestep 2
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  40
+  Number of particles lost this time step:  0
+
+  Now at t=4.14362e-07, dt=2.01715e-07.
+
+Timestep 3
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  60
+  Number of particles lost this time step:  0
+
+  Now at t=5.23066e-07, dt=1.08704e-07.
+
+Timestep 4
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  80
+  Number of particles lost this time step:  0
+
+  Now at t=6.08431e-07, dt=8.53649e-08.
+
+Timestep 5
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  100
+  Number of particles lost this time step:  0
+
+  Now at t=6.81935e-07, dt=7.35039e-08.
+
+Timestep 6
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  120
+  Number of particles lost this time step:  0
+
+  Now at t=7.47864e-07, dt=6.59294e-08.
+
+Timestep 7
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  140
+  Number of particles lost this time step:  0
+
+  Now at t=8.2516e-07, dt=7.72957e-08.
+
+Timestep 8
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  158
+  Number of particles lost this time step:  0
+
+  Now at t=8.95325e-07, dt=7.01652e-08.
+
+Timestep 9
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  172
+  Number of particles lost this time step:  0
+
+  Now at t=9.67852e-07, dt=7.25269e-08.
+
+Timestep 10
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  186
+  Number of particles lost this time step:  0
+
+  Now at t=1.03349e-06, dt=6.56398e-08.
+
+--- Writing checkpoint... ---
+
+Timestep 11
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  198
+  Number of particles lost this time step:  0
+
+  Now at t=1.11482e-06, dt=8.13268e-08.
+
+Timestep 12
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  206
+  Number of particles lost this time step:  0
+
+  Now at t=1.18882e-06, dt=7.39967e-08.
+
+Timestep 13
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  212
+  Number of particles lost this time step:  0
+
+  Now at t=1.26049e-06, dt=7.16705e-08.
+
+[...]
+@endcode
+The only thing that is different is the additional line after the tenth
+output (which also appears after the 20th, 30th, etc., time step) indicating
+that a checkpoint file has been written.
+
+Because we chose to use a text-based format for the checkpoint file (which
+you should not do in production codes because it uses quite a lot of disk space),
+we can actually inspect this file. It will look like this, with many many more
+lines:
+@code
+22 serialization::archive 18 0 0 0 0 0 7 0 0 3 1 0
+4989 -1.00000000000000000e+00 -1.00000000000000000e+00 -1.00000000000000000e+00 -9.96108134982226390e-01 -1.00000000000000000e+00 -9.98361082867431748e-01
+[...]
+@endcode
+What each of these numbers represents is hard to tell in practice, and also
+entirely unimporant for our current purposes -- it's
+a representation of the many objects that make up this program's state, and
+from which one can restore its state. The point simply being that this is what
+serialization produces: A long list (a sequence) of bits that we can put
+into a file, and that we can later read again to recreate the state of the
+program.
+
+Now here's the fun part. Let's say you hit Control-C on the command
+line at the point above (say, during time step 13 or 14). There's
+a set of checkpoint files on disk that saved the state after ten time steps.
+Based on the logic in `main()`, we should be able to restart from
+this point if we run the program with
+```
+  ./step-83 restart
+```
+Indeed, this is what then happens:
+@code
+--- Restarting at t=1.03349e-06, dt=6.56398e-08
+
+Timestep 11
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  198
+  Number of particles lost this time step:  0
+
+  Now at t=1.11482e-06, dt=8.13268e-08.
+
+Timestep 12
+  Field degrees of freedom:                 4989
+  Total number of particles in simulation:  206
+  Number of particles lost this time step:  0
+
+  Now at t=1.18882e-06, dt=7.39967e-08.
+
+[...]
+@endcode
+After the status message that shows that we are restarting, this is
+indeed *the exact same output* for the following time steps we had
+gotten previously -- in other words, saving the complete state
+seems to have worked, and the program has continued as if it had
+never been interrupted!
+
+
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+<h4> Making efficiency a priority </h4>
+
+While the techniques we have shown here are fully sufficient for what
+we want to do, namely checkpoint and restart computations, and are in
+fact also fully sufficient for much larger code bases such as
+[ASPECT](https://aspect.geodynamics.org), one could go beyond what is
+still a relatively simple scheme.
+
+Specifically, among the things we need to recognize is that writing
+large amounts of data to disk is expensive and can take a good long
+time to finish -- for example, for large parallel computations with,
+say, a billion unknowns, checkpoints can run into the hundred gigabyte
+range or beyond. One may ask whether that could be avoided, or at
+least whether we can mitigate the cost.
+
+One way to do that is to first serialize the state of the program into
+a buffer in memory (like the `Archive` objects the `serialize()`
+functions write to and read from), and once that is done, start a *separate*
+thread do the writing while the rest of the program continues with computations.
+This is useful because writing the data to disk often takes a
+long time but not a lot of CPU power: It just takes time to move the
+data through the network to the file server, and from there onto the
+actual disks. This is something that might as well happen while we are
+doing something useful again (namely, solving more time steps). Should
+the machine crash during this phase, nothing is lost: As discussed in
+the introduction, we are writing the checkpoint into a temporary file
+(which will be lost in the case of a machine failure), but we have
+kept the previous checkpoint around until we know that the temporary
+file is complete and can be moved over the old one.
+
+The only thing we have to pay attention in this background-writing
+scheme is that we cannot start with creating a new checkpoint while
+the previous one is still being written in the background.
+
+Doing this all is not technically very difficult; the code just
+requires more explanation than lines of code, and so we omit doing
+this in the program here. But you can take a look at the
+`MainLoop::output()` function of step-69 to see how such a code looks
+like.
+
+A variation of this general approach is that each process writes its
+data immediately, but into files that are held on fast file systems --
+say, a node-local SSD rather than a file server shared by the entire
+cluster. One would then just tell the operating system to move this
+file to the centeral file server in a second step, and this step can
+happen in the background at whatever speed the operating system can
+provide. Or perhaps one leaves *most* of these files on the fast local
+file system in hopes that the restart happens before these files are
+purged (say, by a script that runs nightly) and only moves these files
+to the permanent file system every tenth time we create a checkpoint.
+
+In all of these cases, the logic quickly becomes quite complicated. As
+usual, the solution is not to re-invent the wheel: Libraries such as
+[VeloC](https://www.anl.gov/mcs/veloc-very-low-overhead-transparent-multilevel-checkpointrestart),
+developed by the Exascale Computing Project (ECP) already do all of
+this and more, for codes that are orders of magnitude more complex
+than the little example here.
+
+Separately, one might want to try to reduce the amount of time it
+takes to serialize objects into a buffer in memory. As mentioned
+above, we use the
+[BOOST serialization
+library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html)
+for this task, but it is not the only player in town. One could, for
+example, use the [bitsery](https://github.com/fraillt/bitsery), [cereal](https://github.com/USCiLab/cereal), or
+[zpp](https://github.com/eyalz800/serializer) projects instead, which can be substantially faster than BOOST.
diff --git a/examples/step-83/doc/tooltip b/examples/step-83/doc/tooltip
new file mode 100644 (file)
index 0000000..20eef3c
--- /dev/null
@@ -0,0 +1 @@
+Serialization, checkpoint/restart.
diff --git a/examples/step-83/step-83.cc b/examples/step-83/step-83.cc
new file mode 100644 (file)
index 0000000..ea01a24
--- /dev/null
@@ -0,0 +1,1013 @@
+/* ------------------------------------------------------------------------
+ *
+ * SPDX-License-Identifier: LGPL-2.1-or-later
+ * Copyright (C) 2023 - 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.
+ *
+ * ------------------------------------------------------------------------
+ *
+ * Author: Pasquale Africa, SISSA, 2024,
+ *         Wolfgang Bangerth, Colorado State University, 2024,
+ *         Bruno Blais, Polytechnique Montreal, 2024.
+ */
+
+
+// @sect3{Include files}
+
+// This program, with the exception of the checkpointing component
+// is identical to step-19, and so the following include files are
+// all the same:
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/discrete_time.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.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/particles/particle_handler.h>
+#include <deal.II/particles/data_out.h>
+
+// The only thing new are the following two include files. They are the ones
+// that declare the classes we use as archives for reading (`iarchive` = input
+// archive) and writing (`oarchive` = output archive) serialized data:
+#include <boost/archive/text_iarchive.hpp>
+#include <boost/archive/text_oarchive.hpp>
+
+#include <filesystem>
+#include <fstream>
+#include <string>
+
+// @sect3{Global definitions}
+
+// As is customary, we put everything that corresponds to the details of the
+// program into a namespace of its own.
+namespace Step83
+{
+  using namespace dealii;
+
+  namespace BoundaryIds
+  {
+    constexpr types::boundary_id open          = 101;
+    constexpr types::boundary_id cathode       = 102;
+    constexpr types::boundary_id focus_element = 103;
+    constexpr types::boundary_id anode         = 104;
+  } // namespace BoundaryIds
+
+  namespace Constants
+  {
+    constexpr double electron_mass   = 9.1093837015e-31;
+    constexpr double electron_charge = 1.602176634e-19;
+
+    constexpr double V0 = 1;
+
+    constexpr double E_threshold = 0.05;
+
+    constexpr double electrons_per_particle = 3e15;
+  } // namespace Constants
+
+
+  // @sect3{The main class}
+
+  // The following is then the main class of this program. It is,
+  // fundamentally, identical to step-19 with the exception of
+  // the `checkpoint()` and `restart()` functions, along with the
+  // `serialize()` function we use to serialize and deserialize the
+  // data this class stores. The `serialize()` function is called
+  // by the BOOST serialization framework, and consequently has to
+  // have exactly the set of arguments used here. Furthermore, because
+  // it is called by BOOST functions, it has to be `public`; the other
+  // two new functions are as always made `private`.
+  //
+  // The `run()` function has also been modified to enable simulation restart
+  // via its new argument `do_restart` that indicates whether or not to
+  // start the simulation from a checkpoint.
+  template <int dim>
+  class CathodeRaySimulator
+  {
+  public:
+    CathodeRaySimulator();
+
+    void run(const bool do_restart);
+
+    template <class Archive>
+    void serialize(Archive &ar, const unsigned int version);
+
+  private:
+    void make_grid();
+    void setup_system();
+    void assemble_system();
+    void solve_field();
+    void refine_grid();
+
+    void create_particles();
+    void move_particles();
+    void track_lost_particle(
+      const Particles::ParticleIterator<dim>                  &particle,
+      const typename Triangulation<dim>::active_cell_iterator &cell);
+
+    void update_timestep_size();
+    void output_results() const;
+
+    void checkpoint();
+    void restart();
+
+    Triangulation<dim>        triangulation;
+    const MappingQ<dim>       mapping;
+    const FE_Q<dim>           fe;
+    DoFHandler<dim>           dof_handler;
+    AffineConstraints<double> constraints;
+
+    SparseMatrix<double> system_matrix;
+    SparsityPattern      sparsity_pattern;
+
+    Vector<double> solution;
+    Vector<double> system_rhs;
+
+    Particles::ParticleHandler<dim> particle_handler;
+    types::particle_index           next_unused_particle_id;
+    types::particle_index           n_recently_lost_particles;
+    types::particle_index           n_total_lost_particles;
+    types::particle_index           n_particles_lost_through_anode;
+
+    DiscreteTime time;
+  };
+
+
+
+  // @sect3{The <code>CathodeRaySimulator</code> class implementation}
+
+  // @sect4{The unchanged parts of the class}
+
+  // Let us start with those parts of the class that are all unchanged
+  // from step-19 and about which you can learn there. We will
+  // then pick up with commentary again when we get to the two new
+  // functions, `checkpoint()` and `restart()`, along with how the
+  // `run()` function needs to be modified:
+  template <int dim>
+  CathodeRaySimulator<dim>::CathodeRaySimulator()
+    : mapping(1)
+    , fe(2)
+    , dof_handler(triangulation)
+    , particle_handler(triangulation, mapping, /*n_properties=*/dim)
+    , next_unused_particle_id(0)
+    , n_recently_lost_particles(0)
+    , n_total_lost_particles(0)
+    , n_particles_lost_through_anode(0)
+    , time(0, 1e-4)
+  {
+    particle_handler.signals.particle_lost.connect(
+      [this](const typename Particles::ParticleIterator<dim>         &particle,
+             const typename Triangulation<dim>::active_cell_iterator &cell) {
+        this->track_lost_particle(particle, cell);
+      });
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::make_grid()
+  {
+    static_assert(dim == 2,
+                  "This function is currently only implemented for 2d.");
+
+    const double       delta = 0.5;
+    const unsigned int nx    = 5;
+    const unsigned int ny    = 3;
+
+    const std::vector<Point<dim>> vertices //
+      = {{0, 0},
+         {1, 0},
+         {2, 0},
+         {3, 0},
+         {4, 0},
+         {delta, 1},
+         {1, 1},
+         {2, 1},
+         {3, 1},
+         {4, 1},
+         {0, 2},
+         {1, 2},
+         {2, 2},
+         {3, 2},
+         {4, 2}};
+    AssertDimension(vertices.size(), nx * ny);
+
+    const std::vector<unsigned int> cell_vertices[(nx - 1) * (ny - 1)] = {
+      {0, 1, nx + 0, nx + 1},
+      {1, 2, nx + 1, nx + 2},
+      {2, 3, nx + 2, nx + 3},
+      {3, 4, nx + 3, nx + 4},
+
+      {5, nx + 1, 2 * nx + 0, 2 * nx + 1},
+      {nx + 1, nx + 2, 2 * nx + 1, 2 * nx + 2},
+      {nx + 2, nx + 3, 2 * nx + 2, 2 * nx + 3},
+      {nx + 3, nx + 4, 2 * nx + 3, 2 * nx + 4}};
+
+    std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
+    for (unsigned int i = 0; i < cells.size(); ++i)
+      {
+        cells[i].vertices    = cell_vertices[i];
+        cells[i].material_id = 0;
+      }
+
+    GridTools::consistently_order_cells(cells);
+    triangulation.create_triangulation(
+      vertices,
+      cells,
+      SubCellData()); // No boundary information
+
+    triangulation.refine_global(2);
+
+    for (auto &cell : triangulation.active_cell_iterators())
+      for (auto &face : cell->face_iterators())
+        if (face->at_boundary())
+          {
+            if ((face->center()[0] > 0) && (face->center()[0] < 0.5) &&
+                (face->center()[1] > 0) && (face->center()[1] < 2))
+              face->set_boundary_id(BoundaryIds::cathode);
+            else if ((face->center()[0] > 0) && (face->center()[0] < 2))
+              face->set_boundary_id(BoundaryIds::focus_element);
+            else if ((face->center()[0] > 4 - 1e-12) &&
+                     ((face->center()[1] > 1.5) || (face->center()[1] < 0.5)))
+              face->set_boundary_id(BoundaryIds::anode);
+            else
+              face->set_boundary_id(BoundaryIds::open);
+          }
+
+    triangulation.refine_global(1);
+  }
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    solution.reinit(dof_handler.n_dofs());
+    system_rhs.reinit(dof_handler.n_dofs());
+
+    constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             BoundaryIds::cathode,
+                                             Functions::ConstantFunction<dim>(
+                                               -Constants::V0),
+                                             constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             BoundaryIds::focus_element,
+                                             Functions::ConstantFunction<dim>(
+                                               -Constants::V0),
+                                             constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             BoundaryIds::anode,
+                                             Functions::ConstantFunction<dim>(
+                                               +Constants::V0),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    dsp,
+                                    constraints,
+                                    /* keep_constrained_dofs = */ false);
+    sparsity_pattern.copy_from(dsp);
+
+    system_matrix.reinit(sparsity_pattern);
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::assemble_system()
+  {
+    system_matrix = 0;
+    system_rhs    = 0;
+
+    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.dofs_per_cell;
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(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_rhs    = 0;
+
+        fe_values.reinit(cell);
+
+        for (const unsigned int q_index : 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_index) * // grad phi_i(x_q)
+                   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
+                   fe_values.JxW(q_index));           // dx
+            }
+
+        if (particle_handler.n_particles_in_cell(cell) > 0)
+          for (const auto &particle : particle_handler.particles_in_cell(cell))
+            {
+              const Point<dim> &reference_location =
+                particle.get_reference_location();
+              for (const unsigned int i : fe_values.dof_indices())
+                cell_rhs(i) +=
+                  (fe.shape_value(i, reference_location) * // phi_i(x_p)
+                   (-Constants::electrons_per_particle *   // N
+                    Constants::electron_charge));          // e
+            }
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global(
+          cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+      }
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::solve_field()
+  {
+    SolverControl            solver_control(1000, 1e-12);
+    SolverCG<Vector<double>> solver(solver_control);
+
+    PreconditionSSOR<SparseMatrix<double>> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+
+    solver.solve(system_matrix, solution, system_rhs, preconditioner);
+
+    constraints.distribute(solution);
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::refine_grid()
+  {
+    Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+
+    KellyErrorEstimator<dim>::estimate(dof_handler,
+                                       QGauss<dim - 1>(fe.degree + 1),
+                                       {},
+                                       solution,
+                                       estimated_error_per_cell);
+
+    GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.1,
+                                                    0.03);
+
+    triangulation.execute_coarsening_and_refinement();
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::create_particles()
+  {
+    FEFaceValues<dim> fe_face_values(fe,
+                                     QMidpoint<dim - 1>(),
+                                     update_quadrature_points |
+                                       update_gradients |
+                                       update_normal_vectors);
+
+    std::vector<Tensor<1, dim>> solution_gradients(
+      fe_face_values.n_quadrature_points);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary() &&
+            (face->boundary_id() == BoundaryIds::cathode))
+          {
+            fe_face_values.reinit(cell, face);
+
+            const FEValuesExtractors::Scalar electric_potential(0);
+            fe_face_values[electric_potential].get_function_gradients(
+              solution, solution_gradients);
+            for (const unsigned int q_point :
+                 fe_face_values.quadrature_point_indices())
+              {
+                const Tensor<1, dim> E = solution_gradients[q_point];
+
+                if ((E * fe_face_values.normal_vector(q_point) < 0) &&
+                    (E.norm() > Constants::E_threshold))
+                  {
+                    const Point<dim> &location =
+                      fe_face_values.quadrature_point(q_point);
+
+                    Particles::Particle<dim> new_particle;
+                    new_particle.set_location(location);
+                    new_particle.set_reference_location(
+                      mapping.transform_real_to_unit_cell(cell, location));
+                    new_particle.set_id(next_unused_particle_id);
+                    particle_handler.insert_particle(new_particle, cell);
+
+                    ++next_unused_particle_id;
+                  }
+              }
+          }
+
+    particle_handler.update_cached_numbers();
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::move_particles()
+  {
+    const double dt = time.get_next_step_size();
+
+    Vector<double>            solution_values(fe.n_dofs_per_cell());
+    FEPointEvaluation<1, dim> evaluator(mapping, fe, update_gradients);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (particle_handler.n_particles_in_cell(cell) > 0)
+        {
+          const typename Particles::ParticleHandler<
+            dim>::particle_iterator_range particles_in_cell =
+            particle_handler.particles_in_cell(cell);
+
+          std::vector<Point<dim>> particle_positions;
+          for (const auto &particle : particles_in_cell)
+            particle_positions.push_back(particle.get_reference_location());
+
+          cell->get_dof_values(solution, solution_values);
+
+          evaluator.reinit(cell, particle_positions);
+          evaluator.evaluate(make_array_view(solution_values),
+                             EvaluationFlags::gradients);
+
+          {
+            typename Particles::ParticleHandler<dim>::particle_iterator
+              particle = particles_in_cell.begin();
+            for (unsigned int particle_index = 0;
+                 particle != particles_in_cell.end();
+                 ++particle, ++particle_index)
+              {
+                const Tensor<1, dim> &E =
+                  evaluator.get_gradient(particle_index);
+
+                const Tensor<1, dim> old_velocity(particle->get_properties());
+
+                const Tensor<1, dim> acceleration =
+                  Constants::electron_charge / Constants::electron_mass * E;
+
+                const Tensor<1, dim> new_velocity =
+                  old_velocity + acceleration * dt;
+
+                particle->set_properties(new_velocity);
+
+                const Point<dim> new_location =
+                  particle->get_location() + dt * new_velocity;
+                particle->set_location(new_location);
+              }
+          }
+        }
+
+    particle_handler.sort_particles_into_subdomains_and_cells();
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::track_lost_particle(
+    const typename Particles::ParticleIterator<dim>         &particle,
+    const typename Triangulation<dim>::active_cell_iterator &cell)
+  {
+    ++n_recently_lost_particles;
+    ++n_total_lost_particles;
+
+    const Point<dim> current_location              = particle->get_location();
+    const Point<dim> approximate_previous_location = cell->center();
+
+    if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
+      {
+        const Tensor<1, dim> direction =
+          (current_location - approximate_previous_location) /
+          (current_location[0] - approximate_previous_location[0]);
+
+        const double right_boundary_intercept =
+          approximate_previous_location[1] +
+          (4 - approximate_previous_location[0]) * direction[1];
+        if ((right_boundary_intercept > 0.5) &&
+            (right_boundary_intercept < 1.5))
+          ++n_particles_lost_through_anode;
+      }
+  }
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::update_timestep_size()
+  {
+    if (time.get_step_number() > 0)
+      {
+        double min_cell_size_over_velocity = std::numeric_limits<double>::max();
+
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          if (particle_handler.n_particles_in_cell(cell) > 0)
+            {
+              const double cell_size = cell->minimum_vertex_distance();
+
+              double max_particle_velocity(0.0);
+
+              for (const auto &particle :
+                   particle_handler.particles_in_cell(cell))
+                {
+                  const Tensor<1, dim> velocity(particle.get_properties());
+                  max_particle_velocity =
+                    std::max(max_particle_velocity, velocity.norm());
+                }
+
+              if (max_particle_velocity > 0)
+                min_cell_size_over_velocity =
+                  std::min(min_cell_size_over_velocity,
+                           cell_size / max_particle_velocity);
+            }
+
+        constexpr double c_safety = 0.5;
+        time.set_desired_next_step_size(c_safety * 0.5 *
+                                        min_cell_size_over_velocity);
+      }
+    else
+      {
+        const QTrapezoid<dim> vertex_quadrature;
+        FEValues<dim> fe_values(fe, vertex_quadrature, update_gradients);
+
+        std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
+
+        double min_timestep = std::numeric_limits<double>::max();
+
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          if (particle_handler.n_particles_in_cell(cell) > 0)
+            {
+              const double cell_size = cell->minimum_vertex_distance();
+
+              fe_values.reinit(cell);
+              fe_values.get_function_gradients(solution, field_gradients);
+
+              double max_E = 0;
+              for (const auto q_point : fe_values.quadrature_point_indices())
+                max_E = std::max(max_E, field_gradients[q_point].norm());
+
+              if (max_E > 0)
+                min_timestep =
+                  std::min(min_timestep,
+                           std::sqrt(0.5 * cell_size *
+                                     Constants::electron_mass /
+                                     Constants::electron_charge / max_E));
+            }
+
+        time.set_desired_next_step_size(min_timestep);
+      }
+  }
+
+
+
+  template <int dim>
+  class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
+  {
+  public:
+    ElectricFieldPostprocessor()
+      : DataPostprocessorVector<dim>("electric_field", update_gradients)
+    {}
+
+    virtual void evaluate_scalar_field(
+      const DataPostprocessorInputs::Scalar<dim> &input_data,
+      std::vector<Vector<double>> &computed_quantities) const override
+    {
+      AssertDimension(input_data.solution_gradients.size(),
+                      computed_quantities.size());
+
+      for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
+        {
+          AssertDimension(computed_quantities[p].size(), dim);
+          for (unsigned int d = 0; d < dim; ++d)
+            computed_quantities[p][d] = input_data.solution_gradients[p][d];
+        }
+    }
+  };
+
+
+
+  template <int dim>
+  void CathodeRaySimulator<dim>::output_results() const
+  {
+    {
+      ElectricFieldPostprocessor<dim> electric_field;
+      DataOut<dim>                    data_out;
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(solution, "electric_potential");
+      data_out.add_data_vector(solution, electric_field);
+      data_out.build_patches();
+
+      DataOutBase::VtkFlags output_flags;
+      output_flags.time  = time.get_current_time();
+      output_flags.cycle = time.get_step_number();
+      output_flags.physical_units["electric_potential"] = "V";
+      output_flags.physical_units["electric_field"]     = "V/m";
+
+      data_out.set_flags(output_flags);
+
+      std::ofstream output("solution-" +
+                           Utilities::int_to_string(time.get_step_number(), 4) +
+                           ".vtu");
+      data_out.write_vtu(output);
+    }
+
+    {
+      Particles::DataOut<dim> particle_out;
+      particle_out.build_patches(
+        particle_handler,
+        std::vector<std::string>(dim, "velocity"),
+        std::vector<DataComponentInterpretation::DataComponentInterpretation>(
+          dim, DataComponentInterpretation::component_is_part_of_vector));
+
+      DataOutBase::VtkFlags output_flags;
+      output_flags.time                       = time.get_current_time();
+      output_flags.cycle                      = time.get_step_number();
+      output_flags.physical_units["velocity"] = "m/s";
+
+      particle_out.set_flags(output_flags);
+
+      std::ofstream output("particles-" +
+                           Utilities::int_to_string(time.get_step_number(), 4) +
+                           ".vtu");
+      particle_out.write_vtu(output);
+    }
+  }
+
+
+
+  // @sect4{CathodeRaySimulator::serialize()}
+
+  // The first of the new function is the one that is called by the
+  // BOOST Serialization framework to serialize and deserialize the
+  // data of this class. It has already been discussed in the introduction
+  // to this program and so does not provide any surprises. All it does is
+  // write those member variables of the current class that cannot be
+  // re-created easily into an archive, or read these members from it.
+  // (Whether `operator &` facilitates a write or read operation depends on
+  // whether the `Archive` type is an output or input archive.)
+  //
+  // The function takes a second argument, `version`, that can be used to
+  // create checkpoints that have version numbers. This is useful if one
+  // evolves programs by adding more member variables, but still wants to
+  // retain the ability to read checkpoint files created with earlier
+  // versions of the program. The `version` variable would, in that case,
+  // be used to represent which version of the program wrote the file,
+  // and if necessary to read only those variables that were written with
+  // that past version, finding a different way to initialize the new member
+  // variables that have been added since then. We will not make use of this
+  // ability here.
+  //
+  // Finally, while the program that indents all deal.II source files format
+  // the following code as
+  // @code
+  //   ar &solution;
+  // @endcode
+  // as if we are taking the address of the `triangulation` variable, the
+  // way you *should* read the code is as
+  // @code
+  //   ar & solution;
+  // @endcode
+  // where `operator &` is a binary operator that could either be interpreted
+  // as `operator <<` for output or `operator >>` for input.
+  //
+  // As discussed in the introduction, we do not serialize the `triangulation`
+  // member variable, instead leaving that to separate calls in the
+  // `checkpoint()` and `restart()` functions below.
+  template <int dim>
+  template <class Archive>
+  void CathodeRaySimulator<dim>::serialize(Archive &ar,
+                                           const unsigned int /* version */)
+  {
+    ar &solution;
+    ar &particle_handler;
+    ar &next_unused_particle_id;
+    ar &n_recently_lost_particles;
+    ar &n_total_lost_particles;
+    ar &n_particles_lost_through_anode;
+    ar &time;
+  }
+
+
+
+  // @sect4{CathodeRaySimulator::checkpoint()}
+
+  // The checkpoint function of the principal class of this program is then
+  // quite straightforward: We create an output file (and check that it is
+  // writable), create an output archive, and then move the serialized
+  // contents of the current object (i.e., the `*this` object) into the
+  // archive. The use of `operator<<` here calls the `serialize()` function
+  // above with an output archive as argument. When the destructor of the
+  // `archive` variable is called at the end of the code block within which
+  // it lives, the entire archive is written into the output file stream it
+  // is associated with.
+  //
+  // As mentioned in the introduction, "real" applications would not use
+  // text-based archives as provided by the `boost::archive::text_oarchive`
+  // class, but use binary and potentially compressed versions. This can
+  // easily be achieved by using differently named classes, and the BOOST
+  // documentation explains how to do that.
+  template <int dim>
+  void CathodeRaySimulator<dim>::checkpoint()
+  {
+    std::cout << "--- Writing checkpoint... ---" << std::endl << std::endl;
+
+    {
+      std::ofstream checkpoint_file("tmp.checkpoint_step_83");
+      AssertThrow(checkpoint_file,
+                  ExcMessage(
+                    "Could not write to the <tmp.checkpoint_step_83> file."));
+
+      boost::archive::text_oarchive archive(checkpoint_file);
+
+      archive << *this;
+    }
+
+    // The second part of the serialization is all of the data that we can
+    // attach to cells -- see the discussion about this in the introduction.
+    // Here, the only data we attach to cells are the particles. We then
+    // let the triangulation save these into a set of files that all start
+    // with the same prefix as we chose above, namely "tmp.checkpoint":
+    particle_handler.prepare_for_serialization();
+    triangulation.save("tmp.checkpoint");
+
+
+    // At this point, the serialized data of this file has ended up in a number
+    // of files that all start with `tmp.checkpoint` file. As mentioned in the
+    // introduction, we do not want to directly overwrite the checkpointing
+    // files from the previous checkpoint operation, for fear that the program
+    // may be interrupted *while writing the checkpoint files*. This would
+    // result in corrupted files, and defeat the whole purpose of checkpointing
+    // because one cannot restart from such a file. On the other hand, if we got
+    // here, we know that the "tmp.checkpoint*" files were successfully written,
+    // and we can rename it to "checkpoint*", in the process replacing the old
+    // file.
+    //
+    // We do this move operation by calling the
+    // [C++ function that does the renaming of
+    // files](https://en.cppreference.com/w/cpp/filesystem/rename).
+    // Note that it is documented as being for all practical purposes
+    // "atomic", i.e., we do not need to worry that the program may be
+    // interrupted somewhere within the renaming operation itself. Of
+    // course, it is possible that we get interrupted somewhere between
+    // renaming one file and the next, and that presents problems in
+    // itself -- in essence, we want the entire renaming operation of all of
+    // these files to be atomic. With a couple dozen lines of extra code, one
+    // could address this issue (using strategies that databases use frequently:
+    // if one operation fails, we need to
+    // [rollback](https://en.wikipedia.org/wiki/Rollback_(data_management))
+    // the entire transaction). For the purposes of this program, this is
+    // perhaps too much, and we will simply hope that that doesn't happen,
+    // perhaps based on the belief that renaming files is much faster than
+    // writing them, and that unlike writing checkpoint files, renaming does not
+    // require much memory or disk space and so does not risk running out of
+    // either.
+    //
+    // As a consequence, the following code first loops over all files in
+    // the current directory, picks out those that start with the string
+    // "tmp.checkpoint", and puts them into a list. In a second step,
+    // we loop over the list and rename each of these files to one
+    // whose name consists of the "tmp.checkpoint*" file but stripped off
+    // its first four characters (i.e., only the "checkpoint*" part). We use
+    // this approach, rather than listing the files we want to rename,
+    // because we do not actually know the names of the files written by
+    // the Triangulation::save() function, though we know how each of these
+    // file names starts.
+    std::list<std::string> tmp_checkpoint_files;
+    for (const auto &dir_entry : std::filesystem::directory_iterator("."))
+      if (dir_entry.is_regular_file() &&
+          (dir_entry.path().filename().string().find("tmp.checkpoint") == 0))
+        tmp_checkpoint_files.push_back(dir_entry.path().filename().string());
+
+    for (const std::string &filename : tmp_checkpoint_files)
+      std::filesystem::rename(filename, filename.substr(4, std::string::npos));
+  }
+
+
+  // @sect4{CathodeRaySimulator::restart()}
+
+  // The restart function of this class then simply does the opposite:
+  // It opens an input file (and triggers an error if that file cannot be
+  // opened), associates an input archive with it, and then reads the
+  // contents of the current object from it, again using the
+  // `serialize()` function from above. Clearly, since we have written
+  // data into a text-based archive above, we need to use the corresponding
+  // `boost::archive::text_iarchive` class for reading.
+  //
+  // In a second step, we ask the triangulation to read in cell-attached
+  // data, and then tell the Particles::ParticleHandler object to re-create
+  // its information about all of the particles from the data just read.
+  //
+  // The function ends by printing a status message about having restarted:
+  template <int dim>
+  void CathodeRaySimulator<dim>::restart()
+  {
+    {
+      std::ifstream checkpoint_file("checkpoint_step_83");
+      AssertThrow(checkpoint_file,
+                  ExcMessage(
+                    "Could not read from the <checkpoint_step_83> file."));
+
+      boost::archive::text_iarchive archive(checkpoint_file);
+      archive >> *this;
+    }
+
+    triangulation.load("checkpoint");
+    particle_handler.deserialize();
+
+    std::cout << "--- Restarting at t=" << time.get_current_time()
+              << ", dt=" << time.get_next_step_size() << std::endl
+              << std::endl;
+  }
+
+
+  // @sect4{CathodeRaySimulator::run()}
+
+  // The last member function of the principal class of this program is then the
+  // driver. The driver takes a single argument to indicate if the simulation
+  // is a restart. If it is not a restart, the mesh is set up and the problem is
+  // solved like in step-19. If it is a restart, then we read in everything that
+  // is a history variable from the checkpoint file via the `restart()`
+  // function. Recall that everything that is inside the `if` block at the top
+  // of the function is exactly like in step-19, as is almost everything that
+  // follows:
+  template <int dim>
+  void CathodeRaySimulator<dim>::run(const bool do_restart)
+  {
+    if (do_restart == false)
+      {
+        make_grid();
+
+        const unsigned int n_pre_refinement_cycles = 3;
+        for (unsigned int refinement_cycle = 0;
+             refinement_cycle < n_pre_refinement_cycles;
+             ++refinement_cycle)
+          {
+            setup_system();
+            assemble_system();
+            solve_field();
+            refine_grid();
+          }
+      }
+    else
+      {
+        restart();
+      }
+
+    setup_system();
+    do
+      {
+        std::cout << "Timestep " << time.get_step_number() + 1 << std::endl;
+        std::cout << "  Field degrees of freedom:                 "
+                  << dof_handler.n_dofs() << std::endl;
+
+        assemble_system();
+        solve_field();
+
+        create_particles();
+        std::cout << "  Total number of particles in simulation:  "
+                  << particle_handler.n_global_particles() << std::endl;
+
+        n_recently_lost_particles = 0;
+        update_timestep_size();
+        move_particles();
+
+        time.advance_time();
+
+        output_results();
+
+        std::cout << "  Number of particles lost this time step:  "
+                  << n_recently_lost_particles << std::endl;
+        if (n_total_lost_particles > 0)
+          std::cout << "  Fraction of particles lost through anode: "
+                    << 1. * n_particles_lost_through_anode /
+                         n_total_lost_particles
+                    << std::endl;
+
+        std::cout << std::endl
+                  << "  Now at t=" << time.get_current_time()
+                  << ", dt=" << time.get_previous_step_size() << '.'
+                  << std::endl
+                  << std::endl;
+
+        // The only other difference between this program and step-19 is that
+        // we checkpoint the simulation every ten time steps:
+        if (time.get_step_number() % 10 == 0)
+          checkpoint();
+      }
+    while (time.is_at_end() == false);
+  }
+} // namespace Step83
+
+
+// @sect3{The <code>main</code> function}
+
+// The final function of the program is then again the `main()` function. Its
+// overall structure is unchanged in all tutorial programs since step-6 and
+// so there is nothing new to discuss about this aspect.
+//
+// The only difference is that we need to figure out whether a restart was
+// requested, or whether the program should simply start from scratch when
+// called. We do this using a command line argument: The `argc` argument to
+// `main()` indicates how many command line arguments were provided when
+// the program was called (counting the name of the program as the zeroth
+// argument), and `argv` is an array of strings with as many elements as
+// `argc` that contains these command line arguments. So if you call
+// the program as
+// @code
+//   ./step-83
+// @endcode
+// then `argc` will be 1, and `argv` will be the array with one element
+// and content `[ "./step-83" ]`. On the other hand, if you call the program
+// as
+// @code
+//   ./step-83 restart
+// @endcode
+// then `argc` will be 2, and `argv` will be the array with two elements
+// and content `[ "./step-83", "restart" ]`. Every other choice should be
+// flagged as an error. The following try block then does exactly this:
+int main(int argc, char *argv[])
+{
+  try
+    {
+      Step83::CathodeRaySimulator<2> cathode_ray_simulator;
+
+      if (argc == 1)
+        cathode_ray_simulator.run(false); // no restart
+      else if ((argc == 2) && (std::string(argv[1]) == "restart"))
+        cathode_ray_simulator.run(true); // yes restart
+      else
+        {
+          std::cerr << "Error: The only allowed command line argument to this\n"
+                       "       program is 'restart'."
+                    << std::endl;
+          return 1;
+        }
+    }
+  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.