113 lines
17 KiB
TeX
113 lines
17 KiB
TeX
\subsection*{Numerical Accuracy of Trajectories}
|
|
As there is no closed form solution to the trajectories in a Penning trap containing multiple particles, which interact through Coulomb forces, the case of a single particle is used to verify the numerical accuracy of the implemented algorithms. As a metric for the accuracy, the relative error
|
|
\begin{equation} \label{eq:rel_error}
|
|
\epsilon = \frac{|\vec{r}_{\text{exact}} - \vec{r}_{\text{numerical}}|}{|\vec{r}_{\text{exact}}|}
|
|
\end{equation}
|
|
is used, where $\vec{r}_{\text{exact}}$ as derived in \cref{sec:methods}, and $\vec{r}_{\text{numerical}}$ is the numerical solution obtained from the implemented algorithms. The relative error is calculated at each time step, and the maximum value over the entire simulation is reported as the error for that particular simulation.
|
|
|
|
Apart from the quantitative measure of the accuracy, the qualitative behavior of the trajectories is also examined. The numerical solutions are plotted and compared to the expected behavior. Especially properties like energy drift of the solution and the general shape of the trajectory can give insight into the performance of the algorithms.
|
|
|
|
\begin{figure}
|
|
\centering
|
|
\includegraphics[width=\columnwidth]{../../src/project3/python/plots/xy_trajectory/two_particles_2_no_interactions_4000.000000_steps_rk4_positions_xy_trajectory.pdf}
|
|
\caption{Trajectory of two particles in a Penning trap, simulated using the RK4 method. The particles are initialized as: $\vec r_1 = (\qty{20}{\um}, \qty{0}{um}, \qty{20}{\um})$, $\vec v_1 = (\qty{0}{\um/\us}, \qty{25}{\um/\us}, \qty{0}{\um/\us})$, $\vec r_2 = (\qty{25}{\um}, \qty{25}{\um}, \qty{0}{\um})$, $\vec v_2 = (\qty{0}{\um/\us}, \qty{40}{\um/\us}, \qty{5}{\um/\us})$.}
|
|
\label{fig:two_particle_trajectory_rk4}
|
|
\end{figure}
|
|
The Runge-Kutta 4 method shows a very good qualitative behavior, with no apparent inaccuracies in the trajectory, as shown in \cref{fig:two_particle_trajectory_rk4}. Comparing against the same trajectories calculated using the Velocity-Verlet method, we can immediately see shortcomings of the latter. Trajectories for the Velocity-Verlet algorithm, the Boris algorithm and Euler algorithm are included in \cref{fig:two_particle_trajectory_verlet,fig:two_particle_trajectory_boris,fig:two_particle_trajectory_euler} in \cref{app:further_results}. The trajectory shows a significant drift over time, which is not present in the RK4 solution. The Verlet solution rapidly increases its radius, which is not physical in a system with conserved energy. While one would expect the symplectic Verlet method to perform better in this regard, the velocity-dependant Lorentz force is not suited for computations using the Velocity-Verlet algorithm, as discussed in \cref{sec:methods}.
|
|
A similar behavior of non-physical qualitative behavior is observed for the Euler method, which also shows a significant energy drift. The Boris algorithm, which is specifically designed to handle electromagnetic forces, shows a qualitatively correct behavior, with no apparent energy drift and a stable orbit over time.
|
|
|
|
|
|
As well as the trajectory in the transversal plane, the osciallatory behavior in the $z$ direction and the phase space usage in all three dimensions is examined. As there is no influence of the magnetic field in the $z$ direction, the oscillations should be perfectly harmonic. Furthermore should the Velocity-Verlet method exhibit correct behavior, as there is no velocity-dependant force in this direction. The very good performance of RK4, Boris and Velocity-Verlet in this regard is confirmed, while the Euler method shows energy drift and a non-harmonic behavior with increasing amplitude over time.
|
|
|
|
The numerical analysis of the trajectories shows similar behavior for the two algorithms designed to handle electromagnetic forces, RK4 and Boris. The evolution of the relative error $\epsilon$ as defined in \cref{eq:rel_error} over time is shown in \cref{fig:relative_error}. Both algorithms show a steady osciallation in the relative error around \numrange{1e-3}{2e-3}, with RK4 showing a slightly lower error overall. The Velocity-Verlet and Euler methods show a significantly higher error, with the Euler method reaching a relative error of up to \num{0.1} at certain times. The Velocity-Verlet method shows a relative error of up to \num{0.2}.
|
|
|
|
\begin{figure}
|
|
\centering
|
|
\includegraphics[width=\columnwidth]{../../src/project3/python/plots/error_plot/two_particles_1_no_interactions_32000.000000_steps_analytical_positions_error_plot.pdf}
|
|
\caption{Relative error $\epsilon$ for a single particle in a Penning trap simulated for \num{32000} time steps. The particle is initialized as: $\vec r = (\qty{20}{\um}, \qty{0}{um}, \qty{20}{\um})$, $\vec v = (\qty{0}{\um/\us}, \qty{25}{\um/\us}, \qty{0}{\um/\us})$.}
|
|
\label{fig:relative_error}
|
|
\end{figure}
|
|
|
|
\subsection*{Many-Body Simulations}
|
|
In addition to the case of perfect trajectories without Coulomb interactions, the behavior of a cloud of up to \num{100} particles is examined. This use case is of particular interest, as it is closer to real-world applications of Penning traps, such as in mass spectrometry. The limit of increasing particle number is also of interest, as the computational complexity increases significantly with the number of particles, due to the expected $\mathcal{O}(N^2)$ scaling of the Coulomb interaction calculations. Additionally, the collective behavior of particles, i.e. the statistical properties of the ensemble, and the accuracy thereof is studied.
|
|
|
|
All particles are initialized with random positions and velocities, sampled from a normal distribution. The standard deviation of the position distribution is set to \qty{50}{\um} in all three dimensions, while the standard deviation of the velocity distribution is set to \qty{50}{\um/\us}. To compare the influence of Coulomb interactions on the simulation performance and statistical properties, simulations with and without interactions are performed.
|
|
|
|
\subsubsection*{Resonant Behavior}
|
|
\begin{figure}
|
|
\centering
|
|
\includegraphics[width=\columnwidth]{../../src/project3/python/plots/resonances/Boris_ many_particles_100_particles_40000_steps_enabled_interactions_oscillating_potential.png}
|
|
\caption{Fraction of particles remaining in the trap after \qty{50}{\us} of simulation time using Boris algorithm, as a function of the frequency of the dynamic trapping potential $\vdyn$ with different oscillatory amplitudes. The first scan (dashed) is performed without Coulomb interactions, while the second scan (solid) is performed with Coulomb interactions enabled.}
|
|
\label{fig:resonance_scan}
|
|
\end{figure}
|
|
|
|
A further behavior that can be studied in the many-body simulations is the resonant behavior of the particle cloud when a dynamic trapping potential $\vdyn$ is used. By varying the frequency of the dynamic potential $\omega_V$, we can study how the stability of the trap is influenced by resonances. The stability is measured by the fraction of particles remaining in the trap after a fixed simulation time of \qty{50}{\us}. A particle is considered to be inside the trap, if its distance from the center of the trap is less than $d$, where $d$ is the characteristic dimension of the trap, as defined in \cref{sec:methods}. In a first scan of the frequency range from \qtyrange{0.2}{2.5}{\mega\hertz}, three resonances are found, at \qty{0.7}{\mega\hertz}, \qty{1.4}{\mega\hertz} and \qty{2.1}{\mega\hertz}. The first scan is performed without Coulomb interactions between the particles, as to speed up the simulation. In a second scan, the resonance in the range of \qtyrange{1.0}{1.8}{\mega\hertz} is studied in more detail, with a finer frequency resolution and with Coulomb interactions enabled. The results of both scans are shown in \cref{fig:resonance_scan}.
|
|
|
|
|
|
\subsection*{Performance and Efficiency}
|
|
|
|
|
|
|
|
\begin{table}
|
|
\centering
|
|
\caption{Relative run times of the different solvers in \unit{\um} per particle and time step. The times are averaged over 4 runs with different time steps. Coulomb interactions are enabled and $\vstat$ is used as the trapping potential. Output during the simulation is enabled.}
|
|
\label{tab:solver_scaling}
|
|
\begin{tabular}{lccc}
|
|
\toprule
|
|
Num. Particles & 10 & 50 & 100 \\
|
|
Solver & & & \\
|
|
\midrule
|
|
Euler & \num{1.45(30)e+00} & \num{5.90(41)e+00} & \num{1.244(90)e+01} \\
|
|
Velocity Verlet & \num{1.46(17)e+00} & \num{6.22(52)e+00} & \num{1.271(61)e+01} \\
|
|
Runge-Kutta 4 & \num{3.71(65)e+00} & \num{1.57(13)e+01} & \num{3.58(57)e+01} \\
|
|
Boris & \num{1.43(22)e+00} & \num{5.95(51)e+00} & \num{1.263(60)e+01} \\
|
|
\bottomrule
|
|
\end{tabular}
|
|
\end{table}
|
|
|
|
For different step sizes, i.e. a different number of total time steps, the total run time of a simulation with different numbers of particles is measured. For measuring the run time, the \texttt{ctime} module with the \texttt{std::clock()} function is used. All runs were performed on the same machine, with an AMD Ryzen 7 4750U CPU and \qty{16}{\giga\byte} of RAM, running Fedora 42. For parallelization of the runs the scheduler \texttt{pueue} was used with \num{14} simultaneous jobs. To reduce the influence of other processes on the machine, the runs were performed at night when the machine was otherwise idle. Furthermore, the runs were interweaved in a way, that no solving algorithm could take dvantage from time dependant properties of the machine, such as thermal throttling. The output during the simulation was enabled, as the output is reused for the analysis of the statistical properties of the particle ensemble. The I/O operations were performed outside of the timed section of the code, to get a more accurate measure of the run time of the algorithms themselves. The times were averaged over \num{4} runs with different time steps, and the relative time per particle and time step was calculated. The results are shown in \cref{tab:solver_scaling}.
|
|
|
|
The results show a clear linear scaling of the relative time per particle and time step with increasing number of particles, as expected from the $\mathcal{O}(N^2)$ scaling of the Coulomb interaction calculations. The Euler, Velocity-Verlet and Boris methods show a very similar performance, with the RK4 method being slower by a factor of approximately \num{2.8}. The similar performance of the Euler, Velocity-Verlet and Boris methods is expected, as all three methods require a single force calculation per time step. The RK4 method requires four force calculations per time step, which explains the significantly higher run time. From the scaling behavior for simulations with Coulomb interactions, we can derive that the major contributing factor to the run time is the calculation of the interactions.
|
|
|
|
|
|
\begin{table}
|
|
\centering
|
|
\caption{Total run times of the different solvers in \unit{\s} for a simulation of \num{100} particles and \num{40000} time steps. The dynamic potential $\vdyn$ is used, and Coulomb interactions are both enabled and disabled. The output during the simulation is disabled to get a more accurate measure of the run time.}
|
|
\label{tab:solver_comparison}
|
|
\begin{tabular}{lcc}
|
|
\toprule
|
|
& Interactions Disabled & Interactions Enabled \\
|
|
Solver & & \\
|
|
\midrule
|
|
Velocity Verlet & \num{1.23(26)e+00} & \num{3.14(16)e+01} \\
|
|
Runge-Kutta 4 & \num{8.8(11)e+00} & \num{1.267(48)e+02} \\
|
|
Boris & \num{1.07(24)e+00} & \num{3.11(16)e+01} \\
|
|
\bottomrule
|
|
\end{tabular}
|
|
|
|
\end{table}
|
|
|
|
For a more precise comparison on how the four different algorithms compare, with respect to run time, a simulation with \num{100} particles and \num{40000} time steps is performed. There is no output recording during the simulation, as the simulations for the resonant behavior are reused. The lack of output recording allows for a more accurate measure of the run time of the algorithms themselves. The results are shown in \cref{tab:solver_comparison}. The Boris algorithm will serve as a baseline for the comparison, as it shows the lowest execution time of the three algorithms tested. The Velocity-Verlet algorithm shows a similar performance, with an increase in CPU time of \qty{0.75+-7.33}{\percent} with interactions and \qty{14.9+-34.95}{\percent} without interactions. The RK4 method shows a significantly higher run time, with an increase of \qty{307+-26}{\percent} with interactions and \qty{721+-210}{\percent} without interactions. The stark difference in performance between the RK4 method and the other two methods is expected, as the RK4 method requires four force calculations per time step, while the other two methods only require a single force calculation per time step. The slight increase in run time of the Velocity-Verlet method compared to the Boris method is unexpected, but within the margin of error of the measurements. The error of each run time is estimated via the standard deviation of \num{348} runs without interactions and \num{483} runs with interactions. The error on the relative speed decrease is estimated via Gaussian error propagation assuming uncorrelated errors.
|
|
|
|
The high performance of the Boris method in comparison to RK4 is also an important factor to consider when assessing the numerical accuracy of the trajectories. While both algorithms performed similarly in terms of trajectory accuracy at the same step size, the decreased computation time for Boris algorithm allows for a significantly smaller step size at the same computation time. This in turn allows for a higher numerical accuracy of the trajectories, as the error generally decreases with decreasing step size.
|
|
|
|
\subsection*{Symplectic Properties}
|
|
While a full analysis of the symplectic properties of the different algorithms is outside the scope of this project, a brief qualitative analysis is performed. The symplectic nature of an algorithm is especially important in long-term simulations, as non-symplectic algorithms may produce incorrect ensemble properties over long time scales, even if the short-term accuracy is high. The Velocity-Verlet algorithm is symplectic by design for Hamiltonian systems, while the Boris algorithm is also known to exhibit symplectic properties in electromagnetic systems. The RK4 method and the Euler method are not symplectic.
|
|
|
|
A qualitative analysis of the symplectic properties is performed by examining the systems total energy
|
|
\begin{equation}
|
|
E = \sum_{i=1}^N \left( \frac{m_i v_i^2}{2} + q_i V(\vec r_i) + \frac{1}{2} \sum_{\substack{j=1 \\ j \neq i}}^N \frac{q_i q_j}{4 \pi \epsilon_0|\vec r_i - \vec r_j|} \right)
|
|
\end{equation}
|
|
over the simulation duration. As the system is isolated, the total energy should be conserved. A non-conservation of the total energy indicates a non-symplectic behavior of the algorithm.
|
|
|
|
\Cref{fig:total_energy} shows the total energy development of the system for the different algorithms. A more detailed analysis on the dependance on step size is shown in \cref{app:further_results} with \cref{fig:boris_many_particles_energy,fig:euler_many_particles_energy,fig:verlet_many_particles_energy,fig:rk4_many_particles_energy}. The RK4 and Boris methods show a very stable total energy over the entire simulation duration, with slight oscillations around a mean value. The Velocity-Verlet method shows a significant drift to increasing total energy over time, indicating a non-symplectic behavior in this system. The Euler method shows a similar behavior, with a slightly less pronounced increase in total energy over time.
|
|
\begin{figure}
|
|
\centering
|
|
\includegraphics[width=\columnwidth]{../../src/project3/python/plots/many_particles_100_particles_32000_steps_enabled_interactions_total_energies.pdf}
|
|
\caption{Total energy of a system of \num{100} particles in a Penning trap over \num{32000} time steps, simulated using different algorithms. The kinetic component of the total energy is shown as a dotted line, while the potential component is shown as a dashed line. The solid line shows the total energy of the system.}
|
|
\label{fig:total_energy}
|
|
\end{figure}
|
|
|
|
The non-symplectic behavior of the Velocity-Verlet method in this system is expected due to the velocity-dependant Lorentz force, as discussed in \cref{sec:methods}. The non-symplectic behavior of the Euler method is also expected, as it is not a symplectic algorithm. The good energy conservation of the RK4 method is noteworthy, as it is not a symplectic algorithm either. However, a conservation of total energy in a single simulation for a specific time span does not imply symplectic properties of the algorithm. In general RK4 is not considered a symplectic algorithm, and may show non-symplectic behavior in other systems or over longer time spans.
|
|
|
|
Boris method also shows a very good conservation of total energy, which is expected due to its design for electromagnetic systems and symplectic properties therein. |