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/src/Makefile b/src/Makefile index cb6d014..2ebf10f 100644 --- a/src/Makefile +++ b/src/Makefile @@ -46,6 +46,10 @@ instrument: main: main.o $(LIBOBJS) $(CLASSOBJS) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) +main_no_interface: main_no_interface.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/job.script b/src/job.script new file mode 100644 index 0000000..c87ea6d --- /dev/null +++ b/src/job.script @@ -0,0 +1,13 @@ +#!/bin/bash + +#SBATCH --account=ec54 +#SBATCH --job-name=simple +#SBATCH --time=0-00:05:00 +#SBATCH --mem-per-cpu=4G +#SBATCH --cpus-per-task=8 + +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 + +srun ./main_no_interface diff --git a/src/main_no_interface.cpp b/src/main_no_interface.cpp new file mode 100644 index 0000000..6b00410 --- /dev/null +++ b/src/main_no_interface.cpp @@ -0,0 +1,334 @@ +/** @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 + +// Particles used for testing +Particle p1(vec3{20., 0., 20.}, vec3{0., 25., 0.}); ///< Particle 1 +Particle p2(vec3{25., 25., 0.}, vec3{0., 40., 5.}); ///< Particle 2 + +/** @brief The analytical solution for particle p1 + * + * @param t Time + * + * @return vec3 + * */ +vec3 analytical_solution_particle_1(double t) +{ + double w_0 = T / CA_MASS; + double w_z2 = (50. * V / 1000.) / (CA_MASS * 500. * 500.); + double w_p = (w_0 + std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.; + 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::complex f = + A_p * std::exp(std::complex(0., -w_p * t)) + + A_n * std::exp(std::complex(0., -w_n * t)); + vec3 res{std::real(f), std::imag(f), 20. * std::cos(std::sqrt(w_z2) * t)}; + return res; +} + +/** @brief Simulate a single particle over the period of 50 \f$ \mu s \f$. + * */ +void simulate_single_particle() +{ + // Initialize trap with particle 1 + PenningTrap trap(std::vector{p1}); + + double time = 50.; // microseconds + + // Simulate and write results to file + trap.write_simulation_to_dir("output/simulate_single_particle", time, N, + "rk4", false); +} + +/** @brief Simulate 2 particles over the period of 50 \f$ \mu s \f$ with and + * without particle interactions. + * */ +void simulate_two_particles() +{ + // Initialize traps with particles + PenningTrap trap_no_interaction(std::vector{p1, p2}); + PenningTrap trap_with_interaction(std::vector{p1, p2}); + + double time = 50.; // microseconds + + // Simulate and write results to files + trap_no_interaction.write_simulation_to_dir( + "output/simulate_2_particles/no_interaction", time, N, "rk4", false); + trap_with_interaction.write_simulation_to_dir( + "output/simulate_2_particles/with_interaction", time, N); +} + +/** @brief Simulate a single particle over 50 \f$ \mu s \f$ using different + * amount of steps and different methods. + * */ +void simulate_single_particle_with_different_steps() +{ + double time = 50.; // microseconds + + std::ofstream ofile; + + // Calculate relative error for RK4 + std::string path = "output/relative_error/RK4/"; + mkpath(path); +#pragma omp parallel for private(ofile) + for (int i = 0; i < 4; i++) { + int steps = 4000 * std::pow(2, i); + std::cout << steps << std::endl; + double dt = time / (double)steps; + ofile.open(path + std::to_string(steps) + "_steps.txt"); + PenningTrap trap(std::vector{p1}); + simulation_t res = trap.simulate(time, steps, "rk4", false); + for (int i = 0; i < steps; i++) { + ofile << arma::norm(res.r_vecs[0][i] + - analytical_solution_particle_1(dt * i)) + << '\n'; + } + ofile.close(); + } + + // Calculate relative error for forward Euler + path = "output/relative_error/euler/"; + mkpath(path); +#pragma omp parallel for private(ofile) + for (int i = 0; i < 4; i++) { + int steps = 4000 * std::pow(2, i); + double dt = time / (double)steps; + ofile.open(path + std::to_string(steps) + "_steps.txt"); + PenningTrap trap(std::vector{p1}); + simulation_t res = trap.simulate(time, steps, "euler", false); + for (int i = 0; i < steps; i++) { + ofile << arma::norm(res.r_vecs[0][i] + - analytical_solution_particle_1(dt * i)) + << '\n'; + } + ofile.close(); + } +} + +/** @brief Simulate 100 particles over 50 \f$ \mu s \f$. + * */ +void simulate_100_particles() +{ + PenningTrap trap((unsigned)100); + + double time = 50.; // microseconds + + // trap.write_simulation_to_dir("output/simulate_100_particles", time, N, + //"rk4", false); + trap.simulate(time, N, "rk4", true); +} + +/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time + * dependent potential. + * + * @details The simulation sweeps over different frequencies in [0.2, 2.5] + * MHz. + * + * */ +void potential_resonance_wide_sweep() +{ + double time = 500.; + + double amplitudes[]{.1, .4, .7}; + + double freq_start = .2; + double freq_end = 2.5; + double freq_increment = .02; + 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)100); +#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 + "wide_sweep.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() +{ + double time = 500.; + + double amplitudes[]{.1, .4, .7}; + + double freq_start = 1.; + double freq_end = 1.7; + double freq_increment = .002; + 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)100); +#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.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.; + double freq_end = 1.7; + double freq_increment = .002; + 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)100); +#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.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() +{ + int option = 1; + bool chosen = false; + + double start, end; + + start = omp_get_wtime(); + + simulate_100_particles(); + + potential_resonance_wide_sweep(); + + // 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/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])