diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index 4a131aa..74759e2 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -36,8 +36,8 @@ private: std::vector particles; ///< The particles in the Penning trap sim_arr k_v; sim_arr k_r; - std::function v_funcs[4]; - std::function r_funcs[4]; + std::function* v_funcs; + std::function* r_funcs; public: /** @brief Set B_0, V_0 and d. diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 93c6ed7..4ab38e9 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -29,21 +29,31 @@ PenningTrap::PenningTrap(double B_0, double V_0, double d) this->d = d; // 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) { - 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()); + this->v_funcs = new std::function[] { + [this](int i, double dt) { + DEBUG("Inside v_funcs"); + return .5 * dt * this->k_v[0][i]; + }, + [this](int i, double dt) { return .5 * dt * this->k_v[1][i]; }, + [this](int i, double dt) { return dt * this->k_v[2][i]; }, + [this](int i, double dt) { + 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()); + } + }; + this->r_funcs = new std::function[] { + [this](int i, double dt) { + DEBUG("Inside v_funcs"); + return .5 * dt * this->k_r[0][i]; + }, + [this](int i, double dt) { return .5 * dt * this->k_r[1][i]; }, + [this](int i, double dt) { return dt * this->k_r[2][i]; }, + [this](int i, double dt) { + 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()); + } }; - 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) { - 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) @@ -79,7 +89,7 @@ vec_3d PenningTrap::external_E_field(vec_3d r) vec_3d PenningTrap::external_B_field(vec_3d r) { - return vec_3d{0.,0.,this->B_0}; + return vec_3d{0., 0., this->B_0}; } vec_3d PenningTrap::force_on_particle(int i, int j) @@ -143,9 +153,11 @@ vec_3d PenningTrap::total_force(int i) void PenningTrap::evolve_RK4(double dt, bool particle_interaction) { + DEBUG("Inside RK4"); std::vector original_particles = this->particles; std::vector tmp_particles = this->particles; + DEBUG("Choose force func"); std::function force; if (particle_interaction) { force = [this](int i) { return this->total_force(i); }; @@ -154,22 +166,28 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) force = [this](int i) { return this->total_force_external(i); }; } + DEBUG("Create k_v and k_r"); int size = this->particles.size(); this->k_v = sim_arr(4, sim_cols(size)); this->k_r = sim_arr(4, sim_cols(size)); + DEBUG("Loop"); for (int i = 0; i < 4; i++) { + DEBUG("Inside outer loop"); #pragma omp parallel for for (int j = 0; j < size; j++) { + DEBUG("Set k_v and k_r"); 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[j]; + DEBUG("Update v and r"); 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); } + DEBUG("After inner loop"); this->particles = tmp_particles; } } @@ -201,13 +219,16 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) sim_arr PenningTrap::simulate(double time, unsigned int steps, std::string method, bool particle_interaction) { + DEBUG("Inside simulate"); double dt = time / (double)steps; sim_arr res(this->particles.size(), sim_cols(steps)); + DEBUG("Choose function to use"); std::function func; if (method == "rk4") { func = [this](double dt, bool particle_interaction) { + DEBUG("Inside RK4 lambda"); this->evolve_RK4(dt, particle_interaction); }; } @@ -221,10 +242,13 @@ sim_arr PenningTrap::simulate(double time, unsigned int steps, abort(); } + DEBUG("Loop"); for (size_t j = 0; j < steps; j++) { + DEBUG("Inside outer loop"); for (size_t i = 0; i < this->particles.size(); i++) { res[i][j] = this->particles[i].r_vec; } + DEBUG("After inner loop"); func(dt, particle_interaction); } @@ -235,6 +259,7 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time, int steps, std::string method, bool particle_interaction) { + DEBUG("Inside write to dir"); if (path.back() != '/') { path += '/'; } @@ -243,10 +268,12 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time, return; } + DEBUG("Create sim_arr"); sim_arr res = this->simulate(time, steps, method, particle_interaction); std::ofstream ofile; + DEBUG("Loop"); #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"); diff --git a/src/main.cpp b/src/main.cpp index cd72d0c..35ec896 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -29,10 +29,12 @@ Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.}); void simulate_single_particle() { + DEBUG("Inside single particle sim"); PenningTrap trap(std::vector{p1}); double time = 50.; // microseconds + DEBUG("Write to dir"); trap.write_simulation_to_dir("output/simulate_single_particle", time, N); } @@ -82,14 +84,19 @@ void simulate_100_particles() int main() { + DEBUG("Start"); double start = omp_get_wtime(); + DEBUG("Single particle"); simulate_single_particle(); + DEBUG("Two particles"); simulate_two_particles(); + DEBUG("Single particle with different steps"); simulate_single_particle_with_different_steps(); + DEBUG("100 particles"); simulate_100_particles(); double end = omp_get_wtime();