From 335907d4358e55a8f79eff5385757c2c7ade573a Mon Sep 17 00:00:00 2001 From: Cory Date: Fri, 20 Oct 2023 11:20:57 +0200 Subject: [PATCH] Formatting --- src/PenningTrap.cpp | 21 +++++---- src/main.cpp | 105 +++++++++++++++++++++++++++++++++++--------- 2 files changed, 95 insertions(+), 31 deletions(-) diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 2f97c70..d23f063 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -52,8 +52,8 @@ vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt) return dt * this->k_v[2][j]; case 3: 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])); + * (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(); @@ -71,8 +71,8 @@ vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt) return dt * this->k_r[2][j]; case 3: 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])); + * (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(); @@ -116,8 +116,8 @@ vec_3d PenningTrap::total_force_external(unsigned int i) } return vec_3d(p.q - * (this->external_E_field(p.r_vec) - + arma::cross(p.v_vec, this->external_B_field(p.r_vec)))); + * (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,11 +134,10 @@ vec_3d PenningTrap::total_force_particles(unsigned int i) 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 vec_3d(this->total_force_external(i) - - this->total_force_particles(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)); } void PenningTrap::evolve_RK4(double dt, bool particle_interaction) diff --git a/src/main.cpp b/src/main.cpp index 4b5bb5c..c5f1ab1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -97,6 +97,7 @@ void simulate_single_particle_with_different_steps() // Calculate relative error for RK4 std::string path = "output/relative_error/RK4/"; mkpath(path); +#pragma omp parallel for for (int i = 0; i < 4; i++) { int steps = 4000 * std::pow(2, i); double dt = time / (double)steps; @@ -114,6 +115,7 @@ void simulate_single_particle_with_different_steps() // Calculate relative error for forward Euler path = "output/relative_error/euler/"; mkpath(path); +#pragma omp parallel for for (int i = 0; i < 4; i++) { int steps = 4000 * std::pow(2, i); double dt = time / (double)steps; @@ -137,8 +139,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); } /** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time @@ -148,7 +151,7 @@ void simulate_100_particles() * MHz. * * */ -void simulate_100_particles_with_time_potential() +void simulate_100_particles_with_time_potential_wide_sweep() { double time = 500.; @@ -167,24 +170,85 @@ void simulate_100_particles_with_time_potential() std::ofstream ofile; - #pragma omp parallel for +#pragma omp parallel for for (size_t i = 0; i < freq_iterations; i++) { - res[0][i] = freq_start+ freq_increment*i; + res[0][i] = freq_start + freq_increment * i; } - // Using num_threads() is usually bad practice, but not having it - // causes a SIGKILL. - #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++) { - 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.).fraction_of_particles_left(time, N, "rk4", false); + 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.) + .fraction_of_particles_left(time, N, "rk4", false); + } + } + + ofile.open(path + "res.txt"); + for (size_t i = 0; i < freq_iterations; i++) { + ofile << res[0][i] << "," << res[1][i] << "," << res[2][i] << "," + << res[3][i] << "\n"; + } + ofile.close(); +} + +/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time + * dependent potential. + * + * @details The simulation sweeps over different frequencies in [1., 1.7] + * MHz. + * + * */ +void simulate_100_particles_with_time_potential_narrow_sweep() +{ + double time = 500.; + + double amplitudes[]{.1, .4, .7}; + + double freq_start = 1.; + double freq_end = 1.7; + double freq_increment = .001; + size_t freq_iterations + = (size_t)((freq_end - freq_start) / freq_increment); + + double res[4][freq_iterations]; + + std::string path = "output/time_dependent_potential/"; + mkpath(path); + + std::ofstream ofile; + +#pragma omp parallel for + for (size_t i = 0; i < freq_iterations; i++) { + res[0][i] = freq_start + freq_increment * i; + } + +// 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++) { + 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.) + .fraction_of_particles_left(time, N, "rk4", false); } } @@ -200,15 +264,16 @@ 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(); + //simulate_single_particle_with_different_steps(); simulate_100_particles(); - simulate_100_particles_with_time_potential(); + //simulate_100_particles_with_time_potential_wide_sweep(); + //simulate_100_particles_with_time_potential_narrow_sweep(); double end = omp_get_wtime();