diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index 0581e50..587299a 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -30,6 +30,9 @@ * */ class PenningTrap { + /** @brief Make PenningTrapTest a friend of PenningTrap. + * */ + friend class PenningTrapTest; private: double B_0; ///< Magnetic field strength double V_0; ///< Applied potential @@ -37,10 +40,16 @@ private: double d; ///< Characteristic dimension double t; ///< Current time std::vector particles; ///< The particles in the Penning trap - sim_arr k_v; ///< A 2D vector containing all \f$k_{i,j}\f$ where \f$j\f$ is - ///< the index of a particle - sim_arr k_r; ///< A 2D vector containing all \f$k_{i,j}\f$ where \f$j\f$ is - ///< the index of a particle + + /** @brief A 2D vector containing all \f$k_{v,i,j}\f$ where \f$j\f$ is the + * index of a particle + * */ + sim_arr k_v; + + /** @brief A 2D vector containing all \f$k_{r,i,j}\f$ where \f$j\f$ is the + * index of a particle + * */ + sim_arr k_r; /** @brief Helper for evolve_RK4 when calculating \f$k_{v,i,j}\f$ values * @@ -122,6 +131,14 @@ private: * */ vec3 total_force(uint i); + /** @brief calculate the total force on a particle p_i without interaction + * + * @param i The index of particle p_i + * + * @return vec3 + * */ + vec3 total_force_no_interaction(uint i); + public: /** @brief Constructor for the PenningTrap class * @@ -229,7 +246,6 @@ public: std::string method = "rk4", bool particle_interaction = true); - friend class PenningTrapTest; }; #endif diff --git a/src/Makefile b/src/Makefile index 3718e50..cb6d014 100644 --- a/src/Makefile +++ b/src/Makefile @@ -8,7 +8,7 @@ CLASSOBJS=$(CLASSSRCS:.cpp=.o) INCLUDE=../include -CFLAGS=-Wall -larmadillo -lblas -llapack -std=c++11 -O3 +CFLAGS=-Wall -larmadillo -lblas -llapack -std=c++11 -O3 -fomit-frame-pointer OPENMP=-fopenmp # Add a debug flag when compiling (For the DEBUG macro in utils.hpp) @@ -19,6 +19,7 @@ else DBGFLAG= endif +# Add profiling for serial code PROFILE ?= 0 ifeq ($(PROFILE), 1) PROFFLAG=-pg -fno-inline-functions @@ -30,6 +31,7 @@ endif all: test_suite main +# Instrumentation using scorep for parallel analysis instrument: scorep $(CC) -c PenningTrap.cpp -o PenningTrap.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) scorep $(CC) -c Particle.cpp -o Particle.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP) diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 1239be8..95a03c4 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -10,13 +10,14 @@ * @bug No known bugs * */ -#include "PenningTrap.hpp" -#include "typedefs.hpp" #include #include #include #include +#include "PenningTrap.hpp" +#include "typedefs.hpp" + vec3 PenningTrap::v_func(uint i, uint j, double dt) { switch (i) { @@ -59,8 +60,8 @@ vec3 PenningTrap::external_E_field(vec3 r) { r(2) *= -2.; - return vec3((this->V_0 * this->perturbation(this->t) / (this->d * this->d)) - * r); + return vec3( + ((this->V_0 * this->perturbation(this->t)) / (this->d * this->d)) * r); } vec3 PenningTrap::external_B_field(vec3 r) @@ -83,10 +84,6 @@ vec3 PenningTrap::total_force_external(uint i) { Particle *p = &this->particles[i]; - if (arma::norm(p->r_vec) > this->d) { - return vec3{0., 0., 0.}; - } - return vec3(p->q * (this->external_E_field(p->r_vec) + arma::cross(p->v_vec, this->external_B_field(p->r_vec)))); @@ -112,6 +109,14 @@ vec3 PenningTrap::total_force(uint i) return vec3(this->total_force_external(i) - this->total_force_particles(i)); } +vec3 PenningTrap::total_force_no_interaction(uint i) +{ + if (arma::norm(this->particles[i].r_vec) > this->d) { + return vec3{0., 0., 0.}; + } + return vec3(this->total_force_external(i)); +} + PenningTrap::PenningTrap(double B_0, double V_0, double d, double t) { this->B_0 = B_0; @@ -149,9 +154,12 @@ void PenningTrap::reinitialize(double f, double omega_V, double t) { this->t = t; this->set_pertubation(f, omega_V); + Particle *p; for (size_t i = 0; i < this->particles.size(); i++) { - this->particles[i].r_vec = vec3().randn() * .1 * this->d; + p = &this->particles[i]; + p->v_vec = vec3().randn() * .1 * this->d; + p->v_vec = vec3().randn() * .1 * this->d; } } @@ -162,13 +170,14 @@ void PenningTrap::add_particle(Particle particle) void PenningTrap::evolve_RK4(double dt, bool particle_interaction) { + Particle *p; std::vector original_particles = this->particles; std::vector tmp_particles = this->particles; - vec3 (PenningTrap::*force)(uint) = particle_interaction - ? &PenningTrap::total_force - : &PenningTrap::total_force_external; + vec3 (PenningTrap::*force)(uint) = + particle_interaction ? &PenningTrap::total_force + : &PenningTrap::total_force_no_interaction; size_t size = this->particles.size(); @@ -182,15 +191,15 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) // 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 +#pragma omp parallel for private(p) 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; - 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); + 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); } this->particles = tmp_particles; } @@ -202,10 +211,11 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) { size_t size = this->particles.size(); vec3 force_res[size]; + Particle *p; - vec3 (PenningTrap::*force)(uint) = particle_interaction - ? &PenningTrap::total_force - : &PenningTrap::total_force_external; + vec3 (PenningTrap::*force)(uint) = + particle_interaction ? &PenningTrap::total_force + : &PenningTrap::total_force_no_interaction; // Calculating the force for each particle is independent and therefore // a good candidate for parallel execution @@ -218,8 +228,9 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) // this as well #pragma omp parallel for for (size_t i = 0; i < size; i++) { - this->particles[i].r_vec += dt * this->particles[i].v_vec; - this->particles[i].v_vec += dt * force_res[i] / this->particles[i].m; + p = &this->particles[i]; + p->r_vec += dt * p->v_vec; + p->v_vec += dt * force_res[i] / p->m; } this->t += dt; @@ -228,6 +239,7 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) simulation_t PenningTrap::simulate(double time, uint steps, std::string method, bool particle_interaction) { + Particle *p; double dt = time / (double)steps; uint size = this->particles.size(); @@ -247,8 +259,9 @@ simulation_t PenningTrap::simulate(double time, uint steps, std::string method, for (size_t j = 0; j < steps; j++) { for (size_t i = 0; i < size; i++) { - res.r_vecs[i][j] = this->particles[i].r_vec; - res.v_vecs[i][j] = this->particles[i].v_vec; + p = &this->particles[i]; + res.r_vecs[i][j] = p->r_vec; + res.v_vecs[i][j] = p->v_vec; } (this->*func)(dt, particle_interaction); } diff --git a/src/main.cpp b/src/main.cpp index f6a707c..937c306 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -135,8 +136,9 @@ void simulate_100_particles() double time = 50.; // microseconds - trap.write_simulation_to_dir("output/simulate_100_particles", time, N, - "rk4", false); + // 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 @@ -212,7 +214,8 @@ void potential_resonance_narrow_sweep() 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); + size_t freq_iterations = + (size_t)((freq_end - freq_start) / freq_increment) + 1; double res[4][freq_iterations]; @@ -268,7 +271,8 @@ void potential_resonance_narrow_sweep_interaction() 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); + size_t freq_iterations = + (size_t)((freq_end - freq_start) / freq_increment) + 1; double res[4][freq_iterations]; @@ -308,60 +312,104 @@ void potential_resonance_narrow_sweep_interaction() int main() { - double start, end, t1, t2; + int option = 1; + bool chosen = false; + + system("clear"); + std::cout << "(1) All (default)\n" + << "(2) Simulate single particle\n" + << "(3) simulate 2 particles\n" + << "(4) Simulate single particle with different time steps\n" + << "(5) Simulate 100 particles\n" + << "(6) Potential resonance wide sweep\n" + << "(7) Potential resonance narrow sweep without particle " + "interactions\n" + << "(8) Potential resonance narrow sweep with particle " + "interaction\n" + << "Select what to run: "; + std::cin >> std::noskipws; + do { + if (!(std::cin >> option) || option < 1 || option > 8) { + std::cin.clear(); + std::cin.ignore(std::numeric_limits::max(), '\n'); + system("clear"); + std::cout + << "(1) All (default)\n" + << "(2) Simulate single particle\n" + << "(3) simulate 2 particles\n" + << "(4) Simulate single particle with different time steps\n" + << "(5) Simulate 100 particles\n" + << "(6) Potential resonance wide sweep\n" + << "(7) Potential resonance narrow sweep without particle " + "interactions\n" + << "(8) Potential resonance narrow sweep with particle " + "interaction\n" + << "Not a valid option, please enter a valid number: "; + } else { + chosen = true; + } + } while (!chosen); + + double start, end; + + system("clear"); + start = omp_get_wtime(); + switch (option) { + case 1: + std::cout << "Running simulate_single_particle\n"; + simulate_single_particle(); - simulate_single_particle(); + std::cout << "Running simulate_two_particles\n"; + simulate_two_particles(); - simulate_two_particles(); + std::cout << "Running simulate_single_particle_with_different_steps\n"; + simulate_single_particle_with_different_steps(); - simulate_single_particle_with_different_steps(); + std::cout << "Running simulate_100_particles\n"; + simulate_100_particles(); - t2 = omp_get_wtime(); + std::cout << "Running potential_resonance_wide_sweep\n"; + potential_resonance_wide_sweep(); - std::cout << "Time single and double : " << (t2 - start) - << " seconds" << std::endl; - - t1 = omp_get_wtime(); - - simulate_100_particles(); - - t2 = omp_get_wtime(); - - std::cout << "Time 100 particles : " << (t2 - t1) - << " seconds" << std::endl; - - t1 = omp_get_wtime(); - - potential_resonance_wide_sweep(); - - t2 = omp_get_wtime(); - - std::cout << "Time wide sweep : " << (t2 - t1) - << " seconds" << std::endl; - - t1 = omp_get_wtime(); - - potential_resonance_narrow_sweep(); - - t2 = omp_get_wtime(); - - std::cout << "Time narrow sweep no interaction : " << (t2 - t1) - << " seconds" << std::endl; - - t1 = omp_get_wtime(); - - potential_resonance_narrow_sweep_interaction(); - - t2 = omp_get_wtime(); - - std::cout << "Time narrow sweep with interaction: " << (t2 - t1) - << " seconds" << std::endl; + std::cout << "Running potential_resonance_narrow_sweep\n"; + potential_resonance_narrow_sweep(); + std::cout << "Running potential_resonance_narrow_sweep_interaction\n"; + potential_resonance_narrow_sweep_interaction(); + break; + case 2: + std::cout << "Running simulate_single_particle\n"; + simulate_single_particle(); + break; + case 3: + std::cout << "Running simulate_two_particles\n"; + simulate_two_particles(); + break; + case 4: + std::cout << "Running simulate_single_particle_with_different_steps\n"; + simulate_single_particle_with_different_steps(); + break; + case 5: + std::cout << "Running simulate_100_particles\n"; + simulate_100_particles(); + break; + case 6: + std::cout << "Running potential_resonance_wide_sweep\n"; + potential_resonance_wide_sweep(); + break; + case 7: + std::cout << "Running potential_resonance_narrow_sweep\n"; + potential_resonance_narrow_sweep(); + break; + case 8: + std::cout << "Running potential_resonance_narrow_sweep_interaction\n"; + potential_resonance_narrow_sweep_interaction(); + break; + } end = omp_get_wtime(); - std::cout << "Time : " << (end - start) - << " seconds" << std::endl; + std::cout << "Time: " << end - start << " seconds" << std::endl; return 0; }