From ae7b437ac30d564104e8e36e5e798f1fba86627e Mon Sep 17 00:00:00 2001 From: Cory Date: Fri, 13 Oct 2023 15:36:52 +0200 Subject: [PATCH] Fix compiler warnings --- include/Particle.hpp | 14 ++-- include/PenningTrap.hpp | 31 +++++---- include/typedefs.hpp | 12 ++++ src/PenningTrap.cpp | 142 +++++++++++++++++----------------------- src/main.cpp | 4 +- src/test_suite.cpp | 2 +- src/utils.cpp | 2 +- 7 files changed, 99 insertions(+), 108 deletions(-) create mode 100644 include/typedefs.hpp diff --git a/include/Particle.hpp b/include/Particle.hpp index 29b4ccb..cf3232c 100644 --- a/include/Particle.hpp +++ b/include/Particle.hpp @@ -14,14 +14,16 @@ #include +#include "typedefs.hpp" + /** @brief A class that holds attributes of a particle * */ class Particle { private: - double q; ///< Charge - double m; ///< Mass - arma::vec::fixed<3> r_vec; ///< position - arma::vec::fixed<3> v_vec; ///< velocity + double q; ///< Charge + double m; ///< Mass + vec_3d r_vec; ///< position + vec_3d v_vec; ///< velocity public: /** @brief Initialize the particle. @@ -29,9 +31,7 @@ public: * @details Initialize the particle with a charge, mass, position and * velocity. * */ - Particle(double q, double m, - arma::vec::fixed<3> r_vec, - arma::vec::fixed<3> v_vec); + Particle(double q, double m, vec_3d r_vec, vec_3d v_vec); /** @brief Make private attributes available for PenningTrap. * */ diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index b3a32d1..4a131aa 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -17,12 +17,10 @@ #include "Particle.hpp" #include "constants.hpp" +#include "typedefs.hpp" -typedef std::vector> sim_cols; -typedef std::vector> sim_rows; -typedef std::vector sim_arr; -#pragma omp declare reduction(+ : arma::vec : 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. @@ -36,10 +34,10 @@ private: double V_0; ///< Applied potential double d; ///< Characteristic dimension std::vector particles; ///< The particles in the Penning trap - arma::vec::fixed<3> *k_v; - arma::vec::fixed<3> *k_r; - std::function v_funcs[4]; - std::function r_funcs[4]; + sim_arr k_v; + sim_arr k_r; + std::function v_funcs[4]; + std::function r_funcs[4]; public: /** @brief Set B_0, V_0 and d. @@ -58,33 +56,33 @@ public: /** @brief Calculate E at point r * */ - arma::vec external_E_field(arma::vec r); + vec_3d external_E_field(vec_3d r); /** @brief Calculate B at point r * */ - arma::vec external_B_field(arma::vec r); + vec_3d external_B_field(vec_3d r); /** @brief Calculate the force between 2 particles. * * @details Calculate the force exhibited on particle p_i from * particle p_j. * */ - arma::vec force_on_particle(int i, int j); + vec_3d force_on_particle(int i, int j); /** @brief Calculate the total external force on a particle. * * @details Calculate the total amount of force that E and B exhibits * on particle p_i. * */ - arma::vec total_force_external(int i); + vec_3d total_force_external(int i); /** @brief Calculate the total force on a particle from other particles. * */ - arma::vec total_force_particles(int i); + vec_3d total_force_particles(unsigned int i); /** @brief calculate the total force on a particle. * */ - arma::vec total_force(int i); + vec_3d total_force(int i); /** @brief Go forward one timestep using the RK4 method * */ @@ -94,8 +92,9 @@ public: * */ void evolve_forward_euler(double dt, bool particle_interaction = true); - sim_arr simulate(double time, int steps, std::string method = "rk4", - bool particle_interaction = true); + sim_arr simulate(double time, unsigned int steps, + std::string method = "rk4", + bool particle_interaction = true); void write_simulation_to_dir(std::string path, double time, int steps, std::string method = "rk4", diff --git a/include/typedefs.hpp b/include/typedefs.hpp new file mode 100644 index 0000000..42c00bb --- /dev/null +++ b/include/typedefs.hpp @@ -0,0 +1,12 @@ +#ifndef __TYPEDEFS__ +#define __TYPEDEFS__ + +#include +#include + +typedef std::vector> sim_cols; +typedef std::vector> sim_rows; +typedef std::vector sim_arr; +typedef arma::vec::fixed<3> vec_3d; + +#endif diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 93ed3a6..93c6ed7 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -15,6 +15,7 @@ #include "PenningTrap.hpp" #include "constants.hpp" +#include "typedefs.hpp" #include "utils.hpp" #include #include @@ -27,42 +28,31 @@ PenningTrap::PenningTrap(double B_0, double V_0, double d) this->V_0 = V_0; this->d = d; - v_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_v[i]; }; - v_funcs[1] = [this](int i, double dt) { - return .5 * dt * this->k_v[this->particles.size() + i]; - }; - v_funcs[2] = [this](int i, double dt) { - return dt * this->k_v[2 * this->particles.size() + i]; - }; + // Create lambda functions for each iteration in RK4 + v_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_v[0][i]; }; + v_funcs[1] = [this](int i, double dt) { return .5 * dt * this->k_v[1][i]; }; + v_funcs[2] = [this](int i, double dt) { return dt * this->k_v[2][i]; }; v_funcs[3] = [this](int i, double dt) { - arma::vec res = this->k_v[i] + this->k_v[this->particles.size() + i] + - this->k_v[2 * this->particles.size() + i] + - this->k_v[3 * this->particles.size() + i]; - return (dt / 6) * res; + return (dt / 6) * (this->k_v[0][i].eval() + this->k_v[1][i].eval() + + this->k_v[2][i].eval() + this->k_v[3][i].eval()); }; - r_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_r[i]; }; - r_funcs[1] = [this](int i, double dt) { - return .5 * dt * this->k_r[this->particles.size() + i]; - }; - r_funcs[2] = [this](int i, double dt) { - return dt * this->k_r[2 * this->particles.size() + i]; - }; + r_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_r[0][i]; }; + r_funcs[1] = [this](int i, double dt) { return .5 * dt * this->k_r[0][i]; }; + r_funcs[2] = [this](int i, double dt) { return dt * this->k_r[2][i]; }; r_funcs[3] = [this](int i, double dt) { - arma::vec res = this->k_r[i] + this->k_r[this->particles.size() + i] + - this->k_r[2 * this->particles.size() + i] + - this->k_r[3 * this->particles.size() + i]; - return (dt / 6) * res; + return (dt / 6) * (this->k_r[0][i].eval() + this->k_r[1][i].eval() + + this->k_r[2][i].eval() + this->k_r[3][i].eval()); }; } PenningTrap::PenningTrap(int i, double B_0, double V_0, double d) : PenningTrap::PenningTrap(B_0, V_0, d) { - arma::vec r, v; + vec_3d r, v; for (int j = 0; j < i; j++) { - r = arma::vec::fixed<3>().randn() * .1 * this->d; - v = arma::vec::fixed<3>().randn() * .1 * this->d; + r = vec_3d().randn() * .1 * this->d; + v = vec_3d().randn() * .1 * this->d; this->add_particle(Particle(1., 40., r, v)); } } @@ -79,63 +69,60 @@ void PenningTrap::add_particle(Particle particle) this->particles.push_back(particle); } -arma::vec PenningTrap::external_E_field(arma::vec r) +vec_3d PenningTrap::external_E_field(vec_3d r) { - arma::vec::fixed<3> res{r(0), r(1), -2. * r(2)}; + r(2) *= -2.; double f = this->V_0 / (this->d * this->d); - return f * res; + return f * r; } -arma::vec PenningTrap::external_B_field(arma::vec r) +vec_3d PenningTrap::external_B_field(vec_3d r) { - arma::vec::fixed<3> res{0., 0., this->B_0}; - - return res; + return vec_3d{0.,0.,this->B_0}; } -arma::vec PenningTrap::force_on_particle(int i, int j) +vec_3d PenningTrap::force_on_particle(int i, int j) { // Calculate the difference between the particles' position - arma::vec::fixed<3> res = - this->particles.at(i).r_vec - this->particles.at(j).r_vec; + vec_3d res = this->particles[i].r_vec - this->particles[j].r_vec; // Get the distance between the particles double norm = arma::norm(res); // Multiply res with p_j's charge divided by the norm cubed - res *= this->particles.at(j).q / (norm * norm * norm); + res *= this->particles[j].q / (norm * norm * norm); return res; } -arma::vec PenningTrap::total_force_external(int i) +vec_3d PenningTrap::total_force_external(int i) { - if (arma::norm(this->particles.at(i).r_vec) > this->d) { - return arma::vec::fixed<3>{0., 0., 0.}; + Particle p = this->particles[i]; + + if (arma::norm(p.r_vec) > this->d) { + return vec_3d{0., 0., 0.}; } - Particle p = this->particles.at(i); + vec_3d B = this->external_B_field(p.r_vec); - arma::vec::fixed<3> B = this->external_B_field(p.r_vec); + vec_3d v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1), + p.v_vec(2) * B(0) - p.v_vec(0) * B(2), + p.v_vec(0) * B(1) - p.v_vec(1) * B(0)}; - arma::vec::fixed<3> v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1), - p.v_vec(2) * B(0) - p.v_vec(0) * B(2), - p.v_vec(0) * B(1) - p.v_vec(1) * B(0)}; - - arma::vec::fixed<3> force = - p.q * (this->external_E_field(p.r_vec) + v_cross_B); + vec_3d force = p.q * (this->external_E_field(p.r_vec) + v_cross_B); return force; } -arma::vec PenningTrap::total_force_particles(int i) +vec_3d PenningTrap::total_force_particles(unsigned int i) { - Particle p = this->particles.at(i); + Particle p = this->particles[i]; - arma::vec::fixed<3> res; + vec_3d res; + size_t size = this->particles.size(); - for (int j = 0; j < this->particles.size(); j++) { + for (size_t j = 0; j < size; j++) { if (i == j) { continue; } @@ -148,7 +135,7 @@ arma::vec PenningTrap::total_force_particles(int i) return res; } -arma::vec PenningTrap::total_force(int i) +vec_3d PenningTrap::total_force(int i) { return this->total_force_external(i) - this->total_force_particles(i); } @@ -159,7 +146,7 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) std::vector original_particles = this->particles; std::vector tmp_particles = this->particles; - std::function force; + std::function force; if (particle_interaction) { force = [this](int i) { return this->total_force(i); }; } @@ -169,25 +156,22 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) int size = this->particles.size(); - k_v = new arma::vec::fixed<3>[size * 4]; - k_r = new arma::vec::fixed<3>[size * 4]; + this->k_v = sim_arr(4, sim_cols(size)); + this->k_r = sim_arr(4, sim_cols(size)); for (int i = 0; i < 4; i++) { #pragma omp parallel for for (int j = 0; j < size; j++) { - k_v[i * size + j] = this->total_force(j) / this->particles.at(j).m; - k_r[i * size + j] = this->particles.at(j).v_vec; + this->k_v[i][j] = this->total_force(j) / this->particles[j].m; + this->k_r[i][j] = this->particles[j].v_vec; - Particle *p = &tmp_particles.at(j); + Particle *p = &tmp_particles[j]; - p->v_vec = original_particles.at(j).v_vec + this->v_funcs[i](j, dt); - p->r_vec = original_particles.at(j).r_vec + this->r_funcs[i](j, dt); + p->v_vec = original_particles[j].v_vec + this->v_funcs[i](j, dt); + p->r_vec = original_particles[j].r_vec + this->r_funcs[i](j, dt); } this->particles = tmp_particles; } - - delete[] k_v; - delete[] k_r; } void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) @@ -196,7 +180,7 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) Particle *p; - std::function force; + std::function force; if (particle_interaction) { force = [this](int i) { return this->total_force(i); }; } @@ -205,19 +189,20 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) } #pragma omp parallel for private(p) - for (int i = 0; i < this->particles.size(); i++) { - p = &new_state.at(i); - p->v_vec += dt * force(i) / new_state.at(i).m; - p->r_vec += dt * this->particles.at(i).v_vec; + for (size_t i = 0; i < this->particles.size(); i++) { + p = &new_state[i]; + p->v_vec += dt * force(i) / p->m; + p->r_vec += dt * this->particles[i].v_vec; } this->particles = new_state; } -sim_arr PenningTrap::simulate(double time, int steps, std::string method, - bool particle_interaction) +sim_arr PenningTrap::simulate(double time, unsigned int steps, + std::string method, bool particle_interaction) { double dt = time / (double)steps; + sim_arr res(this->particles.size(), sim_cols(steps)); std::function func; @@ -236,9 +221,8 @@ sim_arr PenningTrap::simulate(double time, int steps, std::string method, abort(); } - for (int j = 0; j < steps; j++) { - #pragma omp parallel for - for (int i = 0; i < this->particles.size(); i++) { + for (size_t j = 0; j < steps; j++) { + for (size_t i = 0; i < this->particles.size(); i++) { res[i][j] = this->particles[i].r_vec; } func(dt, particle_interaction); @@ -261,17 +245,13 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time, sim_arr res = this->simulate(time, steps, method, particle_interaction); - int number_of_particles = this->particles.size(); - std::ofstream ofile; - sim_rows row; -#pragma omp parallel for private(ofile, row) - for (int i = 0; i < number_of_particles; i++) { - row = res[i]; +#pragma omp parallel for private(ofile) + for (size_t i = 0; i < this->particles.size(); i++) { ofile.open(path + "particle_" + std::to_string(i) + ".txt"); - for (arma::vec vector : row) { - ofile << vector(0) << "," << vector(1) << "," << vector(2) << "\n"; + for (vec_3d &vec : res[i]) { + ofile << vec(0) << "," << vec(1) << "," << vec(2) << "\n"; } ofile.close(); } diff --git a/src/main.cpp b/src/main.cpp index 2930784..cd72d0c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -24,8 +24,8 @@ #define CHARGE 1. #define MASS 40. // unit: amu -Particle p1(CHARGE, MASS, arma::vec{20., 0., 20.}, arma::vec{0., 25., 0.}); -Particle p2(CHARGE, MASS, arma::vec{25., 25., 0.}, arma::vec{0., 40., 5.}); +Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.}); +Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.}); void simulate_single_particle() { diff --git a/src/test_suite.cpp b/src/test_suite.cpp index 81ea5c1..cd422c5 100644 --- a/src/test_suite.cpp +++ b/src/test_suite.cpp @@ -44,7 +44,7 @@ public: arma::vec result; arma::vec v; std::stringstream msg; - for (int i = 0; i < tests.size(); i++) { + for (size_t i = 0; i < tests.size(); i++) { v = tests.at(i).first; result = trap.external_E_field(v); diff --git a/src/utils.cpp b/src/utils.cpp index 934e50f..14670f0 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -65,7 +65,7 @@ bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol) return false; } - for (int i = 0; i < a.n_elem; i++) { + for (size_t i = 0; i < a.n_elem; i++) { if (std::abs(a(i) - b(i)) >= tol) { return false; }