Reference documentation for deal.II version 9.6.0
\(\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{
227 ar & values;
228}
229@endcode
230
231In other words, the Tensor class stores an array of `dim` objects of a
232tensor of one rank lower, and in order to serialize this array, the
233function simply uses the overloaded `operator&` which recognizes that
234`values` is a C-style array of fixed length, and what the data type of
235the elements of this array are. It will then in turn call the
236`serialize()` member function of `Tensor<rank-1,dim>`. (The recursion
237has a stopping point which is not of importance for this discussion. The
238actual implementation of class Tensor looks different, but you get the
239general idea here.)
240
241Depending on whether the `Archive` type used for the first argument of
242this function is an input or output archive, `operator&` reads from or
243writes into the `values` member variable. If the class had more member
244variables, one would list them one after the other. The same is true
245if a class has base classes. For example, here is the corresponding
246function for the Quadrature class that stores quadrature points and
247weights, and is derived from the Subscriptor base class:
248@code
249template <int dim>
250template <class Archive>
251inline void
252Quadrature<dim>::serialize(Archive &ar, const unsigned int)
253{
254 // Forward to the (de)serialization function in the base class:
255 ar & static_cast<Subscriptor &>(*this);
256
257 // Then (de)serialize the member variables:
258 ar & quadrature_points & weights;
259}
260@endcode
261The beauty of having to write only one function for both serialization
262and deserialization is that obviously one has to read data from the
263archive in the same order as one writes into it, and that is easy
264to get wrong if one had to write two separate functions -- but that
265is automatically correct if one has the same function for both
266purposes!
267
268The BOOST serialization library also has mechanisms if one wants to
269write separate save and load functions. This is useful if a class
270stores a lot of data in one array and then has a cache of commonly
271accessed data on the side for fast access; one would then only store
272and load the large array, and in the load function re-build the cache
273from the large array. (See also the discussion below on what one
274actually wants to save.)
275
276Based on these principles, let us consider how one would (de)serialize
277a program such as @ref step_19 "step-19". Recall that its principal class looked
278essentially as follows:
279@code
280 template <int dim>
281 class CathodeRaySimulator
282 {
283 public:
284 CathodeRaySimulator();
285
286 void run();
287
288 private:
289 void make_grid();
290 void setup_system();
291 void assemble_system();
292 void solve_field();
293 void refine_grid();
294
295 void create_particles();
296 void move_particles();
297 void track_lost_particle(const typename Particles::ParticleIterator<dim> & particle,
298 const typename Triangulation<dim>::active_cell_iterator &cell);
299
300
301 void update_timestep_size();
302 void output_results() const;
303
305 MappingQGeneric<dim> mapping;
306 FE_Q<dim> fe;
307 DoFHandler<dim> dof_handler;
308 AffineConstraints<double> constraints;
309
310 SparseMatrix<double> system_matrix;
311 SparsityPattern sparsity_pattern;
312
313 Vector<double> solution;
314 Vector<double> system_rhs;
315
316 Particles::ParticleHandler<dim> particle_handler;
317 types::particle_index next_unused_particle_id;
318 types::particle_index n_recently_lost_particles;
319 types::particle_index n_total_lost_particles;
320 types::particle_index n_particles_lost_through_anode;
321
322 DiscreteTime time;
323 };
324@endcode
325One would then first write a
326function that allows (de)serialization of the data of this class
327as follows:
328@code
329 template <int dim>
330 template <class Archive>
331 void CathodeRaySimulator<dim>::serialize(Archive &ar,
332 const unsigned int /* version */)
333 {
334 ar & triangulation;
335 ar & solution;
336 ar & particle_handler;
337 ar & next_unused_particle_id;
338 ar & n_recently_lost_particles;
339 ar & n_total_lost_particles;
340 ar & n_particles_lost_through_anode;
341 ar & time;
342 }
343@endcode
344As discussed below, this is not the entire list of member variables: We
345only want to serialize those that we cannot otherwise re-generate.
346
347Next, we need a function that writes a checkpoint by creating an "output
348archive" and using the usual `operator<<` to put an object into this archive. In
349the code of this program, this will then look as follows (with some
350minor modifications we will discuss later):
351@code
352 template <int dim>
353 void CathodeRaySimulator<dim>::checkpoint()
354 {
355 std::ofstream checkpoint_file("checkpoint");
356 boost::archive::text_oarchive archive(checkpoint_file,
357 boost::archive::no_header);
358
359 archive << *this;
360 }
361@endcode
362
363What `operator<<` does here is to call the serialization functions of the right hand
364operand (here, the `serialize()` function described above, with an output
365archive as template argument), and create a serialized representation of the data. In the
366end, serialization has put everything into the "archive" from which one
367can extract a (sometimes very long) string that one can save in bulk
368into a file, and that is exactly what happens when the destructor of the
369`archive` variable is called.
370
371BOOST serialization offers different archives, including ones that
372store the data in text format (as we do above), in binary format, or
373already compressed with something like gzip to minimize the amount of
374space necessary. The specific type of archive to be used is selected
375by the type of the `archive` variable above (and the corresponding
376variable in the `restart()` function of course). This program uses a
377text archive so that you can look at how a serialized representation
378would actually look at, though a "real" program would of course try to
379be more space efficient by using binary (and possibly compressed)
380representations of the data.
381
382In any case, the data we have thus created is read in very similarly
383in the following function (again with some minor modifications):
384@code
385 template <int dim>
386 void CathodeRaySimulator<dim>::restart()
387 {
388 std::ifstream checkpoint_file("checkpoint");
389
390 boost::archive::text_iarchive archive(checkpoint_file,
391 boost::archive::no_header);
392
393 archive >> *this;
394 }
395@endcode
396
397The magic of this approach is that one doesn't actually have to write very
398much code to checkpoint or restart programs: It is all hidden in the
399serialization functions of classes such as Triangulation,
400Particles::ParticleHandler, etc., which the deal.II library provides for you.
401
402
403<a name="step_83-Serializationofparallelprograms"></a><h4> Serialization of parallel programs </h4>
404
405
406The program as provided here is sequential, i.e., it runs on a single
407processor just like the original @ref step_19 "step-19" did. But what if you had a parallel
408program -- say, something like @ref step_40 "step-40" -- that runs in parallel with MPI?
409In that case, serialization and checkpoint/restart becomes more complicated.
410While parallel execution is not something that is of concern to us
411in this program, it is an issue that has influenced the design of how
412deal.II does serialization; as a consequence, we need to talk through
413what makes serialization of parallel programs more difficult in order
414to understand why this program does things the way it does.
415
416Intuitively, one might simply want to use the same idea as we used
417here except that we let every MPI process serialize its own data, and
418read its own data. This works, but there are some drawbacks:
419- There is a certain subset of data that is replicated across all
420 MPI processes and that would be then written by all processes. An example
421 is the `time` data structure that stores the current time, time step
422 size, and other information, and that better be the same on all processes.
423 Typically, the replicated data isn't a large fraction of a program's
424 overall memory usage, and perhaps writing it more than once isn't
425 going to be too much of a problem, but it is unsatisfactory anyway
426 to have this kind of data on disk more than once.
427- One would have to think in detail how exactly one wants to represent
428 the data on disk. One possibility would be for every MPI process to write
429 its own file. On the other hand, checkpointing is most useful for large
430 programs, using large numbers of processes -- it is not uncommon to use
431 checkpointing on programs that run on 10,000 or more processors in parallel.
432 This would then lead to 10,000 or more files on disk. That's unpleasant, and
433 probably inefficient as well. We could address this by first letting all
434 the processes serialize into a string in memory (using `std::ostringstream`)
435 and then collating all of these strings into one file. The MPI I/O sub-system
436 has facilities to make that happen, for example, but it will require a bit
437 of thought not the least because the serialized data from each process
438 will likely result in strings of different sizes.
439- Perhaps the most important reason to rethink how one does things in parallel
440 is because, with a little bit of thought, it is possible to checkpoint a
441 program running with @f$N@f$ MPI processes and restart it with @f$M\neq N@f$
442 processes. This may, at first, seem like a pointless exercise, but it is
443 useful if one had, for example, a program that repeatedly refines the mesh
444 and where it is inefficient to run the early refinement steps with a
445 coarse mesh on too many processes, whereas it is too slow to run the later
446 refinement steps with a fine mesh on too few processes.
447
448In order to address these issues, in particular the last one, the right approach
449is to deviate a bit from the simple scheme of having a `serialize()` function
450that simply serializes/deserializes *everything* into an archive, and then have two
451functions `checkpoint()` and `restart()` that for all practical purposes defer
452all the work to the `serialize()` function. Instead, one splits all data into
453two categories:
454- Data that is tied to the cells of a triangulation. This includes the mesh itself,
455 but also the particles in the Particles::ParticleHandler class, and most
456 importantly the solution vector(s). The way to serialize such data is to
457 attach the data to cells and then let the Triangulation class serialize the attached
458 data along with its own data. If this is done in a way so that we can re-load
459 a triangulation on a different number of processes than the data was written,
460 then this automatically also ensures that we can restore solution vectors
461 and Particles::ParticleHandler objects (and everything else we can attach
462 to the cells of a triangulation) on a different number of processes.
463- Other data. In finite element programs, this data is almost always replicated
464 across processes, and so it is enough if the "root" process (typically the
465 process with MPI rank zero) writes it to disk. Upon restart, the root
466 process reads the data from disk, sends it to all other processes (however many
467 of them there may be), and these then initialize their own copies of the
468 replicated data structures.
469
470These sorts of considerations have influenced the design of the Triangulation and
471Particles::ParticleHandler classes. In particular, Particles::ParticleHandler's
472`serialize()` function only serializes the "other data" category, but not the
473actual particles; these can instead be attached to the triangulation by calling
475call Triangulation::save() to actually write this information into a set of
476files that become part of the checkpoint. Upon restart, we then first call
477Triangulation::load(), followed by Particles::ParticleHandler::deserialize()
478to retrieve the particles from the cells they are attached to.
479
480(We could, with relatively minimal effort, use the same scheme for the solution
481vector: The SolutionTransfer class can be used to attach the values of degrees
482of freedom to cells, and then Triangulation::save() also writes these into
483checkpoint files. SolutionTransfer would then be able to re-create the solution
484vector upon restart in a similar way. However, in contrast to
485Particles::ParticleHandler, the Vector class we use for the solution vector
486can actually serialize itself completely, and so we will go with this
487approach and save ourselves the dozen or so additional lines of code.)
488
489Finally, even though we wrote the `serialize()` function above in such
490a way that it also serializes the `triangulation` member variable, in practice
491the call to Triangulation::save() we needed to deal with the particles *also*
492saves the same kind of information, and Triangulation::load() reads it.
493In other words, we are saving redundant information; in the actual
494implementation of the program, we therefore skip the call to
495@code
496 ar & triangulation;
497@endcode
498We do still need to say
499@code
500 ar & particle_handler;
501@endcode
502because the information attached to the cells of the triangulation only contains
503information about the particles themselves, whereas the previous line is
504necessary to store information such as how many particles there are, what the
505next unused particle index is, and other internal information about the class.
506
507
508<a name="step_83-Checkpointingstrategies"></a><h3>Checkpointing strategies</h3>
509
510
511Having discussed the general idea of checkpoint/restart, let us turn
512to some more specific questions one has to answer: (i) What do we
513actually want to save/restore? (ii) How often do we want to write
514checkpoints?
515
516
517<a name="step_83-Whattosaverestore"></a><h4>What to save/restore</h4>
518
519
520We will base this tutorial on @ref step_19 "step-19", and so let us use it as an
521example in this section. Recall that that program simulates an
522electric field in which particles move from the electrode on one
523side to the other side of the domain, i.e., we have both field-based
524and particle-based information to store.
525
526Recall the main class of @ref step_19 "step-19", which
527had quite a lot of member variables one might want to
528(de)serialize:
529@code
530 template <int dim>
531 class CathodeRaySimulator
532 {
533 public:
534 CathodeRaySimulator();
535
536 void run();
537
538 private:
539 [... member functions ...]
540
542 MappingQGeneric<dim> mapping;
543 FE_Q<dim> fe;
544 DoFHandler<dim> dof_handler;
545 AffineConstraints<double> constraints;
546
547 SparseMatrix<double> system_matrix;
548 SparsityPattern sparsity_pattern;
549
550 Vector<double> solution;
551 Vector<double> system_rhs;
552
553 Particles::ParticleHandler<dim> particle_handler;
554 types::particle_index next_unused_particle_id;
555 types::particle_index n_recently_lost_particles;
556 types::particle_index n_total_lost_particles;
557 types::particle_index n_particles_lost_through_anode;
558
559 DiscreteTime time;
560 };
561@endcode
562Do we really need to save all of these to disk? That would presumably
563lead to quite a lot of data that needs to be stored and, if necessary,
564re-loaded.
565
566In practice, one does not save all of this information, but only what
567cannot be reasonably re-computed in different ways. What is saved
568should also depend on also *where* in the program's algorithm one
569currently is, and one generally finds a convenient point at which not
570so much data needs to be stored. For the
571current example of @ref step_19 "step-19", a time dependent problem, one could apply
572the following considerations:
573
574- The program runs with the same finite element every time, so there
575 is no need to actually save the element: We know what polynomial
576 degree we want, and can just re-generate the element upon
577 restart. If the polynomial degree was a run-time parameter, then
578 maybe we should serialize the polynomial degree rather than all of
579 the myriad data structures that characterize a FE_Q object, given
580 that we can always re-generate the object by just knowing its
581 polynomial degree. This is the classical trade-off of space vs time:
582 We can achieve the same outcome by saving far less data if we are
583 willing to offer a bit of CPU time to regenerate all of the internal
584 data structures of the FE_Q given the polynomial degree.
585
586- We rebuild the matrix and sparsity pattern in each time step from
587 the DoFHandler and the finite element. These are quite large data
588 structures, but they are conceptually easy to re-create again as
589 necessary. So they do not need to be saved to disk, and this is
590 going to save quite a lot of space. Furthermore, we really only need
591 the matrix for the linear solve; once we are done with the linear
592 solve in the `solve_field()` function, the contents of the matrix
593 are no longer used and are, indeed, overwritten in the next time
594 step. As a consequence, there would really only be a point in saving
595 the matrix if we did the checkpointing between the assembly and the
596 linear solve -- but maybe that is just not a convenient point for
597 this operation, and we should pick a better location. In practice,
598 one generally puts the checkpointing at either the very end or the
599 very beginning of the time stepping loop, given that this is the
600 point where the number of variables whose values are currently
601 active is minimal.
602
603- We also do not need to save the DoFHandler object: If we know the
604 triangulation, we can always just create a DoFHandler object during
605 restart to enumerate degrees of freedom in the same way as we did
606 the last time before a previous program run was checkpointed. In
607 fact, the example implementation of the `checkpoint()` function
608 shown above did not serialize the DoFHandler object for this very
609 reason. On the other hand, we probably do want to save the
610 Triangulation here given that the triangulation is not statically
611 generated once at the beginning of the program and then never
612 changed, but is dynamically adapted every few time steps. In order
613 to re-generate the triangulation, we would therefore have to save
614 which cells were refined/coarsened and when (i.e., the *history* of
615 the triangulation), and this would likely cost substantially more
616 disk space for long-running computations than just saving the
617 triangulation itself.
618
619Similar considerations can be applied to all member variables: Can we
620re-generate their values with relatively little effort (in which case
621they do not have to be saved) or is their state difficult or
622impossible to re-generate if it is not saved to disk (in which case
623the variable needs to be serialized)?
624
625@note If you have carefully read @ref step_19 "step-19", you might now realize that
626 strictly speaking, we do not *need* to checkpoint to solution vector.
627 This is because the solution vector represents the electric field,
628 which the program solves for at the beginning of each timestep and
629 that this solve does not make reference to the electric field at
630 previous time steps -- in other words, the electric field is not
631 a "history variable": If we know the domain, the mesh, the finite
632 element, and the positions of the particles, we can recompute the
633 solution vector, and consequently we would not have to save it
634 into the checkpoint file. However, this is perhaps more work than
635 we want to do for checkpointing (which you will see is otherwise
636 rather little code) and so, for pedagological purposes, we simply
637 save the solution vector along with the other variables that
638 actually do represent the history of the program.
639
640
641<a name="step_83-Howpreciselyshouldwesavethedataofacheckpoint"></a><h4>How precisely should we save the data of a checkpoint</h4>
642
643
644Recall that the goal of checkpointing is to end up with a safe copy of
645where the program currently is in its computations. As a consequence,
646we need to make sure that we do not end up in a situation where, for
647example, we start overwriting the previous checkpoint file and
648somewhere halfway through the serialization process, the machine
649crashes and we end up with an aborted program and no functional
650checkpoint file.
651
652Instead, the procedure one generally follows to guard against this
653kind of scenario is that checkpoints are written into a file that is
654*separate* from the previous checkpoint file; only once we are past
655the writing process and the file is safely on disk can we replace the
656previous checkpoint file by the one just written -- that is, we *move*
657the new file into place of the old one. You will see in the code how
658this two-step process is implemented in practice.
659
660The situation is made slightly more complicated by the fact that
661in the program below, a "checkpoint" actually consists of a number
662of files -- one file into which we write the program's member
663variables, and several into which the triangulation puts its
664information. We then would have to rename several files,
665preferrably as a single, "atomic" operation that cannot be
666interrupted. Implementing this is tricky and non-trivial (though
667possible), and so we will not show this part and instead just
668assume that nothing will happen between renaming the first
669and the last of the files -- maybe not a great strategy in
670general, but good enough for this tutorial program.
671
672
673<a name="step_83-Howoftentosaverestore"></a><h4>How often to save/restore</h4>
674
675
676Now that we know *what* we want to save and how we want to restore it,
677we need to answer the question *how often* we want to checkpoint the
678program. At least theoretically, this question has been answered many
679decades ago already, see @cite Young1974 and @cite Daly2006. In
680practice (as actually also in these theoretical derivations), it comes
681down to (i) how long it takes to checkpoint data, and (ii) how
682frequently we expect that the stored data will have to be used, i.e.,
683how often the system crashes.
684
685For example, if it takes five minutes to save the state of the
686program, then we probably do not want to write a checkpoint every ten
687minutes. On the other hand, if it only takes five seconds, then maybe
688ten minutes is a reasonable frequency if we run on a modest 100 cores
689and the machine does not crash very often, given that in that case the
690overhead is only approximately 1%. Finally, if it takes five seconds
691to save the state, but we are running on 100,000 processes (i.e., a
692very expensive simulation) and the machine frequently crashes, then
693maybe we are willing to offer a 5% penalty in the overall run time and
694write a checkpoint every minute and a half given that we lose far less
695work this way on average if the machine crashes at an unpredictable
696moment in our computations. The papers cited above essentially just
697formalize this sort of consideration.
698
699In the program, we will not dwell on this and simply choose an ad-hoc
700value of saving the state every ten time steps: That's too often in
701practice, but is useful for experiencing how this works in
702practice without having to run the program too long.
703 *
704 *
705 * <a name="step_83-CommProg"></a>
706 * <h1> The commented program</h1>
707 *
708 *
709 * <a name="step_83-Includefiles"></a>
710 * <h3>Include files</h3>
711 *
712
713 *
714 * This program, with the exception of the checkpointing component
715 * is identical to @ref step_19 "step-19", and so the following include files are
716 * all the same:
717 *
718 * @code
719 *   #include <deal.II/base/quadrature_lib.h>
720 *   #include <deal.II/base/discrete_time.h>
721 *  
722 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
723 *   #include <deal.II/lac/full_matrix.h>
724 *   #include <deal.II/lac/precondition.h>
725 *   #include <deal.II/lac/solver_cg.h>
726 *   #include <deal.II/lac/sparse_matrix.h>
727 *   #include <deal.II/lac/vector.h>
728 *   #include <deal.II/lac/affine_constraints.h>
729 *  
730 *   #include <deal.II/grid/tria.h>
731 *   #include <deal.II/grid/grid_refinement.h>
732 *   #include <deal.II/grid/grid_tools.h>
733 *  
734 *   #include <deal.II/fe/mapping_q.h>
735 *   #include <deal.II/matrix_free/fe_point_evaluation.h>
736 *   #include <deal.II/fe/fe_q.h>
737 *   #include <deal.II/fe/fe_values.h>
738 *  
739 *   #include <deal.II/dofs/dof_handler.h>
740 *   #include <deal.II/dofs/dof_tools.h>
741 *  
742 *   #include <deal.II/numerics/data_out.h>
743 *   #include <deal.II/numerics/vector_tools.h>
744 *   #include <deal.II/numerics/error_estimator.h>
745 *  
746 *   #include <deal.II/particles/particle_handler.h>
747 *   #include <deal.II/particles/data_out.h>
748 *  
749 * @endcode
750 *
751 * The only thing new are the following two include files. They are the ones
752 * that declare the classes we use as archives for reading (`iarchive` = input
753 * archive) and writing (`oarchive` = output archive) serialized data:
754 *
755 * @code
756 *   #include <boost/archive/text_iarchive.hpp>
757 *   #include <boost/archive/text_oarchive.hpp>
758 *  
759 *   #include <filesystem>
760 *   #include <fstream>
761 *   #include <string>
762 *  
763 * @endcode
764 *
765 *
766 * <a name="step_83-Globaldefinitions"></a>
767 * <h3>Global definitions</h3>
768 *
769
770 *
771 * As is customary, we put everything that corresponds to the details of the
772 * program into a namespace of its own.
773 *
774 * @code
775 *   namespace Step83
776 *   {
777 *   using namespace dealii;
778 *  
779 *   namespace BoundaryIds
780 *   {
781 *   constexpr types::boundary_id open = 101;
782 *   constexpr types::boundary_id cathode = 102;
783 *   constexpr types::boundary_id focus_element = 103;
784 *   constexpr types::boundary_id anode = 104;
785 *   } // namespace BoundaryIds
786 *  
787 *   namespace Constants
788 *   {
789 *   constexpr double electron_mass = 9.1093837015e-31;
790 *   constexpr double electron_charge = 1.602176634e-19;
791 *  
792 *   constexpr double V0 = 1;
793 *  
794 *   constexpr double E_threshold = 0.05;
795 *  
796 *   constexpr double electrons_per_particle = 3e15;
797 *   } // namespace Constants
798 *  
799 *  
800 * @endcode
801 *
802 *
803 * <a name="step_83-Themainclass"></a>
804 * <h3>The main class</h3>
805 *
806
807 *
808 * The following is then the main class of this program. It is,
809 * fundamentally, identical to @ref step_19 "step-19" with the exception of
810 * the `checkpoint()` and `restart()` functions, along with the
811 * `serialize()` function we use to serialize and deserialize the
812 * data this class stores. The `serialize()` function is called
813 * by the BOOST serialization framework, and consequently has to
814 * have exactly the set of arguments used here. Furthermore, because
815 * it is called by BOOST functions, it has to be `public`; the other
816 * two new functions are as always made `private`.
817 *
818
819 *
820 * The `run()` function has also been modified to enable simulation restart
821 * via its new argument `do_restart` that indicates whether or not to
822 * start the simulation from a checkpoint.
823 *
824 * @code
825 *   template <int dim>
826 *   class CathodeRaySimulator
827 *   {
828 *   public:
829 *   CathodeRaySimulator();
830 *  
831 *   void run(const bool do_restart);
832 *  
833 *   template <class Archive>
834 *   void serialize(Archive &ar, const unsigned int version);
835 *  
836 *   private:
837 *   void make_grid();
838 *   void setup_system();
839 *   void assemble_system();
840 *   void solve_field();
841 *   void refine_grid();
842 *  
843 *   void create_particles();
844 *   void move_particles();
845 *   void track_lost_particle(
846 *   const Particles::ParticleIterator<dim> &particle,
847 *   const typename Triangulation<dim>::active_cell_iterator &cell);
848 *  
849 *   void update_timestep_size();
850 *   void output_results() const;
851 *  
852 *   void checkpoint();
853 *   void restart();
854 *  
855 *   Triangulation<dim> triangulation;
856 *   const MappingQ<dim> mapping;
857 *   const FE_Q<dim> fe;
858 *   DoFHandler<dim> dof_handler;
859 *   AffineConstraints<double> constraints;
860 *  
861 *   SparseMatrix<double> system_matrix;
862 *   SparsityPattern sparsity_pattern;
863 *  
864 *   Vector<double> solution;
865 *   Vector<double> system_rhs;
866 *  
867 *   Particles::ParticleHandler<dim> particle_handler;
868 *   types::particle_index next_unused_particle_id;
869 *   types::particle_index n_recently_lost_particles;
870 *   types::particle_index n_total_lost_particles;
871 *   types::particle_index n_particles_lost_through_anode;
872 *  
873 *   DiscreteTime time;
874 *   };
875 *  
876 *  
877 *  
878 * @endcode
879 *
880 *
881 * <a name="step_83-ThecodeCathodeRaySimulatorcodeclassimplementation"></a>
882 * <h3>The <code>CathodeRaySimulator</code> class implementation</h3>
883 *
884
885 *
886 *
887 * <a name="step_83-Theunchangedpartsoftheclass"></a>
888 * <h4>The unchanged parts of the class</h4>
889 *
890
891 *
892 * Let us start with those parts of the class that are all unchanged
893 * from @ref step_19 "step-19" and about which you can learn there. We will
894 * then pick up with commentary again when we get to the two new
895 * functions, `checkpoint()` and `restart()`, along with how the
896 * `run()` function needs to be modified:
897 *
898 * @code
899 *   template <int dim>
900 *   CathodeRaySimulator<dim>::CathodeRaySimulator()
901 *   : mapping(1)
902 *   , fe(2)
903 *   , dof_handler(triangulation)
904 *   , particle_handler(triangulation, mapping, /*n_properties=*/dim)
905 *   , next_unused_particle_id(0)
906 *   , n_recently_lost_particles(0)
907 *   , n_total_lost_particles(0)
908 *   , n_particles_lost_through_anode(0)
909 *   , time(0, 1e-4)
910 *   {
911 *   particle_handler.signals.particle_lost.connect(
912 *   [this](const typename Particles::ParticleIterator<dim> &particle,
913 *   const typename Triangulation<dim>::active_cell_iterator &cell) {
914 *   this->track_lost_particle(particle, cell);
915 *   });
916 *   }
917 *  
918 *  
919 *  
920 *   template <int dim>
921 *   void CathodeRaySimulator<dim>::make_grid()
922 *   {
923 *   static_assert(dim == 2,
924 *   "This function is currently only implemented for 2d.");
925 *  
926 *   const double delta = 0.5;
927 *   const unsigned int nx = 5;
928 *   const unsigned int ny = 3;
929 *  
930 *   const std::vector<Point<dim>> vertices
931 *   = {{0, 0},
932 *   {1, 0},
933 *   {2, 0},
934 *   {3, 0},
935 *   {4, 0},
936 *   {delta, 1},
937 *   {1, 1},
938 *   {2, 1},
939 *   {3, 1},
940 *   {4, 1},
941 *   {0, 2},
942 *   {1, 2},
943 *   {2, 2},
944 *   {3, 2},
945 *   {4, 2}};
946 *   AssertDimension(vertices.size(), nx * ny);
947 *  
948 *   const std::vector<unsigned int> cell_vertices[(nx - 1) * (ny - 1)] = {
949 *   {0, 1, nx + 0, nx + 1},
950 *   {1, 2, nx + 1, nx + 2},
951 *   {2, 3, nx + 2, nx + 3},
952 *   {3, 4, nx + 3, nx + 4},
953 *  
954 *   {5, nx + 1, 2 * nx + 0, 2 * nx + 1},
955 *   {nx + 1, nx + 2, 2 * nx + 1, 2 * nx + 2},
956 *   {nx + 2, nx + 3, 2 * nx + 2, 2 * nx + 3},
957 *   {nx + 3, nx + 4, 2 * nx + 3, 2 * nx + 4}};
958 *  
959 *   std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
960 *   for (unsigned int i = 0; i < cells.size(); ++i)
961 *   {
962 *   cells[i].vertices = cell_vertices[i];
963 *   cells[i].material_id = 0;
964 *   }
965 *  
966 *   GridTools::consistently_order_cells(cells);
967 *   triangulation.create_triangulation(
968 *   vertices,
969 *   cells,
970 *   SubCellData()); // No boundary information
971 *  
972 *   triangulation.refine_global(2);
973 *  
974 *   for (auto &cell : triangulation.active_cell_iterators())
975 *   for (auto &face : cell->face_iterators())
976 *   if (face->at_boundary())
977 *   {
978 *   if ((face->center()[0] > 0) && (face->center()[0] < 0.5) &&
979 *   (face->center()[1] > 0) && (face->center()[1] < 2))
980 *   face->set_boundary_id(BoundaryIds::cathode);
981 *   else if ((face->center()[0] > 0) && (face->center()[0] < 2))
982 *   face->set_boundary_id(BoundaryIds::focus_element);
983 *   else if ((face->center()[0] > 4 - 1e-12) &&
984 *   ((face->center()[1] > 1.5) || (face->center()[1] < 0.5)))
985 *   face->set_boundary_id(BoundaryIds::anode);
986 *   else
987 *   face->set_boundary_id(BoundaryIds::open);
988 *   }
989 *  
990 *   triangulation.refine_global(1);
991 *   }
992 *  
993 *  
994 *   template <int dim>
995 *   void CathodeRaySimulator<dim>::setup_system()
996 *   {
997 *   dof_handler.distribute_dofs(fe);
998 *  
999 *   solution.reinit(dof_handler.n_dofs());
1000 *   system_rhs.reinit(dof_handler.n_dofs());
1001 *  
1002 *   constraints.clear();
1003 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1004 *  
1005 *   VectorTools::interpolate_boundary_values(dof_handler,
1006 *   BoundaryIds::cathode,
1007 *   Functions::ConstantFunction<dim>(
1008 *   -Constants::V0),
1009 *   constraints);
1010 *   VectorTools::interpolate_boundary_values(dof_handler,
1011 *   BoundaryIds::focus_element,
1012 *   Functions::ConstantFunction<dim>(
1013 *   -Constants::V0),
1014 *   constraints);
1015 *   VectorTools::interpolate_boundary_values(dof_handler,
1016 *   BoundaryIds::anode,
1017 *   Functions::ConstantFunction<dim>(
1018 *   +Constants::V0),
1019 *   constraints);
1020 *   constraints.close();
1021 *  
1022 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
1023 *   DoFTools::make_sparsity_pattern(dof_handler,
1024 *   dsp,
1025 *   constraints,
1026 *   /* keep_constrained_dofs = */ false);
1027 *   sparsity_pattern.copy_from(dsp);
1028 *  
1029 *   system_matrix.reinit(sparsity_pattern);
1030 *   }
1031 *  
1032 *  
1033 *  
1034 *   template <int dim>
1035 *   void CathodeRaySimulator<dim>::assemble_system()
1036 *   {
1037 *   system_matrix = 0;
1038 *   system_rhs = 0;
1039 *  
1040 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
1041 *  
1042 *   FEValues<dim> fe_values(fe,
1043 *   quadrature_formula,
1044 *   update_values | update_gradients |
1045 *   update_quadrature_points | update_JxW_values);
1046 *  
1047 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1048 *  
1049 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1050 *   Vector<double> cell_rhs(dofs_per_cell);
1051 *  
1052 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1053 *  
1054 *   for (const auto &cell : dof_handler.active_cell_iterators())
1055 *   {
1056 *   cell_matrix = 0;
1057 *   cell_rhs = 0;
1058 *  
1059 *   fe_values.reinit(cell);
1060 *  
1061 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
1062 *   for (const unsigned int i : fe_values.dof_indices())
1063 *   {
1064 *   for (const unsigned int j : fe_values.dof_indices())
1065 *   cell_matrix(i, j) +=
1066 *   (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
1067 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
1068 *   fe_values.JxW(q_index)); // dx
1069 *   }
1070 *  
1071 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1072 *   for (const auto &particle : particle_handler.particles_in_cell(cell))
1073 *   {
1074 *   const Point<dim> &reference_location =
1075 *   particle.get_reference_location();
1076 *   for (const unsigned int i : fe_values.dof_indices())
1077 *   cell_rhs(i) +=
1078 *   (fe.shape_value(i, reference_location) * // phi_i(x_p)
1079 *   (-Constants::electrons_per_particle * // N
1080 *   Constants::electron_charge)); // e
1081 *   }
1082 *  
1083 *   cell->get_dof_indices(local_dof_indices);
1084 *   constraints.distribute_local_to_global(
1085 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1086 *   }
1087 *   }
1088 *  
1089 *  
1090 *  
1091 *   template <int dim>
1092 *   void CathodeRaySimulator<dim>::solve_field()
1093 *   {
1094 *   SolverControl solver_control(1000, 1e-12);
1095 *   SolverCG<Vector<double>> solver(solver_control);
1096 *  
1097 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1098 *   preconditioner.initialize(system_matrix, 1.2);
1099 *  
1100 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
1101 *  
1102 *   constraints.distribute(solution);
1103 *   }
1104 *  
1105 *  
1106 *  
1107 *   template <int dim>
1108 *   void CathodeRaySimulator<dim>::refine_grid()
1109 *   {
1110 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1111 *  
1112 *   KellyErrorEstimator<dim>::estimate(dof_handler,
1113 *   QGauss<dim - 1>(fe.degree + 1),
1114 *   {},
1115 *   solution,
1116 *   estimated_error_per_cell);
1117 *  
1118 *   GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1119 *   estimated_error_per_cell,
1120 *   0.1,
1121 *   0.03);
1122 *  
1123 *   triangulation.execute_coarsening_and_refinement();
1124 *   }
1125 *  
1126 *  
1127 *  
1128 *   template <int dim>
1129 *   void CathodeRaySimulator<dim>::create_particles()
1130 *   {
1131 *   FEFaceValues<dim> fe_face_values(fe,
1132 *   QMidpoint<dim - 1>(),
1133 *   update_quadrature_points |
1134 *   update_gradients |
1135 *   update_normal_vectors);
1136 *  
1137 *   std::vector<Tensor<1, dim>> solution_gradients(
1138 *   fe_face_values.n_quadrature_points);
1139 *  
1140 *   for (const auto &cell : dof_handler.active_cell_iterators())
1141 *   for (const auto &face : cell->face_iterators())
1142 *   if (face->at_boundary() &&
1143 *   (face->boundary_id() == BoundaryIds::cathode))
1144 *   {
1145 *   fe_face_values.reinit(cell, face);
1146 *  
1147 *   const FEValuesExtractors::Scalar electric_potential(0);
1148 *   fe_face_values[electric_potential].get_function_gradients(
1149 *   solution, solution_gradients);
1150 *   for (const unsigned int q_point :
1151 *   fe_face_values.quadrature_point_indices())
1152 *   {
1153 *   const Tensor<1, dim> E = solution_gradients[q_point];
1154 *  
1155 *   if ((E * fe_face_values.normal_vector(q_point) < 0) &&
1156 *   (E.norm() > Constants::E_threshold))
1157 *   {
1158 *   const Point<dim> &location =
1159 *   fe_face_values.quadrature_point(q_point);
1160 *  
1161 *   Particles::Particle<dim> new_particle;
1162 *   new_particle.set_location(location);
1163 *   new_particle.set_reference_location(
1164 *   mapping.transform_real_to_unit_cell(cell, location));
1165 *   new_particle.set_id(next_unused_particle_id);
1166 *   particle_handler.insert_particle(new_particle, cell);
1167 *  
1168 *   ++next_unused_particle_id;
1169 *   }
1170 *   }
1171 *   }
1172 *  
1173 *   particle_handler.update_cached_numbers();
1174 *   }
1175 *  
1176 *  
1177 *  
1178 *   template <int dim>
1179 *   void CathodeRaySimulator<dim>::move_particles()
1180 *   {
1181 *   const double dt = time.get_next_step_size();
1182 *  
1183 *   Vector<double> solution_values(fe.n_dofs_per_cell());
1184 *   FEPointEvaluation<1, dim> evaluator(mapping, fe, update_gradients);
1185 *  
1186 *   for (const auto &cell : dof_handler.active_cell_iterators())
1187 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1188 *   {
1189 *   const typename Particles::ParticleHandler<
1190 *   dim>::particle_iterator_range particles_in_cell =
1191 *   particle_handler.particles_in_cell(cell);
1192 *  
1193 *   std::vector<Point<dim>> particle_positions;
1194 *   for (const auto &particle : particles_in_cell)
1195 *   particle_positions.push_back(particle.get_reference_location());
1196 *  
1197 *   cell->get_dof_values(solution, solution_values);
1198 *  
1199 *   evaluator.reinit(cell, particle_positions);
1200 *   evaluator.evaluate(make_array_view(solution_values),
1201 *   EvaluationFlags::gradients);
1202 *  
1203 *   {
1204 *   typename Particles::ParticleHandler<dim>::particle_iterator
1205 *   particle = particles_in_cell.begin();
1206 *   for (unsigned int particle_index = 0;
1207 *   particle != particles_in_cell.end();
1208 *   ++particle, ++particle_index)
1209 *   {
1210 *   const Tensor<1, dim> &E =
1211 *   evaluator.get_gradient(particle_index);
1212 *  
1213 *   const Tensor<1, dim> old_velocity(particle->get_properties());
1214 *  
1215 *   const Tensor<1, dim> acceleration =
1216 *   Constants::electron_charge / Constants::electron_mass * E;
1217 *  
1218 *   const Tensor<1, dim> new_velocity =
1219 *   old_velocity + acceleration * dt;
1220 *  
1221 *   particle->set_properties(new_velocity);
1222 *  
1223 *   const Point<dim> new_location =
1224 *   particle->get_location() + dt * new_velocity;
1225 *   particle->set_location(new_location);
1226 *   }
1227 *   }
1228 *   }
1229 *  
1230 *   particle_handler.sort_particles_into_subdomains_and_cells();
1231 *   }
1232 *  
1233 *  
1234 *  
1235 *   template <int dim>
1236 *   void CathodeRaySimulator<dim>::track_lost_particle(
1237 *   const typename Particles::ParticleIterator<dim> &particle,
1238 *   const typename Triangulation<dim>::active_cell_iterator &cell)
1239 *   {
1240 *   ++n_recently_lost_particles;
1241 *   ++n_total_lost_particles;
1242 *  
1243 *   const Point<dim> current_location = particle->get_location();
1244 *   const Point<dim> approximate_previous_location = cell->center();
1245 *  
1246 *   if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
1247 *   {
1248 *   const Tensor<1, dim> direction =
1249 *   (current_location - approximate_previous_location) /
1250 *   (current_location[0] - approximate_previous_location[0]);
1251 *  
1252 *   const double right_boundary_intercept =
1253 *   approximate_previous_location[1] +
1254 *   (4 - approximate_previous_location[0]) * direction[1];
1255 *   if ((right_boundary_intercept > 0.5) &&
1256 *   (right_boundary_intercept < 1.5))
1257 *   ++n_particles_lost_through_anode;
1258 *   }
1259 *   }
1260 *  
1261 *  
1262 *  
1263 *   template <int dim>
1264 *   void CathodeRaySimulator<dim>::update_timestep_size()
1265 *   {
1266 *   if (time.get_step_number() > 0)
1267 *   {
1268 *   double min_cell_size_over_velocity = std::numeric_limits<double>::max();
1269 *  
1270 *   for (const auto &cell : dof_handler.active_cell_iterators())
1271 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1272 *   {
1273 *   const double cell_size = cell->minimum_vertex_distance();
1274 *  
1275 *   double max_particle_velocity(0.0);
1276 *  
1277 *   for (const auto &particle :
1278 *   particle_handler.particles_in_cell(cell))
1279 *   {
1280 *   const Tensor<1, dim> velocity(particle.get_properties());
1281 *   max_particle_velocity =
1282 *   std::max(max_particle_velocity, velocity.norm());
1283 *   }
1284 *  
1285 *   if (max_particle_velocity > 0)
1286 *   min_cell_size_over_velocity =
1287 *   std::min(min_cell_size_over_velocity,
1288 *   cell_size / max_particle_velocity);
1289 *   }
1290 *  
1291 *   constexpr double c_safety = 0.5;
1292 *   time.set_desired_next_step_size(c_safety * 0.5 *
1293 *   min_cell_size_over_velocity);
1294 *   }
1295 *   else
1296 *   {
1297 *   const QTrapezoid<dim> vertex_quadrature;
1298 *   FEValues<dim> fe_values(fe, vertex_quadrature, update_gradients);
1299 *  
1300 *   std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
1301 *  
1302 *   double min_timestep = std::numeric_limits<double>::max();
1303 *  
1304 *   for (const auto &cell : dof_handler.active_cell_iterators())
1305 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1306 *   {
1307 *   const double cell_size = cell->minimum_vertex_distance();
1308 *  
1309 *   fe_values.reinit(cell);
1310 *   fe_values.get_function_gradients(solution, field_gradients);
1311 *  
1312 *   double max_E = 0;
1313 *   for (const auto q_point : fe_values.quadrature_point_indices())
1314 *   max_E = std::max(max_E, field_gradients[q_point].norm());
1315 *  
1316 *   if (max_E > 0)
1317 *   min_timestep =
1318 *   std::min(min_timestep,
1319 *   std::sqrt(0.5 * cell_size *
1320 *   Constants::electron_mass /
1321 *   Constants::electron_charge / max_E));
1322 *   }
1323 *  
1324 *   time.set_desired_next_step_size(min_timestep);
1325 *   }
1326 *   }
1327 *  
1328 *  
1329 *  
1330 *   template <int dim>
1331 *   class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
1332 *   {
1333 *   public:
1334 *   ElectricFieldPostprocessor()
1335 *   : DataPostprocessorVector<dim>("electric_field", update_gradients)
1336 *   {}
1337 *  
1338 *   virtual void evaluate_scalar_field(
1339 *   const DataPostprocessorInputs::Scalar<dim> &input_data,
1340 *   std::vector<Vector<double>> &computed_quantities) const override
1341 *   {
1342 *   AssertDimension(input_data.solution_gradients.size(),
1343 *   computed_quantities.size());
1344 *  
1345 *   for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
1346 *   {
1347 *   AssertDimension(computed_quantities[p].size(), dim);
1348 *   for (unsigned int d = 0; d < dim; ++d)
1349 *   computed_quantities[p][d] = input_data.solution_gradients[p][d];
1350 *   }
1351 *   }
1352 *   };
1353 *  
1354 *  
1355 *  
1356 *   template <int dim>
1357 *   void CathodeRaySimulator<dim>::output_results() const
1358 *   {
1359 *   {
1360 *   ElectricFieldPostprocessor<dim> electric_field;
1361 *   DataOut<dim> data_out;
1362 *   data_out.attach_dof_handler(dof_handler);
1363 *   data_out.add_data_vector(solution, "electric_potential");
1364 *   data_out.add_data_vector(solution, electric_field);
1365 *   data_out.build_patches();
1366 *  
1367 *   DataOutBase::VtkFlags output_flags;
1368 *   output_flags.time = time.get_current_time();
1369 *   output_flags.cycle = time.get_step_number();
1370 *   output_flags.physical_units["electric_potential"] = "V";
1371 *   output_flags.physical_units["electric_field"] = "V/m";
1372 *  
1373 *   data_out.set_flags(output_flags);
1374 *  
1375 *   std::ofstream output("solution-" +
1376 *   Utilities::int_to_string(time.get_step_number(), 4) +
1377 *   ".vtu");
1378 *   data_out.write_vtu(output);
1379 *   }
1380 *  
1381 *   {
1382 *   Particles::DataOut<dim> particle_out;
1383 *   particle_out.build_patches(
1384 *   particle_handler,
1385 *   std::vector<std::string>(dim, "velocity"),
1386 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>(
1387 *   dim, DataComponentInterpretation::component_is_part_of_vector));
1388 *  
1389 *   DataOutBase::VtkFlags output_flags;
1390 *   output_flags.time = time.get_current_time();
1391 *   output_flags.cycle = time.get_step_number();
1392 *   output_flags.physical_units["velocity"] = "m/s";
1393 *  
1394 *   particle_out.set_flags(output_flags);
1395 *  
1396 *   std::ofstream output("particles-" +
1397 *   Utilities::int_to_string(time.get_step_number(), 4) +
1398 *   ".vtu");
1399 *   particle_out.write_vtu(output);
1400 *   }
1401 *   }
1402 *  
1403 *  
1404 *  
1405 * @endcode
1406 *
1407 *
1408 * <a name="step_83-CathodeRaySimulatorserialize"></a>
1409 * <h4>CathodeRaySimulator::serialize()</h4>
1410 *
1411
1412 *
1413 * The first of the new function is the one that is called by the
1414 * BOOST Serialization framework to serialize and deserialize the
1415 * data of this class. It has already been discussed in the introduction
1416 * to this program and so does not provide any surprises. All it does is
1417 * write those member variables of the current class that cannot be
1418 * re-created easily into an archive, or read these members from it.
1419 * (Whether `operator &` facilitates a write or read operation depends on
1420 * whether the `Archive` type is an output or input archive.)
1421 *
1422
1423 *
1424 * The function takes a second argument, `version`, that can be used to
1425 * create checkpoints that have version numbers. This is useful if one
1426 * evolves programs by adding more member variables, but still wants to
1427 * retain the ability to read checkpoint files created with earlier
1428 * versions of the program. The `version` variable would, in that case,
1429 * be used to represent which version of the program wrote the file,
1430 * and if necessary to read only those variables that were written with
1431 * that past version, finding a different way to initialize the new member
1432 * variables that have been added since then. We will not make use of this
1433 * ability here.
1434 *
1435
1436 *
1437 * Finally, while the program that indents all deal.II source files format
1438 * the following code as
1439 * <div class=CodeFragmentInTutorialComment>
1440 * @code
1441 * ar &solution;
1442 * @endcode
1443 * </div>
1444 * as if we are taking the address of the `triangulation` variable, the
1445 * way you *should* read the code is as
1446 * <div class=CodeFragmentInTutorialComment>
1447 * @code
1448 * ar & solution;
1449 * @endcode
1450 * </div>
1451 * where `operator &` is a binary operator that could either be interpreted
1452 * as `operator <<` for output or `operator >>` for input.
1453 *
1454
1455 *
1456 * As discussed in the introduction, we do not serialize the `triangulation`
1457 * member variable, instead leaving that to separate calls in the
1458 * `checkpoint()` and `restart()` functions below.
1459 *
1460 * @code
1461 *   template <int dim>
1462 *   template <class Archive>
1463 *   void CathodeRaySimulator<dim>::serialize(Archive &ar,
1464 *   const unsigned int /* version */)
1465 *   {
1466 *   ar &solution;
1467 *   ar &particle_handler;
1468 *   ar &next_unused_particle_id;
1469 *   ar &n_recently_lost_particles;
1470 *   ar &n_total_lost_particles;
1471 *   ar &n_particles_lost_through_anode;
1472 *   ar &time;
1473 *   }
1474 *  
1475 *  
1476 *  
1477 * @endcode
1478 *
1479 *
1480 * <a name="step_83-CathodeRaySimulatorcheckpoint"></a>
1481 * <h4>CathodeRaySimulator::checkpoint()</h4>
1482 *
1483
1484 *
1485 * The checkpoint function of the principal class of this program is then
1486 * quite straightforward: We create an output file (and check that it is
1487 * writable), create an output archive, and then move the serialized
1488 * contents of the current object (i.e., the `*this` object) into the
1489 * archive. The use of `operator<<` here calls the `serialize()` function
1490 * above with an output archive as argument. When the destructor of the
1491 * `archive` variable is called at the end of the code block within which
1492 * it lives, the entire archive is written into the output file stream it
1493 * is associated with.
1494 *
1495
1496 *
1497 * As mentioned in the introduction, "real" applications would not use
1498 * text-based archives as provided by the `boost::archive::text_oarchive`
1499 * class, but use binary and potentially compressed versions. This can
1500 * easily be achieved by using differently named classes, and the BOOST
1501 * documentation explains how to do that.
1502 *
1503 * @code
1504 *   template <int dim>
1505 *   void CathodeRaySimulator<dim>::checkpoint()
1506 *   {
1507 *   std::cout << "--- Writing checkpoint... ---" << std::endl << std::endl;
1508 *  
1509 *   {
1510 *   std::ofstream checkpoint_file("tmp.checkpoint_step_83");
1511 *   AssertThrow(checkpoint_file,
1512 *   ExcMessage(
1513 *   "Could not write to the <tmp.checkpoint_step_83> file."));
1514 *  
1515 *   boost::archive::text_oarchive archive(checkpoint_file);
1516 *  
1517 *   archive << *this;
1518 *   }
1519 *  
1520 * @endcode
1521 *
1522 * The second part of the serialization is all of the data that we can
1523 * attach to cells -- see the discussion about this in the introduction.
1524 * Here, the only data we attach to cells are the particles. We then
1525 * let the triangulation save these into a set of files that all start
1526 * with the same prefix as we chose above, namely "tmp.checkpoint":
1527 *
1528 * @code
1529 *   particle_handler.prepare_for_serialization();
1530 *   triangulation.save("tmp.checkpoint");
1531 *  
1532 *  
1533 * @endcode
1534 *
1535 * At this point, the serialized data of this file has ended up in a number
1536 * of files that all start with `tmp.checkpoint` file. As mentioned in the
1537 * introduction, we do not want to directly overwrite the checkpointing
1538 * files from the previous checkpoint operation, for fear that the program
1539 * may be interrupted *while writing the checkpoint files*. This would
1540 * result in corrupted files, and defeat the whole purpose of checkpointing
1541 * because one cannot restart from such a file. On the other hand, if we got
1542 * here, we know that the "tmp.checkpoint*" files were successfully written,
1543 * and we can rename it to "checkpoint*", in the process replacing the old
1544 * file.
1545 *
1546
1547 *
1548 * We do this move operation by calling the
1549 * [C++ function that does the renaming of
1550 * files](https://en.cppreference.com/w/cpp/filesystem/rename).
1551 * Note that it is documented as being for all practical purposes
1552 * "atomic", i.e., we do not need to worry that the program may be
1553 * interrupted somewhere within the renaming operation itself. Of
1554 * course, it is possible that we get interrupted somewhere between
1555 * renaming one file and the next, and that presents problems in
1556 * itself -- in essence, we want the entire renaming operation of all of
1557 * these files to be atomic. With a couple dozen lines of extra code, one
1558 * could address this issue (using strategies that databases use frequently:
1559 * if one operation fails, we need to
1560 * [rollback](https://en.wikipedia.org/wiki/Rollback_(data_management))
1561 * the entire transaction). For the purposes of this program, this is
1562 * perhaps too much, and we will simply hope that that doesn't happen,
1563 * perhaps based on the belief that renaming files is much faster than
1564 * writing them, and that unlike writing checkpoint files, renaming does not
1565 * require much memory or disk space and so does not risk running out of
1566 * either.
1567 *
1568
1569 *
1570 * As a consequence, the following code first loops over all files in
1571 * the current directory, picks out those that start with the string
1572 * "tmp.checkpoint", and puts them into a list. In a second step,
1573 * we loop over the list and rename each of these files to one
1574 * whose name consists of the "tmp.checkpoint*" file but stripped off
1575 * its first four characters (i.e., only the "checkpoint*" part). We use
1576 * this approach, rather than listing the files we want to rename,
1577 * because we do not actually know the names of the files written by
1578 * the Triangulation::save() function, though we know how each of these
1579 * file names starts.
1580 *
1581 * @code
1582 *   std::list<std::string> tmp_checkpoint_files;
1583 *   for (const auto &dir_entry : std::filesystem::directory_iterator("."))
1584 *   if (dir_entry.is_regular_file() &&
1585 *   (dir_entry.path().filename().string().find("tmp.checkpoint") == 0))
1586 *   tmp_checkpoint_files.push_back(dir_entry.path().filename().string());
1587 *  
1588 *   for (const std::string &filename : tmp_checkpoint_files)
1589 *   std::filesystem::rename(filename, filename.substr(4, std::string::npos));
1590 *   }
1591 *  
1592 *  
1593 * @endcode
1594 *
1595 *
1596 * <a name="step_83-CathodeRaySimulatorrestart"></a>
1597 * <h4>CathodeRaySimulator::restart()</h4>
1598 *
1599
1600 *
1601 * The restart function of this class then simply does the opposite:
1602 * It opens an input file (and triggers an error if that file cannot be
1603 * opened), associates an input archive with it, and then reads the
1604 * contents of the current object from it, again using the
1605 * `serialize()` function from above. Clearly, since we have written
1606 * data into a text-based archive above, we need to use the corresponding
1607 * `boost::archive::text_iarchive` class for reading.
1608 *
1609
1610 *
1611 * In a second step, we ask the triangulation to read in cell-attached
1612 * data, and then tell the Particles::ParticleHandler object to re-create
1613 * its information about all of the particles from the data just read.
1614 *
1615
1616 *
1617 * The function ends by printing a status message about having restarted:
1618 *
1619 * @code
1620 *   template <int dim>
1621 *   void CathodeRaySimulator<dim>::restart()
1622 *   {
1623 *   {
1624 *   std::ifstream checkpoint_file("checkpoint_step_83");
1625 *   AssertThrow(checkpoint_file,
1626 *   ExcMessage(
1627 *   "Could not read from the <checkpoint_step_83> file."));
1628 *  
1629 *   boost::archive::text_iarchive archive(checkpoint_file);
1630 *   archive >> *this;
1631 *   }
1632 *  
1633 *   triangulation.load("checkpoint");
1634 *   particle_handler.deserialize();
1635 *  
1636 *   std::cout << "--- Restarting at t=" << time.get_current_time()
1637 *   << ", dt=" << time.get_next_step_size() << std::endl
1638 *   << std::endl;
1639 *   }
1640 *  
1641 *  
1642 * @endcode
1643 *
1644 *
1645 * <a name="step_83-CathodeRaySimulatorrun"></a>
1646 * <h4>CathodeRaySimulator::run()</h4>
1647 *
1648
1649 *
1650 * The last member function of the principal class of this program is then the
1651 * driver. The driver takes a single argument to indicate if the simulation
1652 * is a restart. If it is not a restart, the mesh is set up and the problem is
1653 * solved like in @ref step_19 "step-19". If it is a restart, then we read in everything that
1654 * is a history variable from the checkpoint file via the `restart()`
1655 * function. Recall that everything that is inside the `if` block at the top
1656 * of the function is exactly like in @ref step_19 "step-19", as is almost everything that
1657 * follows:
1658 *
1659 * @code
1660 *   template <int dim>
1661 *   void CathodeRaySimulator<dim>::run(const bool do_restart)
1662 *   {
1663 *   if (do_restart == false)
1664 *   {
1665 *   make_grid();
1666 *  
1667 *   const unsigned int n_pre_refinement_cycles = 3;
1668 *   for (unsigned int refinement_cycle = 0;
1669 *   refinement_cycle < n_pre_refinement_cycles;
1670 *   ++refinement_cycle)
1671 *   {
1672 *   setup_system();
1673 *   assemble_system();
1674 *   solve_field();
1675 *   refine_grid();
1676 *   }
1677 *   }
1678 *   else
1679 *   {
1680 *   restart();
1681 *   }
1682 *  
1683 *   setup_system();
1684 *   do
1685 *   {
1686 *   std::cout << "Timestep " << time.get_step_number() + 1 << std::endl;
1687 *   std::cout << " Field degrees of freedom: "
1688 *   << dof_handler.n_dofs() << std::endl;
1689 *  
1690 *   assemble_system();
1691 *   solve_field();
1692 *  
1693 *   create_particles();
1694 *   std::cout << " Total number of particles in simulation: "
1695 *   << particle_handler.n_global_particles() << std::endl;
1696 *  
1697 *   n_recently_lost_particles = 0;
1698 *   update_timestep_size();
1699 *   move_particles();
1700 *  
1701 *   time.advance_time();
1702 *  
1703 *   output_results();
1704 *  
1705 *   std::cout << " Number of particles lost this time step: "
1706 *   << n_recently_lost_particles << std::endl;
1707 *   if (n_total_lost_particles > 0)
1708 *   std::cout << " Fraction of particles lost through anode: "
1709 *   << 1. * n_particles_lost_through_anode /
1710 *   n_total_lost_particles
1711 *   << std::endl;
1712 *  
1713 *   std::cout << std::endl
1714 *   << " Now at t=" << time.get_current_time()
1715 *   << ", dt=" << time.get_previous_step_size() << '.'
1716 *   << std::endl
1717 *   << std::endl;
1718 *  
1719 * @endcode
1720 *
1721 * The only other difference between this program and @ref step_19 "step-19" is that
1722 * we checkpoint the simulation every ten time steps:
1723 *
1724 * @code
1725 *   if (time.get_step_number() % 10 == 0)
1726 *   checkpoint();
1727 *   }
1728 *   while (time.is_at_end() == false);
1729 *   }
1730 *   } // namespace Step83
1731 *  
1732 *  
1733 * @endcode
1734 *
1735 *
1736 * <a name="step_83-Thecodemaincodefunction"></a>
1737 * <h3>The <code>main</code> function</h3>
1738 *
1739
1740 *
1741 * The final function of the program is then again the `main()` function. Its
1742 * overall structure is unchanged in all tutorial programs since @ref step_6 "step-6" and
1743 * so there is nothing new to discuss about this aspect.
1744 *
1745
1746 *
1747 * The only difference is that we need to figure out whether a restart was
1748 * requested, or whether the program should simply start from scratch when
1749 * called. We do this using a command line argument: The `argc` argument to
1750 * `main()` indicates how many command line arguments were provided when
1751 * the program was called (counting the name of the program as the zeroth
1752 * argument), and `argv` is an array of strings with as many elements as
1753 * `argc` that contains these command line arguments. So if you call
1754 * the program as
1755 * <div class=CodeFragmentInTutorialComment>
1756 * @code
1757 * ./step-83
1758 * @endcode
1759 * </div>
1760 * then `argc` will be 1, and `argv` will be the array with one element
1761 * and content `[ "./step-83" ]`. On the other hand, if you call the program
1762 * as
1763 * <div class=CodeFragmentInTutorialComment>
1764 * @code
1765 * ./step-83 restart
1766 * @endcode
1767 * </div>
1768 * then `argc` will be 2, and `argv` will be the array with two elements
1769 * and content `[ "./step-83", "restart" ]`. Every other choice should be
1770 * flagged as an error. The following try block then does exactly this:
1771 *
1772 * @code
1773 *   int main(int argc, char *argv[])
1774 *   {
1775 *   try
1776 *   {
1777 *   Step83::CathodeRaySimulator<2> cathode_ray_simulator;
1778 *  
1779 *   if (argc == 1)
1780 *   cathode_ray_simulator.run(false); // no restart
1781 *   else if ((argc == 2) && (std::string(argv[1]) == "restart"))
1782 *   cathode_ray_simulator.run(true); // yes restart
1783 *   else
1784 *   {
1785 *   std::cerr << "Error: The only allowed command line argument to this\n"
1786 *   " program is 'restart'."
1787 *   << std::endl;
1788 *   return 1;
1789 *   }
1790 *   }
1791 *   catch (std::exception &exc)
1792 *   {
1793 *   std::cerr << std::endl
1794 *   << std::endl
1795 *   << "----------------------------------------------------"
1796 *   << std::endl;
1797 *   std::cerr << "Exception on processing: " << std::endl
1798 *   << exc.what() << std::endl
1799 *   << "Aborting!" << std::endl
1800 *   << "----------------------------------------------------"
1801 *   << std::endl;
1802 *  
1803 *   return 1;
1804 *   }
1805 *   catch (...)
1806 *   {
1807 *   std::cerr << std::endl
1808 *   << std::endl
1809 *   << "----------------------------------------------------"
1810 *   << std::endl;
1811 *   std::cerr << "Unknown exception!" << std::endl
1812 *   << "Aborting!" << std::endl
1813 *   << "----------------------------------------------------"
1814 *   << std::endl;
1815 *   return 1;
1816 *   }
1817 *   return 0;
1818 *   }
1819 * @endcode
1820<a name="step_83-Results"></a><h1>Results</h1>
1821
1822
1823When you run this program, it produces the following output that is
1824almost exactly identical to what you get from @ref step_19 "step-19":
1825@code
1826Timestep 1
1827 Field degrees of freedom: 4989
1828 Total number of particles in simulation: 20
1829 Number of particles lost this time step: 0
1830
1831 Now at t=2.12647e-07, dt=2.12647e-07.
1832
1833Timestep 2
1834 Field degrees of freedom: 4989
1835 Total number of particles in simulation: 40
1836 Number of particles lost this time step: 0
1837
1838 Now at t=4.14362e-07, dt=2.01715e-07.
1839
1840Timestep 3
1841 Field degrees of freedom: 4989
1842 Total number of particles in simulation: 60
1843 Number of particles lost this time step: 0
1844
1845 Now at t=5.23066e-07, dt=1.08704e-07.
1846
1847Timestep 4
1848 Field degrees of freedom: 4989
1849 Total number of particles in simulation: 80
1850 Number of particles lost this time step: 0
1851
1852 Now at t=6.08431e-07, dt=8.53649e-08.
1853
1854Timestep 5
1855 Field degrees of freedom: 4989
1856 Total number of particles in simulation: 100
1857 Number of particles lost this time step: 0
1858
1859 Now at t=6.81935e-07, dt=7.35039e-08.
1860
1861Timestep 6
1862 Field degrees of freedom: 4989
1863 Total number of particles in simulation: 120
1864 Number of particles lost this time step: 0
1865
1866 Now at t=7.47864e-07, dt=6.59294e-08.
1867
1868Timestep 7
1869 Field degrees of freedom: 4989
1870 Total number of particles in simulation: 140
1871 Number of particles lost this time step: 0
1872
1873 Now at t=8.2516e-07, dt=7.72957e-08.
1874
1875Timestep 8
1876 Field degrees of freedom: 4989
1877 Total number of particles in simulation: 158
1878 Number of particles lost this time step: 0
1879
1880 Now at t=8.95325e-07, dt=7.01652e-08.
1881
1882Timestep 9
1883 Field degrees of freedom: 4989
1884 Total number of particles in simulation: 172
1885 Number of particles lost this time step: 0
1886
1887 Now at t=9.67852e-07, dt=7.25269e-08.
1888
1889Timestep 10
1890 Field degrees of freedom: 4989
1891 Total number of particles in simulation: 186
1892 Number of particles lost this time step: 0
1893
1894 Now at t=1.03349e-06, dt=6.56398e-08.
1895
1896--- Writing checkpoint... ---
1897
1898Timestep 11
1899 Field degrees of freedom: 4989
1900 Total number of particles in simulation: 198
1901 Number of particles lost this time step: 0
1902
1903 Now at t=1.11482e-06, dt=8.13268e-08.
1904
1905Timestep 12
1906 Field degrees of freedom: 4989
1907 Total number of particles in simulation: 206
1908 Number of particles lost this time step: 0
1909
1910 Now at t=1.18882e-06, dt=7.39967e-08.
1911
1912Timestep 13
1913 Field degrees of freedom: 4989
1914 Total number of particles in simulation: 212
1915 Number of particles lost this time step: 0
1916
1917 Now at t=1.26049e-06, dt=7.16705e-08.
1918
1919[...]
1920@endcode
1921The only thing that is different is the additional line after the tenth
1922output (which also appears after the 20th, 30th, etc., time step) indicating
1923that a checkpoint file has been written.
1924
1925Because we chose to use a text-based format for the checkpoint file (which
1926you should not do in production codes because it uses quite a lot of disk space),
1927we can actually inspect this file. It will look like this, with many many more
1928lines:
1929@code
193022 serialization::archive 18 0 0 0 0 0 7 0 0 3 1 0
19314989 -1.00000000000000000e+00 -1.00000000000000000e+00 -1.00000000000000000e+00 -9.96108134982226390e-01 -1.00000000000000000e+00 -9.98361082867431748e-01
1932[...]
1933@endcode
1934What each of these numbers represents is hard to tell in practice, and also
1935entirely unimporant for our current purposes -- it's
1936a representation of the many objects that make up this program's state, and
1937from which one can restore its state. The point simply being that this is what
1938serialization produces: A long list (a sequence) of bits that we can put
1939into a file, and that we can later read again to recreate the state of the
1940program.
1941
1942Now here's the fun part. Let's say you hit Control-C on the command
1943line at the point above (say, during time step 13 or 14). There's
1944a set of checkpoint files on disk that saved the state after ten time steps.
1945Based on the logic in `main()`, we should be able to restart from
1946this point if we run the program with
1947```
1948 ./step-83 restart
1949```
1950Indeed, this is what then happens:
1951@code
1952--- Restarting at t=1.03349e-06, dt=6.56398e-08
1953
1954Timestep 11
1955 Field degrees of freedom: 4989
1956 Total number of particles in simulation: 198
1957 Number of particles lost this time step: 0
1958
1959 Now at t=1.11482e-06, dt=8.13268e-08.
1960
1961Timestep 12
1962 Field degrees of freedom: 4989
1963 Total number of particles in simulation: 206
1964 Number of particles lost this time step: 0
1965
1966 Now at t=1.18882e-06, dt=7.39967e-08.
1967
1968[...]
1969@endcode
1970After the status message that shows that we are restarting, this is
1971indeed *the exact same output* for the following time steps we had
1972gotten previously -- in other words, saving the complete state
1973seems to have worked, and the program has continued as if it had
1974never been interrupted!
1975
1976
1977
1978<a name="extensions"></a>
1979<a name="step_83-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1980
1981
1982<a name="step_83-Makingefficiencyapriority"></a><h4> Making efficiency a priority </h4>
1983
1984
1985While the techniques we have shown here are fully sufficient for what
1986we want to do, namely checkpoint and restart computations, and are in
1987fact also fully sufficient for much larger code bases such as
1988[ASPECT](https://aspect.geodynamics.org), one could go beyond what is
1989still a relatively simple scheme.
1990
1991Specifically, among the things we need to recognize is that writing
1992large amounts of data to disk is expensive and can take a good long
1993time to finish -- for example, for large parallel computations with,
1994say, a billion unknowns, checkpoints can run into the hundred gigabyte
1995range or beyond. One may ask whether that could be avoided, or at
1996least whether we can mitigate the cost.
1997
1998One way to do that is to first serialize the state of the program into
1999a buffer in memory (like the `Archive` objects the `serialize()`
2000functions write to and read from), and once that is done, start a *separate*
2001thread do the writing while the rest of the program continues with computations.
2002This is useful because writing the data to disk often takes a
2003long time but not a lot of CPU power: It just takes time to move the
2004data through the network to the file server, and from there onto the
2005actual disks. This is something that might as well happen while we are
2006doing something useful again (namely, solving more time steps). Should
2007the machine crash during this phase, nothing is lost: As discussed in
2008the introduction, we are writing the checkpoint into a temporary file
2009(which will be lost in the case of a machine failure), but we have
2010kept the previous checkpoint around until we know that the temporary
2011file is complete and can be moved over the old one.
2012
2013The only thing we have to pay attention in this background-writing
2014scheme is that we cannot start with creating a new checkpoint while
2015the previous one is still being written in the background.
2016
2017Doing this all is not technically very difficult; the code just
2018requires more explanation than lines of code, and so we omit doing
2019this in the program here. But you can take a look at the
2020`MainLoop::output()` function of @ref step_69 "step-69" to see how such a code looks
2021like.
2022
2023A variation of this general approach is that each process writes its
2024data immediately, but into files that are held on fast file systems --
2025say, a node-local SSD rather than a file server shared by the entire
2026cluster. One would then just tell the operating system to move this
2027file to the centeral file server in a second step, and this step can
2028happen in the background at whatever speed the operating system can
2029provide. Or perhaps one leaves *most* of these files on the fast local
2030file system in hopes that the restart happens before these files are
2031purged (say, by a script that runs nightly) and only moves these files
2032to the permanent file system every tenth time we create a checkpoint.
2033
2034In all of these cases, the logic quickly becomes quite complicated. As
2035usual, the solution is not to re-invent the wheel: Libraries such as
2036[VeloC](https://www.anl.gov/mcs/veloc-very-low-overhead-transparent-multilevel-checkpointrestart),
2037developed by the Exascale Computing Project (ECP) already do all of
2038this and more, for codes that are orders of magnitude more complex
2039than the little example here.
2040
2041Separately, one might want to try to reduce the amount of time it
2042takes to serialize objects into a buffer in memory. As mentioned
2043above, we use the
2044[BOOST serialization
2045library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html)
2046for this task, but it is not the only player in town. One could, for
2047example, use the [bitsery](https://github.com/fraillt/bitsery), [cereal](https://github.com/USCiLab/cereal), or
2048[zpp](https://github.com/eyalz800/serializer) projects instead, which can be substantially faster than BOOST.
2049 *
2050 *
2051<a name="step_83-PlainProg"></a>
2052<h1> The plain program</h1>
2053@include "step-83.cc"
2054*/
Definition fe_q.h:554
void serialize(Archive &ar, const unsigned int version)
void save(Archive &ar, const unsigned int version) const
virtual void load(const std::string &file_basename) override
Definition tria.cc:2109
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
__global__ void set(Number *val, const Number s, const size_type N)
#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