--- /dev/null
+##
+# 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()
--- /dev/null
+<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.
--- /dev/null
+techniques
--- /dev/null
+<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.
--- /dev/null
+Serialization, checkpoint/restart.
--- /dev/null
+/* ------------------------------------------------------------------------
+ *
+ * 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;
+}
projects / dealii.git / commitdiff