230In other words, the
Tensor class stores an array of `dim` objects of a
231tensor of one rank lower, and in order to serialize this array, the
232function simply uses the overloaded `operator&` which recognizes that
233`
values` is a
C-style array of fixed length, and what the data type of
234the elements of
this array are. It will then in turn call the
235`serialize()` member function of `
Tensor<rank-1,dim>`. (The recursion
236has a stopping
point which is not of importance
for this discussion. The
237actual implementation of
class Tensor looks different, but you get the
240Depending on whether the `Archive` type used for the
first argument of
241this function is an input or output archive, `operator&` reads from or
242writes into the `
values` member variable. If the
class had more member
243variables, one would list them one after the other. The same is true
244if a class has base classes. For example, here is the corresponding
245function for the
Quadrature class that stores quadrature points and
246weights, and is derived from the
Subscriptor base class:
249template <
class Archive>
260The beauty of having to write only one function
for both serialization
261and deserialization is that obviously one has to read data from the
262archive in the same order as one writes into it, and that is easy
263to get wrong
if one had to write two separate
functions -- but that
264is automatically correct
if one has the same function
for both
267The BOOST serialization library also has mechanisms
if one wants to
268write separate save and load
functions. This is useful
if a
class
269stores a lot of data in one array and then has a cache of commonly
270accessed data on the side
for fast access; one would then only store
271and load the large array, and in the load function re-build the cache
272from the large array. (See also the discussion below on what one
273actually wants to save.)
275Based on these principles, let us consider how one would (de)serialize
276a program such as @ref step_19
"step-19". Recall that its principal
class looked
277essentially as follows:
280 class CathodeRaySimulator
283 CathodeRaySimulator();
290 void assemble_system();
294 void create_particles();
295 void move_particles();
300 void update_timestep_size();
301 void output_results()
const;
324One would then
first write a
325function that allows (de)serialization of the data of
this class
329 template <
class Archive>
330 void CathodeRaySimulator<dim>::serialize(Archive &ar,
335 ar & particle_handler;
336 ar & next_unused_particle_id;
337 ar & n_recently_lost_particles;
338 ar & n_total_lost_particles;
339 ar & n_particles_lost_through_anode;
343As discussed below,
this is not the entire list of member variables: We
344only want to serialize those that we cannot otherwise re-generate.
346Next, we need a function that writes a checkpoint by creating an
"output
347archive" and
using the usual `
operator<<` to put an
object into
this archive. In
348the code of
this program,
this will then look as follows (with some
349minor modifications we will discuss later):
352 void CathodeRaySimulator<dim>::checkpoint()
354 std::ofstream checkpoint_file(
"checkpoint");
355 boost::archive::text_oarchive archive(checkpoint_file,
356 boost::archive::no_header);
362What `
operator<<` does here is to call the serialization
functions of the right hand
363operand (here, the `serialize()` function described above, with an output
364archive as
template argument), and create a serialized representation of the data. In the
365end, serialization has put everything into the
"archive" from which one
366can
extract a (sometimes very
long)
string that one can save in bulk
367into a file, and that is exactly what happens when the destructor of the
368`archive` variable is called.
370BOOST serialization offers different archives, including ones that
371store the data in text format (as we
do above), in binary format, or
372already compressed with something like gzip to minimize the amount of
373space necessary. The specific type of archive to be used is selected
374by the type of the `archive` variable above (and the corresponding
375variable in the `restart()` function of course). This program uses a
376text archive so that you can look at how a serialized representation
377would actually look at, though a
"real" program would of course
try to
378be more space efficient by
using binary (and possibly compressed)
379representations of the data.
381In any
case, the data we have thus created is read in very similarly
382in the following function (again with some minor modifications):
385 void CathodeRaySimulator<dim>::restart()
387 std::ifstream checkpoint_file(
"checkpoint");
389 boost::archive::text_iarchive archive(checkpoint_file,
390 boost::archive::no_header);
396The magic of
this approach is that one doesn
't actually have to write very
397much code to checkpoint or restart programs: It is all hidden in the
398serialization functions of classes such as Triangulation,
399Particles::ParticleHandler, etc., which the deal.II library provides for you.
402<a name="step_83-Serializationofparallelprograms"></a><h4> Serialization of parallel programs </h4>
405The program as provided here is sequential, i.e., it runs on a single
406processor just like the original @ref step_19 "step-19" did. But what if you had a parallel
407program -- say, something like @ref step_40 "step-40" -- that runs in parallel with MPI?
408In that case, serialization and checkpoint/restart becomes more complicated.
409While parallel execution is not something that is of concern to us
410in this program, it is an issue that has influenced the design of how
411deal.II does serialization; as a consequence, we need to talk through
412what makes serialization of parallel programs more difficult in order
413to understand why this program does things the way it does.
415Intuitively, one might simply want to use the same idea as we used
416here except that we let every MPI process serialize its own data, and
417read its own data. This works, but there are some drawbacks:
418- There is a certain subset of data that is replicated across all
419 MPI processes and that would be then written by all processes. An example
420 is the `time` data structure that stores the current time, time step
421 size, and other information, and that better be the same on all processes.
422 Typically, the replicated data isn't a large fraction of a program
's
423 overall memory usage, and perhaps writing it more than once isn't
424 going to be too much of a problem, but it is unsatisfactory anyway
425 to have
this kind of data on disk more than once.
426- One would have to think in detail how exactly one wants to represent
427 the data on disk. One possibility would be
for every
MPI process to write
428 its own file. On the other hand, checkpointing is most useful
for large
429 programs,
using large
numbers of processes -- it is not uncommon to use
430 checkpointing on programs that
run on 10,000 or more processors in
parallel.
431 This would then lead to 10,000 or more files on disk. That
's unpleasant, and
432 probably inefficient as well. We could address this by first letting all
433 the processes serialize into a string in memory (using `std::ostringstream`)
434 and then collating all of these strings into one file. The MPI I/O sub-system
435 has facilities to make that happen, for example, but it will require a bit
436 of thought not the least because the serialized data from each process
437 will likely result in strings of different sizes.
438- Perhaps the most important reason to rethink how one does things in parallel
439 is because, with a little bit of thought, it is possible to checkpoint a
440 program running with @f$N@f$ MPI processes and restart it with @f$M\neq N@f$
441 processes. This may, at first, seem like a pointless exercise, but it is
442 useful if one had, for example, a program that repeatedly refines the mesh
443 and where it is inefficient to run the early refinement steps with a
444 coarse mesh on too many processes, whereas it is too slow to run the later
445 refinement steps with a fine mesh on too few processes.
447In order to address these issues, in particular the last one, the right approach
448is to deviate a bit from the simple scheme of having a `serialize()` function
449that simply serializes/deserializes *everything* into an archive, and then have two
450functions `checkpoint()` and `restart()` that for all practical purposes defer
451all the work to the `serialize()` function. Instead, one splits all data into
453- Data that is tied to the cells of a triangulation. This includes the mesh itself,
454 but also the particles in the Particles::ParticleHandler class, and most
455 importantly the solution vector(s). The way to serialize such data is to
456 attach the data to cells and then let the Triangulation class serialize the attached
457 data along with its own data. If this is done in a way so that we can re-load
458 a triangulation on a different number of processes than the data was written,
459 then this automatically also ensures that we can restore solution vectors
460 and Particles::ParticleHandler objects (and everything else we can attach
461 to the cells of a triangulation) on a different number of processes.
462- Other data. In finite element programs, this data is almost always replicated
463 across processes, and so it is enough if the "root" process (typically the
464 process with MPI rank zero) writes it to disk. Upon restart, the root
465 process reads the data from disk, sends it to all other processes (however many
466 of them there may be), and these then initialize their own copies of the
467 replicated data structures.
469These sorts of considerations have influenced the design of the Triangulation and
470Particles::ParticleHandler classes. In particular, Particles::ParticleHandler's
471`serialize()` function only serializes the
"other data" category, but not the
472actual particles; these can instead be attached to the
triangulation by calling
475files that become part of the checkpoint. Upon restart, we then
first call
477to retrieve the particles from the cells they are attached to.
479(We could, with relatively minimal effort, use the same scheme for the solution
480vector: The
SolutionTransfer class can be used to attach the values of degrees
481of freedom to cells, and then
Triangulation::save() also writes these into
482checkpoint files.
SolutionTransfer would then be able to re-create the solution
483vector upon restart in a similar way. However, in contrast to
484Particles::ParticleHandler, the
Vector class we use for the solution vector
485can actually serialize itself completely, and so we will go with this
486approach and save ourselves the dozen or so additional lines of code.)
488Finally, even though we wrote the `serialize()` function above in such
489a way that it also serializes the `
triangulation` member variable, in practice
490the call to
Triangulation::save() we needed to deal with the particles *also*
491saves the same kind of information, and
Triangulation::load() reads it.
492In other words, we are saving redundant information; in the actual
493implementation of the program, we therefore skip the call to
497We do still need to say
499 ar & particle_handler;
501because the information attached to the cells of the
triangulation only contains
502information about the particles themselves, whereas the previous line is
503necessary to store information such as how many particles there are, what the
504next unused particle index is, and other
internal information about the class.
507<a name="step_83-Checkpointingstrategies"></a><h3>Checkpointing strategies</h3>
510Having discussed the general idea of checkpoint/restart, let us turn
511to some more specific questions one has to answer: (i) What do we
512actually want to save/restore? (ii) How often do we want to write
516<a name="step_83-Whattosaverestore"></a><h4>What to save/restore</h4>
519We will base this tutorial on @ref step_19 "step-19", and so let us use it as an
520example in this section. Recall that that program simulates an
521electric field in which particles move from the electrode on one
522side to the other side of the domain, i.e., we have both field-based
523and particle-based information to store.
525Recall the main class of @ref step_19 "step-19", which
526had quite a lot of member variables one might want to
530 class CathodeRaySimulator
533 CathodeRaySimulator();
561Do we really need to save all of these to disk? That would presumably
562lead to quite a lot of data that needs to be stored and,
if necessary,
565In practice, one does not save all of
this information, but only what
566cannot be reasonably re-computed in different ways. What is saved
567should also depend on also *where* in the program
's algorithm one
568currently is, and one generally finds a convenient point at which not
569so much data needs to be stored. For the
570current example of @ref step_19 "step-19", a time dependent problem, one could apply
571the following considerations:
573- The program runs with the same finite element every time, so there
574 is no need to actually save the element: We know what polynomial
575 degree we want, and can just re-generate the element upon
576 restart. If the polynomial degree was a run-time parameter, then
577 maybe we should serialize the polynomial degree rather than all of
578 the myriad data structures that characterize a FE_Q object, given
579 that we can always re-generate the object by just knowing its
580 polynomial degree. This is the classical trade-off of space vs time:
581 We can achieve the same outcome by saving far less data if we are
582 willing to offer a bit of CPU time to regenerate all of the internal
583 data structures of the FE_Q given the polynomial degree.
585- We rebuild the matrix and sparsity pattern in each time step from
586 the DoFHandler and the finite element. These are quite large data
587 structures, but they are conceptually easy to re-create again as
588 necessary. So they do not need to be saved to disk, and this is
589 going to save quite a lot of space. Furthermore, we really only need
590 the matrix for the linear solve; once we are done with the linear
591 solve in the `solve_field()` function, the contents of the matrix
592 are no longer used and are, indeed, overwritten in the next time
593 step. As a consequence, there would really only be a point in saving
594 the matrix if we did the checkpointing between the assembly and the
595 linear solve -- but maybe that is just not a convenient point for
596 this operation, and we should pick a better location. In practice,
597 one generally puts the checkpointing at either the very end or the
598 very beginning of the time stepping loop, given that this is the
599 point where the number of variables whose values are currently
602- We also do not need to save the DoFHandler object: If we know the
603 triangulation, we can always just create a DoFHandler object during
604 restart to enumerate degrees of freedom in the same way as we did
605 the last time before a previous program run was checkpointed. In
606 fact, the example implementation of the `checkpoint()` function
607 shown above did not serialize the DoFHandler object for this very
608 reason. On the other hand, we probably do want to save the
609 Triangulation here given that the triangulation is not statically
610 generated once at the beginning of the program and then never
611 changed, but is dynamically adapted every few time steps. In order
612 to re-generate the triangulation, we would therefore have to save
613 which cells were refined/coarsened and when (i.e., the *history* of
614 the triangulation), and this would likely cost substantially more
615 disk space for long-running computations than just saving the
616 triangulation itself.
618Similar considerations can be applied to all member variables: Can we
619re-generate their values with relatively little effort (in which case
620they do not have to be saved) or is their state difficult or
621impossible to re-generate if it is not saved to disk (in which case
622the variable needs to be serialized)?
624@note If you have carefully read @ref step_19 "step-19", you might now realize that
625 strictly speaking, we do not *need* to checkpoint to solution vector.
626 This is because the solution vector represents the electric field,
627 which the program solves for at the beginning of each timestep and
628 that this solve does not make reference to the electric field at
629 previous time steps -- in other words, the electric field is not
630 a "history variable": If we know the domain, the mesh, the finite
631 element, and the positions of the particles, we can recompute the
632 solution vector, and consequently we would not have to save it
633 into the checkpoint file. However, this is perhaps more work than
634 we want to do for checkpointing (which you will see is otherwise
635 rather little code) and so, for pedagological purposes, we simply
636 save the solution vector along with the other variables that
637 actually do represent the history of the program.
640<a name="step_83-Howpreciselyshouldwesavethedataofacheckpoint"></a><h4>How precisely should we save the data of a checkpoint</h4>
643Recall that the goal of checkpointing is to end up with a safe copy of
644where the program currently is in its computations. As a consequence,
645we need to make sure that we do not end up in a situation where, for
646example, we start overwriting the previous checkpoint file and
647somewhere halfway through the serialization process, the machine
648crashes and we end up with an aborted program and no functional
651Instead, the procedure one generally follows to guard against this
652kind of scenario is that checkpoints are written into a file that is
653*separate* from the previous checkpoint file; only once we are past
654the writing process and the file is safely on disk can we replace the
655previous checkpoint file by the one just written -- that is, we *move*
656the new file into place of the old one. You will see in the code how
657this two-step process is implemented in practice.
659The situation is made slightly more complicated by the fact that
660in the program below, a "checkpoint" actually consists of a number
661of files -- one file into which we write the program's member
663information. We then would have to rename several files,
664preferably as a single,
"atomic" operation that cannot be
665interrupted. Implementing
this is tricky and non-trivial (though
666possible), and so we will not show
this part and instead just
667assume that
nothing will happen between renaming the
first
668and the last of the files -- maybe not a great strategy in
669general, but good enough
for this tutorial program.
672<a name=
"step_83-Howoftentosaverestore"></a><h4>How often to save/restore</h4>
675Now that we know *what* we want to save and how we want to restore it,
676we need to answer the question *how often* we want to checkpoint the
677program. At least theoretically,
this question has been answered many
678decades ago already, see @cite Young1974 and @cite Daly2006. In
679practice (as actually also in these theoretical derivations), it comes
680down to (i) how
long it takes to checkpoint data, and (ii) how
681frequently we expect that the stored data will have to be used, i.e.,
682how often the system crashes.
684For example,
if it takes five minutes to save the state of the
685program, then we probably
do not want to write a checkpoint every ten
686minutes. On the other hand,
if it only takes five seconds, then maybe
687ten minutes is a reasonable frequency
if we
run on a modest 100 cores
688and the machine does not crash very often, given that in that
case the
689overhead is only approximately 1%. Finally,
if it takes five seconds
690to save the state, but we are running on 100,000 processes (i.e., a
691very expensive simulation) and the machine frequently crashes, then
692maybe we are willing to offer a 5% penalty in the overall
run time and
693write a checkpoint every minute and a half given that we lose far less
694work
this way on average
if the machine crashes at an unpredictable
695moment in our computations. The papers cited above essentially just
696formalize
this sort of consideration.
698In the program, we will not dwell on
this and simply choose an ad-hoc
699value of saving the state every ten time steps: That
's too often in
700practice, but is useful for experiencing how this works in
701practice without having to run the program too long.
704 * <a name="step_83-CommProg"></a>
705 * <h1> The commented program</h1>
708 * <a name="step_83-Includefiles"></a>
709 * <h3>Include files</h3>
713 * This program, with the exception of the checkpointing component
714 * is identical to @ref step_19 "step-19", and so the following include files are
718 * #include <deal.II/base/quadrature_lib.h>
719 * #include <deal.II/base/discrete_time.h>
721 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
722 * #include <deal.II/lac/full_matrix.h>
723 * #include <deal.II/lac/precondition.h>
724 * #include <deal.II/lac/solver_cg.h>
725 * #include <deal.II/lac/sparse_matrix.h>
726 * #include <deal.II/lac/vector.h>
727 * #include <deal.II/lac/affine_constraints.h>
729 * #include <deal.II/grid/tria.h>
730 * #include <deal.II/grid/grid_refinement.h>
731 * #include <deal.II/grid/grid_tools.h>
733 * #include <deal.II/fe/mapping_q.h>
734 * #include <deal.II/matrix_free/fe_point_evaluation.h>
735 * #include <deal.II/fe/fe_q.h>
736 * #include <deal.II/fe/fe_values.h>
738 * #include <deal.II/dofs/dof_handler.h>
739 * #include <deal.II/dofs/dof_tools.h>
741 * #include <deal.II/numerics/data_out.h>
742 * #include <deal.II/numerics/vector_tools.h>
743 * #include <deal.II/numerics/error_estimator.h>
745 * #include <deal.II/particles/particle_handler.h>
746 * #include <deal.II/particles/data_out.h>
750 * The only thing new are the following two include files. They are the ones
751 * that declare the classes we use as archives for reading (`iarchive` = input
752 * archive) and writing (`oarchive` = output archive) serialized data:
755 * #include <boost/archive/text_iarchive.hpp>
756 * #include <boost/archive/text_oarchive.hpp>
758 * #include <filesystem>
765 * <a name="step_83-Globaldefinitions"></a>
766 * <h3>Global definitions</h3>
770 * As is customary, we put everything that corresponds to the details of the
771 * program into a namespace of its own.
776 * using namespace dealii;
778 * namespace BoundaryIds
780 * constexpr types::boundary_id open = 101;
781 * constexpr types::boundary_id cathode = 102;
782 * constexpr types::boundary_id focus_element = 103;
783 * constexpr types::boundary_id anode = 104;
784 * } // namespace BoundaryIds
786 * namespace Constants
788 * constexpr double electron_mass = 9.1093837015e-31;
789 * constexpr double electron_charge = 1.602176634e-19;
791 * constexpr double V0 = 1;
793 * constexpr double E_threshold = 0.05;
795 * constexpr double electrons_per_particle = 3e15;
796 * } // namespace Constants
802 * <a name="step_83-Themainclass"></a>
803 * <h3>The main class</h3>
807 * The following is then the main class of this program. It is,
808 * fundamentally, identical to @ref step_19 "step-19" with the exception of
809 * the `checkpoint()` and `restart()` functions, along with the
810 * `serialize()` function we use to serialize and deserialize the
811 * data this class stores. The `serialize()` function is called
812 * by the BOOST serialization framework, and consequently has to
813 * have exactly the set of arguments used here. Furthermore, because
814 * it is called by BOOST functions, it has to be `public`; the other
815 * two new functions are as always made `private`.
819 * The `run()` function has also been modified to enable simulation restart
820 * via its new argument `do_restart` that indicates whether or not to
821 * start the simulation from a checkpoint.
825 * class CathodeRaySimulator
828 * CathodeRaySimulator();
830 * void run(const bool do_restart);
832 * template <class Archive>
833 * void serialize(Archive &ar, const unsigned int version);
837 * void setup_system();
838 * void assemble_system();
839 * void solve_field();
840 * void refine_grid();
842 * void create_particles();
843 * void move_particles();
844 * void track_lost_particle(
845 * const Particles::ParticleIterator<dim> &particle,
846 * const typename Triangulation<dim>::active_cell_iterator &cell);
848 * void update_timestep_size();
849 * void output_results() const;
854 * Triangulation<dim> triangulation;
855 * const MappingQ<dim> mapping;
856 * const FE_Q<dim> fe;
857 * DoFHandler<dim> dof_handler;
858 * AffineConstraints<double> constraints;
860 * SparseMatrix<double> system_matrix;
861 * SparsityPattern sparsity_pattern;
863 * Vector<double> solution;
864 * Vector<double> system_rhs;
866 * Particles::ParticleHandler<dim> particle_handler;
867 * types::particle_index next_unused_particle_id;
868 * types::particle_index n_recently_lost_particles;
869 * types::particle_index n_total_lost_particles;
870 * types::particle_index n_particles_lost_through_anode;
880 * <a name="step_83-ThecodeCathodeRaySimulatorcodeclassimplementation"></a>
881 * <h3>The <code>CathodeRaySimulator</code> class implementation</h3>
886 * <a name="step_83-Theunchangedpartsoftheclass"></a>
887 * <h4>The unchanged parts of the class</h4>
891 * Let us start with those parts of the class that are all unchanged
892 * from @ref step_19 "step-19" and about which you can learn there. We will
893 * then pick up with commentary again when we get to the two new
894 * functions, `checkpoint()` and `restart()`, along with how the
895 * `run()` function needs to be modified:
899 * CathodeRaySimulator<dim>::CathodeRaySimulator()
902 * , dof_handler(triangulation)
903 * , particle_handler(triangulation, mapping, /*n_properties=*/dim)
904 * , next_unused_particle_id(0)
905 * , n_recently_lost_particles(0)
906 * , n_total_lost_particles(0)
907 * , n_particles_lost_through_anode(0)
910 * particle_handler.signals.particle_lost.connect(
911 * [this](const typename Particles::ParticleIterator<dim> &particle,
912 * const typename Triangulation<dim>::active_cell_iterator &cell) {
913 * this->track_lost_particle(particle, cell);
920 * void CathodeRaySimulator<dim>::make_grid()
922 * static_assert(dim == 2,
923 * "This function is currently only implemented for 2d.");
925 * const double delta = 0.5;
926 * const unsigned int nx = 5;
927 * const unsigned int ny = 3;
929 * const std::vector<Point<dim>> vertices
945 * AssertDimension(vertices.size(), nx * ny);
947 * const std::vector<unsigned int> cell_vertices[(nx - 1) * (ny - 1)] = {
948 * {0, 1, nx + 0, nx + 1},
949 * {1, 2, nx + 1, nx + 2},
950 * {2, 3, nx + 2, nx + 3},
951 * {3, 4, nx + 3, nx + 4},
953 * {5, nx + 1, 2 * nx + 0, 2 * nx + 1},
954 * {nx + 1, nx + 2, 2 * nx + 1, 2 * nx + 2},
955 * {nx + 2, nx + 3, 2 * nx + 2, 2 * nx + 3},
956 * {nx + 3, nx + 4, 2 * nx + 3, 2 * nx + 4}};
958 * std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
959 * for (unsigned int i = 0; i < cells.size(); ++i)
961 * cells[i].vertices = cell_vertices[i];
962 * cells[i].material_id = 0;
965 * GridTools::consistently_order_cells(cells);
966 * triangulation.create_triangulation(
969 * SubCellData()); // No boundary information
971 * triangulation.refine_global(2);
973 * for (auto &cell : triangulation.active_cell_iterators())
974 * for (auto &face : cell->face_iterators())
975 * if (face->at_boundary())
977 * if ((face->center()[0] > 0) && (face->center()[0] < 0.5) &&
978 * (face->center()[1] > 0) && (face->center()[1] < 2))
979 * face->set_boundary_id(BoundaryIds::cathode);
980 * else if ((face->center()[0] > 0) && (face->center()[0] < 2))
981 * face->set_boundary_id(BoundaryIds::focus_element);
982 * else if ((face->center()[0] > 4 - 1e-12) &&
983 * ((face->center()[1] > 1.5) || (face->center()[1] < 0.5)))
984 * face->set_boundary_id(BoundaryIds::anode);
986 * face->set_boundary_id(BoundaryIds::open);
989 * triangulation.refine_global(1);
994 * void CathodeRaySimulator<dim>::setup_system()
996 * dof_handler.distribute_dofs(fe);
998 * solution.reinit(dof_handler.n_dofs());
999 * system_rhs.reinit(dof_handler.n_dofs());
1001 * constraints.clear();
1002 * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1004 * VectorTools::interpolate_boundary_values(dof_handler,
1005 * BoundaryIds::cathode,
1006 * Functions::ConstantFunction<dim>(
1009 * VectorTools::interpolate_boundary_values(dof_handler,
1010 * BoundaryIds::focus_element,
1011 * Functions::ConstantFunction<dim>(
1014 * VectorTools::interpolate_boundary_values(dof_handler,
1015 * BoundaryIds::anode,
1016 * Functions::ConstantFunction<dim>(
1019 * constraints.close();
1021 * DynamicSparsityPattern dsp(dof_handler.n_dofs());
1022 * DoFTools::make_sparsity_pattern(dof_handler,
1025 * /* keep_constrained_dofs = */ false);
1026 * sparsity_pattern.copy_from(dsp);
1028 * system_matrix.reinit(sparsity_pattern);
1033 * template <int dim>
1034 * void CathodeRaySimulator<dim>::assemble_system()
1036 * system_matrix = 0;
1039 * const QGauss<dim> quadrature_formula(fe.degree + 1);
1041 * FEValues<dim> fe_values(fe,
1042 * quadrature_formula,
1043 * update_values | update_gradients |
1044 * update_quadrature_points | update_JxW_values);
1046 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1048 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1049 * Vector<double> cell_rhs(dofs_per_cell);
1051 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1053 * for (const auto &cell : dof_handler.active_cell_iterators())
1058 * fe_values.reinit(cell);
1060 * for (const unsigned int q_index : fe_values.quadrature_point_indices())
1061 * for (const unsigned int i : fe_values.dof_indices())
1063 * for (const unsigned int j : fe_values.dof_indices())
1064 * cell_matrix(i, j) +=
1065 * (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
1066 * fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
1067 * fe_values.JxW(q_index)); // dx
1070 * if (particle_handler.n_particles_in_cell(cell) > 0)
1071 * for (const auto &particle : particle_handler.particles_in_cell(cell))
1073 * const Point<dim> &reference_location =
1074 * particle.get_reference_location();
1075 * for (const unsigned int i : fe_values.dof_indices())
1077 * (fe.shape_value(i, reference_location) * // phi_i(x_p)
1078 * (-Constants::electrons_per_particle * // N
1079 * Constants::electron_charge)); // e
1082 * cell->get_dof_indices(local_dof_indices);
1083 * constraints.distribute_local_to_global(
1084 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1090 * template <int dim>
1091 * void CathodeRaySimulator<dim>::solve_field()
1093 * SolverControl solver_control(1000, 1e-12);
1094 * SolverCG<Vector<double>> solver(solver_control);
1096 * PreconditionSSOR<SparseMatrix<double>> preconditioner;
1097 * preconditioner.initialize(system_matrix, 1.2);
1099 * solver.solve(system_matrix, solution, system_rhs, preconditioner);
1101 * constraints.distribute(solution);
1106 * template <int dim>
1107 * void CathodeRaySimulator<dim>::refine_grid()
1109 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1111 * KellyErrorEstimator<dim>::estimate(dof_handler,
1112 * QGauss<dim - 1>(fe.degree + 1),
1115 * estimated_error_per_cell);
1117 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1118 * estimated_error_per_cell,
1122 * triangulation.execute_coarsening_and_refinement();
1127 * template <int dim>
1128 * void CathodeRaySimulator<dim>::create_particles()
1130 * FEFaceValues<dim> fe_face_values(fe,
1131 * QMidpoint<dim - 1>(),
1132 * update_quadrature_points |
1133 * update_gradients |
1134 * update_normal_vectors);
1136 * std::vector<Tensor<1, dim>> solution_gradients(
1137 * fe_face_values.n_quadrature_points);
1139 * for (const auto &cell : dof_handler.active_cell_iterators())
1140 * for (const auto &face : cell->face_iterators())
1141 * if (face->at_boundary() &&
1142 * (face->boundary_id() == BoundaryIds::cathode))
1144 * fe_face_values.reinit(cell, face);
1146 * const FEValuesExtractors::Scalar electric_potential(0);
1147 * fe_face_values[electric_potential].get_function_gradients(
1148 * solution, solution_gradients);
1149 * for (const unsigned int q_point :
1150 * fe_face_values.quadrature_point_indices())
1152 * const Tensor<1, dim> E = solution_gradients[q_point];
1154 * if ((E * fe_face_values.normal_vector(q_point) < 0) &&
1155 * (E.norm() > Constants::E_threshold))
1157 * const Point<dim> &location =
1158 * fe_face_values.quadrature_point(q_point);
1160 * Particles::Particle<dim> new_particle;
1161 * new_particle.set_location(location);
1162 * new_particle.set_reference_location(
1163 * mapping.transform_real_to_unit_cell(cell, location));
1164 * new_particle.set_id(next_unused_particle_id);
1165 * particle_handler.insert_particle(new_particle, cell);
1167 * ++next_unused_particle_id;
1172 * particle_handler.update_cached_numbers();
1177 * template <int dim>
1178 * void CathodeRaySimulator<dim>::move_particles()
1180 * const double dt = time.get_next_step_size();
1182 * Vector<double> solution_values(fe.n_dofs_per_cell());
1183 * FEPointEvaluation<1, dim> evaluator(mapping, fe, update_gradients);
1185 * for (const auto &cell : dof_handler.active_cell_iterators())
1186 * if (particle_handler.n_particles_in_cell(cell) > 0)
1188 * const typename Particles::ParticleHandler<
1189 * dim>::particle_iterator_range particles_in_cell =
1190 * particle_handler.particles_in_cell(cell);
1192 * std::vector<Point<dim>> particle_positions;
1193 * for (const auto &particle : particles_in_cell)
1194 * particle_positions.push_back(particle.get_reference_location());
1196 * cell->get_dof_values(solution, solution_values);
1198 * evaluator.reinit(cell, particle_positions);
1199 * evaluator.evaluate(make_array_view(solution_values),
1200 * EvaluationFlags::gradients);
1203 * typename Particles::ParticleHandler<dim>::particle_iterator
1204 * particle = particles_in_cell.begin();
1205 * for (unsigned int particle_index = 0;
1206 * particle != particles_in_cell.end();
1207 * ++particle, ++particle_index)
1209 * const Tensor<1, dim> &E =
1210 * evaluator.get_gradient(particle_index);
1212 * const Tensor<1, dim> old_velocity(particle->get_properties());
1214 * const Tensor<1, dim> acceleration =
1215 * Constants::electron_charge / Constants::electron_mass * E;
1217 * const Tensor<1, dim> new_velocity =
1218 * old_velocity + acceleration * dt;
1220 * particle->set_properties(new_velocity);
1222 * const Point<dim> new_location =
1223 * particle->get_location() + dt * new_velocity;
1224 * particle->set_location(new_location);
1229 * particle_handler.sort_particles_into_subdomains_and_cells();
1234 * template <int dim>
1235 * void CathodeRaySimulator<dim>::track_lost_particle(
1236 * const typename Particles::ParticleIterator<dim> &particle,
1237 * const typename Triangulation<dim>::active_cell_iterator &cell)
1239 * ++n_recently_lost_particles;
1240 * ++n_total_lost_particles;
1242 * const Point<dim> current_location = particle->get_location();
1243 * const Point<dim> approximate_previous_location = cell->center();
1245 * if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
1247 * const Tensor<1, dim> direction =
1248 * (current_location - approximate_previous_location) /
1249 * (current_location[0] - approximate_previous_location[0]);
1251 * const double right_boundary_intercept =
1252 * approximate_previous_location[1] +
1253 * (4 - approximate_previous_location[0]) * direction[1];
1254 * if ((right_boundary_intercept > 0.5) &&
1255 * (right_boundary_intercept < 1.5))
1256 * ++n_particles_lost_through_anode;
1262 * template <int dim>
1263 * void CathodeRaySimulator<dim>::update_timestep_size()
1265 * if (time.get_step_number() > 0)
1267 * double min_cell_size_over_velocity = std::numeric_limits<double>::max();
1269 * for (const auto &cell : dof_handler.active_cell_iterators())
1270 * if (particle_handler.n_particles_in_cell(cell) > 0)
1272 * const double cell_size = cell->minimum_vertex_distance();
1274 * double max_particle_velocity(0.0);
1276 * for (const auto &particle :
1277 * particle_handler.particles_in_cell(cell))
1279 * const Tensor<1, dim> velocity(particle.get_properties());
1280 * max_particle_velocity =
1281 * std::max(max_particle_velocity, velocity.norm());
1284 * if (max_particle_velocity > 0)
1285 * min_cell_size_over_velocity =
1286 * std::min(min_cell_size_over_velocity,
1287 * cell_size / max_particle_velocity);
1290 * constexpr double c_safety = 0.5;
1291 * time.set_desired_next_step_size(c_safety * 0.5 *
1292 * min_cell_size_over_velocity);
1296 * const QTrapezoid<dim> vertex_quadrature;
1297 * FEValues<dim> fe_values(fe, vertex_quadrature, update_gradients);
1299 * std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
1301 * double min_timestep = std::numeric_limits<double>::max();
1303 * for (const auto &cell : dof_handler.active_cell_iterators())
1304 * if (particle_handler.n_particles_in_cell(cell) > 0)
1306 * const double cell_size = cell->minimum_vertex_distance();
1308 * fe_values.reinit(cell);
1309 * fe_values.get_function_gradients(solution, field_gradients);
1312 * for (const auto q_point : fe_values.quadrature_point_indices())
1313 * max_E = std::max(max_E, field_gradients[q_point].norm());
1317 * std::min(min_timestep,
1318 * std::sqrt(0.5 * cell_size *
1319 * Constants::electron_mass /
1320 * Constants::electron_charge / max_E));
1323 * time.set_desired_next_step_size(min_timestep);
1329 * template <int dim>
1330 * class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
1333 * ElectricFieldPostprocessor()
1334 * : DataPostprocessorVector<dim>("electric_field", update_gradients)
1337 * virtual void evaluate_scalar_field(
1338 * const DataPostprocessorInputs::Scalar<dim> &input_data,
1339 * std::vector<Vector<double>> &computed_quantities) const override
1341 * AssertDimension(input_data.solution_gradients.size(),
1342 * computed_quantities.size());
1344 * for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
1346 * AssertDimension(computed_quantities[p].size(), dim);
1347 * for (unsigned int d = 0; d < dim; ++d)
1348 * computed_quantities[p][d] = input_data.solution_gradients[p][d];
1355 * template <int dim>
1356 * void CathodeRaySimulator<dim>::output_results() const
1359 * ElectricFieldPostprocessor<dim> electric_field;
1360 * DataOut<dim> data_out;
1361 * data_out.attach_dof_handler(dof_handler);
1362 * data_out.add_data_vector(solution, "electric_potential");
1363 * data_out.add_data_vector(solution, electric_field);
1364 * data_out.build_patches();
1366 * DataOutBase::VtkFlags output_flags;
1367 * output_flags.time = time.get_current_time();
1368 * output_flags.cycle = time.get_step_number();
1369 * output_flags.physical_units["electric_potential"] = "V";
1370 * output_flags.physical_units["electric_field"] = "V/m";
1372 * data_out.set_flags(output_flags);
1374 * std::ofstream output("solution-" +
1375 * Utilities::int_to_string(time.get_step_number(), 4) +
1377 * data_out.write_vtu(output);
1381 * Particles::DataOut<dim> particle_out;
1382 * particle_out.build_patches(
1384 * std::vector<std::string>(dim, "velocity"),
1385 * std::vector<DataComponentInterpretation::DataComponentInterpretation>(
1386 * dim, DataComponentInterpretation::component_is_part_of_vector));
1388 * DataOutBase::VtkFlags output_flags;
1389 * output_flags.time = time.get_current_time();
1390 * output_flags.cycle = time.get_step_number();
1391 * output_flags.physical_units["velocity"] = "m/s";
1393 * particle_out.set_flags(output_flags);
1395 * std::ofstream output("particles-" +
1396 * Utilities::int_to_string(time.get_step_number(), 4) +
1398 * particle_out.write_vtu(output);
1407 * <a name="step_83-CathodeRaySimulatorserialize"></a>
1408 * <h4>CathodeRaySimulator::serialize()</h4>
1412 * The first of the new function is the one that is called by the
1413 * BOOST Serialization framework to serialize and deserialize the
1414 * data of this class. It has already been discussed in the introduction
1415 * to this program and so does not provide any surprises. All it does is
1416 * write those member variables of the current class that cannot be
1417 * re-created easily into an archive, or read these members from it.
1418 * (Whether `operator &` facilitates a write or read operation depends on
1419 * whether the `Archive` type is an output or input archive.)
1423 * The function takes a second argument, `version`, that can be used to
1424 * create checkpoints that have version numbers. This is useful if one
1425 * evolves programs by adding more member variables, but still wants to
1426 * retain the ability to read checkpoint files created with earlier
1427 * versions of the program. The `version` variable would, in that case,
1428 * be used to represent which version of the program wrote the file,
1429 * and if necessary to read only those variables that were written with
1430 * that past version, finding a different way to initialize the new member
1431 * variables that have been added since then. We will not make use of this
1436 * Finally, while the program that indents all deal.II source files format
1437 * the following code as
1438 * <div class=CodeFragmentInTutorialComment>
1443 * as if we are taking the address of the `triangulation` variable, the
1444 * way you *should* read the code is as
1445 * <div class=CodeFragmentInTutorialComment>
1450 * where `operator &` is a binary operator that could either be interpreted
1451 * as `operator <<` for output or `operator >>` for input.
1455 * As discussed in the introduction, we do not serialize the `triangulation`
1456 * member variable, instead leaving that to separate calls in the
1457 * `checkpoint()` and `restart()` functions below.
1460 * template <int dim>
1461 * template <class Archive>
1462 * void CathodeRaySimulator<dim>::serialize(Archive &ar,
1463 * const unsigned int /* version */)
1466 * ar &particle_handler;
1467 * ar &next_unused_particle_id;
1468 * ar &n_recently_lost_particles;
1469 * ar &n_total_lost_particles;
1470 * ar &n_particles_lost_through_anode;
1479 * <a name="step_83-CathodeRaySimulatorcheckpoint"></a>
1480 * <h4>CathodeRaySimulator::checkpoint()</h4>
1484 * The checkpoint function of the principal class of this program is then
1485 * quite straightforward: We create an output file (and check that it is
1486 * writable), create an output archive, and then move the serialized
1487 * contents of the current object (i.e., the `*this` object) into the
1488 * archive. The use of `operator<<` here calls the `serialize()` function
1489 * above with an output archive as argument. When the destructor of the
1490 * `archive` variable is called at the end of the code block within which
1491 * it lives, the entire archive is written into the output file stream it
1492 * is associated with.
1496 * As mentioned in the introduction, "real" applications would not use
1497 * text-based archives as provided by the `boost::archive::text_oarchive`
1498 * class, but use binary and potentially compressed versions. This can
1499 * easily be achieved by using differently named classes, and the BOOST
1500 * documentation explains how to do that.
1503 * template <int dim>
1504 * void CathodeRaySimulator<dim>::checkpoint()
1506 * std::cout << "--- Writing checkpoint... ---" << std::endl << std::endl;
1509 * std::ofstream checkpoint_file("tmp.checkpoint_step_83");
1510 * AssertThrow(checkpoint_file,
1512 * "Could not write to the <tmp.checkpoint_step_83> file."));
1514 * boost::archive::text_oarchive archive(checkpoint_file);
1521 * The second part of the serialization is all of the data that we can
1522 * attach to cells -- see the discussion about this in the introduction.
1523 * Here, the only data we attach to cells are the particles. We then
1524 * let the triangulation save these into a set of files that all start
1525 * with the same prefix as we chose above, namely "tmp.checkpoint":
1528 * particle_handler.prepare_for_serialization();
1529 * triangulation.save("tmp.checkpoint");
1534 * At this point, the serialized data of this file has ended up in a number
1535 * of files that all start with `tmp.checkpoint` file. As mentioned in the
1536 * introduction, we do not want to directly overwrite the checkpointing
1537 * files from the previous checkpoint operation, for fear that the program
1538 * may be interrupted *while writing the checkpoint files*. This would
1539 * result in corrupted files, and defeat the whole purpose of checkpointing
1540 * because one cannot restart from such a file. On the other hand, if we got
1541 * here, we know that the "tmp.checkpoint*" files were successfully written,
1542 * and we can rename it to "checkpoint*", in the process replacing the old
1547 * We do this move operation by calling the
1548 * [C++ function that does the renaming of
1549 * files](https://en.cppreference.com/w/cpp/filesystem/rename).
1550 * Note that it is documented as being for all practical purposes
1551 * "atomic", i.e., we do not need to worry that the program may be
1552 * interrupted somewhere within the renaming operation itself. Of
1553 * course, it is possible that we get interrupted somewhere between
1554 * renaming one file and the next, and that presents problems in
1555 * itself -- in essence, we want the entire renaming operation of all of
1556 * these files to be atomic. With a couple dozen lines of extra code, one
1557 * could address this issue (using strategies that databases use frequently:
1558 * if one operation fails, we need to
1559 * [rollback](https://en.wikipedia.org/wiki/Rollback_(data_management))
1560 * the entire transaction). For the purposes of this program, this is
1561 * perhaps too much, and we will simply hope that that doesn't happen,
1562 * perhaps based on the belief that renaming files is much faster than
1563 * writing them, and that unlike writing checkpoint files, renaming does not
1564 * require much memory or disk space and so does not risk running out of
1569 * As a consequence, the following code
first loops over all files in
1570 * the current directory, picks out those that start with the
string
1571 *
"tmp.checkpoint", and puts them into a list. In a
second step,
1572 * we
loop over the list and rename each of these files to one
1573 * whose name consists of the
"tmp.checkpoint*" file but stripped off
1574 * its
first four characters (i.e., only the
"checkpoint*" part). We use
1575 *
this approach, rather than listing the files we want to rename,
1576 * because we
do not actually know the names of the files written by
1578 * file names starts.
1581 *
std::list<
std::
string> tmp_checkpoint_files;
1582 * for (const auto &dir_entry :
std::filesystem::directory_iterator("."))
1583 * if (dir_entry.is_regular_file() &&
1584 * (dir_entry.path().filename().
string().find("tmp.checkpoint") == 0))
1585 * tmp_checkpoint_files.push_back(dir_entry.path().filename().
string());
1587 * for (const
std::
string &filename : tmp_checkpoint_files)
1588 *
std::filesystem::rename(filename, filename.substr(4,
std::
string::npos));
1595 * <a name="step_83-CathodeRaySimulatorrestart"></a>
1596 * <h4>CathodeRaySimulator::restart()</h4>
1600 * The restart function of this class then simply does the opposite:
1601 * It opens an input file (and triggers an error if that file cannot be
1602 * opened), associates an input archive with it, and then reads the
1603 * contents of the current
object from it, again using the
1604 * `serialize()` function from above. Clearly, since we have written
1605 * data into a text-based archive above, we need to use the corresponding
1606 * `
boost::archive::text_iarchive` class for reading.
1611 * data, and then tell the
Particles::ParticleHandler
object to re-create
1612 * its information about all of the particles from the data just read.
1616 * The function ends by printing a status message about having restarted:
1619 * template <
int dim>
1620 *
void CathodeRaySimulator<dim>::restart()
1623 * std::ifstream checkpoint_file(
"checkpoint_step_83");
1626 *
"Could not read from the <checkpoint_step_83> file."));
1628 * boost::archive::text_iarchive archive(checkpoint_file);
1633 * particle_handler.deserialize();
1635 * std::cout <<
"--- Restarting at t=" << time.get_current_time()
1636 * <<
", dt=" << time.get_next_step_size() << std::endl
1644 * <a name=
"step_83-CathodeRaySimulatorrun"></a>
1645 * <h4>CathodeRaySimulator::run()</h4>
1649 * The last member function of the principal
class of this program is then the
1650 * driver. The driver takes a single argument to indicate
if the simulation
1651 * is a restart. If it is not a restart, the mesh is set up and the problem is
1652 * solved like in @ref step_19
"step-19". If it is a restart, then we read in everything that
1653 * is a history variable from the checkpoint file via the `restart()`
1654 * function. Recall that everything that is
inside the `
if` block at the top
1655 * of the function is exactly like in @ref step_19
"step-19", as is almost everything that
1659 *
template <
int dim>
1660 *
void CathodeRaySimulator<dim>::run(
const bool do_restart)
1662 *
if (do_restart ==
false)
1666 *
const unsigned int n_pre_refinement_cycles = 3;
1667 *
for (
unsigned int refinement_cycle = 0;
1668 * refinement_cycle < n_pre_refinement_cycles;
1669 * ++refinement_cycle)
1672 * assemble_system();
1685 * std::cout <<
"Timestep " << time.get_step_number() + 1 << std::endl;
1686 * std::cout <<
" Field degrees of freedom: "
1687 * << dof_handler.n_dofs() << std::endl;
1689 * assemble_system();
1692 * create_particles();
1693 * std::cout <<
" Total number of particles in simulation: "
1694 * << particle_handler.n_global_particles() << std::endl;
1696 * n_recently_lost_particles = 0;
1697 * update_timestep_size();
1700 * time.advance_time();
1704 * std::cout <<
" Number of particles lost this time step: "
1705 * << n_recently_lost_particles << std::endl;
1706 *
if (n_total_lost_particles > 0)
1707 * std::cout <<
" Fraction of particles lost through anode: "
1708 * << 1. * n_particles_lost_through_anode /
1709 * n_total_lost_particles
1712 * std::cout << std::endl
1713 * <<
" Now at t=" << time.get_current_time()
1714 * <<
", dt=" << time.get_previous_step_size() <<
'.'
1720 * The only other difference between
this program and @ref step_19
"step-19" is that
1721 * we checkpoint the simulation every ten time steps:
1724 *
if (time.get_step_number() % 10 == 0)
1727 *
while (time.is_at_end() ==
false);
1735 * <a name=
"step_83-Thecodemaincodefunction"></a>
1736 * <h3>The <code>main</code> function</h3>
1740 * The
final function of the program is then again the `main()` function. Its
1741 * overall structure is unchanged in all tutorial programs since @ref step_6
"step-6" and
1742 * so there is
nothing new to discuss about
this aspect.
1746 * The only difference is that we need to figure out whether a restart was
1747 * requested, or whether the program should simply start from scratch when
1748 * called. We
do this using a command line argument: The `argc` argument to
1749 * `main()` indicates how many command line arguments were provided when
1750 * the program was called (counting the name of the program as the zeroth
1751 * argument), and `argv` is an array of strings with as many elements as
1752 * `argc` that contains these command line arguments. So
if you call
1754 * <div
class=CodeFragmentInTutorialComment>
1759 * then `argc` will be 1, and `argv` will be the array with one element
1760 * and content `[
"./step-83" ]`. On the other hand,
if you call the program
1762 * <div
class=CodeFragmentInTutorialComment>
1767 * then `argc` will be 2, and `argv` will be the array with two elements
1768 * and content `[
"./step-83",
"restart" ]`. Every other choice should be
1769 * flagged as an error. The following
try block then does exactly
this:
1772 *
int main(
int argc,
char *argv[])
1776 * Step83::CathodeRaySimulator<2> cathode_ray_simulator;
1779 * cathode_ray_simulator.run(
false);
1780 *
else if ((argc == 2) && (std::string(argv[1]) ==
"restart"))
1781 * cathode_ray_simulator.run(
true);
1784 * std::cerr <<
"Error: The only allowed command line argument to this\n"
1785 *
" program is 'restart'."
1790 *
catch (std::exception &exc)
1792 * std::cerr << std::endl
1794 * <<
"----------------------------------------------------"
1796 * std::cerr <<
"Exception on processing: " << std::endl
1797 * << exc.what() << std::endl
1798 * <<
"Aborting!" << std::endl
1799 * <<
"----------------------------------------------------"
1806 * std::cerr << std::endl
1808 * <<
"----------------------------------------------------"
1810 * std::cerr <<
"Unknown exception!" << std::endl
1811 * <<
"Aborting!" << std::endl
1812 * <<
"----------------------------------------------------"
1819<a name=
"step_83-Results"></a><h1>Results</h1>
1822When you
run this program, it produces the following output that is
1823almost exactly identical to what you get from @ref step_19
"step-19":
1826 Field degrees of freedom: 4989
1827 Total number of particles in simulation: 20
1828 Number of particles lost
this time step: 0
1830 Now at t=2.12647e-07, dt=2.12647e-07.
1833 Field degrees of freedom: 4989
1834 Total number of particles in simulation: 40
1835 Number of particles lost
this time step: 0
1837 Now at t=4.14362e-07, dt=2.01715e-07.
1840 Field degrees of freedom: 4989
1841 Total number of particles in simulation: 60
1842 Number of particles lost
this time step: 0
1844 Now at t=5.23066e-07, dt=1.08704e-07.
1847 Field degrees of freedom: 4989
1848 Total number of particles in simulation: 80
1849 Number of particles lost
this time step: 0
1851 Now at t=6.08431e-07, dt=8.53649e-08.
1854 Field degrees of freedom: 4989
1855 Total number of particles in simulation: 100
1856 Number of particles lost
this time step: 0
1858 Now at t=6.81935e-07, dt=7.35039e-08.
1861 Field degrees of freedom: 4989
1862 Total number of particles in simulation: 120
1863 Number of particles lost
this time step: 0
1865 Now at t=7.47864e-07, dt=6.59294e-08.
1868 Field degrees of freedom: 4989
1869 Total number of particles in simulation: 140
1870 Number of particles lost
this time step: 0
1872 Now at t=8.2516e-07, dt=7.72957e-08.
1875 Field degrees of freedom: 4989
1876 Total number of particles in simulation: 158
1877 Number of particles lost
this time step: 0
1879 Now at t=8.95325e-07, dt=7.01652e-08.
1882 Field degrees of freedom: 4989
1883 Total number of particles in simulation: 172
1884 Number of particles lost
this time step: 0
1886 Now at t=9.67852e-07, dt=7.25269e-08.
1889 Field degrees of freedom: 4989
1890 Total number of particles in simulation: 186
1891 Number of particles lost
this time step: 0
1893 Now at t=1.03349e-06, dt=6.56398e-08.
1895--- Writing checkpoint... ---
1898 Field degrees of freedom: 4989
1899 Total number of particles in simulation: 198
1900 Number of particles lost
this time step: 0
1902 Now at t=1.11482e-06, dt=8.13268e-08.
1905 Field degrees of freedom: 4989
1906 Total number of particles in simulation: 206
1907 Number of particles lost
this time step: 0
1909 Now at t=1.18882e-06, dt=7.39967e-08.
1912 Field degrees of freedom: 4989
1913 Total number of particles in simulation: 212
1914 Number of particles lost
this time step: 0
1916 Now at t=1.26049e-06, dt=7.16705e-08.
1920The only thing that is different is the additional line after the tenth
1921output (which also appears after the 20th, 30th, etc., time step) indicating
1922that a checkpoint file has been written.
1924Because we chose to use a text-based format
for the checkpoint file (which
1925you should not
do in production codes because it uses quite a lot of disk space),
1926we can actually inspect
this file. It will look like
this, with many many more
192922 serialization::archive 18 0 0 0 0 0 7 0 0 3 1 0
19304989 -1.00000000000000000e+00 -1.00000000000000000e+00 -1.00000000000000000e+00 -9.96108134982226390e-01 -1.00000000000000000e+00 -9.98361082867431748e-01
1933What each of these
numbers represents is hard to tell in practice, and also
1934entirely unimportant
for our current purposes -- it
's
1935a representation of the many objects that make up this program's state, and
1936from which one can restore its state. The
point simply being that
this is what
1937serialization produces: A
long list (a sequence) of bits that we can put
1938into a file, and that we can later read again to recreate the state of the
1941Now here
's the fun part. Let's say you hit Control-
C on the command
1942line at the
point above (say, during time step 13 or 14). There
's
1943a set of checkpoint files on disk that saved the state after ten time steps.
1944Based on the logic in `main()`, we should be able to restart from
1945this point if we run the program with
1949Indeed, this is what then happens:
1951--- Restarting at t=1.03349e-06, dt=6.56398e-08
1954 Field degrees of freedom: 4989
1955 Total number of particles in simulation: 198
1956 Number of particles lost this time step: 0
1958 Now at t=1.11482e-06, dt=8.13268e-08.
1961 Field degrees of freedom: 4989
1962 Total number of particles in simulation: 206
1963 Number of particles lost this time step: 0
1965 Now at t=1.18882e-06, dt=7.39967e-08.
1969After the status message that shows that we are restarting, this is
1970indeed *the exact same output* for the following time steps we had
1971gotten previously -- in other words, saving the complete state
1972seems to have worked, and the program has continued as if it had
1973never been interrupted!
1977<a name="extensions"></a>
1978<a name="step_83-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1981<a name="step_83-Makingefficiencyapriority"></a><h4> Making efficiency a priority </h4>
1984While the techniques we have shown here are fully sufficient for what
1985we want to do, namely checkpoint and restart computations, and are in
1986fact also fully sufficient for much larger code bases such as
1987[ASPECT](https://aspect.geodynamics.org), one could go beyond what is
1988still a relatively simple scheme.
1990Specifically, among the things we need to recognize is that writing
1991large amounts of data to disk is expensive and can take a good long
1992time to finish -- for example, for large parallel computations with,
1993say, a billion unknowns, checkpoints can run into the hundred gigabyte
1994range or beyond. One may ask whether that could be avoided, or at
1995least whether we can mitigate the cost.
1997One way to do that is to first serialize the state of the program into
1998a buffer in memory (like the `Archive` objects the `serialize()`
1999functions write to and read from), and once that is done, start a *separate*
2000thread do the writing while the rest of the program continues with computations.
2001This is useful because writing the data to disk often takes a
2002long time but not a lot of CPU power: It just takes time to move the
2003data through the network to the file server, and from there onto the
2004actual disks. This is something that might as well happen while we are
2005doing something useful again (namely, solving more time steps). Should
2006the machine crash during this phase, nothing is lost: As discussed in
2007the introduction, we are writing the checkpoint into a temporary file
2008(which will be lost in the case of a machine failure), but we have
2009kept the previous checkpoint around until we know that the temporary
2010file is complete and can be moved over the old one.
2012The only thing we have to pay attention in this background-writing
2013scheme is that we cannot start with creating a new checkpoint while
2014the previous one is still being written in the background.
2016Doing this all is not technically very difficult; the code just
2017requires more explanation than lines of code, and so we omit doing
2018this in the program here. But you can take a look at the
2019`MainLoop::output()` function of @ref step_69 "step-69" to see how such a code looks
2022A variation of this general approach is that each process writes its
2023data immediately, but into files that are held on fast file systems --
2024say, a node-local SSD rather than a file server shared by the entire
2025cluster. One would then just tell the operating system to move this
2026file to the centeral file server in a second step, and this step can
2027happen in the background at whatever speed the operating system can
2028provide. Or perhaps one leaves *most* of these files on the fast local
2029file system in hopes that the restart happens before these files are
2030purged (say, by a script that runs nightly) and only moves these files
2031to the permanent file system every tenth time we create a checkpoint.
2033In all of these cases, the logic quickly becomes quite complicated. As
2034usual, the solution is not to re-invent the wheel: Libraries such as
2035[VeloC](https://www.anl.gov/mcs/veloc-very-low-overhead-transparent-multilevel-checkpointrestart),
2036developed by the Exascale Computing Project (ECP) already do all of
2037this and more, for codes that are orders of magnitude more complex
2038than the little example here.
2040Separately, one might want to try to reduce the amount of time it
2041takes to serialize objects into a buffer in memory. As mentioned
2044library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html)
2045for this task, but it is not the only player in town. One could, for
2046example, use the [bitsery](https://github.com/fraillt/bitsery), [cereal](https://github.com/USCiLab/cereal), or
2047[zpp](https://github.com/eyalz800/serializer) projects instead, which can be substantially faster than BOOST.
2050<a name="step_83-PlainProg"></a>
2051<h1> The plain program</h1>
2052@include "step-83.cc"
void prepare_for_serialization()
void serialize(Archive &ar, const unsigned int version)
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
void save(Archive &ar, const unsigned int version) const
#define AssertThrow(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ general
No special properties.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation