to track individual objects ("particles") and how their positions evolve. If
this simply leads to a set of ordinary differential equations, for example
if you want to track the positions of the planets in the solar system over
-time, then deal.II is clearly not your right tool. On the other hand, if
+time, then deal.II is clearly not the right tool. On the other hand, if
this evolution is due to the interaction with the solution of partial differential
equation, or if having a mesh to determine which particles interact
with others (such as in the
@f[
\rho(\mathbf x) = \sum_i e\delta(\mathbf x-\mathbf x_i).
@f]
-
-@note
The issue now is that in reality, a cathode ray tube in an old television
yields a current of somewhere around a few milli-Amperes. In the much higher
energy beams of particle accelerators, the current may only be a few
nano-Ampere. But an Ampere is $6\times 10^{18}$ electrons flowing per
second. Now, as you will see in the results section, we really only simulate
-a few microseconds ($10^{-5}$ seconds), but that still results in very very
+a few microseconds ($10^{-6}$ seconds), but that still results in very very
large numbers of electrons -- far more than we can hope to simulate
with a program as small as the current one. As a consequence, let us
presume that each particle represents $N$ electrons. Then the particle
@f[
(Nm) {\ddot {\mathbf x}}_i = (Ne)\mathbf E,
@f]
-which is of course exactly the same as above. On the other hand, the charge
+which is of course exactly the same as above after dividing both sides by $N$.
+On the other hand, the charge
density for these "clumps" of electrons is given by
@f[
\rho(\mathbf x) = \sum_i (Ne)\delta(\mathbf x-\mathbf x_i).
electric field. But realism is not our goal here.)
-@note One may wonder why the equation for the electric field (or, rather,
+As a final thought about the model, one may wonder why the equation for the electric field (or, rather,
the electric potential) has no time derivative whereas the equations for
the electron positions do. In essence, this is a modeling assumption: We
assume that the particles move so slowly that at any given time the
<h3>Time discretization</h3>
-The equations outlined above form a set of coupled differential equations.
+The equations outlined above then form a set of coupled differential equations.
Let us bring them all together in one place again to make that clear:
@f{align*}{
-\epsilon_0 \Delta V &= \sum_i e\delta(\mathbf x-\mathbf x_i)
\\
\frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}{\Delta t} &= {\mathbf v}_i^{(n)}.
@f}
+This scheme can be understood in the framework of operator splitting methods (specifically,
+the "Lie splitting" method) wherein a coupled system is solved by updating one
+variable at a time, using either the old values of other variables (e.g., using
+$\mathbf x_i^{(n-1)}$ in the first equation) or the values of variables that have
+already been updated in a previous sub-step (e.g., using $V^{(n)}$ in the second
+equation).
There are of course many better ways to do a time discretization (for
-example the simple [leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration))
+example the simple [leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration)
+when updating the velocity, or more general Strang splitting methods for the coupled
+system)
but this isn't the point of the tutorial program, and so we will be content
with what we have here. (We will comment on a piece of this puzzle in the
<a href="#extensions">possibilities for extensions</a> section of this program,
When this program is run, it produces output that looks as follows:
```
Timestep 1
- Field degrees of freedom: 4989
+ Field degrees of freedom: 4989
Total number of particles in simulation: 20
Number of particles lost this time step: 0
<h4> More statistics about electrons </h4>
+At the end of the day, we are rarely interested in the *solution* of an equation,
+but in numbers that can be *extracted* from it -- in other words, we want to
+*postprocess* the solution (see the results section of step-4 for an extensive
+discussion of this concept). Here, what one would likely be most interested in
+assessing are some statistics about the particles.
+
The program already computes the fraction of the electrons that leave the
domain through the hole in the anode. But there are other quantities one might be
interested in. For example, the average velocity of these particles. It would
The problem we are considering in this program is a coupled, multiphysics
problem. But the way we solve it is by first computing the (electric) potential
-field and then update the particle locations. This is what is called an
+field and then updating first the particle velocity and finally the
+particle locations. This is what is called an
"operator-splitting method", a concept we will investigate in more detail
in step-58.
While it is awkward to think of a way to solve this problem that does not involve
-splitting the problem into a PDE piece and a particles piece, one
-*can* (and probably should!) think of a better way to update the particle
+splitting the problem into a PDE piece and a particles piece, one *can*
+(and probably should!) think of a better way to update the particle
locations. Specifically, the equations we use to update the particle location
are
@f{align*}{
\\
\frac{{\mathbf x}_i^{(n)}-{\mathbf x}_i^{(n-1)}}{\Delta t} &= {\mathbf v}_i^{(n)}.
@f}
-This corresponds to a simple forward-Euler time discretization -- a method of
-first order accuracy in the time step size $\Delta t$ that we know we should
-avoid because we can do better. Rather, one might want to consider a scheme such
+This corresponds to a Lie splitting where we first update one variable
+(the velocity) and then, using the already-updated value of the one variable
+to update the other (the position). It is well understood that the Lie
+splitting incurs a first order error accuracy in the time step (i.e.,
+the error introduced by the splitting is ${\cal O}(\Delta t)$) that we know we should
+avoid because we can do better. Independently, of course, we incur an error
+of the same rate a second time because we use an Euler-like scheme for each of
+the two updates (we replace the time derivative by a simple finite difference
+quotient containing the old and new value, and we evaluate the right hand side
+only at the end time of the interval, at $t_n$), and a better scheme would also
+use a better way to do this kind of update.
+
+Better strategies would replace the Lie splitting by something like a
+[Strang splitting](https://en.wikipedia.org/wiki/Strang_splitting), and the
+update of the particle position and velocity by a scheme such
as the
[leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration)
or more generally
In those cases, one would want to run the program not just on a single processor,
but in fact on as many as one has available. This requires parallelization
-both the PDE solution as well as over particles. In practice, while there
+of both the PDE solution as well as over particles. In practice, while there
are substantial challenges to making this efficient and scale well, these
-challenges are all addressed in deal.II itself. For example, step-40 shows
+challenges are all addressed in deal.II itself, and have been demonstrated
+on simulations running on 10,000 or more cores. For example, step-40 shows
how to parallelize the finite element part, and step-70 shows how one can
then also parallelize the particles part.
projects / dealii.git / commitdiff