Further development on project 3
This commit is contained in:
+2
-1
@@ -290,4 +290,5 @@ src/project3/results/*/*.dat
|
||||
src/project3/results/*/*/*.red
|
||||
*.tar.gz
|
||||
*_timing.txt
|
||||
run_stats.json
|
||||
run_stats.json
|
||||
src/project3/python/plots/
|
||||
@@ -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}
|
||||
@@ -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}
|
||||
Binary file not shown.
Executable
+55
@@ -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
|
||||
@@ -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
|
||||
pueue status
|
||||
fi
|
||||
@@ -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<double> w_0;
|
||||
|
||||
@@ -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<bool>("--oscillating-potential");
|
||||
bool use_euler = program.get<bool>("--euler");
|
||||
bool use_velocity_verlet = program.get<bool>("--velocity-verlet");
|
||||
bool use_boris = program.get<bool>("--boris");
|
||||
bool timing = program.get<bool>("--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<VelocityVerletSolver>(trap, dt);
|
||||
solver_name = "VelocityVerlet";
|
||||
} else if (use_boris) {
|
||||
solver = make_unique<BorisSolver>(trap, dt);
|
||||
solver_name = "Boris";
|
||||
} else {
|
||||
solver = make_unique<RK4Solver>(trap, dt);
|
||||
solver_name = "RK4";
|
||||
|
||||
@@ -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__":
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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()
|
||||
@@ -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();
|
||||
|
||||
@@ -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<bool>("--euler");
|
||||
bool use_analytical = program.get<bool>("--analytical");
|
||||
bool use_velocity_verlet = program.get<bool>("--velocity-verlet");
|
||||
bool use_boris = program.get<bool>("--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> 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<arma::mat> 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<EulerSolver>(trap, dt);
|
||||
solver_name = "euler";
|
||||
} else if (use_velocity_verlet) {
|
||||
solver = make_unique<VelocityVerletSolver>(trap, dt);
|
||||
solver_name = "velocity_verlet";
|
||||
} else if (use_boris) {
|
||||
solver = make_unique<BorisSolver>(trap, dt);
|
||||
solver_name = "boris";
|
||||
} else if (use_analytical) {
|
||||
solver = make_unique<AnalyticalSolver>(trap, dt);
|
||||
solver_name = "analytical";
|
||||
} else {
|
||||
solver = make_unique<RK4Solver>(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<arma::mat> 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<arma::mat> 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<arma::mat> 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<arma::mat> positions = solver->get_positions();
|
||||
cout << "Final position of particle 1: " << endl;
|
||||
cout << positions.back().col(0) << endl;
|
||||
if (n_particles > 1) {
|
||||
|
||||
Reference in New Issue
Block a user