763 * , next_unused_particle_id(0)
764 * , n_recently_lost_particles(0)
765 * , n_total_lost_particles(0)
766 * , n_particles_lost_through_anode(0)
769 * particle_handler.signals.particle_lost.connect(
772 * this->track_lost_particle(particle, cell);
781 * <a name=
"step_19-ThecodeCathodeRaySimulatormake_gridcodefunction"></a>
782 * <h4>The <code>CathodeRaySimulator::make_grid</code> function</h4>
786 * The next function is then responsible
for generating the mesh on which
787 * we want to solve. Recall how the domain looks like:
790 * src=
"https://www.dealii.org/images/steps/developer/step-19.geometry.png"
791 * alt=
"The geometry used in this program"
794 * We subdivide
this geometry into a mesh of @f$4\times 2@f$ cells that looks
796 * <div
class=CodeFragmentInTutorialComment>
805 * The way
this is done is by
first defining where the @f$15=5\times 3@f$
806 *
vertices are located -- here, we say that they are on integer points
807 * with the middle one on the left side moved to the right by a
value of
812 * In the following, we then have to say which
vertices together form
813 * the 8 cells. The following code is then entirely equivalent to what
814 * we also
do in @ref step_14
"step-14":
818 *
void CathodeRaySimulator<dim>::make_grid()
820 *
static_assert(dim == 2,
821 *
"This function is currently only implemented for 2d.");
823 *
const double delta = 0.5;
824 *
const unsigned int nx = 5;
825 *
const unsigned int ny = 3;
827 *
const std::vector<Point<dim>>
vertices
845 *
const std::vector<unsigned int> cell_vertices[(nx - 1) * (ny - 1)] = {
846 * {0, 1, nx + 0, nx + 1},
847 * {1, 2, nx + 1, nx + 2},
848 * {2, 3, nx + 2, nx + 3},
849 * {3, 4, nx + 3, nx + 4},
851 * {5, nx + 1, 2 * nx + 0, 2 * nx + 1},
852 * {nx + 1, nx + 2, 2 * nx + 1, 2 * nx + 2},
853 * {nx + 2, nx + 3, 2 * nx + 2, 2 * nx + 3},
854 * {nx + 3, nx + 4, 2 * nx + 3, 2 * nx + 4}};
858 * With these arrays out of the way, we can move to slightly higher
859 * higher-
level data structures. We create a vector of
CellData
860 * objects that store
for each cell to be created the
vertices in
861 * question as well as the @ref GlossMaterialId
"material id" (which
862 * we will here simply set to zero since we don
't use it in the program).
866 * This information is then handed to the
867 * Triangulation::create_triangulation() function, and the mesh is twice
868 * globally refined. As discussed in the corresponding place in @ref step_14 "step-14",
869 * the inputs to Triangulation::create_triangulation() need to be
870 * consistently oriented, which a function in namespace GridTools
874 * std::vector<CellData<dim>> cells((nx - 1) * (ny - 1), CellData<dim>());
875 * for (unsigned int i = 0; i < cells.size(); ++i)
877 * cells[i].vertices = cell_vertices[i];
878 * cells[i].material_id = 0;
881 * GridTools::consistently_order_cells(cells);
882 * triangulation.create_triangulation(
885 * SubCellData()); // No boundary information
887 * triangulation.refine_global(2);
891 * The remaining part of the function loops over all cells and their faces,
892 * and if a face is at the boundary determines which boundary indicator
893 * should be applied to it. The various conditions should make sense if
894 * you compare the code with the picture of the geometry above.
898 * Once done with this step, we refine the mesh once more globally.
901 * for (auto &cell : triangulation.active_cell_iterators())
902 * for (auto &face : cell->face_iterators())
903 * if (face->at_boundary())
905 * if ((face->center()[0] > 0) && (face->center()[0] < 0.5) &&
906 * (face->center()[1] > 0) && (face->center()[1] < 2))
907 * face->set_boundary_id(BoundaryIds::cathode);
908 * else if ((face->center()[0] > 0) && (face->center()[0] < 2))
909 * face->set_boundary_id(BoundaryIds::focus_element);
910 * else if ((face->center()[0] > 4 - 1e-12) &&
911 * ((face->center()[1] > 1.5) || (face->center()[1] < 0.5)))
912 * face->set_boundary_id(BoundaryIds::anode);
914 * face->set_boundary_id(BoundaryIds::open);
917 * triangulation.refine_global(1);
924 * <a name="step_19-ThecodeCathodeRaySimulatorsetup_systemcodefunction"></a>
925 * <h4>The <code>CathodeRaySimulator::setup_system</code> function</h4>
929 * The next function in this program deals with setting up the various
930 * objects related to solving the partial differential equations. It is
931 * in essence a copy of the corresponding function in @ref step_6 "step-6" and requires
932 * no further discussion.
936 * void CathodeRaySimulator<dim>::setup_system()
938 * dof_handler.distribute_dofs(fe);
940 * solution.reinit(dof_handler.n_dofs());
941 * system_rhs.reinit(dof_handler.n_dofs());
943 * constraints.clear();
944 * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
946 * VectorTools::interpolate_boundary_values(dof_handler,
947 * BoundaryIds::cathode,
948 * Functions::ConstantFunction<dim>(
951 * VectorTools::interpolate_boundary_values(dof_handler,
952 * BoundaryIds::focus_element,
953 * Functions::ConstantFunction<dim>(
956 * VectorTools::interpolate_boundary_values(dof_handler,
957 * BoundaryIds::anode,
958 * Functions::ConstantFunction<dim>(
961 * constraints.close();
963 * DynamicSparsityPattern dsp(dof_handler.n_dofs());
964 * DoFTools::make_sparsity_pattern(dof_handler,
967 * /*keep_constrained_dofs = */ false);
968 * sparsity_pattern.copy_from(dsp);
970 * system_matrix.reinit(sparsity_pattern);
977 * <a name="step_19-ThecodeCathodeRaySimulatorassemble_systemcodefunction"></a>
978 * <h4>The <code>CathodeRaySimulator::assemble_system</code> function</h4>
982 * The function that computes
983 * the matrix entries is again in essence a copy of the
984 * corresponding function in @ref step_6 "step-6":
988 * void CathodeRaySimulator<dim>::assemble_system()
993 * const QGauss<dim> quadrature_formula(fe.degree + 1);
995 * FEValues<dim> fe_values(fe,
996 * quadrature_formula,
997 * update_values | update_gradients |
998 * update_quadrature_points | update_JxW_values);
1000 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1002 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1003 * Vector<double> cell_rhs(dofs_per_cell);
1005 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1007 * for (const auto &cell : dof_handler.active_cell_iterators())
1012 * fe_values.reinit(cell);
1014 * for (const unsigned int q_index : fe_values.quadrature_point_indices())
1015 * for (const unsigned int i : fe_values.dof_indices())
1017 * for (const unsigned int j : fe_values.dof_indices())
1018 * cell_matrix(i, j) +=
1019 * (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
1020 * fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
1021 * fe_values.JxW(q_index)); // dx
1026 * The only interesting part of this function is how it forms the right
1027 * hand side of the linear system. Recall that the right hand side
1030 * \sum_p (N e)\delta(\mathbf x-\mathbf x_p),
1032 * where we have used @f$p@f$ to index the particles here to avoid
1033 * confusion with the shape function @f$\varphi_i@f$; @f$\mathbf x_p@f$
1034 * is the position of the @f$p@f$th particle.
1038 * When multiplied by a test function @f$\varphi_i@f$ and integrated over
1039 * the domain results in a right hand side vector
1041 * F_i &= \int_\Omega \varphi_i (\mathbf x)\left[
1042 * \sum_p (N e)\delta(\mathbf x-\mathbf x_p) \right] dx
1043 * \\ &= \sum_p (N e) \varphi_i(\mathbf x_p).
1045 * Note that the final line no longer contains an integral, and
1046 * consequently also no occurrence of @f$dx@f$ which would require the
1047 * appearance of the `JxW` symbol in our code.
1051 * For a given cell @f$K@f$, this cell's contribution to the right hand
1054 * F_i^
K &= \sum_{p, \mathbf x_p\in
K} (N
e) \varphi_i(\mathbf x_p),
1056 * i.e., we only have to worry about those particles that are actually
1057 * located on the current cell @f$K@f$.
1061 * In practice, what we
do here is the following: If there are any
1062 * particles on the current cell, then we
first obtain an iterator range
1063 * pointing to the
first particle of that cell as well as the particle
1064 * past the last one on
this cell (or the end iterator) -- i.e., a
1065 * half-open range as is common
for C++
functions. Knowing now the list
1066 * of particles, we query their reference locations (with respect to
1067 * the reference cell), evaluate the shape
functions in these reference
1068 * locations, and compute the force according to the formula above
1073 * @note It is worth pointing out that calling the
1075 *
Particles::ParticleHandler::n_particles_in_cell() functions is not
1076 * very efficient on problems with a large number of particles. But it
1077 * illustrates the easiest way to write this algorithm, and so we are
1078 * willing to incur this cost for the moment for expository purposes.
1079 * We discuss the issue in more detail in the
1080 * <a href="
#extensions">"possibilities for extensions" section</a>
1081 * below, and use a better approach in @ref step_70 "step-70", for example.
1084 * if (particle_handler.n_particles_in_cell(cell) > 0)
1085 * for (const auto &particle : particle_handler.particles_in_cell(cell))
1087 * const Point<dim> &reference_location =
1088 * particle.get_reference_location();
1089 * for (const unsigned int i : fe_values.dof_indices())
1091 * (fe.shape_value(i, reference_location) * // phi_i(x_p)
1092 * (-Constants::electrons_per_particle * // N
1093 * Constants::electron_charge)); // e
1098 * Finally, we can copy the contributions of this cell into
1099 * the global matrix and right hand side vector:
1102 * cell->get_dof_indices(local_dof_indices);
1103 * constraints.distribute_local_to_global(
1104 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1112 * <a name="step_19-CathodeRaySimulatorsolve"></a>
1113 * <h4>CathodeRaySimulator::solve</h4>
1117 * The function that solves the linear system is then again exactly as in
1118 * @ref step_6 "step-6":
1121 * template <int dim>
1122 * void CathodeRaySimulator<dim>::solve_field()
1124 * SolverControl solver_control(1000, 1e-12);
1125 * SolverCG<Vector<double>> solver(solver_control);
1127 * PreconditionSSOR<SparseMatrix<double>> preconditioner;
1128 * preconditioner.initialize(system_matrix, 1.2);
1130 * solver.solve(system_matrix, solution, system_rhs, preconditioner);
1132 * constraints.distribute(solution);
1139 * <a name="step_19-CathodeRaySimulatorrefine_grid"></a>
1140 * <h4>CathodeRaySimulator::refine_grid</h4>
1144 * The final field-related function is the one that refines the grid. We will
1145 * call it a number of times in the first time step to obtain a mesh that
1146 * is well-adapted to the structure of the solution and, in particular,
1147 * resolves the various singularities in the solution that are due to
1148 * re-entrant corners and places where the boundary condition type
1149 * changes. You might want to refer to @ref step_6 "step-6" again for more details:
1152 * template <int dim>
1153 * void CathodeRaySimulator<dim>::refine_grid()
1155 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1157 * KellyErrorEstimator<dim>::estimate(dof_handler,
1158 * QGauss<dim - 1>(fe.degree + 1),
1161 * estimated_error_per_cell);
1163 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
1164 * estimated_error_per_cell,
1168 * triangulation.execute_coarsening_and_refinement();
1175 * <a name="step_19-CathodeRaySimulatorcreate_particles"></a>
1176 * <h4>CathodeRaySimulator::create_particles</h4>
1180 * Let us now turn to the functions that deal with particles. The first one
1181 * is about the creation of particles. As mentioned in the introduction,
1182 * we want to create a particle at points of the cathode if the electric
1183 * field @f$\mathbf E=\nabla V@f$ exceeds a certain threshold, i.e., if
1184 * @f$|\mathbf E| \ge E_\text{threshold}@f$, and if furthermore the electric field
1185 * points into the domain (i.e., if @f$\mathbf E \cdot \mathbf n < 0@f$). As is
1186 * common in the finite element method, we evaluate fields (and their
1187 * derivatives) at specific evaluation points; typically, these are
1188 * "quadrature points", and so we create a "quadrature formula" that we will
1189 * use to designate the points at which we want to evaluate the solution.
1190 * Here, we will simply take QMidpoint implying that we will only check the
1191 * threshold condition at the midpoints of faces. We then use this to
1192 * initialize an object of type FEFaceValues to evaluate the solution at these
1197 * All of this will then be used in a loop over all cells, their faces, and
1198 * specifically those faces that are at the boundary and, moreover, the
1199 * cathode part of the boundary.
1202 * template <int dim>
1203 * void CathodeRaySimulator<dim>::create_particles()
1205 * FEFaceValues<dim> fe_face_values(fe,
1206 * QMidpoint<dim - 1>(),
1207 * update_quadrature_points |
1208 * update_gradients |
1209 * update_normal_vectors);
1211 * std::vector<Tensor<1, dim>> solution_gradients(
1212 * fe_face_values.n_quadrature_points);
1214 * for (const auto &cell : dof_handler.active_cell_iterators())
1215 * for (const auto &face : cell->face_iterators())
1216 * if (face->at_boundary() &&
1217 * (face->boundary_id() == BoundaryIds::cathode))
1219 * fe_face_values.reinit(cell, face);
1223 * So we have found a face on the cathode. Next, we let the
1224 * FEFaceValues object compute the gradient of the solution at each
1225 * "quadrature" point, and extract the electric field vector from
1226 * the gradient in the form of a Tensor variable through the methods
1228 * @ref vector_valued "vector-valued problems" documentation topic.
1231 * const FEValuesExtractors::Scalar electric_potential(0);
1232 * fe_face_values[electric_potential].get_function_gradients(
1233 * solution, solution_gradients);
1234 * for (const unsigned int q_point :
1235 * fe_face_values.quadrature_point_indices())
1237 * const Tensor<1, dim> E = solution_gradients[q_point];
1241 * Electrons can only escape the cathode if the electric field
1242 * strength exceeds a threshold and,
1243 * crucially, if the electric field points *into* the domain.
1244 * Once we have that checked, we create a new
1245 * Particles::Particle object at this location and insert it
1246 * into the Particles::ParticleHandler object with a unique ID.
1250 * The only thing that may be not obvious here is that we also
1251 * associate with this particle the location in the reference
1252 * coordinates of the cell we are currently on. This is done
1253 * because we will in downstream functions compute quantities
1254 * such as the electric field at the location of the particle
1255 * (e.g., to compute the forces that act on it when updating its
1256 * position in each time step). Evaluating a finite element
1257 * field at arbitrary coordinates is quite an expensive
1258 * operation because shape functions are really only defined on
1259 * the reference cell, and so when asking for the electric field
1260 * at an arbitrary point requires us first to determine what the
1261 * reference coordinates of that point are. To avoid having to
1262 * do this over and over, we determine these coordinates once
1263 * and for all and then store these reference coordinates
1264 * directly with the particle.
1267 * if ((E * fe_face_values.normal_vector(q_point) < 0) &&
1268 * (E.norm() > Constants::E_threshold))
1270 * const Point<dim> &location =
1271 * fe_face_values.quadrature_point(q_point);
1273 * Particles::Particle<dim> new_particle;
1274 * new_particle.set_location(location);
1275 * new_particle.set_reference_location(
1276 * mapping.transform_real_to_unit_cell(cell, location));
1277 * new_particle.set_id(next_unused_particle_id);
1278 * particle_handler.insert_particle(new_particle, cell);
1280 * ++next_unused_particle_id;
1287 * At the end of all of these insertions, we let the `particle_handler`
1288 * update some internal statistics about the particles it stores.
1291 * particle_handler.update_cached_numbers();
1298 * <a name="step_19-CathodeRaySimulatormove_particles"></a>
1299 * <h4>CathodeRaySimulator::move_particles</h4>
1303 * The second particle-related function is the one that moves the particles
1304 * in each time step. To do this, we have to loop over all cells, the
1305 * particles in each cell, and evaluate the electric field at each of the
1306 * particles' positions.
1310 * The approach used here is conceptually the same used in the
1311 * `assemble_system()` function: We loop over all cells, find the particles
1312 * located there (with the same caveat about the inefficiency of the algorithm
1313 * used here to find these particles), and use FEPointEvaluation object to
1314 * evaluate the gradient at these positions:
1317 * template <int dim>
1318 * void CathodeRaySimulator<dim>::move_particles()
1320 * const double dt = time.get_next_step_size();
1322 * Vector<double> solution_values(fe.n_dofs_per_cell());
1323 * FEPointEvaluation<1, dim> evaluator(mapping, fe, update_gradients);
1325 * for (const auto &cell : dof_handler.active_cell_iterators())
1326 * if (particle_handler.n_particles_in_cell(cell) > 0)
1328 * const typename Particles::ParticleHandler<
1329 * dim>::particle_iterator_range particles_in_cell =
1330 * particle_handler.particles_in_cell(cell);
1332 * std::vector<Point<dim>> particle_positions;
1333 * for (const auto &particle : particles_in_cell)
1334 * particle_positions.push_back(particle.get_reference_location());
1336 * cell->get_dof_values(solution, solution_values);
1340 * Then we can ask the FEPointEvaluation object for the gradients of
1341 * the solution (i.e., the electric field @f$\mathbf E@f$) at these
1342 * locations and loop over the individual particles:
1345 * evaluator.reinit(cell, particle_positions);
1346 * evaluator.evaluate(make_array_view(solution_values),
1347 * EvaluationFlags::gradients);
1350 * typename Particles::ParticleHandler<dim>::particle_iterator
1351 * particle = particles_in_cell.begin();
1352 * for (unsigned int particle_index = 0;
1353 * particle != particles_in_cell.end();
1354 * ++particle, ++particle_index)
1356 * const Tensor<1, dim> &E =
1357 * evaluator.get_gradient(particle_index);
1361 * Having now obtained the electric field at the location of one
1362 * of the particles, we use this to update first the velocity
1363 * and then the position. To do so, let us first get the old
1364 * velocity out of the properties stored with the particle,
1365 * compute the acceleration, update the velocity, and store this
1366 * new velocity again in the properties of the particle. Recall
1367 * that this corresponds to the first of the following set of
1368 * update equations discussed in the introduction:
1370 * \frac{{\mathbf v}_i^{(n)}
1371 * -{\mathbf v}_i^{(n-1)}}{\Delta t}
1372 * &= \frac{e\nabla V^{(n)}}{m}
1373 * \\ \frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}
1374 * {\Delta t} &= {\mathbf v}_i^{(n)}.
1378 * const Tensor<1, dim> old_velocity(particle->get_properties());
1380 * const Tensor<1, dim> acceleration =
1381 * Constants::electron_charge / Constants::electron_mass * E;
1383 * const Tensor<1, dim> new_velocity =
1384 * old_velocity + acceleration * dt;
1386 * particle->set_properties(new_velocity);
1390 * With the new velocity, we can then also update the location
1391 * of the particle and tell the particle about it.
1394 * const Point<dim> new_location =
1395 * particle->get_location() + dt * new_velocity;
1396 * particle->set_location(new_location);
1403 * Having updated the locations and properties (i.e., velocities) of all
1404 * particles, we need to make sure that the `particle_handler` again knows
1405 * which cells they are in, and what their locations in the coordinate
1406 * system of the reference cell are. The following function does that. (It
1407 * also makes sure that, in parallel computations, particles are moved from
1408 * one processor to another processor if a particle moves from the subdomain
1409 * owned by the former to the subdomain owned by the latter.)
1412 * particle_handler.sort_particles_into_subdomains_and_cells();
1419 * <a name="step_19-CathodeRaySimulatortrack_lost_particle"></a>
1420 * <h4>CathodeRaySimulator::track_lost_particle</h4>
1424 * The final particle-related function is the one that is called whenever a
1425 * particle is lost from the simulation. This typically happens if it leaves
1426 * the domain. If that happens, this function is called both the cell (which
1427 * we can ask for its new location) and the cell it was previously on. The
1428 * function then keeps track of updating the number of particles lost in this
1429 * time step, the total number of lost particles, and then estimates whether
1430 * the particle left through the hole in the middle of the anode. We do so by
1431 * first checking whether the cell it was in last had an @f$x@f$ coordinate to the
1432 * left of the right boundary (located at @f$x=4@f$) and the particle now has a
1433 * position to the right of the right boundary. If that is so, we compute a
1434 * direction vector of its motion that is normalized so that the @f$x@f$ component
1435 * of the direction vector is equal to @f$1@f$. With this direction vector, we can
1436 * compute where it would have intersected the line @f$x=4@f$. If this intersect
1437 * is between @f$0.5@f$ and @f$1.5@f$, then we claim that the particle left through
1438 * the hole and increment a counter.
1441 * template <int dim>
1442 * void CathodeRaySimulator<dim>::track_lost_particle(
1443 * const typename Particles::ParticleIterator<dim> &particle,
1444 * const typename Triangulation<dim>::active_cell_iterator &cell)
1446 * ++n_recently_lost_particles;
1447 * ++n_total_lost_particles;
1449 * const Point<dim> current_location = particle->get_location();
1450 * const Point<dim> approximate_previous_location = cell->center();
1452 * if ((approximate_previous_location[0] < 4) && (current_location[0] > 4))
1454 * const Tensor<1, dim> direction =
1455 * (current_location - approximate_previous_location) /
1456 * (current_location[0] - approximate_previous_location[0]);
1458 * const double right_boundary_intercept =
1459 * approximate_previous_location[1] +
1460 * (4 - approximate_previous_location[0]) * direction[1];
1461 * if ((right_boundary_intercept > 0.5) &&
1462 * (right_boundary_intercept < 1.5))
1463 * ++n_particles_lost_through_anode;
1472 * <a name="step_19-CathodeRaySimulatorupdate_timestep_size"></a>
1473 * <h4>CathodeRaySimulator::update_timestep_size</h4>
1477 * As discussed at length in the introduction, we need to respect a time step
1478 * condition whereby particles can not move further than one cell in one time
1479 * step. To ensure that this is the case, we again first compute the maximal
1480 * speed of all particles on each cell, and divide the cell size by that
1481 * speed. We then compute the next time step size as the minimum of this
1482 * quantity over all cells, using the safety factor discussed in the
1483 * introduction, and set this as the desired time step size using the
1484 * DiscreteTime::set_desired_time_step_size() function.
1487 * template <int dim>
1488 * void CathodeRaySimulator<dim>::update_timestep_size()
1490 * if (time.get_step_number() > 0)
1492 * double min_cell_size_over_velocity = std::numeric_limits<double>::max();
1494 * for (const auto &cell : dof_handler.active_cell_iterators())
1495 * if (particle_handler.n_particles_in_cell(cell) > 0)
1497 * const double cell_size = cell->minimum_vertex_distance();
1499 * double max_particle_velocity(0.0);
1501 * for (const auto &particle :
1502 * particle_handler.particles_in_cell(cell))
1504 * const Tensor<1, dim> velocity(particle.get_properties());
1505 * max_particle_velocity =
1506 * std::max(max_particle_velocity, velocity.norm());
1509 * if (max_particle_velocity > 0)
1510 * min_cell_size_over_velocity =
1511 * std::min(min_cell_size_over_velocity,
1512 * cell_size / max_particle_velocity);
1515 * constexpr double c_safety = 0.5;
1516 * time.set_desired_next_step_size(c_safety * 0.5 *
1517 * min_cell_size_over_velocity);
1521 * As mentioned in the introduction, we have to treat the very first
1522 * time step differently since there, particles are not available yet or
1523 * do not yet have the information associated that we need for the
1524 * computation of a reasonable step length. The formulas below follow the
1525 * discussion in the introduction.
1530 * const QTrapezoid<dim> vertex_quadrature;
1531 * FEValues<dim> fe_values(fe, vertex_quadrature, update_gradients);
1533 * std::vector<Tensor<1, dim>> field_gradients(vertex_quadrature.size());
1535 * double min_timestep = std::numeric_limits<double>::max();
1537 * for (const auto &cell : dof_handler.active_cell_iterators())
1538 * if (particle_handler.n_particles_in_cell(cell) > 0)
1540 * const double cell_size = cell->minimum_vertex_distance();
1542 * fe_values.reinit(cell);
1543 * fe_values.get_function_gradients(solution, field_gradients);
1546 * for (const auto q_point : fe_values.quadrature_point_indices())
1547 * max_E = std::max(max_E, field_gradients[q_point].norm());
1551 * std::min(min_timestep,
1552 * std::sqrt(0.5 * cell_size *
1553 * Constants::electron_mass /
1554 * Constants::electron_charge / max_E));
1557 * time.set_desired_next_step_size(min_timestep);
1566 * <a name="step_19-ThecodeCathodeRaySimulatoroutput_resultscodefunction"></a>
1567 * <h4>The <code>CathodeRaySimulator::output_results()</code> function</h4>
1571 * The final function implementing pieces of the overall algorithm is the one
1572 * that generates graphical output. In the current context, we want to output
1573 * both the electric potential field as well as the particle locations and
1574 * velocities. But we also want to output the electric field, i.e., the
1575 * gradient of the solution.
1579 * deal.II has a general way how one can compute derived quantities from
1580 * the solution and output those as well. Here, this is the electric
1581 * field, but it could also be some other quantity -- say, the norm of the
1582 * electric field, or in fact anything else one could want to compute from
1583 * the solution @f$V_h(\mathbf x)@f$ or its derivatives. This general solution
1584 * uses the DataPostprocessor class and, in cases like the one here where we
1585 * want to output a quantity that represents a vector field, the
1586 * DataPostprocessorVector class.
1590 * Rather than try and explain how this class works, let us simply refer to
1591 * the documentation of the DataPostprocessorVector class that has essentially
1592 * this case as a well-documented example.
1595 * template <int dim>
1596 * class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
1599 * ElectricFieldPostprocessor()
1600 * : DataPostprocessorVector<dim>("electric_field", update_gradients)
1603 * virtual void evaluate_scalar_field(
1604 * const DataPostprocessorInputs::Scalar<dim> &input_data,
1605 * std::vector<Vector<double>> &computed_quantities) const override
1607 * AssertDimension(input_data.solution_gradients.size(),
1608 * computed_quantities.size());
1610 * for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
1612 * AssertDimension(computed_quantities[p].size(), dim);
1613 * for (unsigned int d = 0; d < dim; ++d)
1614 * computed_quantities[p][d] = input_data.solution_gradients[p][d];
1623 * With this, the `output_results()` function becomes relatively
1624 * straightforward: We use the DataOut class as we have in almost
1625 * every one of the previous tutorial programs to output the
1626 * solution (the "electric potential") and we use the postprocessor
1627 * defined above to also output its gradient (the "electric
1628 * field
"). This all is then written into a file in VTU format after
1629 * also associating the current time and time step number with this
1630 * file, and providing physical units for the two fields to be
1634 * template <int dim>
1635 * void CathodeRaySimulator<dim>::output_results() const
1638 * ElectricFieldPostprocessor<dim> electric_field;
1639 * DataOut<dim> data_out;
1640 * data_out.attach_dof_handler(dof_handler);
1641 * data_out.add_data_vector(solution, "electric_potential
");
1642 * data_out.add_data_vector(solution, electric_field);
1643 * data_out.build_patches();
1645 * DataOutBase::VtkFlags output_flags;
1646 * output_flags.time = time.get_current_time();
1647 * output_flags.cycle = time.get_step_number();
1648 * output_flags.physical_units["electric_potential
"] = "V
";
1649 * output_flags.physical_units["electric_field
"] = "V/m
";
1651 * data_out.set_flags(output_flags);
1653 * std::ofstream output("solution-
" +
1654 * Utilities::int_to_string(time.get_step_number(), 4) +
1656 * data_out.write_vtu(output);
1661 * Output the particle positions and properties is not more complicated. The
1662 * Particles::DataOut class plays the role of the DataOut class for
1663 * particles, and all we have to do is tell that class where to take
1664 * particles from and how to interpret the `dim` components of the
1665 * properties -- namely, as a single vector indicating the velocity, rather
1666 * than as `dim` scalar properties. The rest is then the same as above:
1670 * Particles::DataOut<dim> particle_out;
1671 * particle_out.build_patches(
1673 * std::vector<std::string>(dim, "velocity
"),
1674 * std::vector<DataComponentInterpretation::DataComponentInterpretation>(
1675 * dim, DataComponentInterpretation::component_is_part_of_vector));
1677 * DataOutBase::VtkFlags output_flags;
1678 * output_flags.time = time.get_current_time();
1679 * output_flags.cycle = time.get_step_number();
1680 * output_flags.physical_units["velocity
"] = "m/s
";
1682 * particle_out.set_flags(output_flags);
1684 * std::ofstream output("particles-
" +
1685 * Utilities::int_to_string(time.get_step_number(), 4) +
1687 * particle_out.write_vtu(output);
1695 * <a name="step_19-CathodeRaySimulatorrun
"></a>
1696 * <h4>CathodeRaySimulator::run</h4>
1700 * The last member function of the principal class of this program is then the
1701 * driver. At the top, it refines the mesh a number of times by solving the
1702 * problem (with no particles yet created) on a sequence of finer and finer
1706 * template <int dim>
1707 * void CathodeRaySimulator<dim>::run()
1711 * const unsigned int n_pre_refinement_cycles = 3;
1712 * for (unsigned int refinement_cycle = 0;
1713 * refinement_cycle < n_pre_refinement_cycles;
1714 * ++refinement_cycle)
1717 * assemble_system();
1725 * Now do the loop over time. The sequence of steps follows closely the
1726 * outline of the algorithm discussed in the introduction. As discussed in
1727 * great detail in the documentation of the DiscreteTime class, while we
1728 * move the field and particle information forward by one time step, the
1729 * time stored in the `time` variable is not consistent with where (some of)
1730 * these quantities are (in the diction of DiscreteTime, this is the "update
1731 * stage
"). The call to `time.advance_time()` makes everything consistent
1732 * again by setting the `time` variable to the time at which the field and
1733 * particles already are, and once we are in this "consistent stage
", we can
1734 * generate graphical output and write information about the current state
1735 * of the simulation to screen.
1741 * std::cout << "Timestep
" << time.get_step_number() + 1 << std::endl;
1742 * std::cout << " Field degrees of freedom:
"
1743 * << dof_handler.n_dofs() << std::endl;
1745 * assemble_system();
1748 * create_particles();
1749 * std::cout << " Total number of particles in simulation:
"
1750 * << particle_handler.n_global_particles() << std::endl;
1752 * n_recently_lost_particles = 0;
1753 * update_timestep_size();
1756 * time.advance_time();
1760 * std::cout << " Number of particles lost
this time step:
"
1761 * << n_recently_lost_particles << std::endl;
1762 * if (n_total_lost_particles > 0)
1763 * std::cout << " Fraction of particles lost through anode:
"
1764 * << 1. * n_particles_lost_through_anode /
1765 * n_total_lost_particles
1768 * std::cout << std::endl
1769 * << " Now at t=
" << time.get_current_time()
1770 * << ", dt=
" << time.get_previous_step_size() << '.'
1774 * while (time.is_at_end() == false);
1776 * } // namespace Step19
1783 * <a name="step_19-Thecodemaincodefunction
"></a>
1784 * <h3>The <code>main</code> function</h3>
1788 * The final function of the program is then again the `main()` function. It is
1789 * unchanged in all tutorial programs since @ref step_6 "step-6
" and so there is nothing new
1797 * Step19::CathodeRaySimulator<2> cathode_ray_simulator_2d;
1798 * cathode_ray_simulator_2d.run();
1800 * catch (std::exception &exc)
1802 * std::cerr << std::endl
1804 * << "----------------------------------------------------
"
1806 * std::cerr << "Exception on processing:
" << std::endl
1807 * << exc.what() << std::endl
1808 * << "Aborting!
" << std::endl
1809 * << "----------------------------------------------------
"
1816 * std::cerr << std::endl
1818 * << "----------------------------------------------------
"
1820 * std::cerr << "Unknown exception!
" << std::endl
1821 * << "Aborting!
" << std::endl
1822 * << "----------------------------------------------------
"
1829<a name="step_19-Results
"></a><h1>Results</h1>
1832When this program is run, it produces output that looks as follows:
1835 Field degrees of freedom: 4989
1836 Total number of particles in simulation: 20
1837 Number of particles lost this time step: 0
1839 Now at t=2.12647e-07, dt=2.12647e-07.
1842 Field degrees of freedom: 4989
1843 Total number of particles in simulation: 40
1844 Number of particles lost this time step: 0
1846 Now at t=4.14362e-07, dt=2.01715e-07.
1849 Field degrees of freedom: 4989
1850 Total number of particles in simulation: 60
1851 Number of particles lost this time step: 0
1853 Now at t=5.23066e-07, dt=1.08704e-07.
1856 Field degrees of freedom: 4989
1857 Total number of particles in simulation: 80
1858 Number of particles lost this time step: 0
1860 Now at t=6.08431e-07, dt=8.53649e-08.
1867 Field degrees of freedom: 4989
1868 Total number of particles in simulation: 394
1869 Number of particles lost this time step: 2
1870 Fraction of particles lost through anode: 0.823188
1872 Now at t=3.83822e-05, dt=4.18519e-08.
1875 Field degrees of freedom: 4989
1876 Total number of particles in simulation: 392
1877 Number of particles lost this time step: 0
1878 Fraction of particles lost through anode: 0.823188
1880 Now at t=3.84244e-05, dt=4.22178e-08.
1887 Field degrees of freedom: 4989
1888 Total number of particles in simulation: 364
1889 Number of particles lost this time step: 0
1890 Fraction of particles lost through anode: 0.867745
1892 Now at t=9.99509e-05, dt=4.20455e-08.
1895 Field degrees of freedom: 4989
1896 Total number of particles in simulation: 364
1897 Number of particles lost this time step: 2
1898 Fraction of particles lost through anode: 0.867871
1900 Now at t=9.99929e-05, dt=4.20105e-08.
1903 Field degrees of freedom: 4989
1904 Total number of particles in simulation: 362
1905 Number of particles lost this time step: 2
1906 Fraction of particles lost through anode: 0.867047
1908 Now at t=0.0001, dt=7.11401e-09.
1911Picking a random few time steps, we can visualize the solution in the
1912form of streamlines for the electric field and dots for the electrons
1913(the pictures were created with an earlier version of the program and
1914are no longer exact representations of the solution as currently
1915computed -- but they are conceptually still correct):
1916<div class="twocolumn
" style="width: 80%
">
1919 alt=
"The solution at time step 0 (t=0 seconds)."
1922 Solution at time step 0 (t=0 seconds).
1926 <img src=
"https://www.dealii.org/images/steps/developer/step-19.solution.1400.png"
1927 alt=
"The solution at time step 1400 (t=0.000068 seconds)."
1930 Solution at time step 1400 (t=0.000068 seconds).
1934 <img src=
"https://www.dealii.org/images/steps/developer/step-19.solution.0700.png"
1935 alt=
"The solution at time step 700 (t=0.000035 seconds)."
1938 Solution at time step 700 (t=0.000035 seconds).
1942 <img src=
"https://www.dealii.org/images/steps/developer/step-19.solution.2092.png"
1943 alt=
"The solution at time step 2092 (t=0.0001 seconds)."
1946 Solution at time step 2092 (t=0.0001 seconds).
1951That said, a more appropriate way to visualize the results of
this
1952program are by creating a video that shows how these electrons move, and how
1953the electric field changes in response to their motion (again showing
1954results of an earlier version of the program):
1958 <iframe width=
"560" height=
"315" src=
"https://www.youtube.com/embed/HwUtE7xuteE"
1960 allow=
"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
1961 allowfullscreen></iframe>
1965What you can see here is how the
"focus element" of the boundary with its
negative
1966voltage repels the electrons and makes sure that they do not just fly away
1967perpendicular from the cathode (as they do in the
initial part of their
1968trajectories). It also shows how the electric field lines
1969move around over time, in response to the charges flying by -- in other words,
1970the feedback the particles have on the electric field that itself drives the
1971motion of the electrons.
1973The movie suggests that electrons move in
"bunches" or
"bursts". One element of
1974this appearance is an artifact of how the movie was created: Every frame of the
1975movie corresponds to one time step, but the time step length varies. More specifically,
1976the fastest particle moving through the smallest cell determines the length of the
1977time step (see the discussion in the introduction), and consequently time steps
1978are small whenever a (fast) particle moves through the small cells at the right
1979edge of the domain; time steps are longer again once the particle has left
1980the domain. This slowing-accelerating effect can easily be visualized by plotting
1981the time step length shown in the screen output.
1983The
second part of
this is real, however: The simulation creates a large
group
1984of particles in the beginning, and fewer after about the 300th time step. This
1985is probably because of the
negative charge of the particles in the simulation:
1986They
reduce the magnitude of the electric field at the (also negatively charged
1987electrode) and consequently
reduce the number of points on the cathode at which
1988the magnitude exceeds the threshold necessary to draw an electron out of the
1992<a name=
"step-19-extensions"></a>
1993<a name=
"step_19-Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
1997<a name=
"step_19-Morestatisticsaboutelectrons"></a><h4> More statistics about electrons </h4>
2000At the
end of the day, we are rarely interested in the *solution* of an equation,
2001but in
numbers that can be *extracted* from it -- in other words, we want to
2002*postprocess* the solution (see the results section of @ref step_4
"step-4" for an extensive
2003discussion of
this concept). Here, what one would likely be most interested in
2004assessing are some statistics about the particles.
2006The program already computes the fraction of the electrons that leave the
2007domain through the hole in the anode. But there are other quantities one might be
2008interested in. For example, the average velocity of these particles. It would
2009not be very difficult to obtain each particle
's velocity from its properties,
2010in the same way as we do in the `move_particles()` function, and compute
2014<a name="step_19-Abettersynchronizedvisualization"></a><h4> A better-synchronized visualization </h4>
2017As discussed above, there is a varying time difference between different frames
2018of the video because we create output for every time step. A better way to
2019create movies would be to generate a new output file in fixed time intervals,
2020regardless of how many time steps lie between each such point.
2023<a name="step_19-Abettertimeintegrator"></a><h4> A better time integrator </h4>
2026The problem we are considering in this program is a coupled, multiphysics
2027problem. But the way we solve it is by first computing the (electric) potential
2028field and then updating first the particle velocity and finally the
2029particle locations. This is what is called an
2030"operator-splitting method", a concept we will investigate in more detail
2031in @ref step_58 "step-58".
2033While it is awkward to think of a way to solve this problem that does not involve
2034splitting the problem into a PDE piece and a particles piece, one *can*
2035(and probably should!) think of a better way to update the particle
2036locations. Specifically, the equations we use to update the particle location
2039 \frac{{\mathbf v}_i^{(n)}-{\mathbf v}_i^{(n-1)}}{\Delta t} &= \frac{e\nabla V^{(n)}}{m}
2041 \frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}{\Delta t} &= {\mathbf v}_i^{(n)}.
2043This corresponds to a Lie splitting where we first update one variable
2044(the velocity) and then, using the already-updated value of the one variable
2045to update the other (the position). It is well understood that the Lie
2046splitting incurs a first order error accuracy in the time step (i.e.,
2047the error introduced by the splitting is @f${\cal O}(\Delta t)@f$) that we know we should
2048avoid because we can do better. Independently, of course, we incur an error
2049of the same rate a second time because we use an Euler-like scheme for each of
2050the two updates (we replace the time derivative by a simple finite difference
2051quotient containing the old and new value, and we evaluate the right hand side
2052only at the end time of the interval, at @f$t_n@f$), and a better scheme would also
2053use a better way to do this kind of update.
2055Better strategies would replace the Lie splitting by something like a
2056[Strang splitting](https://en.wikipedia.org/wiki/Strang_splitting), and the
2057update of the particle position and velocity by a scheme such
2059[leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration)
2061[symplectic integrators](https://en.wikipedia.org/wiki/Symplectic_integrator)
2063[Verlet scheme](https://en.wikipedia.org/wiki/Verlet_integration).
2066<a name="step_19-Parallelization"></a><h4> Parallelization </h4>
2069In release mode, the program runs in about 3.5 minutes on one of the author's
2070laptops at the time of writing this. That
's acceptable. But what if we wanted
2071to make the simulation three-dimensional? If we wanted to not use a maximum
2072of around 100 particles at any given time (as happens with the parameters
2073used here) but 100,000? If we needed a substantially finer mesh?
2075In those cases, one would want to run the program not just on a single processor,
2076but in fact on as many as one has available. This requires parallelization
2077of both the PDE solution as well as over particles. In practice, while there
2078are substantial challenges to making this efficient and scale well, these
2079challenges are all addressed in deal.II itself, and have been demonstrated
2080on simulations running on 10,000 or more cores. For example, @ref step_40 "step-40" shows
2081how to parallelize the finite element part, and @ref step_70 "step-70" shows how one can
2082then also parallelize the particles part.
2085<a name="step_19-PlainProg"></a>
2086<h1> The plain program</h1>
2087@include "step-19.cc"
double JxW(const unsigned int q_point) const
particle_iterator_range particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
#define AssertDimension(dim1, dim2)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
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)