diff --git a/src/Makefile b/src/Makefile index 2ebf10f..47ec673 100644 --- a/src/Makefile +++ b/src/Makefile @@ -46,7 +46,7 @@ instrument: main: main.o $(LIBOBJS) $(CLASSOBJS) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) -main_no_interface: main_no_interface.o $(LIBOBJS) $(CLASSOBJS) +frequency_narrow_sweeps_long: frequency_narrow_sweeps_long.o $(LIBOBJS) $(CLASSOBJS) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) 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 index c87ea6d..1d36bd7 100644 --- a/src/job.script +++ b/src/job.script @@ -1,13 +1,15 @@ #!/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 +#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 -srun ./main_no_interface +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_no_interface.cpp b/src/main_no_interface.cpp deleted file mode 100644 index 6b00410..0000000 --- a/src/main_no_interface.cpp +++ /dev/null @@ -1,334 +0,0 @@ -/** @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; -}