diff --git a/.gitignore b/.gitignore index 17c5519..b3816d1 100644 --- a/.gitignore +++ b/.gitignore @@ -48,3 +48,6 @@ src/* !src/*.py !src/Doxyfile !src/scripts + +# Job +!src/job.script diff --git a/images/100_particles.gif b/images/100_particles.gif index 507ff1b..ba16c67 100644 Binary files a/images/100_particles.gif and b/images/100_particles.gif differ diff --git a/latex/images/particles_left_narrow_sweep.pdf b/latex/images/particles_left_narrow_sweep.pdf index f049bdb..d5c2bbc 100644 Binary files a/latex/images/particles_left_narrow_sweep.pdf and b/latex/images/particles_left_narrow_sweep.pdf differ diff --git a/latex/images/particles_left_narrow_sweep_fine.pdf b/latex/images/particles_left_narrow_sweep_fine.pdf new file mode 100644 index 0000000..31da0eb Binary files /dev/null and b/latex/images/particles_left_narrow_sweep_fine.pdf differ diff --git a/latex/images/particles_left_narrow_sweep_interactions.pdf b/latex/images/particles_left_narrow_sweep_interactions.pdf index cb06e28..f68274c 100644 Binary files a/latex/images/particles_left_narrow_sweep_interactions.pdf and b/latex/images/particles_left_narrow_sweep_interactions.pdf differ diff --git a/latex/images/particles_left_narrow_sweep_interactions_fine.pdf b/latex/images/particles_left_narrow_sweep_interactions_fine.pdf new file mode 100644 index 0000000..40d09e3 Binary files /dev/null and b/latex/images/particles_left_narrow_sweep_interactions_fine.pdf differ diff --git a/latex/images/particles_left_wide_sweep.pdf b/latex/images/particles_left_wide_sweep.pdf index 5cf0bff..c757fc8 100644 Binary files a/latex/images/particles_left_wide_sweep.pdf and b/latex/images/particles_left_wide_sweep.pdf differ diff --git a/src/Makefile b/src/Makefile index cb6d014..22514b4 100644 --- a/src/Makefile +++ b/src/Makefile @@ -8,7 +8,8 @@ CLASSOBJS=$(CLASSSRCS:.cpp=.o) INCLUDE=../include -CFLAGS=-Wall -larmadillo -lblas -llapack -std=c++11 -O3 -fomit-frame-pointer +CFLAGS=-Wall -larmadillo -std=c++11 -O3 -fomit-frame-pointer +#CFLAGS=-Wall -larmadillo -lblas -llapack -std=c++11 -O3 -fomit-frame-pointer OPENMP=-fopenmp # Add a debug flag when compiling (For the DEBUG macro in utils.hpp) @@ -46,6 +47,10 @@ instrument: main: main.o $(LIBOBJS) $(CLASSOBJS) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) +frequency_narrow_sweeps_long: frequency_narrow_sweeps_long.o $(LIBOBJS) $(CLASSOBJS) + $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) + + test_suite: test_suite.o $(LIBOBJS) $(CLASSOBJS) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 95a03c4..8f7996f 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -158,7 +158,7 @@ void PenningTrap::reinitialize(double f, double omega_V, double t) for (size_t i = 0; i < this->particles.size(); i++) { p = &this->particles[i]; - p->v_vec = vec3().randn() * .1 * this->d; + p->r_vec = vec3().randn() * .1 * this->d; p->v_vec = vec3().randn() * .1 * this->d; } } @@ -172,6 +172,7 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) { Particle *p; + // Keep original particles std::vector original_particles = this->particles; std::vector tmp_particles = this->particles; diff --git a/src/frequency_narrow_sweeps_long.cpp b/src/frequency_narrow_sweeps_long.cpp new file mode 100644 index 0000000..d22db4b --- /dev/null +++ b/src/frequency_narrow_sweeps_long.cpp @@ -0,0 +1,154 @@ +/** @file main.cpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 1.0 + * + * @brief The main program for this project + * + * @bug No known bugs + * */ + +#include +#include +#include +#include +#include +#include + +#include "PenningTrap.hpp" +#include "constants.hpp" +#include "utils.hpp" + +#define PARTICLES 100 +#define N 40000 + +/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time + * dependent potential. + * + * @details The simulation sweeps over different frequencies in [1., 1.7] + * MHz. + * + * */ +void potential_resonance_narrow_sweep() +{ + double time = 500.; + + double amplitudes[]{.1, .4, .7}; + + double freq_start = 1.1; + double freq_end = 1.7; + double freq_increment = .0005; + size_t freq_iterations = + (size_t)((freq_end - freq_start) / freq_increment) + 1; + + double res[4][freq_iterations]; + + std::string path = "output/time_dependent_potential/"; + mkpath(path); + + std::ofstream ofile; + +#pragma omp parallel for + // Insert frequencies + for (size_t i = 0; i < freq_iterations; i++) { + res[0][i] = freq_start + freq_increment * i; + } + +#pragma omp parallel + { + // Each thread creates a PenningTrap instance and reuses it throughout + // the sweep. + PenningTrap trap((uint)PARTICLES); +#pragma omp for collapse(2) + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < freq_iterations; j++) { + // Reset particles and give new time dependent potential. + trap.reinitialize(amplitudes[i], res[0][j]); + res[i + 1][j] = + trap.fraction_of_particles_left(time, N, "rk4", false); + } + } + } + + // Write results to file + ofile.open(path + "narrow_sweep_fine.txt"); + for (size_t i = 0; i < freq_iterations; i++) { + ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ',' + << res[3][i] << '\n'; + } + ofile.close(); +} + +/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time + * dependent potential. + * + * @details The simulation sweeps over different frequencies in [1., 1.7] + * MHz. + * + * */ +void potential_resonance_narrow_sweep_interaction() +{ + double time = 500.; + + double amplitudes[]{.1, .4, .7}; + + double freq_start = 1.1; + double freq_end = 1.7; + double freq_increment = .0005; + size_t freq_iterations = + (size_t)((freq_end - freq_start) / freq_increment) + 1; + + double res[4][freq_iterations]; + + std::string path = "output/time_dependent_potential/"; + mkpath(path); + + std::ofstream ofile; + +#pragma omp parallel for + for (size_t i = 0; i < freq_iterations; i++) { + res[0][i] = freq_start + freq_increment * i; + } + +#pragma omp parallel + { + // Each thread creates a PenningTrap instance and reuses it throughout + // the sweep. + PenningTrap trap((uint)PARTICLES); +#pragma omp for collapse(2) + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < freq_iterations; j++) { + // Reset particles and give new time dependent potential. + trap.reinitialize(amplitudes[i], res[0][j]); + res[i + 1][j] = trap.fraction_of_particles_left(time, N); + } + } + } + + // Write results to file + ofile.open(path + "narrow_sweep_interactions_fine.txt"); + for (size_t i = 0; i < freq_iterations; i++) { + ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ',' + << res[3][i] << '\n'; + } + ofile.close(); +} + +int main() +{ + double start, end; + + start = omp_get_wtime(); + + potential_resonance_narrow_sweep(); + + potential_resonance_narrow_sweep_interaction(); + + end = omp_get_wtime(); + + std::cout << "Time: " << end - start << " seconds" << std::endl; + + return 0; +} diff --git a/src/job.script b/src/job.script new file mode 100644 index 0000000..1d36bd7 --- /dev/null +++ b/src/job.script @@ -0,0 +1,15 @@ +#!/bin/bash + +#SBATCH --account=ec54 +#SBATCH --job-name=particle-sim +#SBATCH --time=0-01:30:00 +#SBATCH --mem-per-cpu=2G +#SBATCH --cpus-per-task=16 + +set -o errexit # Exit the script on any error +set -o nounset # Treat any unset variables as an error + +module --quiet purge # Reset the modules to the system default +module load Armadillo/11.4.3-foss-2022b + +srun ./frequency_narrow_sweeps_long diff --git a/src/main.cpp b/src/main.cpp index 0537470..e7e982b 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -42,6 +42,7 @@ vec3 analytical_solution_particle_1(double t) double w_n = (w_0 - std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.; double A_p = (25. + w_n * 20.) / (w_n - w_p); double A_n = -(25. + w_p * 20.) / (w_n - w_p); + std::cout << A_p << "," << A_n << std::endl; std::complex f = A_p * std::exp(std::complex(0., -w_p * t)) + A_n * std::exp(std::complex(0., -w_n * t)); @@ -211,7 +212,7 @@ void potential_resonance_narrow_sweep() double amplitudes[]{.1, .4, .7}; - double freq_start = 1.; + double freq_start = 1.1; double freq_end = 1.7; double freq_increment = .002; size_t freq_iterations = @@ -268,7 +269,7 @@ void potential_resonance_narrow_sweep_interaction() double amplitudes[]{.1, .4, .7}; - double freq_start = 1.; + double freq_start = 1.1; double freq_end = 1.7; double freq_increment = .002; size_t freq_iterations = diff --git a/src/scripts/animate_100_particles.py b/src/scripts/animate_100_particles.py index 5508b28..717e31b 100644 --- a/src/scripts/animate_100_particles.py +++ b/src/scripts/animate_100_particles.py @@ -34,7 +34,7 @@ def animate(): arr = get_data([f"output/simulate_100_particles/particle_{i}_r.txt" for i in range(100)]) - arr = arr[:, :, ::10] + arr = arr[:, :, ::40] N = len(arr[0][0]) diff --git a/src/scripts/plot_particles_left.py b/src/scripts/plot_particles_left.py index 125b73d..5116be7 100644 --- a/src/scripts/plot_particles_left.py +++ b/src/scripts/plot_particles_left.py @@ -23,12 +23,16 @@ def main(): files = [ "output/time_dependent_potential/wide_sweep.txt", "output/time_dependent_potential/narrow_sweep.txt", + "output/time_dependent_potential/narrow_sweep_fine.txt", "output/time_dependent_potential/narrow_sweep_interactions.txt", + "output/time_dependent_potential/narrow_sweep_interactions_fine.txt", ] outputs = [ "../latex/images/particles_left_wide_sweep.pdf", "../latex/images/particles_left_narrow_sweep.pdf", + "../latex/images/particles_left_narrow_sweep_fine.pdf", "../latex/images/particles_left_narrow_sweep_interactions.pdf", + "../latex/images/particles_left_narrow_sweep_interactions_fine.pdf", ] for file, output in zip(files, outputs): with open(file) as f: