From eb2a985b3b23a3f40842517c10f84312ceb73067 Mon Sep 17 00:00:00 2001 From: Cory Date: Sun, 22 Oct 2023 03:47:23 +0200 Subject: [PATCH] Bunch of stuff - Add reinitialize to reuse PenningTrap instance and its particles - Reformatting - Refactor fraction_of_particles_left to not use simulate --- include/PenningTrap.hpp | 31 +++++++---- src/PenningTrap.cpp | 117 ++++++++++++++++++++++++---------------- 2 files changed, 92 insertions(+), 56 deletions(-) diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index 3afcde8..41af4ab 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -20,7 +20,7 @@ #include "typedefs.hpp" #include "utils.hpp" -#pragma omp declare reduction(+ : vec_3d : omp_out += omp_in) \ +#pragma omp declare reduction(+ : vec_3d : omp_out += omp_in) \ initializer(omp_priv = omp_orig) /** @brief A class that simulates a Penning trap. @@ -74,8 +74,8 @@ public: * */ PenningTrap( double B_0 = T, - std::function V_0 = - [](double t) { return 25. * V / 1000.; }, + std::function V_0 + = [](double t) { return 25. * V / 1000.; }, double d = 500., double t = 0.); /** @brief Constructor for the PenningTrap class @@ -88,8 +88,8 @@ public: * */ PenningTrap( unsigned int i, double B_0 = T, - std::function V_0 = - [](double t) { return 25. * V / 1000.; }, + std::function V_0 + = [](double t) { return 25. * V / 1000.; }, double d = 500., double t = 0.); /** @brief Constructor for the PenningTrap class @@ -102,10 +102,21 @@ public: * */ PenningTrap( std::vector particles, double B_0 = T, - std::function V_0 = - [](double t) { return 25. * V / 1000.; }, + std::function V_0 + = [](double t) { return 25. * V / 1000.; }, double d = 500., double t = 0.); + /** @brief Give all particles new positions and velocities, and change t + * and V_0. + * + * @param V_0 The tiome dependent applied potential + * @param t The starting time + * */ + void reinitialize( + std::function V_0 + = [](double t) { return 25. * V / 1000.; }, + double t = 0.); + /** @brief Add a particle to the system * * @param particle The particle to add to the Penning trap @@ -151,7 +162,8 @@ public: * */ vec_3d total_force_external(unsigned int i); - /** @brief Calculate the total force on a particle p_i from other particles. + /** @brief Calculate the total force on a particle p_i from other + * particles. * * @param i The index of particle p_i * @@ -204,7 +216,8 @@ public: * @param particle_interaction Turn particle interactions on/off * */ void write_simulation_to_dir(std::string path, double time, - unsigned int steps, std::string method = "rk4", + unsigned int steps, + std::string method = "rk4", bool particle_interaction = true); /** @brief Simulate and calculate what fraction of particles are still diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index d23f063..d018778 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -11,7 +11,10 @@ * */ #include "PenningTrap.hpp" +#include "typedefs.hpp" +#include #include +#include PenningTrap::PenningTrap(double B_0, std::function V_0, double d, double t) @@ -26,11 +29,9 @@ PenningTrap::PenningTrap(unsigned int i, double B_0, std::function V_0, double d, double t) : PenningTrap::PenningTrap(B_0, V_0, d) { - vec_3d r, v; for (size_t j = 0; j < i; j++) { - r = vec_3d().randn() * .1 * this->d; - v = vec_3d().randn() * .1 * this->d; - this->add_particle(Particle(1., 40., r, v)); + this->particles.push_back(Particle(1., 40., vec_3d(vec_3d().randn() * .1 * this->d), + vec_3d(vec_3d().randn() * .1 * this->d))); } } @@ -41,6 +42,16 @@ PenningTrap::PenningTrap(std::vector particles, double B_0, this->particles = particles; } +void PenningTrap::reinitialize(std::function V_0, double t) +{ + this->V_0 = V_0; + this->t = t; + + for (size_t i = 0; i < this->particles.size(); i++) { + this->particles[i].r_vec = vec_3d().randn() * .1 * this->d; + } +} + vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt) { switch (i) { @@ -109,15 +120,16 @@ vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j) vec_3d PenningTrap::total_force_external(unsigned int i) { - Particle p = this->particles[i]; + Particle *p = &this->particles[i]; - if (arma::norm(p.r_vec) > this->d) { + if (arma::norm(p->r_vec) > this->d) { return vec_3d{0., 0., 0.}; } - return vec_3d(p.q - * (this->external_E_field(p.r_vec) - + arma::cross(p.v_vec, this->external_B_field(p.r_vec)))); + return vec_3d( + p->q + * (this->external_E_field(p->r_vec) + + arma::cross(p->v_vec, this->external_B_field(p->r_vec)))); } vec_3d PenningTrap::total_force_particles(unsigned int i) @@ -134,10 +146,11 @@ vec_3d PenningTrap::total_force_particles(unsigned int i) vec_3d PenningTrap::total_force(unsigned int i) { - return arma::norm(this->particles[i].r_vec) > this->d - ? vec_3d{0., 0., 0.} - : vec_3d(this->total_force_external(i) - - this->total_force_particles(i)); + if (arma::norm(this->particles[i].r_vec) > this->d) { + return vec_3d{0., 0., 0.}; + } + return vec_3d(this->total_force_external(i) + - this->total_force_particles(i)); } void PenningTrap::evolve_RK4(double dt, bool particle_interaction) @@ -146,13 +159,9 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) std::vector original_particles = this->particles; std::vector tmp_particles = this->particles; - vec_3d (PenningTrap::*force)(unsigned int); - if (particle_interaction) { - force = &PenningTrap::total_force; - } - else { - force = &PenningTrap::total_force_external; - } + vec_3d (PenningTrap::*force)(unsigned int) + = particle_interaction ? &PenningTrap::total_force + : &PenningTrap::total_force_external; size_t size = this->particles.size(); @@ -163,7 +172,9 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) this->k_r = sim_arr(4, sim_cols(size)); } + // Each k_{i+1} is dependent on k_i, so outer loop is not parallelizable for (size_t i = 0; i < 4; i++) { + // Inner loop is able to be parallelized #pragma omp parallel for for (size_t j = 0; j < size; j++) { this->k_v[i][j] = (this->*force)(j) / this->particles[j].m; @@ -176,6 +187,7 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) } this->particles = tmp_particles; } + this->t += dt; } @@ -185,19 +197,19 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) vec_3d force_res[size]; Particle *p; - vec_3d (PenningTrap::*force)(unsigned int); - if (particle_interaction) { - force = &PenningTrap::total_force_external; - } - else { - force = &PenningTrap::total_force_external; - } + vec_3d (PenningTrap::*force)(unsigned int) + = particle_interaction ? &PenningTrap::total_force + : &PenningTrap::total_force_external; + // Calculating the force for each particle is independent and therefore + // a good candidate for parallel execution #pragma omp parallel for for (size_t i = 0; i < size; i++) { force_res[i] = (this->*force)(i); } + // Updating the particles is also independent, so we can parallelize + // this as well #pragma omp parallel for private(p) for (size_t i = 0; i < size; i++) { p = &this->particles[i]; @@ -251,8 +263,8 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time, path += '/'; } if (mkpath(path, 0777) != 0) { - std::cout << "Hello" << std::endl; - return; + std::cout << "Failed to make path" << std::endl; + abort(); } simulation_t res @@ -260,18 +272,22 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time, std::ofstream ofile; + // Writing each particle to its own file is independent and can be run in + // parallel. #pragma omp parallel for private(ofile) for (size_t i = 0; i < this->particles.size(); i++) { ofile.open(path + "particle_" + std::to_string(i) + "_r.txt"); for (vec_3d &vec : res.r_vecs[i]) { - ofile << vec(0) << "," << vec(1) << "," << vec(2) << "\n"; + ofile << scientific_format(vec(0), 10, 8) << ',' + << scientific_format(vec(1), 10, 8) << ',' + << scientific_format(vec(2), 10, 8) << '\n'; } ofile.close(); ofile.open(path + "particle_" + std::to_string(i) + "_v.txt"); for (vec_3d &vec : res.v_vecs[i]) { - ofile << scientific_format(vec(0), 10, 8) << "," - << scientific_format(vec(1), 8, 10) << "," - << scientific_format(vec(2), 8, 10) << "\n"; + ofile << scientific_format(vec(0), 10, 8) << ',' + << scientific_format(vec(1), 10, 8) << ',' + << scientific_format(vec(2), 10, 8) << '\n'; } ofile.close(); } @@ -281,26 +297,33 @@ double PenningTrap::fraction_of_particles_left(double time, unsigned int steps, std::string method, bool particle_interaction) { - simulation_t res - = this->simulate(time, steps, method, particle_interaction); + double dt = time / (double)steps; + + void (PenningTrap::*func)(double, bool); + if (method == "rk4") { + func = &PenningTrap::evolve_RK4; + } + else if (method == "euler") { + func = &PenningTrap::evolve_forward_euler; + } + else { + std::cout << "Not a valid method!" << std::endl; + abort(); + } + + for (size_t j = 0; j < steps; j++) { + (this->*func)(dt, particle_interaction); + } int particles_left = 0; - for (Particle p : this->particles) { - if (arma::norm(p.r_vec) < this->d) { + // A reduction is perfect here + #pragma omp parallel for reduction(+:particles_left) + for (size_t i=0; i < this->particles.size(); i++) { + if (arma::norm(this->particles[i].r_vec) < this->d) { particles_left++; } } return (double)particles_left / (double)this->particles.size(); } - -vec_3d PenningTrap::get_r(int i) -{ - return this->particles[i].r_vec; -} - -double PenningTrap::get_t() -{ - return this->t; -}