diff --git a/.gitignore b/.gitignore index ed13ab5..068d5c3 100644 --- a/.gitignore +++ b/.gitignore @@ -290,4 +290,5 @@ src/project3/results/*/*.dat src/project3/results/*/*/*.red *.tar.gz *_timing.txt -run_stats.json \ No newline at end of file +run_stats.json +src/project3/python/plots/ \ No newline at end of file diff --git a/projects/project3/chapters/methods.tex b/projects/project3/chapters/methods.tex index e47f9eb..8223735 100644 --- a/projects/project3/chapters/methods.tex +++ b/projects/project3/chapters/methods.tex @@ -91,11 +91,14 @@ The method is time-reversible and conserves energy better over long time periods \end{align} It is important to note that the Velocity-Verlet integrator requires forces to be independent of velocity. In the case, that the algorithm is applicable, it is a very efficient method, as it only requires a single force evaluation per time step, while still being a second order method with a local truncation error per step on the order of $\mathcal{O}(h^3)$ and a global error after $N$ steps on the order of $\mathcal{O}(h^2)$. +\subsubsection*{Boris algorithm} +\textcolor{red}{TODO} + \subsection*{Code Structure} The basic framework for the numerical analysis is based on a \texttt{PenningTrap} class, which contains all the particles present in the trap, as well as a parametrization of the electric and magnetic fields. The particles are represented by a \texttt{Particle} class, which contains the physical properties of the particle, as well as its current position and velocity. With this information, the \texttt{PenningTrap} class can calculate the forces acting on each particle, including the external fields and the particle-particle interactions. The particle-particle interactions can be toggled on and off, allowing for a simulation of both scenarios. The external fields can be modified by supplying a field-method of the form \texttt{external\_field(const arma::vec\& r, double t, const PenningTrap\& trap)}. The reference to the \texttt{PenningTrap} allows for the parameters of the field to be stored in the trap object. The implementation of the \texttt{Particle} and \texttt{PenningTrap} classes can be found in \texttt{/src/project3/include/classes.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/classes.cpp}. \subsection*{Numerical Methods Implementation} -All numerical methods are implemented as classes inheriting from a base class \texttt{Solver}. The base class contains a reference to the \texttt{PenningTrap} object, as well as the time step size. The recording of particle properties over time, like postion and velocity is part of the general \texttt{Solver} class. Each derived class implements the \texttt{step()} method, which updates the state of the system by one time step using the respective numerical method. The implementation of the \texttt{Solver} class and its derived classes can be found in \texttt{/src/project3/include/solvers.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/solvers.cpp}. The following solvers are implemented:\texttt{Forward\-Euler\-Solver}, \texttt{RK4\-Solver}, \texttt{Velocity\-Verlet\-Solver} and \texttt{Analytical\-Solver}, which implements the analytical solution for a special case of initial conditions (see \cref{app:special_case_analytical_solution}). +All numerical methods are implemented as classes inheriting from a base class \texttt{Solver}. The base class contains a reference to the \texttt{PenningTrap} object, as well as the time step size. The recording of particle properties over time, like position and velocity is part of the general \texttt{Solver} class. Each derived class implements the \texttt{step()} method, which updates the state of the system by one time step using the respective numerical method. The implementation of the \texttt{Solver} class and its derived classes can be found in \texttt{/src/project3/include/solvers.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/solvers.cpp}. The following solvers are implemented:\texttt{Forward\-Euler\-Solver}, \texttt{RK4\-Solver}, \texttt{Velocity\-Verlet\-Solver} and \texttt{Analytical\-Solver}, which implements the analytical solution for a special case of initial conditions (see \cref{app:special_case_analytical_solution}). \subsection*{Tools} \textcolor{red}{TODO} \ No newline at end of file diff --git a/projects/project3/chapters/results.tex b/projects/project3/chapters/results.tex index 7ea86a7..1c4a505 100644 --- a/projects/project3/chapters/results.tex +++ b/projects/project3/chapters/results.tex @@ -1,9 +1,7 @@ -\subsection{Numerical Accuracy of Trajectories} +\subsection*{Numerical Accuracy of Trajectories} -\subsection{Energy Conservation} +\subsection*{Performance and Efficiency} +\subsection*{Many-Body Simulations} -\subsection{Performance and Efficiency} - -\subsection{Many-Body Simulations} - +\subsection*{Symplectic Properties} \ No newline at end of file diff --git a/projects/project3/main.pdf b/projects/project3/main.pdf index f84eb78..3999fe9 100644 Binary files a/projects/project3/main.pdf and b/projects/project3/main.pdf differ diff --git a/src/project3/generate_plots.sh b/src/project3/generate_plots.sh new file mode 100755 index 0000000..23e51ac --- /dev/null +++ b/src/project3/generate_plots.sh @@ -0,0 +1,55 @@ +#!/bin/bash +PUEUE=true + +echo "Generating plots..." +rm -rf python/plots + +if $PUEUE; then + echo "Using pueue to generate results in parallel..." + systemctl --user start pueued + pueue parallel 8 + pueue clean + COMMAND_PREFIX="pueue add " +else + echo "Generating results sequentially..." + COMMAND_PREFIX="" +fi + +export SHOW_PLOTS=0 + +# Create simple trajectory plots for the case of few particles +for file in results/two_particles_*_positions.dat + do + $COMMAND_PREFIX uv run python/plot_trajectories.py plot-z-evolution "$file" + $COMMAND_PREFIX uv run python/plot_trajectories.py plot-xy-trajectory "$file" + + velocity_file="${file/positions/velocities}" + for axis in x y z + do + $COMMAND_PREFIX uv run python/plot_trajectories.py phase-space-plot --axis "$axis" "$file" "$velocity_file" + done +done + +# Create error plots for a single particle +for analytical_file in results/two_particles_1_no_interactions_*_analytical_positions.dat +do + euler_file="${analytical_file/analytical/euler}" + # verlet_file="${analytical_file/analytical/velocity_verlet}" + # boris_file="${analytical_file/analytical/boris}" + rk4_file="${analytical_file/analytical/rk4}" + $COMMAND_PREFIX uv run python/plot_trajectories.py plot-error "$euler_file" "$rk4_file" "$analytical_file" +done + +# Create resonance plots for the many particle case +for folder in results/*/many_particles_*_oscillating_potential/ +do + $COMMAND_PREFIX uv run python/plot_resonances.py "$folder" +done + +# Create system energy plots for the many particle case +$COMMAND_PREFIX uv run python/plot_energies.py results/*/many_particles_100_particles_32000_steps_enabled_interactions_total_energies.dat python/plots/many_particles_100_particles_32000_steps_enabled_interactions_total_energies.pdf +for solver in results/*/ +do + name=$(basename "$solver") + $COMMAND_PREFIX uv run python/plot_energies.py "$solver"many_particles_100_particles_*_steps_enabled_interactions_total_energies.dat python/plots/"$name"_many_particles_100_particles_enabled_interactions_total_energies.pdf +done \ No newline at end of file diff --git a/src/project3/generate_results.sh b/src/project3/generate_results.sh index 4821636..b66a535 100755 --- a/src/project3/generate_results.sh +++ b/src/project3/generate_results.sh @@ -2,11 +2,12 @@ PUEUE=true NORMAL_MANY_PARTICLE_SIMULATION=true -OSCILLATING_FIELD_SIMULATION=false +OSCILLATING_FIELD_SIMULATION=true echo "Removing old results..." mkdir -p results rm results/*.dat +rm results/*.txt echo "Compiling binaries..." rm -f two_particles many_particles @@ -31,9 +32,14 @@ do do for n in 1 2 do + if [ $n -eq 1] && [ "$interactions" == "" ]; then + # Skip the case of one particle with interactions, as it is meaningless + continue + fi $COMMAND_PREFIX ./two_particles -N $N -n $n $interactions $COMMAND_PREFIX ./two_particles -N $N -n $n $interactions -E $COMMAND_PREFIX ./two_particles -N $N -n $n $interactions -V + $COMMAND_PREFIX ./two_particles -N $N -n $n $interactions -B done done # The analytical solution is only valid for one particle (second particle has v_z) without interactions @@ -47,9 +53,10 @@ if $NORMAL_MANY_PARTICLE_SIMULATION; then do for n in 10 50 100 do - $COMMAND_PREFIX ./many_particles -N $N -n $n $interactions - $COMMAND_PREFIX ./many_particles -N $N -n $n $interactions -E - $COMMAND_PREFIX ./many_particles -N $N -n $n $interactions -V + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -E + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -V + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -B done done done @@ -66,6 +73,7 @@ if $OSCILLATING_FIELD_SIMULATION; then do $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -O -f $amplitude -w $frequency $reduced_output $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -V -O -f $amplitude -w $frequency $reduced_output + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -B -O -f $amplitude -w $frequency $reduced_output done done interactions="" @@ -75,13 +83,15 @@ if $OSCILLATING_FIELD_SIMULATION; then do $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -O -f $amplitude -w $frequency $reduced_output $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -V -O -f $amplitude -w $frequency $reduced_output + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -B -O -f $amplitude -w $frequency $reduced_output done done fi +if $PUEUE; then + echo "All tasks added to pueue." -echo "All tasks added to pueue." + echo "Use 'pueue status' to check the status of the tasks." -echo "Use 'pueue status' to check the status of the tasks." - -pueue status \ No newline at end of file + pueue status +fi \ No newline at end of file diff --git a/src/project3/include/solvers.hpp b/src/project3/include/solvers.hpp index 32c11cc..c37adfe 100644 --- a/src/project3/include/solvers.hpp +++ b/src/project3/include/solvers.hpp @@ -75,6 +75,13 @@ class VelocityVerletSolver : public Solver { void step() override; }; + +class BorisSolver : public Solver { + public: + BorisSolver(PenningTrap& trap_in, double dt_in); + void step() override; +}; + class AnalyticalSolver : public Solver { private: std::vector w_0; diff --git a/src/project3/many_particles.cpp b/src/project3/many_particles.cpp index 19efa67..0b2ce75 100644 --- a/src/project3/many_particles.cpp +++ b/src/project3/many_particles.cpp @@ -63,6 +63,9 @@ int main(int argc, char* argv[]) { program.add_argument("-V", "--velocity-verlet") .help("Use Velocity Verlet solver instead of RK4") .flag(); + program.add_argument("-B", "--boris") + .help("Use Boris solver instead of RK4") + .flag(); program.add_argument("-t", "--timing") .help("Enable timing of the simulation") .flag(); @@ -83,6 +86,7 @@ int main(int argc, char* argv[]) { bool oscillating_potential = program.get("--oscillating-potential"); bool use_euler = program.get("--euler"); bool use_velocity_verlet = program.get("--velocity-verlet"); + bool use_boris = program.get("--boris"); bool timing = program.get("--timing"); double B_0 = 1.0; // Tesla @@ -125,6 +129,9 @@ int main(int argc, char* argv[]) { } else if (use_velocity_verlet) { solver = make_unique(trap, dt); solver_name = "VelocityVerlet"; + } else if (use_boris) { + solver = make_unique(trap, dt); + solver_name = "Boris"; } else { solver = make_unique(trap, dt); solver_name = "RK4"; diff --git a/src/project3/python/plot_energies.py b/src/project3/python/plot_energies.py index 0b7b959..311a21e 100644 --- a/src/project3/python/plot_energies.py +++ b/src/project3/python/plot_energies.py @@ -2,8 +2,41 @@ import matplotlib.pyplot as plt import numpy as np import os import typer +from pathlib import Path +from typing import List app = typer.Typer() +SHOW_PLOTS = os.getenv("SHOW_PLOTS", "1") == "1" + + +def get_rc_params(): + colors = ["FF220C", "70D6FF", "8AAA79", "666370", "1C1F33"] + rcParams = plt.rcParams + # Use LaTeX for rendering + # Setup fonts + rcParams["text.usetex"] = True + rcParams["font.family"] = "serif" + rcParams["font.size"] = 12 + rcParams["axes.labelsize"] = 12 + rcParams["axes.titlesize"] = 12 + rcParams["legend.fontsize"] = 10 + rcParams["xtick.labelsize"] = 10 + rcParams["ytick.labelsize"] = 10 + # Figure size and resolution + rcParams["figure.figsize"] = (6, 4) + rcParams["figure.dpi"] = 300 + # Use colors from the palette + rcParams["axes.prop_cycle"] = plt.cycler(color=[f"#{color}" for color in colors]) + # Grid + rcParams["axes.grid"] = True + rcParams["grid.alpha"] = 0.5 + rcParams["grid.linestyle"] = "--" + # Point ticks to the inside of the axes + rcParams["xtick.direction"] = "in" + rcParams["ytick.direction"] = "in" + rcParams["xtick.top"] = True + rcParams["ytick.right"] = True + return rcParams def load_data(total_energy_file): total_energy = np.loadtxt(total_energy_file) @@ -15,20 +48,25 @@ def load_data(total_energy_file): return total_energy, kinetic_energy, potential_energy, step_count @app.command() -def plot_kinetic_energy(energy_files: list[str]=[], labels: list[str]=[], mass: float=1.0): +def plot_kinetic_energy(total_energy_files: List[Path], output_file: Path): """Plot kinetic energy from velocity files.""" + plt.rcParams.update(get_rc_params()) + fig, ax = plt.subplots() + energy_files = [str(fn) for fn in total_energy_files if fn.is_file()] + labels = [f"{fn.split('/')[-2]} ({fn.split('/')[-1].split('_')[4]})" for fn in energy_files] for i, energy_file in enumerate(energy_files): total_energy, kinetic_energy, potential_energy, step_count = load_data(energy_file) time = np.linspace(0, 50e-6, step_count) color = f"C{i % 10}" - plt.plot(time, total_energy, label=labels[i], color=color, linestyle="-") - plt.plot(time, potential_energy, color=color, linestyle="--") - plt.plot(time, kinetic_energy, color=color, linestyle=":") - plt.xlabel('Time (s)') - plt.ylabel('Energy (J)') - plt.title('Energy Evolution') - plt.legend() - plt.show() + ax.plot(time, total_energy, label=labels[i], color=color, linestyle="-") + ax.plot(time, potential_energy, color=color, linestyle="--") + ax.plot(time, kinetic_energy, color=color, linestyle=":") + ax.set_xlabel("Time (s)") + ax.set_ylabel(r"Energy (Da $\cdot$ m$^2$/s$^2$)") + ax.legend() + fig.savefig(output_file) + if SHOW_PLOTS: + fig.show() if __name__ == "__main__": diff --git a/src/project3/python/plot_resonances.py b/src/project3/python/plot_resonances.py index 359acbb..5e1d93f 100644 --- a/src/project3/python/plot_resonances.py +++ b/src/project3/python/plot_resonances.py @@ -2,12 +2,52 @@ import matplotlib.pyplot as plt import numpy as np import os import argparse +from pathlib import Path argparser = argparse.ArgumentParser() argparser.add_argument("directory", type=str, help="Directory containing the data files") args = argparser.parse_args() directory = args.directory +def get_rc_params(): + colors = ["FF220C", "70D6FF", "8AAA79", "666370", "1C1F33"] + rcParams = plt.rcParams + # Use LaTeX for rendering + # Setup fonts + rcParams["text.usetex"] = True + rcParams["font.family"] = "serif" + rcParams["font.size"] = 12 + rcParams["axes.labelsize"] = 12 + rcParams["axes.titlesize"] = 12 + rcParams["legend.fontsize"] = 10 + rcParams["xtick.labelsize"] = 10 + rcParams["ytick.labelsize"] = 10 + # Figure size and resolution + rcParams["figure.figsize"] = (6, 4) + rcParams["figure.dpi"] = 300 + # Use colors from the palette + rcParams["axes.prop_cycle"] = plt.cycler(color=[f"#{color}" for color in colors]) + # Grid + rcParams["axes.grid"] = True + rcParams["grid.alpha"] = 0.5 + rcParams["grid.linestyle"] = "--" + # Point ticks to the inside of the axes + rcParams["xtick.direction"] = "in" + rcParams["ytick.direction"] = "in" + rcParams["xtick.top"] = True + rcParams["ytick.right"] = True + return rcParams + +SHOW_PLOTS = os.getenv("SHOW_PLOTS", "1") == "1" + +def get_filename(path: Path, plottername: str) -> str: + general_path = Path(__file__).parent / "plots" + plotter_dir = general_path / plottername + plotter_dir.mkdir(parents=True, exist_ok=True) + filename = plotter_dir / Path(Path(path).parent.name + "_ " + Path(path).name.replace(".dat", f"_{plottername}.pdf")) + return str(filename) + + def get_filenames(directory): filenames = [os.path.join(directory, f) for f in os.listdir(directory) if f.endswith(".red")] return filenames @@ -23,18 +63,21 @@ def load_data(filename): return x def plot_resonances(amplitudes, omegas, counts): - plt.figure(figsize=(8, 6)) + plt.rcParams.update(get_rc_params()) + fig, ax = plt.subplots() unique_amps = sorted(set(amplitudes)) for amp in unique_amps: mask = np.array(amplitudes) == amp omegas_subset = np.array(omegas)[mask] counts_subset = np.array(counts)[mask] omegas_sorted, counts_sorted = zip(*sorted(zip(omegas_subset, counts_subset))) - plt.plot(omegas_sorted, counts_sorted, label=f'Amplitude {amp}') - plt.xlabel('Omega') - plt.ylabel('Count inside trap') - plt.legend() - plt.show() + ax.plot(omegas_sorted, counts_sorted, label=f'Amplitude {amp}') + ax.set_xlabel('Omega') + ax.set_ylabel('Count inside trap') + ax.legend() + fig.savefig(get_filename(Path(directory), "resonances")) + if SHOW_PLOTS: + fig.show() if __name__ == "__main__": filenames = get_filenames(directory) diff --git a/src/project3/python/plot_trajectories.py b/src/project3/python/plot_trajectories.py index 620924d..681db46 100644 --- a/src/project3/python/plot_trajectories.py +++ b/src/project3/python/plot_trajectories.py @@ -2,14 +2,53 @@ import matplotlib.pyplot as plt import numpy as np import os import typer -from typing import Literal +from typing import Literal, List +from pathlib import Path app = typer.Typer() -def load_trajectory_data(filename): +SHOW_PLOTS = os.getenv("SHOW_PLOTS", "1") == "1" + +def get_rc_params(): + colors = ["FF220C", "70D6FF", "8AAA79", "666370", "1C1F33"] + rcParams = plt.rcParams + # Use LaTeX for rendering + # Setup fonts + rcParams["text.usetex"] = True + rcParams["font.family"] = "serif" + rcParams["font.size"] = 12 + rcParams["axes.labelsize"] = 12 + rcParams["axes.titlesize"] = 12 + rcParams["legend.fontsize"] = 10 + rcParams["xtick.labelsize"] = 10 + rcParams["ytick.labelsize"] = 10 + # Figure size and resolution + rcParams["figure.figsize"] = (6, 4) + rcParams["figure.dpi"] = 300 + # Use colors from the palette + rcParams["axes.prop_cycle"] = plt.cycler(color=[f"#{color}" for color in colors]) + # Grid + rcParams["axes.grid"] = True + rcParams["grid.alpha"] = 0.5 + rcParams["grid.linestyle"] = "--" + # Point ticks to the inside of the axes + rcParams["xtick.direction"] = "in" + rcParams["ytick.direction"] = "in" + rcParams["xtick.top"] = True + rcParams["ytick.right"] = True + return rcParams + +def get_filename(path: Path, plottername: str) -> str: + general_path = Path(__file__).parent / "plots" + plotter_dir = general_path / plottername + plotter_dir.mkdir(parents=True, exist_ok=True) + filename = plotter_dir / Path(path).name.replace(".dat", f"_{plottername}.pdf") + return str(filename) + +def load_trajectory_data(filename: Path): data = np.loadtxt(filename) - step_count = int(float(filename.split("/")[-1].split("_")[5])) - particle_count = int(float(filename.split("/")[-1].split("_")[2])) + step_count = int(float(filename.name.split("_")[5])) + particle_count = int(float(filename.name.split("_")[2])) single_particle = particle_count == 1 if single_particle: data = data.reshape((data.shape[0], 1)) @@ -17,77 +56,93 @@ def load_trajectory_data(filename): return data, step_count, particle_count @app.command() -def plot_z_evolution(trajectory_file: str): +def plot_z_evolution(trajectory_file: Path): """Plot the z-coordinate evolution from a trajectory file.""" + plt.rcParams.update(get_rc_params()) + fig, ax = plt.subplots() data, step_count, particle_count = load_trajectory_data(trajectory_file) time = np.linspace(0, 50e-6, step_count) for i in range(particle_count): - plt.plot(time, data[2, i, :] * 1e-6, label=f'Particle {i+1}') - plt.xlabel('Time (s)') - plt.ylabel('Z Position (m)') - plt.title('Z Coordinate Evolution') - plt.legend() - plt.show() + ax.plot(time, data[2, i, :] * 1e-6, label=f'Particle {i+1}') + ax.set_xlabel('Time (s)') + ax.set_ylabel('Z Position (m)') + ax.set_title('Z Coordinate Evolution') + ax.legend() + fig.savefig(get_filename(trajectory_file, "z_evolution")) + if SHOW_PLOTS: + fig.show() @app.command() -def plot_xy_trajectory(trajectory_file: str, plot_start_end: bool = True, plot_time: bool = False): +def plot_xy_trajectory(trajectory_file: Path, plot_start_end: bool = True, plot_time: bool = False): """Plot the XY trajectory from a trajectory file.""" + plt.rcParams.update(get_rc_params()) + fig, ax = plt.subplots() data, step_count, particle_count = load_trajectory_data(trajectory_file) for i in range(particle_count): if plot_start_end: - plt.scatter(data[0, i, 0] * 1e-6, data[1, i, 0] * 1e-6, marker="o", color="black", label="Start" if i == 0 else "") - plt.scatter(data[0, i, -1] * 1e-6, data[1, i, -1] * 1e-6, marker="x", color="red", label="End" if i == 0 else "") + ax.scatter(data[0, i, 0] * 1e-6, data[1, i, 0] * 1e-6, marker="o", color="black", label="Start" if i == 0 else "") + ax.scatter(data[0, i, -1] * 1e-6, data[1, i, -1] * 1e-6, marker="x", color="red", label="End" if i == 0 else "") if plot_time: marker = [".", ",", "o", "v", "^", "<", ">", "1", "2", "3", "4", "8"][i % 12] cm = plt.get_cmap('viridis') c = cm(np.linspace(0, 1, step_count)) - plt.scatter(data[0, i, :] * 1e-6, data[1, i, :] * 1e-6, label=f'Particle {i+1}', marker=marker, c=c) + ax.scatter(data[0, i, :] * 1e-6, data[1, i, :] * 1e-6, label=f'Particle {i+1}', marker=marker, c=c) else: ls = "-" c = f"C{i % 10}" - plt.plot(data[0, i, :] * 1e-6, data[1, i, :] * 1e-6, label=f'Particle {i+1}', linestyle=ls, color=c) - plt.xlabel('X Position (m)') - plt.ylabel('Y Position (m)') - plt.title('XY Trajectory') - plt.axis('equal') + ax.plot(data[0, i, :] * 1e-6, data[1, i, :] * 1e-6, label=f'Particle {i+1}', linestyle=ls, color=c) + ax.set_xlabel('X Position (m)') + ax.set_ylabel('Y Position (m)') + ax.set_title('XY Trajectory') + ax.axis('equal') if plot_time: #sm = plt.cm.ScalarMappable(cmap=cm, norm=plt.Normalize(vmin=0, vmax=50e-6)) #sm.set_array([]) #cbar = plt.colorbar(sm) #cbar.set_label('Time') pass - plt.legend() - plt.show() + ax.legend() + fig.savefig(get_filename(trajectory_file, "xy_trajectory")) + if SHOW_PLOTS: + fig.show() @app.command() -def phase_space_plot(trajectory_file: str, velocity_file: str, axis: Literal['x', 'y', 'z'] = 'z'): +def phase_space_plot(trajectory_file: Path, velocity_file: Path, axis: Literal['x', 'y', 'z'] = 'z'): """Plot phase space (Z vs Vz) from trajectory and velocity files.""" + plt.rcParams.update(get_rc_params()) + fig, ax = plt.subplots() axis_index = {'x': 0, 'y': 1, 'z': 2}[axis] traj_data, step_count, particle_count = load_trajectory_data(trajectory_file) vel_data, _, _ = load_trajectory_data(velocity_file) for i in range(particle_count): - plt.plot(traj_data[axis_index, i, :] * 1e-6, vel_data[axis_index, i, :] * 1e-6, label=f'Particle {i+1}') - plt.xlabel(f'{axis.upper()} Position (m)') - plt.ylabel(f'{axis.upper()} Velocity (m/s)') - plt.title(f'Phase Space Plot ({axis.upper()} vs V{axis.upper()})') - plt.legend() - plt.show() + ax.plot(traj_data[axis_index, i, :] * 1e-6, vel_data[axis_index, i, :] * 1e-6, label=f'Particle {i+1}') + ax.set_xlabel(f'{axis.upper()} Position (m)') + ax.set_ylabel(f'{axis.upper()} Velocity (m/s)') + ax.set_title(f'Phase Space Plot ({axis.upper()} vs V{axis.upper()})') + ax.legend() + fig.savefig(get_filename(trajectory_file, f"phase_space_{axis}")) + if SHOW_PLOTS: + fig.show() @app.command() -def plot_error(trajectory_files: list[str], analytical_file: str): +def plot_error(trajectory_files: list[Path], analytical_file: Path): """Plot error between numerical and analytical trajectories.""" + plt.rcParams.update(get_rc_params()) + fig, ax = plt.subplots() analytical_data, step_count, particle_count = load_trajectory_data(analytical_file) time = np.linspace(0, 50e-6, step_count) for traj_file in trajectory_files: numerical_data, _, _ = load_trajectory_data(traj_file) error = numerical_data - analytical_data relative_error = np.mean(np.linalg.norm(error, axis=0) / np.linalg.norm(analytical_data, axis=0), axis=0) - plt.plot(time, relative_error, label=os.path.basename(traj_file)) - plt.xlabel('Time (s)') - plt.ylabel('Relative Error') - plt.title('Relative Error between Numerical and Analytical Trajectories') - plt.legend() - plt.show() + ax.plot(time, relative_error, label=os.path.basename(traj_file)) + ax.set_xlabel('Time (s)') + ax.set_ylabel('Relative Error') + ax.set_title('Relative Error between Numerical and Analytical Trajectories') + ax.legend() + fig.savefig(get_filename(analytical_file, "error_plot")) + if SHOW_PLOTS: + fig.show() if __name__ == "__main__": app() \ No newline at end of file diff --git a/src/project3/src/solvers.cpp b/src/project3/src/solvers.cpp index 5a88cdc..524fa59 100644 --- a/src/project3/src/solvers.cpp +++ b/src/project3/src/solvers.cpp @@ -294,6 +294,26 @@ AnalyticalSolver::AnalyticalSolver(PenningTrap& trap_in, double dt_in) } } +BorisSolver::BorisSolver(PenningTrap& trap_in, double dt_in) + : Solver(trap_in, dt_in) {} + +void BorisSolver::step() { + int n = trap.num_particles(); + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + arma::vec E = trap.external_E_field(p.position, steps * dt); + E += trap.total_force_particles(i) / p.charge; // Add Coulomb force as effective E-field + arma::vec B = trap.external_B_field(p.position, steps * dt); + arma:: vec t = (p.charge * B / p.mass) * (dt / 2.0); + arma::vec s = 2.0 * t / (1.0 + arma::dot(t, t)); + arma::vec v_minus = p.velocity + (p.charge * E / p.mass) * (dt / 2.0); + arma::vec v_prime = v_minus + arma::cross(v_minus, t); + arma::vec v_plus = v_minus + arma::cross(v_prime, s); + p.velocity = v_plus + (p.charge * E / p.mass) * (dt / 2.0); + p.position += p.velocity * dt; + } +} + // Perform a single Analytical "integration" step void AnalyticalSolver::step() { int n = trap.num_particles(); diff --git a/src/project3/two_particles.cpp b/src/project3/two_particles.cpp index f247d2d..f06ac23 100644 --- a/src/project3/two_particles.cpp +++ b/src/project3/two_particles.cpp @@ -31,6 +31,9 @@ int main(int argc, char* argv[]) { program.add_argument("-V", "--velocity-verlet") .help("Use Velocity Verlet solver instead of RK4") .flag(); + program.add_argument("-B", "--boris") + .help("Use Boris solver instead of RK4") + .flag(); try { program.parse_args(argc, argv); @@ -46,6 +49,7 @@ int main(int argc, char* argv[]) { bool use_euler = program.get("--euler"); bool use_analytical = program.get("--analytical"); bool use_velocity_verlet = program.get("--velocity-verlet"); + bool use_boris = program.get("--boris"); double B_0 = 1.0; // Tesla double V_0 = 0.025; // Volt @@ -71,52 +75,29 @@ int main(int argc, char* argv[]) { double dt = 50e-6 / N; // seconds + unique_ptr solver; + string solver_name; + if (use_euler) { - EulerSolver solver(trap, dt); - solver.simulate(N); - solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_euler"); - vector positions = solver.get_positions(); - cout << "Final position of particle 1: " << endl; - cout << positions.back().col(0) << endl; - if (n_particles > 1) { - cout << "Final position of particle 2: " << endl; - cout << positions.back().col(1) << endl; - } - return 0; + solver = make_unique(trap, dt); + solver_name = "euler"; + } else if (use_velocity_verlet) { + solver = make_unique(trap, dt); + solver_name = "velocity_verlet"; + } else if (use_boris) { + solver = make_unique(trap, dt); + solver_name = "boris"; + } else if (use_analytical) { + solver = make_unique(trap, dt); + solver_name = "analytical"; + } else { + solver = make_unique(trap, dt); + solver_name = "rk4"; } - if (use_analytical) { - AnalyticalSolver solver(trap, dt); - solver.simulate(N); - solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_analytical"); - vector positions = solver.get_positions(); - cout << "Final position of particle 1: " << endl; - cout << positions.back().col(0) << endl; - if (n_particles > 1) { - cout << "Final position of particle 2: " << endl; - cout << positions.back().col(1) << endl; - } - return 0; - } - - if (use_velocity_verlet) { - VelocityVerletSolver solver(trap, dt); - solver.simulate(N); - solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_velocity_verlet"); - vector positions = solver.get_positions(); - cout << "Final position of particle 1: " << endl; - cout << positions.back().col(0) << endl; - if (n_particles > 1) { - cout << "Final position of particle 2: " << endl; - cout << positions.back().col(1) << endl; - } - return 0; - } - - RK4Solver solver(trap, dt); - solver.simulate(N); - solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_rk4"); - vector positions = solver.get_positions(); + solver->simulate(N); + solver->save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_" + solver_name); + vector positions = solver->get_positions(); cout << "Final position of particle 1: " << endl; cout << positions.back().col(0) << endl; if (n_particles > 1) {