FYS3150/Project-3
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Develop #14

Merged
coryab merged 124 commits from develop into main 2023-10-24 20:43:56 +00:00
Conversation 0 Commits 124 Files Changed 193 +66 -63
2 changed files with 66 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 7269429a37 - Show all commits

62
src/PenningTrap.cpp
View File

@ -11,6 +11,7 @@
* */ * */
#include "PenningTrap.hpp" #include "PenningTrap.hpp"
#include <functional>
PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0, PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0,
double d, double t) double d, double t)
@ -50,8 +51,9 @@ vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
case 2: case 2:
return dt * this->k_v[2][j]; return dt * this->k_v[2][j];
case 3: case 3:
return (dt / 6.) * (this->k_v[0][j] + 2. * this->k_v[1][j] + return vec_3d((dt / 6.)
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: default:
std::cout << "Not valid!" << std::endl; std::cout << "Not valid!" << std::endl;
abort(); abort();
@ -68,8 +70,9 @@ vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt)
case 2: case 2:
return dt * this->k_r[2][j]; return dt * this->k_r[2][j];
case 3: case 3:
return (dt / 6.) * (this->k_r[0][j] + 2. * this->k_r[1][j] + return vec_3d((dt / 6.)
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: default:
std::cout << "Not valid!" << std::endl; std::cout << "Not valid!" << std::endl;
abort(); abort();
@ -84,9 +87,8 @@ void PenningTrap::add_particle(Particle particle)
vec_3d PenningTrap::external_E_field(vec_3d r) vec_3d PenningTrap::external_E_field(vec_3d r)
{ {
r(2) *= -2.; 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) 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 // Get the distance between the particles
double norm = arma::norm(res, 2); 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) 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.}; return vec_3d{0., 0., 0.};
} }
vec_3d force = return vec_3d(p.q
p.q * (this->external_E_field(p.r_vec) + * (this->external_E_field(p.r_vec)
arma::cross(p.v_vec, this->external_B_field(p.r_vec))); + arma::cross(p.v_vec, this->external_B_field(p.r_vec))));
return force;
} }
vec_3d PenningTrap::total_force_particles(unsigned int i) vec_3d PenningTrap::total_force_particles(unsigned int i)
{ {
Particle p = this->particles[i];
vec_3d res; vec_3d res;
for (size_t j = 0; j < this->particles.size(); j++) { for (size_t j = 0; j < this->particles.size(); j++) {
if (i == j) { if (i != j)
continue;
}
res += this->force_on_particle(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) 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) { if (arma::norm(this->particles[i].r_vec) > this->d) {
return vec_3d{0., 0., 0.}; 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) 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(); size_t size = this->particles.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_v = sim_arr(4, sim_cols(size));
this->k_r = sim_arr(4, sim_cols(size)); this->k_r = sim_arr(4, sim_cols(size));
}
for (size_t i = 0; i < 4; i++) { for (size_t i = 0; i < 4; i++) {
#pragma omp parallel for #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_v[i][j] = (this->*force)(j) / this->particles[j].m;
this->k_r[i][j] = this->particles[j].v_vec; this->k_r[i][j] = this->particles[j].v_vec;
Particle *p = &tmp_particles[j]; tmp_particles[j].v_vec
= original_particles[j].v_vec + this->v_func(i, j, dt);
p->v_vec = original_particles[j].v_vec + this->v_func(i, j, dt); tmp_particles[j].r_vec
p->r_vec = original_particles[j].r_vec + this->r_func(i, j, dt); = original_particles[j].r_vec + this->r_func(i, j, dt);
} }
this->particles.swap(tmp_particles); this->particles = tmp_particles;
} }
this->t += dt; this->t += dt;
} }
@ -188,7 +188,7 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
vec_3d (PenningTrap::*force)(unsigned int); vec_3d (PenningTrap::*force)(unsigned int);
if (particle_interaction) { if (particle_interaction) {
force = &PenningTrap::total_force; force = &PenningTrap::total_force_external;
} }
else { else {
force = &PenningTrap::total_force_external; force = &PenningTrap::total_force_external;
@ -256,8 +256,8 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
return; return;
} }
simulation_t res = simulation_t res
this->simulate(time, steps, method, particle_interaction); = this->simulate(time, steps, method, particle_interaction);
std::ofstream ofile; std::ofstream ofile;
@ -282,8 +282,8 @@ double PenningTrap::fraction_of_particles_left(double time, unsigned int steps,
std::string method, std::string method,
bool particle_interaction) bool particle_interaction)
{ {
simulation_t res = simulation_t res
this->simulate(time, steps, method, particle_interaction); = this->simulate(time, steps, method, particle_interaction);
int particles_left = 0; int particles_left = 0;

57
src/main.cpp
View File

@ -26,8 +26,10 @@
#define MASS 40. // unit: amu #define MASS 40. // unit: amu
// Particles used for testing // Particles used for testing
Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.}); ///< Particle 1 Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.},
Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.}); ///< Particle 2 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 /** @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 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_p = (25. + w_n * 20.) / (w_n - w_p);
double A_n = -(25. + w_p * 20.) / (w_n - w_p); double A_n = -(25. + w_p * 20.) / (w_n - w_p);
std::complex<double> f = std::complex<double> f
A_p * std::exp(std::complex<double>(0., -w_p * t)) + = A_p * std::exp(std::complex<double>(0., -w_p * t))
A_n * std::exp(std::complex<double>(0., -w_n * t)); + A_n * std::exp(std::complex<double>(0., -w_n * t));
vec_3d res{std::real(f), std::imag(f), 20. * std::cos(std::sqrt(w_z2) * t)}; vec_3d res{std::real(f), std::imag(f),
20. * std::cos(std::sqrt(w_z2) * t)};
return res; return res;
} }
@ -101,8 +104,8 @@ void simulate_single_particle_with_different_steps()
PenningTrap trap(std::vector<Particle>{p1}); PenningTrap trap(std::vector<Particle>{p1});
simulation_t res = trap.simulate(time, steps, "rk4", false); simulation_t res = trap.simulate(time, steps, "rk4", false);
for (int i = 0; i < steps; i++) { for (int i = 0; i < steps; i++) {
ofile << arma::norm(res.r_vecs[0][i] - ofile << arma::norm(res.r_vecs[0][i]
analytical_solution_particle_1(dt*i)) - analytical_solution_particle_1(dt * i))
<< "\n"; << "\n";
} }
ofile.close(); ofile.close();
@ -118,8 +121,8 @@ void simulate_single_particle_with_different_steps()
PenningTrap trap(std::vector<Particle>{p1}); PenningTrap trap(std::vector<Particle>{p1});
simulation_t res = trap.simulate(time, steps, "euler", false); simulation_t res = trap.simulate(time, steps, "euler", false);
for (int i = 0; i < steps; i++) { for (int i = 0; i < steps; i++) {
ofile << arma::norm(res.r_vecs[0][i] - ofile << arma::norm(res.r_vecs[0][i]
analytical_solution_particle_1(dt*i)) - analytical_solution_particle_1(dt * i))
<< "\n"; << "\n";
} }
ofile.close(); ofile.close();
@ -134,13 +137,15 @@ void simulate_100_particles()
double time = 50.; // microseconds 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. * 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() void simulate_100_particles_with_time_potential()
@ -152,7 +157,8 @@ void simulate_100_particles_with_time_potential()
double freq_start = .2; double freq_start = .2;
double freq_end = 2.5; double freq_end = 2.5;
double freq_increment = .02; 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]; double res[4][freq_iterations];
@ -161,25 +167,24 @@ void simulate_100_particles_with_time_potential()
std::ofstream ofile; std::ofstream ofile;
double freq = freq_start; #pragma omp parallel for
for (size_t i = 0; i < freq_iterations; i++) { for (size_t i = 0; i < freq_iterations; i++) {
res[0][i] = freq; res[0][i] = freq_start+ freq_increment*i;
freq += freq_increment;
} }
#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 i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) { for (size_t j = 0; j < freq_iterations; j++) {
PenningTrap trap( res[i+1][j] = PenningTrap(
(unsigned)100, T, (unsigned)100, T,
std::bind( std::bind(
[](double f, double r, double t) { [](double f, double r, double t) {
return (25. * V / 1000.) * (1. + f * std::cos(r * t)); return (25. * V / 1000.) * (1. + f * std::cos(r * t));
}, },
amplitudes[i], res[0][j], std::placeholders::_1), amplitudes[i], res[0][j], std::placeholders::_1),
500., 0.); 500., 0.).fraction_of_particles_left(time, N, "rk4", false);
res[i + 1][j] =
trap.fraction_of_particles_left(time, N, "rk4", false);
} }
} }
@ -195,21 +200,19 @@ int main()
{ {
double t0 = omp_get_wtime(); 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();
double t1 = omp_get_wtime();
simulate_100_particles(); simulate_100_particles();
//simulate_100_particles_with_time_potential(); simulate_100_particles_with_time_potential();
double end = omp_get_wtime(); double end = omp_get_wtime();
std::cout << "Time: " << (end - t1) << " seconds" << std::endl; std::cout << "Time: " << (end - t0) << " seconds" << std::endl;
return 0; return 0;
} }