diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index a66ca51..2f97c70 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -11,6 +11,7 @@ * */ #include "PenningTrap.hpp" +#include PenningTrap::PenningTrap(double B_0, std::function V_0, double d, double t) @@ -50,8 +51,9 @@ vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt) case 2: return dt * this->k_v[2][j]; case 3: - return (dt / 6.) * (this->k_v[0][j] + 2. * this->k_v[1][j] + - 2. * this->k_v[2][j] + this->k_v[3][j]); + return vec_3d((dt / 6.) + * (this->k_v[0][j] + 2. * this->k_v[1][j] + 2. * this->k_v[2][j] + + this->k_v[3][j])); default: std::cout << "Not valid!" << std::endl; abort(); @@ -68,8 +70,9 @@ vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt) case 2: return dt * this->k_r[2][j]; case 3: - return (dt / 6.) * (this->k_r[0][j] + 2. * this->k_r[1][j] + - 2. * this->k_r[2][j] + this->k_r[3][j]); + return vec_3d((dt / 6.) + * (this->k_r[0][j] + 2. * this->k_r[1][j] + 2. * this->k_r[2][j] + + this->k_r[3][j])); default: std::cout << "Not valid!" << std::endl; abort(); @@ -84,9 +87,8 @@ void PenningTrap::add_particle(Particle particle) vec_3d PenningTrap::external_E_field(vec_3d r) { r(2) *= -2.; - double f = this->V_0(this->t) / (this->d * this->d); - return f * r; + return vec_3d(this->V_0(this->t) / (this->d * this->d) * r); } vec_3d PenningTrap::external_B_field(vec_3d r) @@ -102,7 +104,7 @@ vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j) // Get the distance between the particles double norm = arma::norm(res, 2); - return vec_3d(res * this->particles[j].q / (norm * norm * norm)); + return vec_3d((this->particles[j].q / (norm * norm * norm)) * res); } vec_3d PenningTrap::total_force_external(unsigned int i) @@ -113,28 +115,21 @@ vec_3d PenningTrap::total_force_external(unsigned int i) return vec_3d{0., 0., 0.}; } - vec_3d force = - p.q * (this->external_E_field(p.r_vec) + - arma::cross(p.v_vec, this->external_B_field(p.r_vec))); - - return force; + 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) { - Particle p = this->particles[i]; - vec_3d res; for (size_t j = 0; j < this->particles.size(); j++) { - if (i == j) { - continue; - } - - res += this->force_on_particle(i, j); + if (i != j) + res += this->force_on_particle(i, j); } - return vec_3d(res * K_E * (p.q / p.m)); + return vec_3d(res * K_E * (this->particles[i].q / this->particles[i].m)); } vec_3d PenningTrap::total_force(unsigned int i) @@ -142,7 +137,8 @@ vec_3d PenningTrap::total_force(unsigned int i) if (arma::norm(this->particles[i].r_vec) > this->d) { return vec_3d{0., 0., 0.}; } - return this->total_force_external(i) - this->total_force_particles(i); + return vec_3d(this->total_force_external(i) + - this->total_force_particles(i)); } void PenningTrap::evolve_RK4(double dt, bool particle_interaction) @@ -161,21 +157,25 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) size_t size = this->particles.size(); - this->k_v = sim_arr(4, sim_cols(size)); - this->k_r = sim_arr(4, sim_cols(size)); + // Allocating takes a long time, so reuse sim_arr if possible + if (this->k_v.size() != 4 || this->k_r.size() != 4 + || this->k_v[0].size() != size || this->k_r[0].size() != size) { + this->k_v = sim_arr(4, sim_cols(size)); + this->k_r = sim_arr(4, sim_cols(size)); + } for (size_t i = 0; i < 4; i++) { #pragma omp parallel for - for (size_t j = 0; j < this->particles.size(); j++) { + for (size_t j = 0; j < size; j++) { this->k_v[i][j] = (this->*force)(j) / this->particles[j].m; this->k_r[i][j] = this->particles[j].v_vec; - Particle *p = &tmp_particles[j]; - - p->v_vec = original_particles[j].v_vec + this->v_func(i, j, dt); - p->r_vec = original_particles[j].r_vec + this->r_func(i, j, dt); + tmp_particles[j].v_vec + = original_particles[j].v_vec + this->v_func(i, j, dt); + tmp_particles[j].r_vec + = original_particles[j].r_vec + this->r_func(i, j, dt); } - this->particles.swap(tmp_particles); + this->particles = tmp_particles; } this->t += dt; } @@ -188,7 +188,7 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) vec_3d (PenningTrap::*force)(unsigned int); if (particle_interaction) { - force = &PenningTrap::total_force; + force = &PenningTrap::total_force_external; } else { force = &PenningTrap::total_force_external; @@ -256,8 +256,8 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time, return; } - simulation_t res = - this->simulate(time, steps, method, particle_interaction); + simulation_t res + = this->simulate(time, steps, method, particle_interaction); std::ofstream ofile; @@ -282,8 +282,8 @@ 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); + simulation_t res + = this->simulate(time, steps, method, particle_interaction); int particles_left = 0; diff --git a/src/main.cpp b/src/main.cpp index 604af46..4b5bb5c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -26,8 +26,10 @@ #define MASS 40. // unit: amu // Particles used for testing -Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.}); ///< Particle 1 -Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.}); ///< Particle 2 +Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, + vec_3d{0., 25., 0.}); ///< Particle 1 +Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, + vec_3d{0., 40., 5.}); ///< Particle 2 /** @brief The analytical solution for particle p1 * @@ -43,10 +45,11 @@ vec_3d 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::complex f = - A_p * std::exp(std::complex(0., -w_p * t)) + - A_n * std::exp(std::complex(0., -w_n * t)); - vec_3d res{std::real(f), std::imag(f), 20. * std::cos(std::sqrt(w_z2) * t)}; + std::complex f + = A_p * std::exp(std::complex(0., -w_p * t)) + + A_n * std::exp(std::complex(0., -w_n * t)); + vec_3d res{std::real(f), std::imag(f), + 20. * std::cos(std::sqrt(w_z2) * t)}; return res; } @@ -101,8 +104,8 @@ void simulate_single_particle_with_different_steps() 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)) + ofile << arma::norm(res.r_vecs[0][i] + - analytical_solution_particle_1(dt * i)) << "\n"; } ofile.close(); @@ -118,8 +121,8 @@ void simulate_single_particle_with_different_steps() 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)) + ofile << arma::norm(res.r_vecs[0][i] + - analytical_solution_particle_1(dt * i)) << "\n"; } ofile.close(); @@ -134,13 +137,15 @@ void simulate_100_particles() double time = 50.; // microseconds - trap.write_simulation_to_dir("output/simulate_100_particles", time, N); + trap.write_simulation_to_dir("output/simulate_100_particles", time, N, + "rk4", false); } -/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time +/** @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. + * @details The simulation sweeps over different frequencies in [0.2, 2.5] + * MHz. * * */ void simulate_100_particles_with_time_potential() @@ -152,7 +157,8 @@ void simulate_100_particles_with_time_potential() 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); + size_t freq_iterations + = (size_t)((freq_end - freq_start) / freq_increment); double res[4][freq_iterations]; @@ -161,25 +167,24 @@ void simulate_100_particles_with_time_potential() std::ofstream ofile; - double freq = freq_start; + #pragma omp parallel for for (size_t i = 0; i < freq_iterations; i++) { - res[0][i] = freq; - freq += freq_increment; + res[0][i] = freq_start+ freq_increment*i; } -#pragma omp parallel for collapse(2) num_threads(4) + // Using num_threads() is usually bad practice, but not having it + // causes a SIGKILL. + #pragma omp parallel for collapse(2) num_threads(4) for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < freq_iterations; j++) { - PenningTrap trap( + res[i+1][j] = PenningTrap( (unsigned)100, T, std::bind( [](double f, double r, double t) { return (25. * V / 1000.) * (1. + f * std::cos(r * t)); }, amplitudes[i], res[0][j], std::placeholders::_1), - 500., 0.); - res[i + 1][j] = - trap.fraction_of_particles_left(time, N, "rk4", false); + 500., 0.).fraction_of_particles_left(time, N, "rk4", false); } } @@ -195,21 +200,19 @@ int main() { double t0 = omp_get_wtime(); - // simulate_single_particle(); + simulate_single_particle(); - // simulate_two_particles(); + simulate_two_particles(); - simulate_single_particle_with_different_steps(); - - double t1 = omp_get_wtime(); + simulate_single_particle_with_different_steps(); simulate_100_particles(); - //simulate_100_particles_with_time_potential(); + simulate_100_particles_with_time_potential(); double end = omp_get_wtime(); - std::cout << "Time: " << (end - t1) << " seconds" << std::endl; + std::cout << "Time: " << (end - t0) << " seconds" << std::endl; return 0; }