Reference documentation for deal.II version GIT relicensing-224-gc660c0d696 2024-03-28 18:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-19.h
Go to the documentation of this file.
1dim)
765 *   , next_unused_particle_id(0)
766 *   , n_recently_lost_particles(0)
767 *   , n_total_lost_particles(0)
768 *   , n_particles_lost_through_anode(0)
769 *   , time(0, 1e-4)
770 *   {
771 *   particle_handler.signals.particle_lost.connect(
772 *   [this](const typename Particles::ParticleIterator<dim> &particle,
773 *   const typename Triangulation<dim>::active_cell_iterator &cell) {
774 *   this->track_lost_particle(particle, cell);
775 *   });
776 *   }
777 *  
778 *  
779 *  
780 * @endcode
781 *
782 *
783 * <a name="step_19-ThecodeCathodeRaySimulatormake_gridcodefunction"></a>
784 * <h4>The <code>CathodeRaySimulator::make_grid</code> function</h4>
785 *
786
787 *
788 * The next function is then responsible for generating the mesh on which
789 * we want to solve. Recall how the domain looks like:
790 * <p align="center">
791 * <img
792 * src="https://www.dealii.org/images/steps/developer/step-19.geometry.png"
793 * alt="The geometry used in this program"
794 * width="600">
795 * </p>
796 * We subdivide this geometry into a mesh of @f$4\times 2@f$ cells that looks
797 * like this:
798 * <div class=CodeFragmentInTutorialComment>
799 * @code
800 * *---*---*---*---*
801 * \ | | | |
802 * *--*---*---*---*
803 * / | | | |
804 * *---*---*---*---*
805 * @endcode
806 * </div>
807 * The way this is done is by first defining where the @f$15=5\times 3@f$
808 * vertices are located -- here, we say that they are on integer points
809 * with the middle one on the left side moved to the right by a value of
810 * `delta=0.5`.
811 *
812
813 *
814 * In the following, we then have to say which vertices together form
815 * the 8 cells. The following code is then entirely equivalent to what
816 * we also do in @ref step_14 "step-14":
817 *
818 * @code
819 *   template <int dim>
820 *   void CathodeRaySimulator<dim>::make_grid()
821 *   {
822 *   static_assert(dim == 2,
823 *   "This function is currently only implemented for 2d.");
824 *  
825 *   const double delta = 0.5;
826 *   const unsigned int nx = 5;
827 *   const unsigned int ny = 3;
828 *  
829 *   const std::vector<Point<dim>> vertices
830 *   = {{0, 0},
831 *   {1, 0},
832 *   {2, 0},
833 *   {3, 0},
834 *   {4, 0},
835 *   {delta, 1},
836 *   {1, 1},
837 *   {2, 1},
838 *   {3, 1},
839 *   {4, 1},
840 *   {0, 2},
841 *   {1, 2},
842 *   {2, 2},
843 *   {3, 2},
844 *   {4, 2}};
845 *   AssertDimension(vertices.size(), nx * ny);
846 *  
847 *   const std::vector<unsigned int> cell_vertices[(nx - 1) * (ny - 1)] = {
848 *   {0, 1, nx + 0, nx + 1},
849 *   {1, 2, nx + 1, nx + 2},
850 *   {2, 3, nx + 2, nx + 3},
851 *   {3, 4, nx + 3, nx + 4},
852 *  
853 *   {5, nx + 1, 2 * nx + 0, 2 * nx + 1},
854 *   {nx + 1, nx + 2, 2 * nx + 1, 2 * nx + 2},
855 *   {nx + 2, nx + 3, 2 * nx + 2, 2 * nx + 3},
856 *   {nx + 3, nx + 4, 2 * nx + 3, 2 * nx + 4}};
857 *  
858 * @endcode
859 *
860 * With these arrays out of the way, we can move to slightly higher
861 * higher-level data structures. We create a vector of CellData
862 * objects that store for each cell to be created the vertices in
863 * question as well as the @ref GlossMaterialId "material id" (which
864 * we will here simply set to zero since we don't use it in the program).
865 *
866
867 *
868 * This information is then handed to the
869 * Triangulation::create_triangulation() function, and the mesh is twice
870 * globally refined. As discussed in the corresponding place in @ref step_14 "step-14",
871 * the inputs to Triangulation::create_triangulation() need to be
872 * consistently oriented, which a function in namespace GridTools
873 * does for us.
874 *
875 * @code
876 *   std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
877 *   for (unsigned int i = 0; i < cells.size(); ++i)
878 *   {
879 *   cells[i].vertices = cell_vertices[i];
880 *   cells[i].material_id = 0;
881 *   }
882 *  
883 *   GridTools::consistently_order_cells(cells);
884 *   triangulation.create_triangulation(
885 *   vertices,
886 *   cells,
887 *   SubCellData()); // No boundary information
888 *  
889 *   triangulation.refine_global(2);
890 *  
891 * @endcode
892 *
893 * The remaining part of the function loops over all cells and their faces,
894 * and if a face is at the boundary determines which boundary indicator
895 * should be applied to it. The various conditions should make sense if
896 * you compare the code with the picture of the geometry above.
897 *
898
899 *
900 * Once done with this step, we refine the mesh once more globally.
901 *
902 * @code
903 *   for (auto &cell : triangulation.active_cell_iterators())
904 *   for (auto &face : cell->face_iterators())
905 *   if (face->at_boundary())
906 *   {
907 *   if ((face->center()[0] > 0) && (face->center()[0] < 0.5) &&
908 *   (face->center()[1] > 0) && (face->center()[1] < 2))
909 *   face->set_boundary_id(BoundaryIds::cathode);
910 *   else if ((face->center()[0] > 0) && (face->center()[0] < 2))
911 *   face->set_boundary_id(BoundaryIds::focus_element);
912 *   else if ((face->center()[0] > 4 - 1e-12) &&
913 *   ((face->center()[1] > 1.5) || (face->center()[1] < 0.5)))
914 *   face->set_boundary_id(BoundaryIds::anode);
915 *   else
916 *   face->set_boundary_id(BoundaryIds::open);
917 *   }
918 *  
919 *   triangulation.refine_global(1);
920 *   }
921 *  
922 *  
923 * @endcode
924 *
925 *
926 * <a name="step_19-ThecodeCathodeRaySimulatorsetup_systemcodefunction"></a>
927 * <h4>The <code>CathodeRaySimulator::setup_system</code> function</h4>
928 *
929
930 *
931 * The next function in this program deals with setting up the various
932 * objects related to solving the partial differential equations. It is
933 * in essence a copy of the corresponding function in @ref step_6 "step-6" and requires
934 * no further discussion.
935 *
936 * @code
937 *   template <int dim>
938 *   void CathodeRaySimulator<dim>::setup_system()
939 *   {
940 *   dof_handler.distribute_dofs(fe);
941 *  
942 *   solution.reinit(dof_handler.n_dofs());
943 *   system_rhs.reinit(dof_handler.n_dofs());
944 *  
945 *   constraints.clear();
946 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
947 *  
948 *   VectorTools::interpolate_boundary_values(dof_handler,
949 *   BoundaryIds::cathode,
950 *   Functions::ConstantFunction<dim>(
951 *   -Constants::V0),
952 *   constraints);
953 *   VectorTools::interpolate_boundary_values(dof_handler,
954 *   BoundaryIds::focus_element,
955 *   Functions::ConstantFunction<dim>(
956 *   -Constants::V0),
957 *   constraints);
958 *   VectorTools::interpolate_boundary_values(dof_handler,
959 *   BoundaryIds::anode,
960 *   Functions::ConstantFunction<dim>(
961 *   +Constants::V0),
962 *   constraints);
963 *   constraints.close();
964 *  
965 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
966 *   DoFTools::make_sparsity_pattern(dof_handler,
967 *   dsp,
968 *   constraints,
969 *   /*keep_constrained_dofs = */ false);
970 *   sparsity_pattern.copy_from(dsp);
971 *  
972 *   system_matrix.reinit(sparsity_pattern);
973 *   }
974 *  
975 *  
976 * @endcode
977 *
978 *
979 * <a name="step_19-ThecodeCathodeRaySimulatorassemble_systemcodefunction"></a>
980 * <h4>The <code>CathodeRaySimulator::assemble_system</code> function</h4>
981 *
982
983 *
984 * The function that computes
985 * the matrix entries is again in essence a copy of the
986 * corresponding function in @ref step_6 "step-6":
987 *
988 * @code
989 *   template <int dim>
990 *   void CathodeRaySimulator<dim>::assemble_system()
991 *   {
992 *   system_matrix = 0;
993 *   system_rhs = 0;
994 *  
995 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
996 *  
997 *   FEValues<dim> fe_values(fe,
998 *   quadrature_formula,
999 *   update_values | update_gradients |
1000 *   update_quadrature_points | update_JxW_values);
1001 *  
1002 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1003 *  
1004 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1005 *   Vector<double> cell_rhs(dofs_per_cell);
1006 *  
1007 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1008 *  
1009 *   for (const auto &cell : dof_handler.active_cell_iterators())
1010 *   {
1011 *   cell_matrix = 0;
1012 *   cell_rhs = 0;
1013 *  
1014 *   fe_values.reinit(cell);
1015 *  
1016 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
1017 *   for (const unsigned int i : fe_values.dof_indices())
1018 *   {
1019 *   for (const unsigned int j : fe_values.dof_indices())
1020 *   cell_matrix(i, j) +=
1021 *   (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
1022 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
1023 *   fe_values.JxW(q_index)); // dx
1024 *   }
1025 *  
1026 * @endcode
1027 *
1028 * The only interesting part of this function is how it forms the right
1029 * hand side of the linear system. Recall that the right hand side
1030 * of the PDE is
1031 * @f[
1032 * \sum_p (N e)\delta(\mathbf x-\mathbf x_p),
1033 * @f]
1034 * where we have used @f$p@f$ to index the particles here to avoid
1035 * confusion with the shape function @f$\varphi_i@f$; @f$\mathbf x_p@f$
1036 * is the position of the @f$p@f$th particle.
1037 *
1038
1039 *
1040 * When multiplied by a test function @f$\varphi_i@f$ and integrated over
1041 * the domain results in a right hand side vector
1042 * @f{align*}{
1043 * F_i &= \int_\Omega \varphi_i (\mathbf x)\left[
1044 * \sum_p (N e)\delta(\mathbf x-\mathbf x_p) \right] dx
1045 * \\ &= \sum_p (N e) \varphi_i(\mathbf x_p).
1046 * @f}
1047 * Note that the final line no longer contains an integral, and
1048 * consequently also no occurrence of @f$dx@f$ which would require the
1049 * appearance of the `JxW` symbol in our code.
1050 *
1051
1052 *
1053 * For a given cell @f$K@f$, this cell's contribution to the right hand
1054 * side is then
1055 * @f{align*}{
1056 * F_i^K &= \sum_{p, \mathbf x_p\in K} (N e) \varphi_i(\mathbf x_p),
1057 * @f}
1058 * i.e., we only have to worry about those particles that are actually
1059 * located on the current cell @f$K@f$.
1060 *
1061
1062 *
1063 * In practice, what we do here is the following: If there are any
1064 * particles on the current cell, then we first obtain an iterator range
1065 * pointing to the first particle of that cell as well as the particle
1066 * past the last one on this cell (or the end iterator) -- i.e., a
1067 * half-open range as is common for C++ functions. Knowing now the list
1068 * of particles, we query their reference locations (with respect to
1069 * the reference cell), evaluate the shape functions in these reference
1070 * locations, and compute the force according to the formula above
1071 * (without any FEValues::JxW).
1072 *
1073
1074 *
1075 * @note It is worth pointing out that calling the
1077 * Particles::ParticleHandler::n_particles_in_cell() functions is not
1078 * very efficient on problems with a large number of particles. But it
1079 * illustrates the easiest way to write this algorithm, and so we are
1080 * willing to incur this cost for the moment for expository purposes.
1081 * We discuss the issue in more detail in the
1082 * <a href="#extensions">"possibilities for extensions" section</a>
1083 * below, and use a better approach in @ref step_70 "step-70", for example.
1084 *
1085 * @code
1086 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1087 *   for (const auto &particle : particle_handler.particles_in_cell(cell))
1088 *   {
1089 *   const Point<dim> &reference_location =
1090 *   particle.get_reference_location();
1091 *   for (const unsigned int i : fe_values.dof_indices())
1092 *   cell_rhs(i) +=
1093 *   (fe.shape_value(i, reference_location) * // phi_i(x_p)
1094 *   (-Constants::electrons_per_particle * // N
1095 *   Constants::electron_charge)); // e
1096 *   }
1097 *  
1098 * @endcode
1099 *
1100 * Finally, we can copy the contributions of this cell into
1101 * the global matrix and right hand side vector:
1102 *
1103 * @code
1104 *   cell->get_dof_indices(local_dof_indices);
1105 *   constraints.distribute_local_to_global(
1106 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1107 *   }
1108 *   }
1109 *  
1110 *  
1111 * @endcode
1112 *
1113 *
1114 * <a name="step_19-CathodeRaySimulatorsolve"></a>
1115 * <h4>CathodeRaySimulator::solve</h4>
1116 *
1117
1118 *
1119 * The function that solves the linear system is then again exactly as in
1120 * @ref step_6 "step-6":
1121 *
1122 * @code
1123 *   template <int dim>
1124 *   void CathodeRaySimulator<dim>::solve_field()
1125 *   {
1126 *   SolverControl solver_control(1000, 1e-12);
1127 *   SolverCG<Vector<double>> solver(solver_control);
1128 *  
1129 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1130 *   preconditioner.initialize(system_matrix, 1.2);
1131 *  
1132 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
1133 *  
1134 *   constraints.distribute(solution);
1135 *   }
1136 *  
1137 *  
1138 * @endcode
1139 *
1140 *
1141 * <a name="step_19-CathodeRaySimulatorrefine_grid"></a>
1142 * <h4>CathodeRaySimulator::refine_grid</h4>
1143 *
1144
1145 *
1146 * The final field-related function is the one that refines the grid. We will
1147 * call it a number of times in the first time step to obtain a mesh that
1148 * is well-adapted to the structure of the solution and, in particular,
1149 * resolves the various singularities in the solution that are due to
1150 * re-entrant corners and places where the boundary condition type
1151 * changes. You might want to refer to @ref step_6 "step-6" again for more details:
1152 *
1153 * @code
1154 *   template <int dim>
1155 *   void CathodeRaySimulator<dim>::refine_grid()
1156 *   {
1157 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1158 *  
1159 *   KellyErrorEstimator<dim>::estimate(dof_handler,
1160 *   QGauss<dim - 1>(fe.degree + 1),
1161 *   {},
1162 *   solution,
1163 *   estimated_error_per_cell);
1164 *  
1165 *   GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1166 *   estimated_error_per_cell,
1167 *   0.1,
1168 *   0.03);
1169 *  
1170 *   triangulation.execute_coarsening_and_refinement();
1171 *   }
1172 *  
1173 *  
1174 * @endcode
1175 *
1176 *
1177 * <a name="step_19-CathodeRaySimulatorcreate_particles"></a>
1178 * <h4>CathodeRaySimulator::create_particles</h4>
1179 *
1180
1181 *
1182 * Let us now turn to the functions that deal with particles. The first one
1183 * is about the creation of particles. As mentioned in the introduction,
1184 * we want to create a particle at points of the cathode if the electric
1185 * field @f$\mathbf E=\nabla V@f$ exceeds a certain threshold, i.e., if
1186 * @f$|\mathbf E| \ge E_\text{threshold}@f$, and if furthermore the electric field
1187 * points into the domain (i.e., if @f$\mathbf E \cdot \mathbf n < 0@f$). As is
1188 * common in the finite element method, we evaluate fields (and their
1189 * derivatives) at specific evaluation points; typically, these are
1190 * "quadrature points", and so we create a "quadrature formula" that we will
1191 * use to designate the points at which we want to evaluate the solution.
1192 * Here, we will simply take QMidpoint implying that we will only check the
1193 * threshold condition at the midpoints of faces. We then use this to
1194 * initialize an object of type FEFaceValues to evaluate the solution at these
1195 * points.
1196 *
1197
1198 *
1199 * All of this will then be used in a loop over all cells, their faces, and
1200 * specifically those faces that are at the boundary and, moreover, the
1201 * cathode part of the boundary.
1202 *
1203 * @code
1204 *   template <int dim>
1205 *   void CathodeRaySimulator<dim>::create_particles()
1206 *   {
1207 *   FEFaceValues<dim> fe_face_values(fe,
1208 *   QMidpoint<dim - 1>(),
1209 *   update_quadrature_points |
1210 *   update_gradients |
1211 *   update_normal_vectors);
1212 *  
1213 *   std::vector<Tensor<1, dim>> solution_gradients(
1214 *   fe_face_values.n_quadrature_points);
1215 *  
1216 *   for (const auto &cell : dof_handler.active_cell_iterators())
1217 *   for (const auto &face : cell->face_iterators())
1218 *   if (face->at_boundary() &&
1219 *   (face->boundary_id() == BoundaryIds::cathode))
1220 *   {
1221 *   fe_face_values.reinit(cell, face);
1222 *  
1223 * @endcode
1224 *
1225 * So we have found a face on the cathode. Next, we let the
1226 * FEFaceValues object compute the gradient of the solution at each
1227 * "quadrature" point, and extract the electric field vector from
1228 * the gradient in the form of a Tensor variable through the methods
1229 * discussed in the
1230 * @ref vector_valued "vector-valued problems" documentation module.
1231 *
1232 * @code
1233 *   const FEValuesExtractors::Scalar electric_potential(0);
1234 *   fe_face_values[electric_potential].get_function_gradients(
1235 *   solution, solution_gradients);
1236 *   for (const unsigned int q_point :
1237 *   fe_face_values.quadrature_point_indices())
1238 *   {
1239 *   const Tensor<1, dim> E = solution_gradients[q_point];
1240 *  
1241 * @endcode
1242 *
1243 * Electrons can only escape the cathode if the electric field
1244 * strength exceeds a threshold and,
1245 * crucially, if the electric field points *into* the domain.
1246 * Once we have that checked, we create a new
1247 * Particles::Particle object at this location and insert it
1248 * into the Particles::ParticleHandler object with a unique ID.
1249 *
1250
1251 *
1252 * The only thing that may be not obvious here is that we also
1253 * associate with this particle the location in the reference
1254 * coordinates of the cell we are currently on. This is done
1255 * because we will in downstream functions compute quantities
1256 * such as the electric field at the location of the particle
1257 * (e.g., to compute the forces that act on it when updating its
1258 * position in each time step). Evaluating a finite element
1259 * field at arbitrary coordinates is quite an expensive
1260 * operation because shape functions are really only defined on
1261 * the reference cell, and so when asking for the electric field
1262 * at an arbitrary point requires us first to determine what the
1263 * reference coordinates of that point are. To avoid having to
1264 * do this over and over, we determine these coordinates once
1265 * and for all and then store these reference coordinates
1266 * directly with the particle.
1267 *
1268 * @code
1269 *   if ((E * fe_face_values.normal_vector(q_point) < 0) &&
1270 *   (E.norm() > Constants::E_threshold))
1271 *   {
1272 *   const Point<dim> &location =
1273 *   fe_face_values.quadrature_point(q_point);
1274 *  
1275 *   Particles::Particle<dim> new_particle;
1276 *   new_particle.set_location(location);
1277 *   new_particle.set_reference_location(
1278 *   mapping.transform_real_to_unit_cell(cell, location));
1279 *   new_particle.set_id(next_unused_particle_id);
1280 *   particle_handler.insert_particle(new_particle, cell);
1281 *  
1282 *   ++next_unused_particle_id;
1283 *   }
1284 *   }
1285 *   }
1286 *  
1287 * @endcode
1288 *
1289 * At the end of all of these insertions, we let the `particle_handler`
1290 * update some internal statistics about the particles it stores.
1291 *
1292 * @code
1293 *   particle_handler.update_cached_numbers();
1294 *   }
1295 *  
1296 *  
1297 * @endcode
1298 *
1299 *
1300 * <a name="step_19-CathodeRaySimulatormove_particles"></a>
1301 * <h4>CathodeRaySimulator::move_particles</h4>
1302 *
1303
1304 *
1305 * The second particle-related function is the one that moves the particles
1306 * in each time step. To do this, we have to loop over all cells, the
1307 * particles in each cell, and evaluate the electric field at each of the
1308 * particles' positions.
1309 *
1310
1311 *
1312 * The approach used here is conceptually the same used in the
1313 * `assemble_system()` function: We loop over all cells, find the particles
1314 * located there (with the same caveat about the inefficiency of the algorithm
1315 * used here to find these particles), and use FEPointEvaluation object to
1316 * evaluate the gradient at these positions:
1317 *
1318 * @code
1319 *   template <int dim>
1320 *   void CathodeRaySimulator<dim>::move_particles()
1321 *   {
1322 *   const double dt = time.get_next_step_size();
1323 *  
1324 *   Vector<double> solution_values(fe.n_dofs_per_cell());
1325 *   FEPointEvaluation<1, dim> evaluator(mapping, fe, update_gradients);
1326 *  
1327 *   for (const auto &cell : dof_handler.active_cell_iterators())
1328 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1329 *   {
1330 *   const typename Particles::ParticleHandler<
1331 *   dim>::particle_iterator_range particles_in_cell =
1332 *   particle_handler.particles_in_cell(cell);
1333 *  
1334 *   std::vector<Point<dim>> particle_positions;
1335 *   for (const auto &particle : particles_in_cell)
1336 *   particle_positions.push_back(particle.get_reference_location());
1337 *  
1338 *   cell->get_dof_values(solution, solution_values);
1339 *  
1340 * @endcode
1341 *
1342 * Then we can ask the FEPointEvaluation object for the gradients of
1343 * the solution (i.e., the electric field @f$\mathbf E@f$) at these
1344 * locations and loop over the individual particles:
1345 *
1346 * @code
1347 *   evaluator.reinit(cell, particle_positions);
1348 *   evaluator.evaluate(make_array_view(solution_values),
1349 *   EvaluationFlags::gradients);
1350 *  
1351 *   {
1352 *   typename Particles::ParticleHandler<dim>::particle_iterator
1353 *   particle = particles_in_cell.begin();
1354 *   for (unsigned int particle_index = 0;
1355 *   particle != particles_in_cell.end();
1356 *   ++particle, ++particle_index)
1357 *   {
1358 *   const Tensor<1, dim> &E =
1359 *   evaluator.get_gradient(particle_index);
1360 *  
1361 * @endcode
1362 *
1363 * Having now obtained the electric field at the location of one
1364 * of the particles, we use this to update first the velocity
1365 * and then the position. To do so, let us first get the old
1366 * velocity out of the properties stored with the particle,
1367 * compute the acceleration, update the velocity, and store this
1368 * new velocity again in the properties of the particle. Recall
1369 * that this corresponds to the first of the following set of
1370 * update equations discussed in the introduction:
1371 * @f{align*}{
1372 * \frac{{\mathbf v}_i^{(n)}
1373 * -{\mathbf v}_i^{(n-1)}}{\Delta t}
1374 * &= \frac{e\nabla V^{(n)}}{m}
1375 * \\ \frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}
1376 * {\Delta t} &= {\mathbf v}_i^{(n)}.
1377 * @f}
1378 *
1379 * @code
1380 *   const Tensor<1, dim> old_velocity(particle->get_properties());
1381 *  
1382 *   const Tensor<1, dim> acceleration =
1383 *   Constants::electron_charge / Constants::electron_mass * E;
1384 *  
1385 *   const Tensor<1, dim> new_velocity =
1386 *   old_velocity + acceleration * dt;
1387 *  
1388 *   particle->set_properties(new_velocity);
1389 *  
1390 * @endcode
1391 *
1392 * With the new velocity, we can then also update the location
1393 * of the particle and tell the particle about it.
1394 *
1395 * @code
1396 *   const Point<dim> new_location =
1397 *   particle->get_location() + dt * new_velocity;
1398 *   particle->set_location(new_location);
1399 *   }
1400 *   }
1401 *   }
1402 *  
1403 * @endcode
1404 *
1405 * Having updated the locations and properties (i.e., velocities) of all
1406 * particles, we need to make sure that the `particle_handler` again knows
1407 * which cells they are in, and what their locations in the coordinate
1408 * system of the reference cell are. The following function does that. (It
1409 * also makes sure that, in parallel computations, particles are moved from
1410 * one processor to another processor if a particle moves from the subdomain
1411 * owned by the former to the subdomain owned by the latter.)
1412 *
1413 * @code
1414 *   particle_handler.sort_particles_into_subdomains_and_cells();
1415 *   }
1416 *  
1417 *  
1418 * @endcode
1419 *
1420 *
1421 * <a name="step_19-CathodeRaySimulatortrack_lost_particle"></a>
1422 * <h4>CathodeRaySimulator::track_lost_particle</h4>
1423 *
1424
1425 *
1426 * The final particle-related function is the one that is called whenever a
1427 * particle is lost from the simulation. This typically happens if it leaves
1428 * the domain. If that happens, this function is called both the cell (which
1429 * we can ask for its new location) and the cell it was previously on. The
1430 * function then keeps track of updating the number of particles lost in this
1431 * time step, the total number of lost particles, and then estimates whether
1432 * the particle left through the hole in the middle of the anode. We do so by
1433 * first checking whether the cell it was in last had an @f$x@f$ coordinate to the
1434 * left of the right boundary (located at @f$x=4@f$) and the particle now has a
1435 * position to the right of the right boundary. If that is so, we compute a
1436 * direction vector of its motion that is normalized so that the @f$x@f$ component
1437 * of the direction vector is equal to @f$1@f$. With this direction vector, we can
1438 * compute where it would have intersected the line @f$x=4@f$. If this intersect
1439 * is between @f$0.5@f$ and @f$1.5@f$, then we claim that the particle left through
1440 * the hole and increment a counter.
1441 *
1442 * @code
1443 *   template <int dim>
1444 *   void CathodeRaySimulator<dim>::track_lost_particle(
1445 *   const typename Particles::ParticleIterator<dim> &particle,
1446 *   const typename Triangulation<dim>::active_cell_iterator &cell)
1447 *   {
1448 *   ++n_recently_lost_particles;
1449 *   ++n_total_lost_particles;
1450 *  
1451 *   const Point<dim> current_location = particle->get_location();
1452 *   const Point<dim> approximate_previous_location = cell->center();
1453 *  
1454 *   if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
1455 *   {
1456 *   const Tensor<1, dim> direction =
1457 *   (current_location - approximate_previous_location) /
1458 *   (current_location[0] - approximate_previous_location[0]);
1459 *  
1460 *   const double right_boundary_intercept =
1461 *   approximate_previous_location[1] +
1462 *   (4 - approximate_previous_location[0]) * direction[1];
1463 *   if ((right_boundary_intercept > 0.5) &&
1464 *   (right_boundary_intercept < 1.5))
1465 *   ++n_particles_lost_through_anode;
1466 *   }
1467 *   }
1468 *  
1469 *  
1470 *  
1471 * @endcode
1472 *
1473 *
1474 * <a name="step_19-CathodeRaySimulatorupdate_timestep_size"></a>
1475 * <h4>CathodeRaySimulator::update_timestep_size</h4>
1476 *
1477
1478 *
1479 * As discussed at length in the introduction, we need to respect a time step
1480 * condition whereby particles can not move further than one cell in one time
1481 * step. To ensure that this is the case, we again first compute the maximal
1482 * speed of all particles on each cell, and divide the cell size by that
1483 * speed. We then compute the next time step size as the minimum of this
1484 * quantity over all cells, using the safety factor discussed in the
1485 * introduction, and set this as the desired time step size using the
1486 * DiscreteTime::set_desired_time_step_size() function.
1487 *
1488 * @code
1489 *   template <int dim>
1490 *   void CathodeRaySimulator<dim>::update_timestep_size()
1491 *   {
1492 *   if (time.get_step_number() > 0)
1493 *   {
1494 *   double min_cell_size_over_velocity = std::numeric_limits<double>::max();
1495 *  
1496 *   for (const auto &cell : dof_handler.active_cell_iterators())
1497 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1498 *   {
1499 *   const double cell_size = cell->minimum_vertex_distance();
1500 *  
1501 *   double max_particle_velocity(0.0);
1502 *  
1503 *   for (const auto &particle :
1504 *   particle_handler.particles_in_cell(cell))
1505 *   {
1506 *   const Tensor<1, dim> velocity(particle.get_properties());
1507 *   max_particle_velocity =
1508 *   std::max(max_particle_velocity, velocity.norm());
1509 *   }
1510 *  
1511 *   if (max_particle_velocity > 0)
1512 *   min_cell_size_over_velocity =
1513 *   std::min(min_cell_size_over_velocity,
1514 *   cell_size / max_particle_velocity);
1515 *   }
1516 *  
1517 *   constexpr double c_safety = 0.5;
1518 *   time.set_desired_next_step_size(c_safety * 0.5 *
1519 *   min_cell_size_over_velocity);
1520 *   }
1521 * @endcode
1522 *
1523 * As mentioned in the introduction, we have to treat the very first
1524 * time step differently since there, particles are not available yet or
1525 * do not yet have the information associated that we need for the
1526 * computation of a reasonable step length. The formulas below follow the
1527 * discussion in the introduction.
1528 *
1529 * @code
1530 *   else
1531 *   {
1532 *   const QTrapezoid<dim> vertex_quadrature;
1533 *   FEValues<dim> fe_values(fe, vertex_quadrature, update_gradients);
1534 *  
1535 *   std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
1536 *  
1537 *   double min_timestep = std::numeric_limits<double>::max();
1538 *  
1539 *   for (const auto &cell : dof_handler.active_cell_iterators())
1540 *   if (particle_handler.n_particles_in_cell(cell) > 0)
1541 *   {
1542 *   const double cell_size = cell->minimum_vertex_distance();
1543 *  
1544 *   fe_values.reinit(cell);
1545 *   fe_values.get_function_gradients(solution, field_gradients);
1546 *  
1547 *   double max_E = 0;
1548 *   for (const auto q_point : fe_values.quadrature_point_indices())
1549 *   max_E = std::max(max_E, field_gradients[q_point].norm());
1550 *  
1551 *   if (max_E > 0)
1552 *   min_timestep =
1553 *   std::min(min_timestep,
1554 *   std::sqrt(0.5 * cell_size *
1555 *   Constants::electron_mass /
1556 *   Constants::electron_charge / max_E));
1557 *   }
1558 *  
1559 *   time.set_desired_next_step_size(min_timestep);
1560 *   }
1561 *   }
1562 *  
1563 *  
1564 *  
1565 * @endcode
1566 *
1567 *
1568 * <a name="step_19-ThecodeCathodeRaySimulatoroutput_resultscodefunction"></a>
1569 * <h4>The <code>CathodeRaySimulator::output_results()</code> function</h4>
1570 *
1571
1572 *
1573 * The final function implementing pieces of the overall algorithm is the one
1574 * that generates graphical output. In the current context, we want to output
1575 * both the electric potential field as well as the particle locations and
1576 * velocities. But we also want to output the electric field, i.e., the
1577 * gradient of the solution.
1578 *
1579
1580 *
1581 * deal.II has a general way how one can compute derived quantities from
1582 * the solution and output those as well. Here, this is the electric
1583 * field, but it could also be some other quantity -- say, the norm of the
1584 * electric field, or in fact anything else one could want to compute from
1585 * the solution @f$V_h(\mathbf x)@f$ or its derivatives. This general solution
1586 * uses the DataPostprocessor class and, in cases like the one here where we
1587 * want to output a quantity that represents a vector field, the
1588 * DataPostprocessorVector class.
1589 *
1590
1591 *
1592 * Rather than try and explain how this class works, let us simply refer to
1593 * the documentation of the DataPostprocessorVector class that has essentially
1594 * this case as a well-documented example.
1595 *
1596 * @code
1597 *   template <int dim>
1598 *   class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
1599 *   {
1600 *   public:
1601 *   ElectricFieldPostprocessor()
1602 *   : DataPostprocessorVector<dim>("electric_field", update_gradients)
1603 *   {}
1604 *  
1605 *   virtual void evaluate_scalar_field(
1606 *   const DataPostprocessorInputs::Scalar<dim> &input_data,
1607 *   std::vector<Vector<double>> &computed_quantities) const override
1608 *   {
1609 *   AssertDimension(input_data.solution_gradients.size(),
1610 *   computed_quantities.size());
1611 *  
1612 *   for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
1613 *   {
1614 *   AssertDimension(computed_quantities[p].size(), dim);
1615 *   for (unsigned int d = 0; d < dim; ++d)
1616 *   computed_quantities[p][d] = input_data.solution_gradients[p][d];
1617 *   }
1618 *   }
1619 *   };
1620 *  
1621 *  
1622 *  
1623 * @endcode
1624 *
1625 * With this, the `output_results()` function becomes relatively
1626 * straightforward: We use the DataOut class as we have in almost
1627 * every one of the previous tutorial programs to output the
1628 * solution (the "electric potential") and we use the postprocessor
1629 * defined above to also output its gradient (the "electric
1630 * field"). This all is then written into a file in VTU format after
1631 * also associating the current time and time step number with this
1632 * file, and providing physical units for the two fields to be
1633 * output.
1634 *
1635 * @code
1636 *   template <int dim>
1637 *   void CathodeRaySimulator<dim>::output_results() const
1638 *   {
1639 *   {
1640 *   ElectricFieldPostprocessor<dim> electric_field;
1641 *   DataOut<dim> data_out;
1642 *   data_out.attach_dof_handler(dof_handler);
1643 *   data_out.add_data_vector(solution, "electric_potential");
1644 *   data_out.add_data_vector(solution, electric_field);
1645 *   data_out.build_patches();
1646 *  
1647 *   DataOutBase::VtkFlags output_flags;
1648 *   output_flags.time = time.get_current_time();
1649 *   output_flags.cycle = time.get_step_number();
1650 *   output_flags.physical_units["electric_potential"] = "V";
1651 *   output_flags.physical_units["electric_field"] = "V/m";
1652 *  
1653 *   data_out.set_flags(output_flags);
1654 *  
1655 *   std::ofstream output("solution-" +
1656 *   Utilities::int_to_string(time.get_step_number(), 4) +
1657 *   ".vtu");
1658 *   data_out.write_vtu(output);
1659 *   }
1660 *  
1661 * @endcode
1662 *
1663 * Output the particle positions and properties is not more complicated. The
1664 * Particles::DataOut class plays the role of the DataOut class for
1665 * particles, and all we have to do is tell that class where to take
1666 * particles from and how to interpret the `dim` components of the
1667 * properties -- namely, as a single vector indicating the velocity, rather
1668 * than as `dim` scalar properties. The rest is then the same as above:
1669 *
1670 * @code
1671 *   {
1672 *   Particles::DataOut<dim> particle_out;
1673 *   particle_out.build_patches(
1674 *   particle_handler,
1675 *   std::vector<std::string>(dim, "velocity"),
1676 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>(
1677 *   dim, DataComponentInterpretation::component_is_part_of_vector));
1678 *  
1679 *   DataOutBase::VtkFlags output_flags;
1680 *   output_flags.time = time.get_current_time();
1681 *   output_flags.cycle = time.get_step_number();
1682 *   output_flags.physical_units["velocity"] = "m/s";
1683 *  
1684 *   particle_out.set_flags(output_flags);
1685 *  
1686 *   std::ofstream output("particles-" +
1687 *   Utilities::int_to_string(time.get_step_number(), 4) +
1688 *   ".vtu");
1689 *   particle_out.write_vtu(output);
1690 *   }
1691 *   }
1692 *  
1693 *  
1694 * @endcode
1695 *
1696 *
1697 * <a name="step_19-CathodeRaySimulatorrun"></a>
1698 * <h4>CathodeRaySimulator::run</h4>
1699 *
1700
1701 *
1702 * The last member function of the principal class of this program is then the
1703 * driver. At the top, it refines the mesh a number of times by solving the
1704 * problem (with no particles yet created) on a sequence of finer and finer
1705 * meshes.
1706 *
1707 * @code
1708 *   template <int dim>
1709 *   void CathodeRaySimulator<dim>::run()
1710 *   {
1711 *   make_grid();
1712 *  
1713 *   const unsigned int n_pre_refinement_cycles = 3;
1714 *   for (unsigned int refinement_cycle = 0;
1715 *   refinement_cycle < n_pre_refinement_cycles;
1716 *   ++refinement_cycle)
1717 *   {
1718 *   setup_system();
1719 *   assemble_system();
1720 *   solve_field();
1721 *   refine_grid();
1722 *   }
1723 *  
1724 *  
1725 * @endcode
1726 *
1727 * Now do the loop over time. The sequence of steps follows closely the
1728 * outline of the algorithm discussed in the introduction. As discussed in
1729 * great detail in the documentation of the DiscreteTime class, while we
1730 * move the field and particle information forward by one time step, the
1731 * time stored in the `time` variable is not consistent with where (some of)
1732 * these quantities are (in the diction of DiscreteTime, this is the "update
1733 * stage"). The call to `time.advance_time()` makes everything consistent
1734 * again by setting the `time` variable to the time at which the field and
1735 * particles already are, and once we are in this "consistent stage", we can
1736 * generate graphical output and write information about the current state
1737 * of the simulation to screen.
1738 *
1739 * @code
1740 *   setup_system();
1741 *   do
1742 *   {
1743 *   std::cout << "Timestep " << time.get_step_number() + 1 << std::endl;
1744 *   std::cout << " Field degrees of freedom: "
1745 *   << dof_handler.n_dofs() << std::endl;
1746 *  
1747 *   assemble_system();
1748 *   solve_field();
1749 *  
1750 *   create_particles();
1751 *   std::cout << " Total number of particles in simulation: "
1752 *   << particle_handler.n_global_particles() << std::endl;
1753 *  
1754 *   n_recently_lost_particles = 0;
1755 *   update_timestep_size();
1756 *   move_particles();
1757 *  
1758 *   time.advance_time();
1759 *  
1760 *   output_results();
1761 *  
1762 *   std::cout << " Number of particles lost this time step: "
1763 *   << n_recently_lost_particles << std::endl;
1764 *   if (n_total_lost_particles > 0)
1765 *   std::cout << " Fraction of particles lost through anode: "
1766 *   << 1. * n_particles_lost_through_anode /
1767 *   n_total_lost_particles
1768 *   << std::endl;
1769 *  
1770 *   std::cout << std::endl
1771 *   << " Now at t=" << time.get_current_time()
1772 *   << ", dt=" << time.get_previous_step_size() << '.'
1773 *   << std::endl
1774 *   << std::endl;
1775 *   }
1776 *   while (time.is_at_end() == false);
1777 *   }
1778 *   } // namespace Step19
1779 *  
1780 *  
1781 *  
1782 * @endcode
1783 *
1784 *
1785 * <a name="step_19-Thecodemaincodefunction"></a>
1786 * <h3>The <code>main</code> function</h3>
1787 *
1788
1789 *
1790 * The final function of the program is then again the `main()` function. It is
1791 * unchanged in all tutorial programs since @ref step_6 "step-6" and so there is nothing new
1792 * to discuss:
1793 *
1794 * @code
1795 *   int main()
1796 *   {
1797 *   try
1798 *   {
1799 *   Step19::CathodeRaySimulator<2> cathode_ray_simulator_2d;
1800 *   cathode_ray_simulator_2d.run();
1801 *   }
1802 *   catch (std::exception &exc)
1803 *   {
1804 *   std::cerr << std::endl
1805 *   << std::endl
1806 *   << "----------------------------------------------------"
1807 *   << std::endl;
1808 *   std::cerr << "Exception on processing: " << std::endl
1809 *   << exc.what() << std::endl
1810 *   << "Aborting!" << std::endl
1811 *   << "----------------------------------------------------"
1812 *   << std::endl;
1813 *  
1814 *   return 1;
1815 *   }
1816 *   catch (...)
1817 *   {
1818 *   std::cerr << std::endl
1819 *   << std::endl
1820 *   << "----------------------------------------------------"
1821 *   << std::endl;
1822 *   std::cerr << "Unknown exception!" << std::endl
1823 *   << "Aborting!" << std::endl
1824 *   << "----------------------------------------------------"
1825 *   << std::endl;
1826 *   return 1;
1827 *   }
1828 *   return 0;
1829 *   }
1830 * @endcode
1831<a name="step_19-Results"></a><h1>Results</h1>
1832
1833
1834When this program is run, it produces output that looks as follows:
1835```
1836Timestep 1
1837 Field degrees of freedom: 4989
1838 Total number of particles in simulation: 20
1839 Number of particles lost this time step: 0
1840
1841 Now at t=2.12647e-07, dt=2.12647e-07.
1842
1843Timestep 2
1844 Field degrees of freedom: 4989
1845 Total number of particles in simulation: 40
1846 Number of particles lost this time step: 0
1847
1848 Now at t=4.14362e-07, dt=2.01715e-07.
1849
1850Timestep 3
1851 Field degrees of freedom: 4989
1852 Total number of particles in simulation: 60
1853 Number of particles lost this time step: 0
1854
1855 Now at t=5.23066e-07, dt=1.08704e-07.
1856
1857Timestep 4
1858 Field degrees of freedom: 4989
1859 Total number of particles in simulation: 80
1860 Number of particles lost this time step: 0
1861
1862 Now at t=6.08431e-07, dt=8.53649e-08.
1863
1864
1865...
1866
1867
1868Timestep 1000
1869 Field degrees of freedom: 4989
1870 Total number of particles in simulation: 394
1871 Number of particles lost this time step: 2
1872 Fraction of particles lost through anode: 0.823188
1873
1874 Now at t=3.83822e-05, dt=4.18519e-08.
1875
1876Timestep 1001
1877 Field degrees of freedom: 4989
1878 Total number of particles in simulation: 392
1879 Number of particles lost this time step: 0
1880 Fraction of particles lost through anode: 0.823188
1881
1882 Now at t=3.84244e-05, dt=4.22178e-08.
1883
1884
1885...
1886
1887
1888Timestep 2478
1889 Field degrees of freedom: 4989
1890 Total number of particles in simulation: 364
1891 Number of particles lost this time step: 0
1892 Fraction of particles lost through anode: 0.867745
1893
1894 Now at t=9.99509e-05, dt=4.20455e-08.
1895
1896Timestep 2479
1897 Field degrees of freedom: 4989
1898 Total number of particles in simulation: 364
1899 Number of particles lost this time step: 2
1900 Fraction of particles lost through anode: 0.867871
1901
1902 Now at t=9.99929e-05, dt=4.20105e-08.
1903
1904Timestep 2480
1905 Field degrees of freedom: 4989
1906 Total number of particles in simulation: 362
1907 Number of particles lost this time step: 2
1908 Fraction of particles lost through anode: 0.867047
1909
1910 Now at t=0.0001, dt=7.11401e-09.
1911```
1912
1913Picking a random few time steps, we can visualize the solution in the
1914form of streamlines for the electric field and dots for the electrons
1915(the pictures were created with an earlier version of the program and
1916are no longer exact representations of the solution as currently
1917computed -- but they are conceptually still correct):
1918<div class="twocolumn" style="width: 80%">
1919 <div>
1920 <img src="https://www.dealii.org/images/steps/developer/step-19.solution.0000.png"
1921 alt="The solution at time step 0 (t=0 seconds)."
1922 width="500">
1923 <br>
1924 Solution at time step 0 (t=0 seconds).
1925 <br>
1926 </div>
1927 <div>
1928 <img src="https://www.dealii.org/images/steps/developer/step-19.solution.1400.png"
1929 alt="The solution at time step 1400 (t=0.000068 seconds)."
1930 width="500">
1931 <br>
1932 Solution at time step 1400 (t=0.000068 seconds).
1933 <br>
1934 </div>
1935 <div>
1936 <img src="https://www.dealii.org/images/steps/developer/step-19.solution.0700.png"
1937 alt="The solution at time step 700 (t=0.000035 seconds)."
1938 width="500">
1939 <br>
1940 Solution at time step 700 (t=0.000035 seconds).
1941 <br>
1942 </div>
1943 <div>
1944 <img src="https://www.dealii.org/images/steps/developer/step-19.solution.2092.png"
1945 alt="The solution at time step 2092 (t=0.0001 seconds)."
1946 width="500">
1947 <br>
1948 Solution at time step 2092 (t=0.0001 seconds).
1949 <br>
1950 </div>
1951</div>
1952
1953That said, a more appropriate way to visualize the results of this
1954program are by creating a video that shows how these electrons move, and how
1955the electric field changes in response to their motion (again showing
1956results of an earlier version of the program):
1957
1958@htmlonly
1959<p align="center">
1960 <iframe width="560" height="315" src="https://www.youtube.com/embed/HwUtE7xuteE"
1961 frameborder="0"
1962 allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
1963 allowfullscreen></iframe>
1964 </p>
1965@endhtmlonly
1966
1967What you can see here is how the "focus element" of the boundary with its negative
1968voltage repels the electrons and makes sure that they do not just fly away
1969perpendicular from the cathode (as they do in the initial part of their
1970trajectories). It also shows how the electric field lines
1971move around over time, in response to the charges flying by -- in other words,
1972the feedback the particles have on the electric field that itself drives the
1973motion of the electrons.
1974
1975The movie suggests that electrons move in "bunches" or "bursts". One element of
1976this appearance is an artifact of how the movie was created: Every frame of the
1977movie corresponds to one time step, but the time step length varies. More specifically,
1978the fastest particle moving through the smallest cell determines the length of the
1979time step (see the discussion in the introduction), and consequently time steps
1980are small whenever a (fast) particle moves through the small cells at the right
1981edge of the domain; time steps are longer again once the particle has left
1982the domain. This slowing-accelerating effect can easily be visualized by plotting
1983the time step length shown in the screen output.
1984
1985The second part of this is real, however: The simulation creates a large group
1986of particles in the beginning, and fewer after about the 300th time step. This
1987is probably because of the negative charge of the particles in the simulation:
1988They reduce the magnitude of the electric field at the (also negatively charged
1989electrode) and consequently reduce the number of points on the cathode at which
1990the magnitude exceeds the threshold necessary to draw an electron out of the
1991electrode.
1992
1993
1994<a name="step-19-extensions"></a>
1995<a name="step_19-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1996
1997
1998
1999<a name="step_19-Morestatisticsaboutelectrons"></a><h4> More statistics about electrons </h4>
2000
2001
2002At the end of the day, we are rarely interested in the *solution* of an equation,
2003but in numbers that can be *extracted* from it -- in other words, we want to
2004*postprocess* the solution (see the results section of @ref step_4 "step-4" for an extensive
2005discussion of this concept). Here, what one would likely be most interested in
2006assessing are some statistics about the particles.
2007
2008The program already computes the fraction of the electrons that leave the
2009domain through the hole in the anode. But there are other quantities one might be
2010interested in. For example, the average velocity of these particles. It would
2011not be very difficult to obtain each particle's velocity from its properties,
2012in the same way as we do in the `move_particles()` function, and compute
2013statistics from it.
2014
2015
2016<a name="step_19-Abettersynchronizedvisualization"></a><h4> A better-synchronized visualization </h4>
2017
2018
2019As discussed above, there is a varying time difference between different frames
2020of the video because we create output for every time step. A better way to
2021create movies would be to generate a new output file in fixed time intervals,
2022regardless of how many time steps lie between each such point.
2023
2024
2025<a name="step_19-Abettertimeintegrator"></a><h4> A better time integrator </h4>
2026
2027
2028The problem we are considering in this program is a coupled, multiphysics
2029problem. But the way we solve it is by first computing the (electric) potential
2030field and then updating first the particle velocity and finally the
2031particle locations. This is what is called an
2032"operator-splitting method", a concept we will investigate in more detail
2033in @ref step_58 "step-58".
2034
2035While it is awkward to think of a way to solve this problem that does not involve
2036splitting the problem into a PDE piece and a particles piece, one *can*
2037(and probably should!) think of a better way to update the particle
2038locations. Specifically, the equations we use to update the particle location
2039are
2040@f{align*}{
2041 \frac{{\mathbf v}_i^{(n)}-{\mathbf v}_i^{(n-1)}}{\Delta t} &= \frac{e\nabla V^{(n)}}{m}
2042 \\
2043 \frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}{\Delta t} &= {\mathbf v}_i^{(n)}.
2044@f}
2045This corresponds to a Lie splitting where we first update one variable
2046(the velocity) and then, using the already-updated value of the one variable
2047to update the other (the position). It is well understood that the Lie
2048splitting incurs a first order error accuracy in the time step (i.e.,
2049the error introduced by the splitting is @f${\cal O}(\Delta t)@f$) that we know we should
2050avoid because we can do better. Independently, of course, we incur an error
2051of the same rate a second time because we use an Euler-like scheme for each of
2052the two updates (we replace the time derivative by a simple finite difference
2053quotient containing the old and new value, and we evaluate the right hand side
2054only at the end time of the interval, at @f$t_n@f$), and a better scheme would also
2055use a better way to do this kind of update.
2056
2057Better strategies would replace the Lie splitting by something like a
2058[Strang splitting](https://en.wikipedia.org/wiki/Strang_splitting), and the
2059update of the particle position and velocity by a scheme such
2060as the
2061[leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration)
2062or more generally
2063[symplectic integrators](https://en.wikipedia.org/wiki/Symplectic_integrator)
2064such as the
2065[Verlet scheme](https://en.wikipedia.org/wiki/Verlet_integration).
2066
2067
2068<a name="step_19-Parallelization"></a><h4> Parallelization </h4>
2069
2070
2071In release mode, the program runs in about 3.5 minutes on one of the author's
2072laptops at the time of writing this. That's acceptable. But what if we wanted
2073to make the simulation three-dimensional? If we wanted to not use a maximum
2074of around 100 particles at any given time (as happens with the parameters
2075used here) but 100,000? If we needed a substantially finer mesh?
2076
2077In those cases, one would want to run the program not just on a single processor,
2078but in fact on as many as one has available. This requires parallelization
2079of both the PDE solution as well as over particles. In practice, while there
2080are substantial challenges to making this efficient and scale well, these
2081challenges are all addressed in deal.II itself, and have been demonstrated
2082on simulations running on 10,000 or more cores. For example, @ref step_40 "step-40" shows
2083how to parallelize the finite element part, and @ref step_70 "step-70" shows how one can
2084then also parallelize the particles part.
2085 *
2086 *
2087<a name="step_19-PlainProg"></a>
2088<h1> The plain program</h1>
2089@include "step-19.cc"
2090*/
double JxW(const unsigned int q_point) const
particle_iterator_range particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4614
Point< 2 > first
Definition grid_out.cc:4613
unsigned int level
Definition grid_out.cc:4616
__global__ void set(Number *val, const Number s, const size_type N)
#define AssertDimension(dim1, dim2)
const Event initial
Definition event.cc:64
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
spacedim const Point< spacedim > & p
Definition grid_tools.h:981
spacedim & mesh
Definition grid_tools.h:980
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
int(&) functions(const void *v1, const void *v2)
const InputIterator OutputIterator out
Definition parallel.h:167
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const InputIterator OutputIterator const Function & function
Definition parallel.h:168