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