]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-19 code base.
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 29 Jun 2020 01:02:02 +0000 (19:02 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 29 Jun 2020 01:04:01 +0000 (19:04 -0600)
examples/step-19/step-19.cc [new file with mode: 0644]

diff --git a/examples/step-19/step-19.cc b/examples/step-19/step-19.cc
new file mode 100644 (file)
index 0000000..dd31838
--- /dev/null
@@ -0,0 +1,1112 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2020 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Authors: Wolfgang Bangerth, Rene Gassmoeller, Peter Munch, 2020.
+ */
+
+
+// @sect3{Include files}
+
+// The majority of the include files used in this program are
+// well known from step-6 and similar programs:
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_refinement.h>
+
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+
+
+// The ones that are new are only the following three: The first declared the
+// DiscreteTime class that helps us keep track of time in a time-dependent
+// simulation. The latter two provide all of the particle functionality,
+// namely a way to keep track of particles located on a mesh (the
+// Particles::ParticleHandler class) and the ability to output these
+// particles' locations and their properties for the purposes of
+// visualization (the Particles::DataOut class).
+#include <deal.II/base/discrete_time.h>
+#include <deal.II/particles/particle_handler.h>
+#include <deal.II/particles/data_out.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+
+// @sect3{Global definitions}
+
+// As is customary, we put everything that corresponds to the details of the
+// program into a namespace of its own. At the top, we define a few constants
+// for which we would rather use symbolic names than hard-coded numbers.
+//
+// Specifically, we define numbers for
+// @ref GlossBoundaryIndicator "boundary indicators"
+// for the various parts of the geometry, as well as the physical properties
+// of electrons and other specifics of the setup we use here.
+namespace Step19
+{
+  namespace BoundaryIds
+  {
+    constexpr types::boundary_id open          = 101;
+    constexpr types::boundary_id cathode       = 102;
+    constexpr types::boundary_id focus_element = 103;
+    constexpr types::boundary_id anode         = 104;
+  } // namespace BoundaryIds
+
+  namespace Constants
+  {
+    constexpr double electron_mass   = 9.1093837015e-31;
+    constexpr double electron_charge = 1.602176634e-19;
+
+    constexpr double V0 = 1;
+
+    constexpr double E_threshold = 0.05;
+
+    constexpr double electrons_per_particle = 3e15;
+  } // namespace Constants
+
+
+  // @sect3{The main class}
+
+  // The following is then the main class of this program. It has,
+  // fundamentally, the same structure as step-6 and many other
+  // tutorial programs. This includes the majority of the member
+  // functions (with the purpose of the rest probably self-explanatory
+  // from their names) as well as only a small number of member
+  // variables beyond those of step-6, all of which are related to
+  // dealing with particles.
+  template <int dim>
+  class CathodeRaySimulator
+  {
+  public:
+    CathodeRaySimulator();
+
+    void run();
+
+  private:
+    void make_grid();
+    void setup_system();
+    void assemble_system();
+    void solve_field();
+    void refine_grid();
+
+    void create_particles();
+    void move_particles();
+    void track_lost_particle(
+      const typename Particles::ParticleIterator<dim> &        particle,
+      const typename Triangulation<dim>::active_cell_iterator &cell);
+
+
+    void update_timestep_size();
+    void output_results() const;
+
+    Triangulation<dim>        triangulation;
+    MappingQ<dim>             mapping;
+    FE_Q<dim>                 fe;
+    DoFHandler<dim>           dof_handler;
+    AffineConstraints<double> constraints;
+
+    SparseMatrix<double> system_matrix;
+    SparsityPattern      sparsity_pattern;
+
+    Vector<double> solution;
+    Vector<double> system_rhs;
+
+    Particles::ParticleHandler<dim> particle_handler;
+    types::particle_index           next_unused_particle_id;
+    types::particle_index           n_recently_lost_particles;
+    types::particle_index           n_total_lost_particles;
+    types::particle_index           n_particles_lost_through_anode;
+
+    DiscreteTime time;
+  };
+
+
+
+  // @sect3{The <code>CathodeRaySimulator</code> class implementation}
+
+  // @sect4{The <code>CathodeRaySimulator</code> constructor}
+
+  // So then let us get started on the implementation. What the constructor
+  // does is really only a straight-forward initialization of all of the member
+  // variables at the top. The only two worth mentioning are the
+  // `particle_handler`, which is handed a reference to the triangulation
+  // on which the particles will live (currently of course still empty,
+  // but the particle handler stores the reference and will use it once
+  // particles are added -- which happens after the triangulation is built).
+  // The other piece of information it gets is how many "properties"
+  // each particle needs to store. Here, all we need each particle to
+  // remember is its current velocity, i.e., a vector with `dim`
+  // components. There are, however, other intrinsic properties that
+  // each particle has and that the Particles::ParticleHandler class
+  // automatically and always makes sure are available; in particular,
+  // these are the current location of a particle, the cell it is on,
+  // and a particle's ID.
+  //
+  // The only other variable of interest is `time`, an object of type
+  // DiscreteTime. It keeps track of the current time we are in a
+  // time-dependent simulation, and is initialized with the start time
+  // (zero) and end time ($10^{-4}$). We will later set the time step
+  // size in `update_timestep_size()`.
+  //
+  // The body of the constructor consists of a piece of code we have
+  // already discussed in the introduction. Namely, we make sure that the
+  // `track_lost_particle()` function is called by the `particle_handler`
+  // object every time a particle leaves the domain.
+  template <int dim>
+  CathodeRaySimulator<dim>::CathodeRaySimulator()
+    : mapping(1)
+    , fe(2)
+    , dof_handler(triangulation)
+    , particle_handler(triangulation, mapping, /*n_properties=*/dim)
+    , next_unused_particle_id(0)
+    , n_recently_lost_particles(0)
+    , n_total_lost_particles(0)
+    , n_particles_lost_through_anode(0)
+    , time(0, 1e-4)
+  {
+    particle_handler.signals.particle_lost.connect(
+      [this](const typename Particles::ParticleIterator<dim> &        particle,
+             const typename Triangulation<dim>::active_cell_iterator &cell) {
+        this->track_lost_particle(particle, cell);
+      });
+  }
+
+
+
+  // @sect4{The <code>CathodeRaySimulator::make_grid</code> function}
+
+  // The next function is then responsible for generating the mesh on which
+  // we want to solve. Recall how the domain looks like:
+  //   <p align="center">
+  //     <img
+  //     src="https://www.dealii.org/images/steps/developer/step-19.geometry.png"
+  //          alt="The geometry used in this program"
+  //          width="600">
+  //   </p>
+  // We subdivide this geometry into a mesh of $4\times 2$ cells that looks
+  // like this:
+  // @code
+  //   *---*---*---*---*
+  //   \   |   |   |   |
+  //    *--*---*---*---*
+  //   /   |   |   |   |
+  //   *---*---*---*---*
+  // @endcode
+  // The way this is done is by first defining where the $15=5\times 3$
+  // vertices are located -- here, we say that they are on integer points
+  // with the middle one on the left side moved to the right by a value of
+  // `delta=0.5`.
+  //
+  // In the following, we then have to say which vertices together form
+  // the 8 cells. The following code is then entirely equivalent to what
+  // we also do in step-14:
+  template <int dim>
+  void CathodeRaySimulator<dim>::make_grid()
+  {
+    static_assert(dim == 2,
+                  "This function is currently only implemented for 2d.");
+
+    const double       delta = 0.5;
+    const unsigned int nx    = 5;
+    const unsigned int ny    = 3;
+
+    const std::vector<Point<dim>> vertices //
+      = {{0, 0},
+         {1, 0},
+         {2, 0},
+         {3, 0},
+         {4, 0},
+         {delta, 1},
+         {1, 1},
+         {2, 1},
+         {3, 1},
+         {4, 1},
+         {0, 2},
+         {1, 2},
+         {2, 2},
+         {3, 2},
+         {4, 2}};
+    AssertDimension(vertices.size(), nx * ny);
+
+    const std::vector<unsigned int> cell_vertices[(nx - 1) * (ny - 1)] = {
+      {0, 1, nx + 0, nx + 1},
+      {1, 2, nx + 1, nx + 2},
+      {2, 3, nx + 2, nx + 3},
+      {3, 4, nx + 3, nx + 4},
+
+      {5, nx + 1, 2 * nx + 0, 2 * nx + 1},
+      {nx + 1, nx + 2, 2 * nx + 1, 2 * nx + 2},
+      {nx + 2, nx + 3, 2 * nx + 2, 2 * nx + 3},
+      {nx + 3, nx + 4, 2 * nx + 3, 2 * nx + 4}};
+
+    // With these arrays out of the way, we can move to slightly higher
+    // higher-level data structures. We create a vector of CellData
+    // objects that store for each cell to be created the vertices in
+    // question as well as the @ref GlossMaterialId "material id" (which
+    // we will here simply set to zero since we don't use it in the program).
+    //
+    // This information is then handed to the
+    // Triangulation::create_triangulation() function, and the mesh is twice
+    // globally refined.
+    std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
+    for (unsigned int i = 0; i < cells.size(); ++i)
+      {
+        cells[i].vertices    = cell_vertices[i];
+        cells[i].material_id = 0;
+      }
+
+    triangulation.create_triangulation(
+      {std::begin(vertices), std::end(vertices)},
+      cells,
+      SubCellData()); // No boundary information
+
+    triangulation.refine_global(2);
+
+    // The remaining part of the function loops over all cells and their faces,
+    // and if a face is at the boundary determines which boundary indicator
+    // should be applied to it. The various conditions should make sense if
+    // you compare the code with the picture of the geometry above.
+    //
+    // Once done with this step, we refine the mesh once more globally.
+    for (auto &cell : triangulation.active_cell_iterators())
+      for (auto &face : cell->face_iterators())
+        if (face->at_boundary())
+          {
+            if ((face->center()[0] > 0) && (face->center()[0] < 0.5) &&
+                (face->center()[1] > 0) && (face->center()[1] < 2))
+              face->set_boundary_id(BoundaryIds::cathode);
+            else if ((face->center()[0] > 0) && (face->center()[0] < 2))
+              face->set_boundary_id(BoundaryIds::focus_element);
+            else if ((face->center()[0] > 4 - 1e-12) &&
+                     ((face->center()[1] > 1.5) || (face->center()[1] < 0.5)))
+              face->set_boundary_id(BoundaryIds::anode);
+            else
+              face->set_boundary_id(BoundaryIds::open);
+          }
+
+    triangulation.refine_global(1);
+  }
+
+
+  // @sect4{The <code>CathodeRaySimulator::setup_system</code> function}
+
+  // The next function in this program deals with setting up the various
+  // objects related to solving the partial differential equations. It is
+  // in essence a copy of the corresponding function in step-6 and requires
+  // no further discussion.
+  template <int dim>
+  void CathodeRaySimulator<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    solution.reinit(dof_handler.n_dofs());
+    system_rhs.reinit(dof_handler.n_dofs());
+
+    constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             BoundaryIds::cathode,
+                                             Functions::ConstantFunction<dim>(
+                                               -Constants::V0),
+                                             constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             BoundaryIds::focus_element,
+                                             Functions::ConstantFunction<dim>(
+                                               -Constants::V0),
+                                             constraints);
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             BoundaryIds::anode,
+                                             Functions::ConstantFunction<dim>(
+                                               +Constants::V0),
+                                             constraints);
+    constraints.close();
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    dsp,
+                                    constraints,
+                                    /*keep_constrained_dofs = */ false);
+    sparsity_pattern.copy_from(dsp);
+
+    system_matrix.reinit(sparsity_pattern);
+  }
+
+
+  // @sect4{The <code>CathodeRaySimulator::assemble_system</code> function}
+
+  // The same is true for the function that assembles the linear system to be
+  // solved in each time step. At least that is true for the computation
+  // of the matrix entries, which is again in essence a copy of the
+  // corresponding function in step-6:
+  template <int dim>
+  void CathodeRaySimulator<dim>::assemble_system()
+  {
+    system_matrix = 0;
+    system_rhs    = 0;
+
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+    FEValues<dim> fe_values(fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        cell_matrix = 0;
+        cell_rhs    = 0;
+
+        fe_values.reinit(cell);
+
+        for (const unsigned int q_index : fe_values.quadrature_point_indices())
+          for (const unsigned int i : fe_values.dof_indices())
+            {
+              for (const unsigned int j : fe_values.dof_indices())
+                cell_matrix(i, j) +=
+                  (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
+                   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
+                   fe_values.JxW(q_index));           // dx
+            }
+
+        // The only interesting part of this function is how it forms the right
+        // hand side of the linear system. Recall that the right hand side
+        // of the PDE is
+        // @f[
+        //   \sum_p (N e)\delta(\mathbf x-\mathbf x_p),
+        // @f]
+        // where we have used $p$ to index the particles here to avoid
+        // confusion with the shape function $\varphi_i$; $\mathbf x_p$
+        // is the position of the $p$th particle.
+        //
+        // When multiplied by a test function $\varphi_i$ and integrated over
+        // the domain results in a right hand side vector
+        // @f{align*}{
+        //   F_i &= \int_\Omega \varphi_i (\mathbf x)\left[
+        //                \sum_p (N e)\delta(\mathbf x-\mathbf x_p) \right] dx
+        //   \\  &=  \sum_p (N e) \varphi_i(\mathbf x_p).
+        // @f}
+        // Note that the final line no longer contains an integral, and
+        // consequently also no occurrence of $dx$ which would require the
+        // appearance of the `JxW` symbol in our code.
+        //
+        // For a given cell $K$, this cell's contribution to the right hand
+        // side is then
+        // @f{align*}{
+        //   F_i^K &= \sum_{p, \mathbf x_p\in K} (N e) \varphi_i(\mathbf x_p),
+        // @f}
+        // i.e., we only have to worry about those particles that are actually
+        // located on the current cell $K$.
+        //
+        // In practice, what we do here is the following: If there are any
+        // particles on the current cell, then we first obtain an iterator range
+        // pointing to the first particle of that cell as well as the particle
+        // past the last one on this cell (or the end iterator) -- i.e., a
+        // half-open range as is common for C++ functions. Knowing now the
+        // number of particles, we can start to collect their reference
+        // locations (with respect to the reference cell), which we have stored
+        // with each particle when they were created or when it was last moved
+        // (see below).
+        //
+        // @note It is worth pointing out that calling the
+        //   Particles::ParticleHandler::particles_in_cell() and
+        //   Particles::ParticleHandler::n_particles_in_cell() functions is not
+        //   very effective on problems with a large number of particles. But it
+        //   illustrates the easiest way to write this algorithm, and so we are
+        //   willing to incur this cost for the moment for expository purposes.
+        //   We discuss the issue in more detail in the
+        //   <a href="#extensions">"possibilities for extensions" section</a>
+        //   below, and use a better approach in step-70, for example.
+        if (particle_handler.n_particles_in_cell(cell) > 0)
+          {
+            std::vector<Point<dim>> particle_reference_locations;
+
+            const typename Particles::ParticleHandler<
+              dim>::particle_iterator_range particles_in_cell =
+              particle_handler.particles_in_cell(cell);
+
+            const unsigned int n_particles_in_cell =
+              particle_handler.n_particles_in_cell(cell);
+
+            particle_reference_locations.resize(n_particles_in_cell);
+
+            {
+              typename Particles::ParticleHandler<dim>::particle_iterator
+                particle = particles_in_cell.begin();
+              for (unsigned int particle_index = 0;
+                   particle != particles_in_cell.end();
+                   ++particle, ++particle_index)
+                particle_reference_locations[particle_index] =
+                  particle->get_reference_location();
+            }
+
+            // Now that we know where the particles on the current cell are
+            // located with regard to the reference cell's coordinate system, we
+            // can create a Quadrature object with these locations and then an
+            // FEValues object that we will use to evaluate the shape functions
+            // at these locations. The contribution to the right hand side
+            // vector then immediately follows from the formula shown above.
+            // Note again the absence of the call to FEValues::JxW that would
+            // have to be present if we were evaluating an integral.
+            const Quadrature<dim> quadrature_formula(
+              particle_reference_locations);
+            FEValues<dim> fe_value(mapping,
+                                   fe,
+                                   quadrature_formula,
+                                   update_values);
+
+            for (const unsigned int q_index :
+                 fe_values.quadrature_point_indices())
+              for (const unsigned int i : fe_values.dof_indices())
+                cell_rhs(i) +=
+                  (fe_values.shape_value(i, q_index) *   // phi_i(x_p)
+                   (-Constants::electrons_per_particle * // N
+                    Constants::electron_charge));        // e
+          }
+
+
+        // Finally, we can copy the contributions of this cell into
+        // the global matrix and right hand side vector:
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global(
+          cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+      }
+  }
+
+
+  // @sect4{CathodeRaySimulator::solve}
+
+  // The function that solves the linear system is then again exactly as in
+  // step-6:
+  template <int dim>
+  void CathodeRaySimulator<dim>::solve_field()
+  {
+    SolverControl            solver_control(1000, 1e-12);
+    SolverCG<Vector<double>> solver(solver_control);
+
+    PreconditionSSOR<SparseMatrix<double>> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+
+    solver.solve(system_matrix, solution, system_rhs, preconditioner);
+
+    constraints.distribute(solution);
+  }
+
+
+  // @sect4{CathodeRaySimulator::refine_grid}
+
+  // The final field-related function is the one that refines the grid. We will
+  // call it a number of times in the first time step to obtain a mesh that
+  // is well-adapted to the structure of the solution and, in particular,
+  // resolves the various singularities in the solution that are due to
+  // re-entrant corners and places where the boundary condition type
+  // changes. You might want to refer to step-6 again for more details:
+  template <int dim>
+  void CathodeRaySimulator<dim>::refine_grid()
+  {
+    Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
+
+    KellyErrorEstimator<dim>::estimate(dof_handler,
+                                       QGauss<dim - 1>(fe.degree + 1),
+                                       {},
+                                       solution,
+                                       estimated_error_per_cell);
+
+    GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.1,
+                                                    0.03);
+
+    triangulation.execute_coarsening_and_refinement();
+  }
+
+
+  // @sect4{CathodeRaySimulator::create_particles}
+
+  // Let us now turn to the functions that deal with particles. The first one
+  // is about the creation of particles. As mentioned in the introduction,
+  // we want to create a particle at points of the cathode if the the electric
+  // field $\mathbf E=\nabla V$ exceeds a certain threshold, i.e., if
+  // $|\mathbf E| \ge E_\text{threshold}$, and if furthermore the electric field
+  // points into the domain (i.e., if $\mathbf E \cdot \mathbf n < 0$). As is
+  // common in the finite element method, we evaluate fields (and their
+  // derivatives) at specific evaluation points; typically, these are
+  // "quadrature points", and so we create a "quadrature formula" that we will
+  // use to designate the points at which we want to evaluate the solution.
+  // Here, we will simply take QMidpoint implying that we will only check the
+  // threshold condition at the midpoints of faces. We then use this to
+  // initialize an object of type FEFaceValues to evaluate the solution at these
+  // points.
+  //
+  // All of this will then be used in a loop over all cells, their faces, and
+  // specifically those faces that are at the boundary and, moreover, the
+  // cathode part of the boundary.
+  template <int dim>
+  void CathodeRaySimulator<dim>::create_particles()
+  {
+    FEFaceValues<dim> fe_face_values(fe,
+                                     QMidpoint<dim - 1>(),
+                                     update_quadrature_points |
+                                       update_gradients |
+                                       update_normal_vectors);
+
+    std::vector<Tensor<1, dim>> solution_gradients(
+      fe_face_values.n_quadrature_points);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary() &&
+            (face->boundary_id() == BoundaryIds::cathode))
+          {
+            fe_face_values.reinit(cell, face);
+
+            // So we have found a face on the cathode. Next, we let the
+            // FEFaceValues object compute the gradient of the solution at each
+            // "quadrature" point, and extract the electric field vector from
+            // the gradient in the form of a Tensor variable through the methods
+            // discussed in the
+            // @ref vector_valued "vector-valued problems" documentation module.
+            const FEValuesExtractors::Scalar electric_potential(0);
+            fe_face_values[electric_potential].get_function_gradients(
+              solution, solution_gradients);
+            for (const unsigned int q_point :
+                 fe_face_values.quadrature_point_indices())
+              {
+                const Tensor<1, dim> E = solution_gradients[q_point];
+
+                // Electrons can only escape the cathode if the electric field
+                // strength exceeds a threshold and,
+                // crucially, if the electric field points *into* the domain.
+                // Once we have that checked, we create a new
+                // Particles::Particle object at this location and insert it
+                // into the Particles::ParticleHandler object with a unique ID.
+                //
+                // The only thing that may be not obvious here is that we also
+                // associate with this particle the location in the reference
+                // coordinates of the cell we are currently on. This is done
+                // because we will in downstream functions compute quantities
+                // such as the electric field at the location of the particle
+                // (e.g., to compute the forces that act on it when updating its
+                // position in each time step). Evaluating a finite element
+                // field at arbitrary coordinates is quite an expensive
+                // operation because shape functions are really only defined on
+                // the reference cell, and so when asking for the electric field
+                // at an arbitrary point requires us first to determine what the
+                // reference coordinates of that point are. To avoid having to
+                // do this over and over, we determine these coordinates once
+                // and for all and then store these reference coordinates
+                // directly with the particle.
+                if ((E * fe_face_values.normal_vector(q_point) < 0) &&
+                    (E.norm() > Constants::E_threshold))
+                  {
+                    const Point<dim> location =
+                      fe_face_values.quadrature_point(q_point);
+
+                    Particles::Particle<dim> new_particle;
+                    new_particle.set_location(location);
+                    new_particle.set_reference_location(
+                      mapping.transform_real_to_unit_cell(cell, location));
+                    new_particle.set_id(next_unused_particle_id);
+                    particle_handler.insert_particle(new_particle, cell);
+
+                    ++next_unused_particle_id;
+                  }
+              }
+          }
+
+    // At the end of all of these insertions, we let the `particle_handler`
+    // update some internal statistics about the particles it stores.
+    particle_handler.update_cached_numbers();
+  }
+
+
+  // @sect4{CathodeRaySimulator::move_particles}
+
+  // The second particle-related function is the one that moves the particles
+  // in each time step. To do this, we have to loop over all cells, the
+  // particles in each cell, and evaluate the electric field at each of the
+  // particles' positions.
+  //
+  // The approach used here is conceptually the same used in the
+  // `assemble_system()` function: We loop over all cells, find the particles
+  // located there (with the same caveat about the inefficiency of the algorithm
+  // used here to find these particles), and create an FEValues object with
+  // these positions:
+  template <int dim>
+  void CathodeRaySimulator<dim>::move_particles()
+  {
+    const double dt = time.get_next_step_size();
+
+    std::vector<Point<dim>>     particle_positions;
+    std::vector<Tensor<1, dim>> field_gradients;
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (particle_handler.n_particles_in_cell(cell) > 0)
+        {
+          const typename Particles::ParticleHandler<
+            dim>::particle_iterator_range particles_in_cell =
+            particle_handler.particles_in_cell(cell);
+
+          const unsigned int n_particles_in_cell =
+            particle_handler.n_particles_in_cell(cell);
+
+          particle_positions.resize(n_particles_in_cell);
+          {
+            typename Particles::ParticleHandler<dim>::particle_iterator
+              particle = particles_in_cell.begin();
+            for (unsigned int particle_index = 0;
+                 particle != particles_in_cell.end();
+                 ++particle, ++particle_index)
+              particle_positions[particle_index] =
+                particle->get_reference_location();
+          }
+
+          const Quadrature<dim> quadrature_formula(particle_positions);
+          FEValues<dim>         fe_value(mapping,
+                                 fe,
+                                 quadrature_formula,
+                                 update_gradients);
+
+          fe_value.reinit(cell);
+
+          // Then we can ask the FEValues object for the gradients of the
+          // solution (i.e., the electric field $\mathbf E$) at these locations
+          // and loop over the individual particles:
+          field_gradients.resize(n_particles_in_cell);
+          fe_value.get_function_gradients(solution, field_gradients);
+
+          {
+            typename Particles::ParticleHandler<dim>::particle_iterator
+              particle = particles_in_cell.begin();
+            for (unsigned int particle_index = 0;
+                 particle != particles_in_cell.end();
+                 ++particle, ++particle_index)
+              {
+                const Tensor<1, dim> E = field_gradients[particle_index];
+
+                // Having now obtained the electric field at the location of one
+                // of the particles, we use this to update first the velocity
+                // and then the position. To do so, let us first get the old
+                // velocity out of the properties stored with the particle,
+                // compute the acceleration, update the velocity, and store this
+                // new velocity again in the properties of the particle. Recall
+                // that this corresponds to the first of the following set of
+                // update equations discussed in the introduction:
+                // @f{align*}{
+                //     \frac{{\mathbf v}_i^{(n)}
+                //           -{\mathbf v}_i^{(n-1)}}{\Delta t}
+                //     &= \frac{e\nabla V^{(n)}}{m}
+                //  \\ \frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}
+                //          {\Delta t} &= {\mathbf v}_i^{(n)}.
+                // @f}
+                const Tensor<1, dim> old_velocity(particle->get_properties());
+
+                const Tensor<1, dim> acceleration =
+                  Constants::electron_charge / Constants::electron_mass * E;
+
+                const Tensor<1, dim> new_velocity =
+                  old_velocity + acceleration * dt;
+
+                particle->set_properties(make_array_view(new_velocity));
+
+                // With the new velocity, we can then also update the location
+                // of the particle and set tell the particle about it.
+                const Point<dim> new_location =
+                  particle->get_location() + dt * old_velocity;
+                particle->set_location(new_location);
+              }
+          }
+        }
+
+    // Having updated the locations and properties (i.e., velocities) of all
+    // particles, we need to make sure that the `particle_handler` again knows
+    // which cells they are in, and what their locations in the coordinate
+    // system of the reference cell are. The following function does that. (It
+    // also makes sure that, in parallel computations, particles are moved from
+    // one processor to another processor if a particle moves from the subdomain
+    // owned by the former to the subdomain owned by the latter.)
+    particle_handler.sort_particles_into_subdomains_and_cells();
+  }
+
+
+  // @sect4{CathodeRaySimulator::track_lost_particle}
+
+  // The final particle-related function is the one that is called whenever a
+  // particle is lost from the simulation. This typically happens if it leaves
+  // the domain. If that happens, this function is called both the cell (which
+  // we can ask for its new location) and the cell it was previously on. The
+  // function then keeps track of updating the number of particles lost in this
+  // time step, the total number of lost particles, and then estimates whether
+  // the particle left through the hole in the middle of the anode. We do so by
+  // first checking whether the cell it was in last had an $x$ coordinate to the
+  // left of the right boundary (located at $x=4$) and the particle now has a
+  // position to the right of the right boundary. If that is so, we compute a
+  // direction vector of its motion that is normalized so that the $x$ component
+  // of the direction vector is equal to $1$. With this direction vector, we can
+  // compute where it would have intersected the line $x=4$. If this intersect
+  // is between $0.5$ and $1.5$, then we claim that the particle left through
+  // the hole and increment a counter.
+  template <int dim>
+  void CathodeRaySimulator<dim>::track_lost_particle(
+    const typename Particles::ParticleIterator<dim> &        particle,
+    const typename Triangulation<dim>::active_cell_iterator &cell)
+  {
+    ++n_recently_lost_particles;
+    ++n_total_lost_particles;
+
+    const Point<dim> current_location              = particle->get_location();
+    const Point<dim> approximate_previous_location = cell->center();
+
+    if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
+      {
+        const Tensor<1, dim> direction =
+          (current_location - approximate_previous_location) /
+          (current_location[0] - approximate_previous_location[0]);
+
+        const double right_boundary_intercept =
+          approximate_previous_location[1] +
+          (4 - approximate_previous_location[0]) * direction[1];
+        if ((right_boundary_intercept > 0.5) &&
+            (right_boundary_intercept < 1.5))
+          ++n_particles_lost_through_anode;
+      }
+  }
+
+
+
+  // @sect4{CathodeRaySimulator::update_timestep_size}
+
+  // As discussed at length in the introduction, we need to respect a CFL
+  // condition whereby particles can not move further than one cell in one time
+  // step. To ensure that this is the case, we again first compute the maximal
+  // speed of all particles on each cell, and divide the cell size by that
+  // speed. We then compute the next time step size as the minimum of this
+  // quantity over all cells, using the safety factor discussed in the
+  // introduction, and set this as the desired time step size using the
+  // DiscreteTime::set_desired_time_step_size() function.
+  template <int dim>
+  void CathodeRaySimulator<dim>::update_timestep_size()
+  {
+    if (time.get_step_number() > 0)
+      {
+        double min_cell_size_over_velocity = std::numeric_limits<double>::max();
+
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          if (particle_handler.n_particles_in_cell(cell) > 0)
+            {
+              const double cell_size = cell->minimum_vertex_distance();
+
+              double max_particle_velocity(0.0);
+
+              for (const auto &particle :
+                   particle_handler.particles_in_cell(cell))
+                {
+                  const Tensor<1, dim> velocity(particle.get_properties());
+                  max_particle_velocity =
+                    std::max(max_particle_velocity, velocity.norm());
+                }
+
+              if (max_particle_velocity > 0)
+                min_cell_size_over_velocity =
+                  std::min(min_cell_size_over_velocity,
+                           cell_size / max_particle_velocity);
+            }
+
+        constexpr double c_safety = 0.5;
+        time.set_desired_next_step_size(c_safety * 0.5 *
+                                        min_cell_size_over_velocity);
+      }
+    // As mentioned in the introduction, we have to treat the very first
+    // time step differently since there, particles are not available yet or
+    // do not yet have the information associated that we need for the
+    // computation of a reasonable step length. The formulas below follow the
+    // discussion in the introduction.
+    else
+      {
+        const QTrapez<dim> vertex_quadrature;
+        FEValues<dim>      fe_values(fe, vertex_quadrature, update_gradients);
+
+        std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
+
+        double min_timestep = std::numeric_limits<double>::max();
+
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          if (particle_handler.n_particles_in_cell(cell) > 0)
+            {
+              const double cell_size = cell->minimum_vertex_distance();
+
+              fe_values.reinit(cell);
+              fe_values.get_function_gradients(solution, field_gradients);
+
+              double max_E = 0;
+              for (const auto q_point : fe_values.quadrature_point_indices())
+                max_E = std::max(max_E, field_gradients[q_point].norm());
+
+              if (max_E > 0)
+                min_timestep =
+                  std::min(min_timestep,
+                           std::sqrt(0.5 * cell_size *
+                                     Constants::electron_mass /
+                                     Constants::electron_charge / max_E));
+            }
+
+        time.set_desired_next_step_size(min_timestep);
+      }
+  }
+
+
+
+  // @sect4{The <code>CathodeRaySimulator::output_results()</code> function}
+
+  // The final function implementing pieces of the overall algorithm is the one
+  // that generates graphical output. In the current context, we want to output
+  // both the electric potential field as well as the particle locations and
+  // velocities. But we also want to output the electric field, i.e., the
+  // gradient of the solution.
+  //
+  // deal.II has a general way how one can compute derived quantities from
+  // the solution and output those as well. Here, this is the electric
+  // field, but it could also be some other quantity -- say, the norm of the
+  // electric field, or in fact anything else one could want to compute from
+  // the solution $V_h(\mathbf x)$ or its derivatives. This general solution
+  // uses the DataPostprocessor class and, in cases like the one here where we
+  // want to output a quantity that represents a vector field, the
+  // DataPostprocessorVector class.
+  //
+  // Rather than try and explain how this class works, let us simply refer to
+  // the documentation of the DataPostprocessorVector class that has essentially
+  // this case as a well-documented example.
+  template <int dim>
+  class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
+  {
+  public:
+    ElectricFieldPostprocessor()
+      : DataPostprocessorVector<dim>("electric_field", update_gradients)
+    {}
+
+    virtual void evaluate_scalar_field(
+      const DataPostprocessorInputs::Scalar<dim> &input_data,
+      std::vector<Vector<double>> &computed_quantities) const override
+    {
+      AssertDimension(input_data.solution_gradients.size(),
+                      computed_quantities.size());
+
+      for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
+        {
+          AssertDimension(computed_quantities[p].size(), dim);
+          for (unsigned int d = 0; d < dim; ++d)
+            computed_quantities[p][d] = input_data.solution_gradients[p][d];
+        }
+    }
+  };
+
+
+
+  // With this, the `output_results()` function becomes relatively
+  // straightforward: We use the DataOut class as we have in almost every one of
+  // the previous tutorial programs to output the solution (the "electric
+  // potential") and we use the postprocessor defined above to also output its
+  // gradient (the "electric field"). This all is then written into a file in
+  // VTU format after also associating the current time and time step number
+  // with this file.
+  template <int dim>
+  void CathodeRaySimulator<dim>::output_results() const
+  {
+    {
+      ElectricFieldPostprocessor<dim> electric_field;
+      DataOut<dim>                    data_out;
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(solution, "electric_potential");
+      data_out.add_data_vector(solution, electric_field);
+      data_out.build_patches();
+
+      data_out.set_flags(
+        DataOutBase::VtkFlags(time.get_current_time(), time.get_step_number()));
+
+      std::ofstream output("solution-" +
+                           Utilities::int_to_string(time.get_step_number(), 4) +
+                           ".vtu");
+      data_out.write_vtu(output);
+    }
+
+    // Output the particle positions and properties is not more complicated. The
+    // Particles::DataOut class plays the role of the DataOut class for
+    // particles, and all we have to do is tell that class where to take
+    // particles from and how to interpret the `dim` components of the
+    // properties -- namely, as a single vector indicating the velocity, rather
+    // than as `dim` scalar properties. The rest is then the same as above:
+    {
+      Particles::DataOut<dim, dim> particle_out;
+      particle_out.build_patches(
+        particle_handler,
+        std::vector<std::string>(dim, "velocity"),
+        std::vector<DataComponentInterpretation::DataComponentInterpretation>(
+          dim, DataComponentInterpretation::component_is_part_of_vector));
+
+      particle_out.set_flags(
+        DataOutBase::VtkFlags(time.get_current_time(), time.get_step_number()));
+
+      std::ofstream output("particles-" +
+                           Utilities::int_to_string(time.get_step_number(), 4) +
+                           ".vtu");
+      particle_out.write_vtu(output);
+    }
+  }
+
+
+  // @sect4{CathodeRaySimulator::run}
+
+  // The last member function of the principal class of this program is then the
+  // driver. At the top, it refines the mesh a number of times by solving the
+  // problem (with not particles yet created) on a sequence of finer and finer
+  // meshes.
+  template <int dim>
+  void CathodeRaySimulator<dim>::run()
+  {
+    make_grid();
+
+    // do a few refinement cycles up front
+    const unsigned int n_pre_refinement_cycles = 3;
+    for (unsigned int refinement_cycle = 0;
+         refinement_cycle < n_pre_refinement_cycles;
+         ++refinement_cycle)
+      {
+        setup_system();
+        assemble_system();
+        solve_field();
+        refine_grid();
+      }
+
+
+    // Now do the loop over time. The sequence of steps follows closely the
+    // outline of the algorithm discussed in the introduction. As discussed in
+    // great detail in the documentation of the DiscreteTime class, while we
+    // move the field and particle information forward by one time step, the
+    // time stored in the `time` variable is not consistent with where (some of)
+    // these quantities are (in the diction of DiscreteTime, this is the "update
+    // stage"). The call to `time.advance_time()` makes everything consistent
+    // again by setting the `time` variable to the time at which the field and
+    // particles already are, and once we are in this "consistent stage", we can
+    // generate graphical output and write information about the current state
+    // of the simulation to screen.
+    setup_system();
+    do
+      {
+        std::cout << "Timestep " << time.get_step_number() + 1 << std::endl;
+        std::cout << "  Field degrees of freedom:                 "
+                  << dof_handler.n_dofs() << std::endl;
+
+        assemble_system();
+        solve_field();
+
+        create_particles();
+        std::cout << "  Total number of particles in simulation:  "
+                  << particle_handler.n_global_particles() << std::endl;
+
+        n_recently_lost_particles = 0;
+        update_timestep_size();
+        move_particles();
+
+        time.advance_time();
+
+        output_results();
+
+        std::cout << "  Number of particles lost this time step:  "
+                  << n_recently_lost_particles << std::endl;
+        if (n_total_lost_particles > 0)
+          std::cout << "  Fraction of particles lost through anode: "
+                    << 1. * n_particles_lost_through_anode /
+                         n_total_lost_particles
+                    << std::endl;
+
+        std::cout << std::endl
+                  << "  Now at t=" << time.get_current_time()
+                  << ", dt=" << time.get_previous_step_size() << '.'
+                  << std::endl
+                  << std::endl;
+      }
+    while (time.is_at_end() == false);
+  }
+} // namespace Step19
+
+
+
+// @sect3{The <code>main</code> function}
+
+// The final function of the program is then again the `main()` function. It is
+// unchanged in all tutorial programs since step-6 and so there is nothing new
+// to discuss:
+int main()
+{
+  try
+    {
+      Step19::CathodeRaySimulator<2> cathode_ray_simulator_2d;
+      cathode_ray_simulator_2d.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.