deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-83.h
Go to the documentation of this file.
1)
225{
226 ar & values;
227}
228@endcode
229
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
238general idea here.)
239
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:
247@code
248template <int dim>
249template <class Archive>
250inline void
251Quadrature<dim>::serialize(Archive &ar, const unsigned int)
252{
253 // Forward to the (de)serialization function in the base class:
254 ar & static_cast<Subscriptor &>(*this);
255
256 // Then (de)serialize the member variables:
257 ar & quadrature_points & weights;
258}
259@endcode
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
265purposes!
266
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.)
274
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:
278@code
279 template <int dim>
280 class CathodeRaySimulator
281 {
282 public:
283 CathodeRaySimulator();
284
285 void run();
286
287 private:
288 void make_grid();
289 void setup_system();
290 void assemble_system();
291 void solve_field();
292 void refine_grid();
293
294 void create_particles();
295 void move_particles();
296 void track_lost_particle(const typename Particles::ParticleIterator<dim> & particle,
297 const typename Triangulation<dim>::active_cell_iterator &cell);
298
299
300 void update_timestep_size();
301 void output_results() const;
302
304 MappingQGeneric<dim> mapping;
305 FE_Q<dim> fe;
306 DoFHandler<dim> dof_handler;
307 AffineConstraints<double> constraints;
308
309 SparseMatrix<double> system_matrix;
310 SparsityPattern sparsity_pattern;
311
312 Vector<double> solution;
313 Vector<double> system_rhs;
314
315 Particles::ParticleHandler<dim> particle_handler;
316 types::particle_index next_unused_particle_id;
317 types::particle_index n_recently_lost_particles;
318 types::particle_index n_total_lost_particles;
319 types::particle_index n_particles_lost_through_anode;
320
321 DiscreteTime time;
322 };
323@endcode
324One would then first write a
325function that allows (de)serialization of the data of this class
326as follows:
327@code
328 template <int dim>
329 template <class Archive>
330 void CathodeRaySimulator<dim>::serialize(Archive &ar,
331 const unsigned int /* version */)
332 {
333 ar & triangulation;
334 ar & solution;
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;
340 ar & time;
341 }
342@endcode
343As discussed below, this is not the entire list of member variables: We
344only want to serialize those that we cannot otherwise re-generate.
345
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):
350@code
351 template <int dim>
352 void CathodeRaySimulator<dim>::checkpoint()
353 {
354 std::ofstream checkpoint_file("checkpoint");
355 boost::archive::text_oarchive archive(checkpoint_file,
356 boost::archive::no_header);
357
358 archive << *this;
359 }
360@endcode
361
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.
369
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.
380
381In any case, the data we have thus created is read in very similarly
382in the following function (again with some minor modifications):
383@code
384 template <int dim>
385 void CathodeRaySimulator<dim>::restart()
386 {
387 std::ifstream checkpoint_file("checkpoint");
388
389 boost::archive::text_iarchive archive(checkpoint_file,
390 boost::archive::no_header);
391
392 archive >> *this;
393 }
394@endcode
395
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.
400
401
402<a name="step_83-Serializationofparallelprograms"></a><h4> Serialization of parallel programs </h4>
403
404
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.
414
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.
446
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
452two categories:
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.
468
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
474call Triangulation::save() to actually write this information into a set of
475files that become part of the checkpoint. Upon restart, we then first call
476Triangulation::load(), followed by Particles::ParticleHandler::deserialize()
477to retrieve the particles from the cells they are attached to.
478
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.)
487
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
494@code
495 ar & triangulation;
496@endcode
497We do still need to say
498@code
499 ar & particle_handler;
500@endcode
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.
505
506
507<a name="step_83-Checkpointingstrategies"></a><h3>Checkpointing strategies</h3>
508
509
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
513checkpoints?
514
515
516<a name="step_83-Whattosaverestore"></a><h4>What to save/restore</h4>
517
518
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.
524
525Recall the main class of @ref step_19 "step-19", which
526had quite a lot of member variables one might want to
527(de)serialize:
528@code
529 template <int dim>
530 class CathodeRaySimulator
531 {
532 public:
533 CathodeRaySimulator();
534
535 void run();
536
537 private:
538 [... member functions ...]
539
541 MappingQGeneric<dim> mapping;
542 FE_Q<dim> fe;
543 DoFHandler<dim> dof_handler;
544 AffineConstraints<double> constraints;
545
546 SparseMatrix<double> system_matrix;
547 SparsityPattern sparsity_pattern;
548
549 Vector<double> solution;
550 Vector<double> system_rhs;
551
552 Particles::ParticleHandler<dim> particle_handler;
553 types::particle_index next_unused_particle_id;
554 types::particle_index n_recently_lost_particles;
555 types::particle_index n_total_lost_particles;
556 types::particle_index n_particles_lost_through_anode;
557
558 DiscreteTime time;
559 };
560@endcode
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,
563re-loaded.
564
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:
572
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.
584
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
600 active is minimal.
601
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.
617
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)?
623
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.
638
639
640<a name="step_83-Howpreciselyshouldwesavethedataofacheckpoint"></a><h4>How precisely should we save the data of a checkpoint</h4>
641
642
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
649checkpoint file.
650
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.
658
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
662variables, and several into which the triangulation puts its
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.
670
671
672<a name="step_83-Howoftentosaverestore"></a><h4>How often to save/restore</h4>
673
674
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.
683
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.
697
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.
702 *
703 *
704 * <a name="step_83-CommProg"></a>
705 * <h1> The commented program</h1>
706 *
707 *
708 * <a name="step_83-Includefiles"></a>
709 * <h3>Include files</h3>
710 *
711
712 *
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
715 * all the same:
716 *
717 * @code
718 *   #include <deal.II/base/quadrature_lib.h>
719 *   #include <deal.II/base/discrete_time.h>
720 *  
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>
728 *  
729 *   #include <deal.II/grid/tria.h>
730 *   #include <deal.II/grid/grid_refinement.h>
731 *   #include <deal.II/grid/grid_tools.h>
732 *  
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>
737 *  
738 *   #include <deal.II/dofs/dof_handler.h>
739 *   #include <deal.II/dofs/dof_tools.h>
740 *  
741 *   #include <deal.II/numerics/data_out.h>
742 *   #include <deal.II/numerics/vector_tools.h>
743 *   #include <deal.II/numerics/error_estimator.h>
744 *  
745 *   #include <deal.II/particles/particle_handler.h>
746 *   #include <deal.II/particles/data_out.h>
747 *  
748 * @endcode
749 *
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:
753 *
754 * @code
755 *   #include <boost/archive/text_iarchive.hpp>
756 *   #include <boost/archive/text_oarchive.hpp>
757 *  
758 *   #include <filesystem>
759 *   #include <fstream>
760 *   #include <string>
761 *  
762 * @endcode
763 *
764 *
765 * <a name="step_83-Globaldefinitions"></a>
766 * <h3>Global definitions</h3>
767 *
768
769 *
770 * As is customary, we put everything that corresponds to the details of the
771 * program into a namespace of its own.
772 *
773 * @code
774 *   namespace Step83
775 *   {
776 *   using namespace dealii;
777 *  
778 *   namespace BoundaryIds
779 *   {
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
785 *  
786 *   namespace Constants
787 *   {
788 *   constexpr double electron_mass = 9.1093837015e-31;
789 *   constexpr double electron_charge = 1.602176634e-19;
790 *  
791 *   constexpr double V0 = 1;
792 *  
793 *   constexpr double E_threshold = 0.05;
794 *  
795 *   constexpr double electrons_per_particle = 3e15;
796 *   } // namespace Constants
797 *  
798 *  
799 * @endcode
800 *
801 *
802 * <a name="step_83-Themainclass"></a>
803 * <h3>The main class</h3>
804 *
805
806 *
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`.
816 *
817
818 *
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.
822 *
823 * @code
824 *   template <int dim>
825 *   class CathodeRaySimulator
826 *   {
827 *   public:
828 *   CathodeRaySimulator();
829 *  
830 *   void run(const bool do_restart);
831 *  
832 *   template <class Archive>
833 *   void serialize(Archive &ar, const unsigned int version);
834 *  
835 *   private:
836 *   void make_grid();
837 *   void setup_system();
838 *   void assemble_system();
839 *   void solve_field();
840 *   void refine_grid();
841 *  
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);
847 *  
848 *   void update_timestep_size();
849 *   void output_results() const;
850 *  
851 *   void checkpoint();
852 *   void restart();
853 *  
854 *   Triangulation<dim> triangulation;
855 *   const MappingQ<dim> mapping;
856 *   const FE_Q<dim> fe;
857 *   DoFHandler<dim> dof_handler;
858 *   AffineConstraints<double> constraints;
859 *  
860 *   SparseMatrix<double> system_matrix;
861 *   SparsityPattern sparsity_pattern;
862 *  
863 *   Vector<double> solution;
864 *   Vector<double> system_rhs;
865 *  
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;
871 *  
872 *   DiscreteTime time;
873 *   };
874 *  
875 *  
876 *  
877 * @endcode
878 *
879 *
880 * <a name="step_83-ThecodeCathodeRaySimulatorcodeclassimplementation"></a>
881 * <h3>The <code>CathodeRaySimulator</code> class implementation</h3>
882 *
883
884 *
885 *
886 * <a name="step_83-Theunchangedpartsoftheclass"></a>
887 * <h4>The unchanged parts of the class</h4>
888 *
889
890 *
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:
896 *
897 * @code
898 *   template <int dim>
899 *   CathodeRaySimulator<dim>::CathodeRaySimulator()
900 *   : mapping(1)
901 *   , fe(2)
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)
908 *   , time(0, 1e-4)
909 *   {
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);
914 *   });
915 *   }
916 *  
917 *  
918 *  
919 *   template <int dim>
920 *   void CathodeRaySimulator<dim>::make_grid()
921 *   {
922 *   static_assert(dim == 2,
923 *   "This function is currently only implemented for 2d.");
924 *  
925 *   const double delta = 0.5;
926 *   const unsigned int nx = 5;
927 *   const unsigned int ny = 3;
928 *  
929 *   const std::vector<Point<dim>> vertices
930 *   = {{0, 0},
931 *   {1, 0},
932 *   {2, 0},
933 *   {3, 0},
934 *   {4, 0},
935 *   {delta, 1},
936 *   {1, 1},
937 *   {2, 1},
938 *   {3, 1},
939 *   {4, 1},
940 *   {0, 2},
941 *   {1, 2},
942 *   {2, 2},
943 *   {3, 2},
944 *   {4, 2}};
945 *   AssertDimension(vertices.size(), nx * ny);
946 *  
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},
952 *  
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}};
957 *  
958 *   std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
959 *   for (unsigned int i = 0; i < cells.size(); ++i)
960 *   {
961 *   cells[i].vertices = cell_vertices[i];
962 *   cells[i].material_id = 0;
963 *   }
964 *  
965 *   GridTools::consistently_order_cells(cells);
966 *   triangulation.create_triangulation(
967 *   vertices,
968 *   cells,
969 *   SubCellData()); // No boundary information
970 *  
971 *   triangulation.refine_global(2);
972 *  
973 *   for (auto &cell : triangulation.active_cell_iterators())
974 *   for (auto &face : cell->face_iterators())
975 *   if (face->at_boundary())
976 *   {
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);
985 *   else
986 *   face->set_boundary_id(BoundaryIds::open);
987 *   }
988 *  
989 *   triangulation.refine_global(1);
990 *   }
991 *  
992 *  
993 *   template <int dim>
994 *   void CathodeRaySimulator<dim>::setup_system()
995 *   {
996 *   dof_handler.distribute_dofs(fe);
997 *  
998 *   solution.reinit(dof_handler.n_dofs());
999 *   system_rhs.reinit(dof_handler.n_dofs());
1000 *  
1001 *   constraints.clear();
1002 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1003 *  
1004 *   VectorTools::interpolate_boundary_values(dof_handler,
1005 *   BoundaryIds::cathode,
1006 *   Functions::ConstantFunction<dim>(
1007 *   -Constants::V0),
1008 *   constraints);
1009 *   VectorTools::interpolate_boundary_values(dof_handler,
1010 *   BoundaryIds::focus_element,
1011 *   Functions::ConstantFunction<dim>(
1012 *   -Constants::V0),
1013 *   constraints);
1014 *   VectorTools::interpolate_boundary_values(dof_handler,
1015 *   BoundaryIds::anode,
1016 *   Functions::ConstantFunction<dim>(
1017 *   +Constants::V0),
1018 *   constraints);
1019 *   constraints.close();
1020 *  
1021 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
1022 *   DoFTools::make_sparsity_pattern(dof_handler,
1023 *   dsp,
1024 *   constraints,
1025 *   /* keep_constrained_dofs = */ false);
1026 *   sparsity_pattern.copy_from(dsp);
1027 *  
1028 *   system_matrix.reinit(sparsity_pattern);
1029 *   }
1030 *  
1031 *  
1032 *  
1033 *   template <int dim>
1034 *   void CathodeRaySimulator<dim>::assemble_system()
1035 *   {
1036 *   system_matrix = 0;
1037 *   system_rhs = 0;
1038 *  
1039 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
1040 *  
1041 *   FEValues<dim> fe_values(fe,
1042 *   quadrature_formula,
1043 *   update_values | update_gradients |
1044 *   update_quadrature_points | update_JxW_values);
1045 *  
1046 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1047 *  
1048 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1049 *   Vector<double> cell_rhs(dofs_per_cell);
1050 *  
1051 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1052 *  
1053 *   for (const auto &cell : dof_handler.active_cell_iterators())
1054 *   {
1055 *   cell_matrix = 0;
1056 *   cell_rhs = 0;
1057 *  
1058 *   fe_values.reinit(cell);
1059 *  
1060 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
1061 *   for (const unsigned int i : fe_values.dof_indices())
1062 *   {
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
1068 *   }
1069 *  
1070 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1071 *   for (const auto &particle : particle_handler.particles_in_cell(cell))
1072 *   {
1073 *   const Point<dim> &reference_location =
1074 *   particle.get_reference_location();
1075 *   for (const unsigned int i : fe_values.dof_indices())
1076 *   cell_rhs(i) +=
1077 *   (fe.shape_value(i, reference_location) * // phi_i(x_p)
1078 *   (-Constants::electrons_per_particle * // N
1079 *   Constants::electron_charge)); // e
1080 *   }
1081 *  
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);
1085 *   }
1086 *   }
1087 *  
1088 *  
1089 *  
1090 *   template <int dim>
1091 *   void CathodeRaySimulator<dim>::solve_field()
1092 *   {
1093 *   SolverControl solver_control(1000, 1e-12);
1094 *   SolverCG<Vector<double>> solver(solver_control);
1095 *  
1096 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1097 *   preconditioner.initialize(system_matrix, 1.2);
1098 *  
1099 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
1100 *  
1101 *   constraints.distribute(solution);
1102 *   }
1103 *  
1104 *  
1105 *  
1106 *   template <int dim>
1107 *   void CathodeRaySimulator<dim>::refine_grid()
1108 *   {
1109 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1110 *  
1111 *   KellyErrorEstimator<dim>::estimate(dof_handler,
1112 *   QGauss<dim - 1>(fe.degree + 1),
1113 *   {},
1114 *   solution,
1115 *   estimated_error_per_cell);
1116 *  
1117 *   GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1118 *   estimated_error_per_cell,
1119 *   0.1,
1120 *   0.03);
1121 *  
1122 *   triangulation.execute_coarsening_and_refinement();
1123 *   }
1124 *  
1125 *  
1126 *  
1127 *   template <int dim>
1128 *   void CathodeRaySimulator<dim>::create_particles()
1129 *   {
1130 *   FEFaceValues<dim> fe_face_values(fe,
1131 *   QMidpoint<dim - 1>(),
1132 *   update_quadrature_points |
1133 *   update_gradients |
1134 *   update_normal_vectors);
1135 *  
1136 *   std::vector<Tensor<1, dim>> solution_gradients(
1137 *   fe_face_values.n_quadrature_points);
1138 *  
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))
1143 *   {
1144 *   fe_face_values.reinit(cell, face);
1145 *  
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())
1151 *   {
1152 *   const Tensor<1, dim> E = solution_gradients[q_point];
1153 *  
1154 *   if ((E * fe_face_values.normal_vector(q_point) < 0) &&
1155 *   (E.norm() > Constants::E_threshold))
1156 *   {
1157 *   const Point<dim> &location =
1158 *   fe_face_values.quadrature_point(q_point);
1159 *  
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);
1166 *  
1167 *   ++next_unused_particle_id;
1168 *   }
1169 *   }
1170 *   }
1171 *  
1172 *   particle_handler.update_cached_numbers();
1173 *   }
1174 *  
1175 *  
1176 *  
1177 *   template <int dim>
1178 *   void CathodeRaySimulator<dim>::move_particles()
1179 *   {
1180 *   const double dt = time.get_next_step_size();
1181 *  
1182 *   Vector<double> solution_values(fe.n_dofs_per_cell());
1183 *   FEPointEvaluation<1, dim> evaluator(mapping, fe, update_gradients);
1184 *  
1185 *   for (const auto &cell : dof_handler.active_cell_iterators())
1186 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1187 *   {
1188 *   const typename Particles::ParticleHandler<
1189 *   dim>::particle_iterator_range particles_in_cell =
1190 *   particle_handler.particles_in_cell(cell);
1191 *  
1192 *   std::vector<Point<dim>> particle_positions;
1193 *   for (const auto &particle : particles_in_cell)
1194 *   particle_positions.push_back(particle.get_reference_location());
1195 *  
1196 *   cell->get_dof_values(solution, solution_values);
1197 *  
1198 *   evaluator.reinit(cell, particle_positions);
1199 *   evaluator.evaluate(make_array_view(solution_values),
1200 *   EvaluationFlags::gradients);
1201 *  
1202 *   {
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)
1208 *   {
1209 *   const Tensor<1, dim> &E =
1210 *   evaluator.get_gradient(particle_index);
1211 *  
1212 *   const Tensor<1, dim> old_velocity(particle->get_properties());
1213 *  
1214 *   const Tensor<1, dim> acceleration =
1215 *   Constants::electron_charge / Constants::electron_mass * E;
1216 *  
1217 *   const Tensor<1, dim> new_velocity =
1218 *   old_velocity + acceleration * dt;
1219 *  
1220 *   particle->set_properties(new_velocity);
1221 *  
1222 *   const Point<dim> new_location =
1223 *   particle->get_location() + dt * new_velocity;
1224 *   particle->set_location(new_location);
1225 *   }
1226 *   }
1227 *   }
1228 *  
1229 *   particle_handler.sort_particles_into_subdomains_and_cells();
1230 *   }
1231 *  
1232 *  
1233 *  
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)
1238 *   {
1239 *   ++n_recently_lost_particles;
1240 *   ++n_total_lost_particles;
1241 *  
1242 *   const Point<dim> current_location = particle->get_location();
1243 *   const Point<dim> approximate_previous_location = cell->center();
1244 *  
1245 *   if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
1246 *   {
1247 *   const Tensor<1, dim> direction =
1248 *   (current_location - approximate_previous_location) /
1249 *   (current_location[0] - approximate_previous_location[0]);
1250 *  
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;
1257 *   }
1258 *   }
1259 *  
1260 *  
1261 *  
1262 *   template <int dim>
1263 *   void CathodeRaySimulator<dim>::update_timestep_size()
1264 *   {
1265 *   if (time.get_step_number() > 0)
1266 *   {
1267 *   double min_cell_size_over_velocity = std::numeric_limits<double>::max();
1268 *  
1269 *   for (const auto &cell : dof_handler.active_cell_iterators())
1270 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1271 *   {
1272 *   const double cell_size = cell->minimum_vertex_distance();
1273 *  
1274 *   double max_particle_velocity(0.0);
1275 *  
1276 *   for (const auto &particle :
1277 *   particle_handler.particles_in_cell(cell))
1278 *   {
1279 *   const Tensor<1, dim> velocity(particle.get_properties());
1280 *   max_particle_velocity =
1281 *   std::max(max_particle_velocity, velocity.norm());
1282 *   }
1283 *  
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);
1288 *   }
1289 *  
1290 *   constexpr double c_safety = 0.5;
1291 *   time.set_desired_next_step_size(c_safety * 0.5 *
1292 *   min_cell_size_over_velocity);
1293 *   }
1294 *   else
1295 *   {
1296 *   const QTrapezoid<dim> vertex_quadrature;
1297 *   FEValues<dim> fe_values(fe, vertex_quadrature, update_gradients);
1298 *  
1299 *   std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
1300 *  
1301 *   double min_timestep = std::numeric_limits<double>::max();
1302 *  
1303 *   for (const auto &cell : dof_handler.active_cell_iterators())
1304 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1305 *   {
1306 *   const double cell_size = cell->minimum_vertex_distance();
1307 *  
1308 *   fe_values.reinit(cell);
1309 *   fe_values.get_function_gradients(solution, field_gradients);
1310 *  
1311 *   double max_E = 0;
1312 *   for (const auto q_point : fe_values.quadrature_point_indices())
1313 *   max_E = std::max(max_E, field_gradients[q_point].norm());
1314 *  
1315 *   if (max_E > 0)
1316 *   min_timestep =
1317 *   std::min(min_timestep,
1318 *   std::sqrt(0.5 * cell_size *
1319 *   Constants::electron_mass /
1320 *   Constants::electron_charge / max_E));
1321 *   }
1322 *  
1323 *   time.set_desired_next_step_size(min_timestep);
1324 *   }
1325 *   }
1326 *  
1327 *  
1328 *  
1329 *   template <int dim>
1330 *   class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
1331 *   {
1332 *   public:
1333 *   ElectricFieldPostprocessor()
1334 *   : DataPostprocessorVector<dim>("electric_field", update_gradients)
1335 *   {}
1336 *  
1337 *   virtual void evaluate_scalar_field(
1338 *   const DataPostprocessorInputs::Scalar<dim> &input_data,
1339 *   std::vector<Vector<double>> &computed_quantities) const override
1340 *   {
1341 *   AssertDimension(input_data.solution_gradients.size(),
1342 *   computed_quantities.size());
1343 *  
1344 *   for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
1345 *   {
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];
1349 *   }
1350 *   }
1351 *   };
1352 *  
1353 *  
1354 *  
1355 *   template <int dim>
1356 *   void CathodeRaySimulator<dim>::output_results() const
1357 *   {
1358 *   {
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();
1365 *  
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";
1371 *  
1372 *   data_out.set_flags(output_flags);
1373 *  
1374 *   std::ofstream output("solution-" +
1375 *   Utilities::int_to_string(time.get_step_number(), 4) +
1376 *   ".vtu");
1377 *   data_out.write_vtu(output);
1378 *   }
1379 *  
1380 *   {
1381 *   Particles::DataOut<dim> particle_out;
1382 *   particle_out.build_patches(
1383 *   particle_handler,
1384 *   std::vector<std::string>(dim, "velocity"),
1385 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>(
1386 *   dim, DataComponentInterpretation::component_is_part_of_vector));
1387 *  
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";
1392 *  
1393 *   particle_out.set_flags(output_flags);
1394 *  
1395 *   std::ofstream output("particles-" +
1396 *   Utilities::int_to_string(time.get_step_number(), 4) +
1397 *   ".vtu");
1398 *   particle_out.write_vtu(output);
1399 *   }
1400 *   }
1401 *  
1402 *  
1403 *  
1404 * @endcode
1405 *
1406 *
1407 * <a name="step_83-CathodeRaySimulatorserialize"></a>
1408 * <h4>CathodeRaySimulator::serialize()</h4>
1409 *
1410
1411 *
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.)
1420 *
1421
1422 *
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
1432 * ability here.
1433 *
1434
1435 *
1436 * Finally, while the program that indents all deal.II source files format
1437 * the following code as
1438 * <div class=CodeFragmentInTutorialComment>
1439 * @code
1440 * ar &solution;
1441 * @endcode
1442 * </div>
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>
1446 * @code
1447 * ar & solution;
1448 * @endcode
1449 * </div>
1450 * where `operator &` is a binary operator that could either be interpreted
1451 * as `operator <<` for output or `operator >>` for input.
1452 *
1453
1454 *
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.
1458 *
1459 * @code
1460 *   template <int dim>
1461 *   template <class Archive>
1462 *   void CathodeRaySimulator<dim>::serialize(Archive &ar,
1463 *   const unsigned int /* version */)
1464 *   {
1465 *   ar &solution;
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;
1471 *   ar &time;
1472 *   }
1473 *  
1474 *  
1475 *  
1476 * @endcode
1477 *
1478 *
1479 * <a name="step_83-CathodeRaySimulatorcheckpoint"></a>
1480 * <h4>CathodeRaySimulator::checkpoint()</h4>
1481 *
1482
1483 *
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.
1493 *
1494
1495 *
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.
1501 *
1502 * @code
1503 *   template <int dim>
1504 *   void CathodeRaySimulator<dim>::checkpoint()
1505 *   {
1506 *   std::cout << "--- Writing checkpoint... ---" << std::endl << std::endl;
1507 *  
1508 *   {
1509 *   std::ofstream checkpoint_file("tmp.checkpoint_step_83");
1510 *   AssertThrow(checkpoint_file,
1511 *   ExcMessage(
1512 *   "Could not write to the <tmp.checkpoint_step_83> file."));
1513 *  
1514 *   boost::archive::text_oarchive archive(checkpoint_file);
1515 *  
1516 *   archive << *this;
1517 *   }
1518 *  
1519 * @endcode
1520 *
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":
1526 *
1527 * @code
1528 *   particle_handler.prepare_for_serialization();
1529 *   triangulation.save("tmp.checkpoint");
1530 *  
1531 *  
1532 * @endcode
1533 *
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
1543 * file.
1544 *
1545
1546 *
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
1565 * either.
1566 *
1567
1568 *
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
1577 * the Triangulation::save() function, though we know how each of these
1578 * file names starts.
1579 *
1580 * @code
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());
1586 *  
1587 *   for (const std::string &filename : tmp_checkpoint_files)
1588 *   std::filesystem::rename(filename, filename.substr(4, std::string::npos));
1589 *   }
1590 *  
1591 *  
1592 * @endcode
1593 *
1594 *
1595 * <a name="step_83-CathodeRaySimulatorrestart"></a>
1596 * <h4>CathodeRaySimulator::restart()</h4>
1597 *
1598
1599 *
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.
1607 *
1608
1609 *
1610 * In a second step, we ask the triangulation to read in cell-attached
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.
1613 *
1614
1615 *
1616 * The function ends by printing a status message about having restarted:
1617 *
1618 * @code
1619 *   template <int dim>
1620 *   void CathodeRaySimulator<dim>::restart()
1621 *   {
1622 *   {
1623 *   std::ifstream checkpoint_file("checkpoint_step_83");
1624 *   AssertThrow(checkpoint_file,
1625 *   ExcMessage(
1626 *   "Could not read from the <checkpoint_step_83> file."));
1627 *  
1628 *   boost::archive::text_iarchive archive(checkpoint_file);
1629 *   archive >> *this;
1630 *   }
1631 *  
1632 *   triangulation.load("checkpoint");
1633 *   particle_handler.deserialize();
1634 *  
1635 *   std::cout << "--- Restarting at t=" << time.get_current_time()
1636 *   << ", dt=" << time.get_next_step_size() << std::endl
1637 *   << std::endl;
1638 *   }
1639 *  
1640 *  
1641 * @endcode
1642 *
1643 *
1644 * <a name="step_83-CathodeRaySimulatorrun"></a>
1645 * <h4>CathodeRaySimulator::run()</h4>
1646 *
1647
1648 *
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
1656 * follows:
1657 *
1658 * @code
1659 *   template <int dim>
1660 *   void CathodeRaySimulator<dim>::run(const bool do_restart)
1661 *   {
1662 *   if (do_restart == false)
1663 *   {
1664 *   make_grid();
1665 *  
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)
1670 *   {
1671 *   setup_system();
1672 *   assemble_system();
1673 *   solve_field();
1674 *   refine_grid();
1675 *   }
1676 *   }
1677 *   else
1678 *   {
1679 *   restart();
1680 *   }
1681 *  
1682 *   setup_system();
1683 *   do
1684 *   {
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;
1688 *  
1689 *   assemble_system();
1690 *   solve_field();
1691 *  
1692 *   create_particles();
1693 *   std::cout << " Total number of particles in simulation: "
1694 *   << particle_handler.n_global_particles() << std::endl;
1695 *  
1696 *   n_recently_lost_particles = 0;
1697 *   update_timestep_size();
1698 *   move_particles();
1699 *  
1700 *   time.advance_time();
1701 *  
1702 *   output_results();
1703 *  
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
1710 *   << std::endl;
1711 *  
1712 *   std::cout << std::endl
1713 *   << " Now at t=" << time.get_current_time()
1714 *   << ", dt=" << time.get_previous_step_size() << '.'
1715 *   << std::endl
1716 *   << std::endl;
1717 *  
1718 * @endcode
1719 *
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:
1722 *
1723 * @code
1724 *   if (time.get_step_number() % 10 == 0)
1725 *   checkpoint();
1726 *   }
1727 *   while (time.is_at_end() == false);
1728 *   }
1729 *   } // namespace Step83
1730 *  
1731 *  
1732 * @endcode
1733 *
1734 *
1735 * <a name="step_83-Thecodemaincodefunction"></a>
1736 * <h3>The <code>main</code> function</h3>
1737 *
1738
1739 *
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.
1743 *
1744
1745 *
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
1753 * the program as
1754 * <div class=CodeFragmentInTutorialComment>
1755 * @code
1756 * ./step-83
1757 * @endcode
1758 * </div>
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
1761 * as
1762 * <div class=CodeFragmentInTutorialComment>
1763 * @code
1764 * ./step-83 restart
1765 * @endcode
1766 * </div>
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:
1770 *
1771 * @code
1772 *   int main(int argc, char *argv[])
1773 *   {
1774 *   try
1775 *   {
1776 *   Step83::CathodeRaySimulator<2> cathode_ray_simulator;
1777 *  
1778 *   if (argc == 1)
1779 *   cathode_ray_simulator.run(false); // no restart
1780 *   else if ((argc == 2) && (std::string(argv[1]) == "restart"))
1781 *   cathode_ray_simulator.run(true); // yes restart
1782 *   else
1783 *   {
1784 *   std::cerr << "Error: The only allowed command line argument to this\n"
1785 *   " program is 'restart'."
1786 *   << std::endl;
1787 *   return 1;
1788 *   }
1789 *   }
1790 *   catch (std::exception &exc)
1791 *   {
1792 *   std::cerr << std::endl
1793 *   << std::endl
1794 *   << "----------------------------------------------------"
1795 *   << std::endl;
1796 *   std::cerr << "Exception on processing: " << std::endl
1797 *   << exc.what() << std::endl
1798 *   << "Aborting!" << std::endl
1799 *   << "----------------------------------------------------"
1800 *   << std::endl;
1801 *  
1802 *   return 1;
1803 *   }
1804 *   catch (...)
1805 *   {
1806 *   std::cerr << std::endl
1807 *   << std::endl
1808 *   << "----------------------------------------------------"
1809 *   << std::endl;
1810 *   std::cerr << "Unknown exception!" << std::endl
1811 *   << "Aborting!" << std::endl
1812 *   << "----------------------------------------------------"
1813 *   << std::endl;
1814 *   return 1;
1815 *   }
1816 *   return 0;
1817 *   }
1818 * @endcode
1819<a name="step_83-Results"></a><h1>Results</h1>
1820
1821
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":
1824@code
1825Timestep 1
1826 Field degrees of freedom: 4989
1827 Total number of particles in simulation: 20
1828 Number of particles lost this time step: 0
1829
1830 Now at t=2.12647e-07, dt=2.12647e-07.
1831
1832Timestep 2
1833 Field degrees of freedom: 4989
1834 Total number of particles in simulation: 40
1835 Number of particles lost this time step: 0
1836
1837 Now at t=4.14362e-07, dt=2.01715e-07.
1838
1839Timestep 3
1840 Field degrees of freedom: 4989
1841 Total number of particles in simulation: 60
1842 Number of particles lost this time step: 0
1843
1844 Now at t=5.23066e-07, dt=1.08704e-07.
1845
1846Timestep 4
1847 Field degrees of freedom: 4989
1848 Total number of particles in simulation: 80
1849 Number of particles lost this time step: 0
1850
1851 Now at t=6.08431e-07, dt=8.53649e-08.
1852
1853Timestep 5
1854 Field degrees of freedom: 4989
1855 Total number of particles in simulation: 100
1856 Number of particles lost this time step: 0
1857
1858 Now at t=6.81935e-07, dt=7.35039e-08.
1859
1860Timestep 6
1861 Field degrees of freedom: 4989
1862 Total number of particles in simulation: 120
1863 Number of particles lost this time step: 0
1864
1865 Now at t=7.47864e-07, dt=6.59294e-08.
1866
1867Timestep 7
1868 Field degrees of freedom: 4989
1869 Total number of particles in simulation: 140
1870 Number of particles lost this time step: 0
1871
1872 Now at t=8.2516e-07, dt=7.72957e-08.
1873
1874Timestep 8
1875 Field degrees of freedom: 4989
1876 Total number of particles in simulation: 158
1877 Number of particles lost this time step: 0
1878
1879 Now at t=8.95325e-07, dt=7.01652e-08.
1880
1881Timestep 9
1882 Field degrees of freedom: 4989
1883 Total number of particles in simulation: 172
1884 Number of particles lost this time step: 0
1885
1886 Now at t=9.67852e-07, dt=7.25269e-08.
1887
1888Timestep 10
1889 Field degrees of freedom: 4989
1890 Total number of particles in simulation: 186
1891 Number of particles lost this time step: 0
1892
1893 Now at t=1.03349e-06, dt=6.56398e-08.
1894
1895--- Writing checkpoint... ---
1896
1897Timestep 11
1898 Field degrees of freedom: 4989
1899 Total number of particles in simulation: 198
1900 Number of particles lost this time step: 0
1901
1902 Now at t=1.11482e-06, dt=8.13268e-08.
1903
1904Timestep 12
1905 Field degrees of freedom: 4989
1906 Total number of particles in simulation: 206
1907 Number of particles lost this time step: 0
1908
1909 Now at t=1.18882e-06, dt=7.39967e-08.
1910
1911Timestep 13
1912 Field degrees of freedom: 4989
1913 Total number of particles in simulation: 212
1914 Number of particles lost this time step: 0
1915
1916 Now at t=1.26049e-06, dt=7.16705e-08.
1917
1918[...]
1919@endcode
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.
1923
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
1927lines:
1928@code
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
1931[...]
1932@endcode
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
1939program.
1940
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
1946```
1947 ./step-83 restart
1948```
1949Indeed, this is what then happens:
1950@code
1951--- Restarting at t=1.03349e-06, dt=6.56398e-08
1952
1953Timestep 11
1954 Field degrees of freedom: 4989
1955 Total number of particles in simulation: 198
1956 Number of particles lost this time step: 0
1957
1958 Now at t=1.11482e-06, dt=8.13268e-08.
1959
1960Timestep 12
1961 Field degrees of freedom: 4989
1962 Total number of particles in simulation: 206
1963 Number of particles lost this time step: 0
1964
1965 Now at t=1.18882e-06, dt=7.39967e-08.
1966
1967[...]
1968@endcode
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!
1974
1975
1976
1977<a name="extensions"></a>
1978<a name="step_83-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1979
1980
1981<a name="step_83-Makingefficiencyapriority"></a><h4> Making efficiency a priority </h4>
1982
1983
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.
1989
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.
1996
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.
2011
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.
2015
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
2020like.
2021
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.
2032
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.
2039
2040Separately, one might want to try to reduce the amount of time it
2041takes to serialize objects into a buffer in memory. As mentioned
2042above, we use the
2043[BOOST serialization
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.
2048 *
2049 *
2050<a name="step_83-PlainProg"></a>
2051<h1> The plain program</h1>
2052@include "step-83.cc"
2053*/
Definition fe_q.h:554
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
Definition tensor.h:854
void save(Archive &ar, const unsigned int version) const
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#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())
Definition loop.h:564
@ general
No special properties.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
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)
STL namespace.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation