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

Merge branch 'develop' into janitaws/latex-setup

Browse Source
This commit is contained in:
Janita Willumsen
2023-10-22 17:51:18 +02:00
parent 937e681195 bbb0203d35
commit 346a800428
88 changed files with 3404 additions and 1274 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/Makefile
View File

@@ -21,7 +21,7 @@ endif
PROFILE ?= 0
ifeq ($(PROFILE), 1)
PROFFLAG=-pg
PROFFLAG=-pg -fno-inline-functions
else
PROFFLAG=
endif
@@ -30,6 +30,16 @@ endif
all: test_suite main
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)
scorep $(CC) -c utils.cpp -o utils.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
scorep $(CC) -c main.cpp -o main.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
scorep $(CC) $(LIBOBJS) $(CLASSOBJS) main.o -o main $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
# Rules for executables
main: main.o $(LIBOBJS) $(CLASSOBJS)
$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)

4
src/Particle.cpp
View File

@@ -13,8 +13,8 @@
#include "Particle.hpp"
Particle::Particle(double q, double m,
arma::vec::fixed<3> r_vec,
arma::vec::fixed<3> v_vec)
vec_3d r_vec,
vec_3d v_vec)
{
// Giving the particle its properties
this->q = q;

173
src/PenningTrap.cpp
View File

@@ -11,9 +11,10 @@
* */
#include "PenningTrap.hpp"
#include "constants.hpp"
#include "typedefs.hpp"
#include "utils.hpp"
#include <algorithm>
#include <functional>
#include <vector>
PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0,
double d, double t)
@@ -28,11 +29,9 @@ PenningTrap::PenningTrap(unsigned int i, double B_0,
std::function<double(double)> V_0, double d, double t)
: PenningTrap::PenningTrap(B_0, V_0, d)
{
vec_3d r, v;
for (size_t j = 0; j < i; j++) {
r = vec_3d().randn() * .1 * this->d;
v = vec_3d().randn() * .1 * this->d;
this->add_particle(Particle(1., 40., r, v));
this->particles.push_back(Particle(1., 40., vec_3d(vec_3d().randn() * .1 * this->d),
vec_3d(vec_3d().randn() * .1 * this->d)));
}
}
@@ -43,6 +42,16 @@ PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0,
this->particles = particles;
}
void PenningTrap::reinitialize(std::function<double(double)> V_0, double t)
{
this->V_0 = V_0;
this->t = t;
for (size_t i = 0; i < this->particles.size(); i++) {
this->particles[i].r_vec = vec_3d().randn() * .1 * this->d;
}
}
vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
{
switch (i) {
@@ -53,8 +62,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();
@@ -71,8 +81,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();
@@ -87,9 +98,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)
@@ -99,51 +109,48 @@ vec_3d PenningTrap::external_B_field(vec_3d r)
vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
{
Particle p_j = this->particles[j];
// Calculate the difference between the particles' position
vec_3d res = this->particles[i].r_vec - p_j.r_vec;
vec_3d res = this->particles[i].r_vec - this->particles[j].r_vec;
// Get the distance between the particles
double norm = arma::norm(res, 2);
return vec_3d(res * p_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)
{
Particle p = this->particles[i];
Particle *p = &this->particles[i];
if (arma::norm(p.r_vec) > this->d) {
if (arma::norm(p->r_vec) > this->d) {
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));
}
vec_3d PenningTrap::total_force(unsigned int i)
{
return this->total_force_external(i) - this->total_force_particles(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));
}
void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
@@ -152,32 +159,35 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
std::vector<Particle> original_particles = this->particles;
std::vector<Particle> tmp_particles = this->particles;
vec_3d (PenningTrap::*force)(unsigned int);
if (particle_interaction) {
force = &PenningTrap::total_force;
}
else {
force = &PenningTrap::total_force_external;
}
vec_3d (PenningTrap::*force)(unsigned int)
= particle_interaction ? &PenningTrap::total_force
: &PenningTrap::total_force_external;
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));
}
// 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
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;
}
@@ -187,19 +197,19 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
vec_3d force_res[size];
Particle *p;
vec_3d (PenningTrap::*force)(unsigned int);
if (particle_interaction) {
force = &PenningTrap::total_force;
}
else {
force = &PenningTrap::total_force_external;
}
vec_3d (PenningTrap::*force)(unsigned int)
= particle_interaction ? &PenningTrap::total_force
: &PenningTrap::total_force_external;
// Calculating the force for each particle is independent and therefore
// a good candidate for parallel execution
#pragma omp parallel for
for (size_t i = 0; i < size; i++) {
force_res[i] = (this->*force)(i);
}
// Updating the particles is also independent, so we can parallelize
// this as well
#pragma omp parallel for private(p)
for (size_t i = 0; i < size; i++) {
p = &this->particles[i];
@@ -217,7 +227,7 @@ simulation_t PenningTrap::simulate(double time, unsigned int steps,
double dt = time / (double)steps;
unsigned int size = this->particles.size();
// sim_arr res(this->particles.size(), sim_cols(steps));
simulation_t res{sim_arr(size, sim_cols(steps)),
sim_arr(size, sim_cols(steps))};
@@ -253,27 +263,31 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
path += '/';
}
if (mkpath(path, 0777) != 0) {
std::cout << "Hello" << std::endl;
return;
std::cout << "Failed to make path" << std::endl;
abort();
}
simulation_t res =
this->simulate(time, steps, method, particle_interaction);
simulation_t res
= this->simulate(time, steps, method, particle_interaction);
std::ofstream ofile;
// Writing each particle to its own file is independent and can be run in
// parallel.
#pragma omp parallel for private(ofile)
for (size_t i = 0; i < this->particles.size(); i++) {
ofile.open(path + "particle_" + std::to_string(i) + "_r.txt");
for (vec_3d &vec : res.r_vecs[i]) {
ofile << vec(0) << "," << vec(1) << "," << vec(2) << "\n";
ofile << scientific_format(vec(0), 10, 8) << ','
<< scientific_format(vec(1), 10, 8) << ','
<< scientific_format(vec(2), 10, 8) << '\n';
}
ofile.close();
ofile.open(path + "particle_" + std::to_string(i) + "_v.txt");
for (vec_3d &vec : res.v_vecs[i]) {
ofile << scientific_format(vec(0), 10, 8) << ","
<< scientific_format(vec(1), 8, 10) << ","
<< scientific_format(vec(2), 8, 10) << "\n";
ofile << scientific_format(vec(0), 10, 8) << ','
<< scientific_format(vec(1), 10, 8) << ','
<< scientific_format(vec(2), 10, 8) << '\n';
}
ofile.close();
}
@@ -283,26 +297,33 @@ 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);
double dt = time / (double)steps;
void (PenningTrap::*func)(double, bool);
if (method == "rk4") {
func = &PenningTrap::evolve_RK4;
}
else if (method == "euler") {
func = &PenningTrap::evolve_forward_euler;
}
else {
std::cout << "Not a valid method!" << std::endl;
abort();
}
for (size_t j = 0; j < steps; j++) {
(this->*func)(dt, particle_interaction);
}
int particles_left = 0;
for (Particle p : this->particles) {
if (arma::norm(p.r_vec) < this->d) {
// A reduction is perfect here
#pragma omp parallel for reduction(+:particles_left)
for (size_t i=0; i < this->particles.size(); i++) {
if (arma::norm(this->particles[i].r_vec) < this->d) {
particles_left++;
}
}
return (double)particles_left / (double)this->particles.size();
}
vec_3d PenningTrap::get_r(int i)
{
return this->particles[i].r_vec;
}
double PenningTrap::get_t()
{
return this->t;
}

295
src/main.cpp
View File

@@ -13,10 +13,8 @@
#include <cmath>
#include <complex>
#include <fstream>
#include <functional>
#include <omp.h>
#include <string>
#include <sys/stat.h>
#include <vector>
#include "PenningTrap.hpp"
@@ -24,12 +22,21 @@
#define PARTICLES 100
#define N 40000
#define CHARGE 1.
#define MASS 40. // unit: amu
#define CHARGE 1. // unit: e
#define MASS 40.078 // unit: amu
Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.});
Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.});
// 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
/** @brief The analytical solution for particle p1
*
* @param t Time
*
* @return vec_3d
* */
vec_3d analytical_solution_particle_1(double t)
{
double w_0 = T / MASS;
@@ -38,98 +45,122 @@ 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<double> f = A_p * std::exp(std::complex<double>(0., -w_p * 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)};
std::complex<double> f
= A_p * std::exp(std::complex<double>(0., -w_p * 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)};
return res;
}
/** @brief Simulate a single particle over the period of 50 \f$ \mu s \f$.
* */
void simulate_single_particle()
{
DEBUG("Inside single particle sim");
// Initialize trap with particle 1
PenningTrap trap(std::vector<Particle>{p1});
double time = 50.; // microseconds
DEBUG("Write to dir");
// Simulate and write results to file
trap.write_simulation_to_dir("output/simulate_single_particle", time, N,
"rk4", false);
}
/** @brief Simulate 2 particles over the period of 50 \f$ \mu s \f$ with and
* without particle interactions.
* */
void simulate_two_particles()
{
// Initialize traps with particles
PenningTrap trap_no_interaction(std::vector<Particle>{p1, p2});
PenningTrap trap_with_interaction(std::vector<Particle>{p1, p2});
double time = 50.; // microseconds
// Simulate and write results to files
trap_no_interaction.write_simulation_to_dir(
"output/simulate_2_particles/no_interaction", time, N, "rk4", false);
trap_with_interaction.write_simulation_to_dir(
"output/simulate_2_particles/with_interaction", time, N);
}
/** @brief Simulate a single particle over 50 \f$ \mu s \f$ using different
* amount of steps and different methods.
* */
void simulate_single_particle_with_different_steps()
{
double time = 50.; // microseconds
std::ofstream ofile;
// 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;
std::string path = "output/relative_error/RK4/";
mkpath(path);
ofile.open(path + std::to_string(steps) + "_steps.txt");
PenningTrap trap(std::vector<Particle>{p1});
simulation_t res = trap.simulate(time, steps, "rk4", false);
for (int i = 0; i < steps; i++) {
trap.evolve_RK4(dt);
ofile << arma::norm(trap.get_r(0) -
analytical_solution_particle_1(trap.get_t()))
<< "\n";
ofile << arma::norm(res.r_vecs[0][i]
- analytical_solution_particle_1(dt * i))
<< '\n';
}
ofile.close();
}
// 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;
std::string path = "output/relative_error/euler/";
mkpath(path);
ofile.open(path + std::to_string(steps) + "_steps.txt");
PenningTrap trap(std::vector<Particle>{p1});
simulation_t res = trap.simulate(time, steps, "euler", false);
for (int i = 0; i < steps; i++) {
trap.evolve_forward_euler(dt);
ofile << arma::norm(trap.get_r(0) -
analytical_solution_particle_1(trap.get_t()))
<< "\n";
ofile << arma::norm(res.r_vecs[0][i]
- analytical_solution_particle_1(dt * i))
<< '\n';
}
ofile.close();
}
}
/** @brief Simulate 100 particles over 50 \f$ \mu s \f$.
* */
void simulate_100_particles()
{
PenningTrap trap((unsigned)100, T,
[](double t) {
return 25. * V / 1000. * (1. + .4 * std::cos(1.5 * t));
},
500., 0);
PenningTrap trap((unsigned)100);
double time = 500.; // microseconds
double time = 50.; // microseconds
trap.write_simulation_to_dir("output/simulate_100_particles", time, N * 4);
trap.write_simulation_to_dir("output/simulate_100_particles", time, N,
"rk4", false);
}
void simulate_100_particles_with_time_potential()
/** @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.
*
* */
void potential_resonance_wide_sweep()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
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) + 1;
double res[4][freq_iterations];
@@ -138,39 +169,159 @@ void simulate_100_particles_with_time_potential()
std::ofstream ofile;
double freq = freq_start;
#pragma omp parallel for
// Insert frequencies
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)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
PenningTrap trap(
(unsigned)100, T,
std::bind(
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((unsigned int)100);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(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(500., 40000, "rk4", false);
amplitudes[i], res[0][j], std::placeholders::_1));
res[i + 1][j]
= trap.fraction_of_particles_left(time, N, "rk4", false);
}
}
}
ofile.open(path + "res.txt");
// Write results to file
ofile.open(path + "wide_sweep.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 << 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 potential_resonance_no_interaction_narrow_sweep()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
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);
double res[4][freq_iterations];
std::string path = "output/time_dependent_potential/";
mkpath(path);
std::ofstream ofile;
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((unsigned int)100);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(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));
res[i + 1][j]
= trap.fraction_of_particles_left(time, N, "rk4", false);
}
}
}
// Write results to file
ofile.open(path + "narrow_sweep.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 potential_resonance_with_interaction_narrow_sweep()
{
double time = 500.;
double amplitudes[]{.1, .4, .7};
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);
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;
}
#pragma omp parallel
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((unsigned int)100);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {
// Reset particles and give new time dependent potential.
trap.reinitialize(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));
res[i + 1][j] = trap.fraction_of_particles_left(time, N);
}
}
}
// Write results to file
ofile.open(path + "narrow_sweep_interactions.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();
}
int main()
{
double start = omp_get_wtime();
double start, end, t1, t2;
start = omp_get_wtime();
simulate_single_particle();
@@ -178,13 +329,51 @@ int main()
simulate_single_particle_with_different_steps();
t2 = omp_get_wtime();
std::cout << "Time single and double : " << (t2 - start)
<< " seconds" << std::endl;
t1 = omp_get_wtime();
simulate_100_particles();
simulate_100_particles_with_time_potential();
t2 = omp_get_wtime();
double end = omp_get_wtime();
std::cout << "Time 100 particles : " << (t2 - t1)
<< " seconds" << std::endl;
std::cout << "Time: " << (end - start) << " 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_no_interaction_narrow_sweep();
t2 = omp_get_wtime();
std::cout << "Time narrow sweep no interaction : " << (t2 - t1)
<< " seconds" << std::endl;
t1 = omp_get_wtime();
potential_resonance_with_interaction_narrow_sweep();
t2 = omp_get_wtime();
std::cout << "Time narrow sweep with interaction: " << (t2 - t1)
<< " seconds" << std::endl;
end = omp_get_wtime();
std::cout << "Time : " << (end - start)
<< " seconds" << std::endl;
return 0;
}

27
src/scripts/plot_relative_error.py
View File

@@ -16,10 +16,10 @@ params = {
plt.rcParams.update(params)
def main():
directories = {
directories = [
"output/relative_error/RK4/",
"output/relative_error/euler/",
}
]
files = [
"4000_steps.txt",
"8000_steps.txt",
@@ -36,8 +36,12 @@ def main():
"Relative error for the RK4 method",
"Relative error for the forward Euler method"
]
fig1, axs1 = plt.subplots(2,1, sharex=True)
for i, (dir, title) in enumerate(zip(directories, titles)):
methods = [
"rk4",
"euler"
]
fig1, axs1 = plt.subplots(2,1)
for i, (dir, title) in enumerate(list(zip(directories, titles))):
max_err = []
for label, file in zip(labels, files):
with open(dir+file) as f:
@@ -47,17 +51,14 @@ def main():
max_err.append(max(r))
axs1[i].plot(t, r, label=label)
axs1[i].set(ylabel = r"Relative error $(\mu m)$") # xlabel=r"t $(\mu s)$",
axs1[i].set(xlabel=r"t $(\mu s)$", ylabel = r"relative_error $(\mu m)$")
axs1[i].legend()
axs1[i].set_title(title)
axs1[1].set_xlabel(r"t $(\mu s)$")
axs1[i].legend(loc="upper right")
# axs1[i].set_title(title)
conv_rate = 1/3 * sum([np.log2(max_err[i+1]/max_err[i])/np.log2(.5) for i in range(3)])
print(f"{methods[i]}: {conv_rate}")
conv_rate = 1/3 * sum([np.log10(max_err[i+1]/max_err[i])/np.log10(.5) for i in range(3)])
print(conv_rate)
fig1.savefig("../latex/images/relative_error.pdf")
# plt.show()
fig1.savefig("../latex/images/phase_space_2_particles_x.pdf")
if __name__ == "__main__":

95
src/test_suite.cpp
View File

@@ -10,39 +10,43 @@
* @bug No known bugs
* */
#include "PenningTrap.hpp"
#include "utils.hpp"
#include <iomanip>
#include <sstream>
#include <string>
#include "PenningTrap.hpp"
#include "utils.hpp"
/** @brief Test class for the Penning trap.
* */
class PenningTrapTest {
public:
/** @brief Test that the external E field gives correct values.
* */
static void test_external_E_field()
{
PenningTrap trap;
// Vector containing inputs and expected results
std::vector<std::pair<arma::vec, arma::vec>> tests;
std::vector<std::pair<vec_3d, vec_3d>> tests;
tests.push_back(
std::make_pair(arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.}));
std::make_pair(vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.}));
tests.push_back(std::make_pair(arma::vec{10., 0., 0.},
arma::vec{96.4852558, 0., 0.}));
tests.push_back(std::make_pair(vec_3d{10., 0., 0.},
vec_3d{96.4852558, 0., 0.}));
tests.push_back(std::make_pair(arma::vec{10., 0., 0.},
arma::vec{96.4852558, 0., 0.}));
tests.push_back(std::make_pair(vec_3d{10., 0., 0.},
vec_3d{96.4852558, 0., 0.}));
tests.push_back(std::make_pair(arma::vec{0., 10., 0.},
arma::vec{0., 96.4852558, 0.}));
tests.push_back(std::make_pair(vec_3d{0., 10., 0.},
vec_3d{0., 96.4852558, 0.}));
tests.push_back(std::make_pair(arma::vec{0., 0., 10.},
arma::vec{0., 0., -192.9705116}));
tests.push_back(std::make_pair(vec_3d{0., 0., 10.},
vec_3d{0., 0., -192.9705116}));
arma::vec result;
arma::vec v;
vec_3d result;
vec_3d v;
std::stringstream msg;
for (size_t i = 0; i < tests.size(); i++) {
v = tests.at(i).first;
@@ -52,81 +56,90 @@ public:
msg << "Testing the external E field at (" << std::setprecision(2)
<< v(0) << "," << v(1) << "," << v(2) << ").";
ASSERT(arma_vector_close_to(result, tests.at(i).second), msg.str());
ASSERT(close_to(result, tests.at(i).second), msg.str());
}
}
/** @brief Test that the external B field gives correct values.
* */
static void test_external_B_field()
{
// No point in testing at different points since it's not dependent
// on position.
PenningTrap trap;
arma::vec expected{0., 0., T};
arma::vec result = trap.external_B_field(arma::vec{0., 0., 0.});
ASSERT(arma_vector_close_to(expected, result),
vec_3d expected{0., 0., T};
vec_3d result = trap.external_B_field(vec_3d{0., 0., 0.});
ASSERT(close_to(expected, result),
"Testing the external B field at (0,0,0)");
}
/** @brief Test that the force between particles gives expected results.
* */
static void test_force_on_particle()
{
PenningTrap trap;
arma::vec v{0., 0., 0.};
vec_3d v{0., 0., 0.};
// Add particles to test
trap.add_particle(Particle(1., 40., arma::vec{0., 0., 0.}, v));
trap.add_particle(Particle(1., 40., arma::vec{1., 0., 0.}, v));
trap.add_particle(Particle(1., 40., arma::vec{0., 3., 4.}, v));
trap.add_particle(Particle(1., 40., vec_3d{0., 0., 0.}, v));
trap.add_particle(Particle(1., 40., vec_3d{1., 0., 0.}, v));
trap.add_particle(Particle(1., 40., vec_3d{0., 3., 4.}, v));
// Test p0 and p1
arma::vec expected{-1., 0., 0.};
arma::vec result = trap.force_on_particle(0, 1);
ASSERT(arma_vector_close_to(expected, result),
vec_3d expected{-1., 0., 0.};
vec_3d result = trap.force_on_particle(0, 1);
ASSERT(close_to(expected, result),
"Testing the force on a particle at (0,0,0) from a "
"particle at (1,0,0).");
// Test p0 and p2
expected = arma::vec{0, -.024, -.032};
expected = vec_3d{0, -.024, -.032};
result = trap.force_on_particle(0, 2);
ASSERT(arma_vector_close_to(expected, result),
ASSERT(close_to(expected, result),
"Testing the force on a particle at (0,0,0) from a "
"particle at (0,3,4).");
}
/** @brief Test that the total external force returns expected results
* */
static void test_total_force_external()
{
PenningTrap trap;
trap.add_particle(
Particle(1., 40., arma::vec{1., 2., 3.}, arma::vec{3., 4., 5.}));
Particle(1., 40., vec_3d{1., 2., 3.}, vec_3d{3., 4., 5.}));
arma::vec expected{395.58954878, -270.15871624, -57.89115348};
arma::vec result = trap.total_force_external(0);
ASSERT(arma_vector_close_to(expected, result),
vec_3d expected{395.58954878, -270.15871624, -57.89115348};
vec_3d result = trap.total_force_external(0);
ASSERT(close_to(expected, result),
"Testing the total external force on a particle at "
"(1,2,3) with velocity (3,4,5)");
}
/** @brief Test that the total force of all particles on a single particle
* returns expected results.
* */
static void test_total_force_particles()
{
PenningTrap trap;
trap.add_particle(
Particle(1., 40., arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.}));
Particle(1., 40., vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.}));
arma::vec expected{0., 0., 0.};
arma::vec result = trap.total_force_particles(0);
ASSERT(arma_vector_close_to(expected, result),
vec_3d expected{0., 0., 0.};
vec_3d result = trap.total_force_particles(0);
ASSERT(close_to(expected, result),
"Testing the total force of all particles on particle 0 "
"with only a single particle");
trap.add_particle(
Particle(1., 40., arma::vec{1., 0., 0.}, arma::vec{0., 0., 0.}));
Particle(1., 40., vec_3d{1., 0., 0.}, vec_3d{0., 0., 0.}));
trap.add_particle(
Particle(1., 40., arma::vec{0., 1., 0.}, arma::vec{0., 0., 0.}));
Particle(1., 40., vec_3d{0., 1., 0.}, vec_3d{0., 0., 0.}));
trap.add_particle(
Particle(1., 40., arma::vec{0., 0., 1.}, arma::vec{0., 0., 0.}));
Particle(1., 40., vec_3d{0., 0., 1.}, vec_3d{0., 0., 0.}));
expected = arma::vec(3, arma::fill::value(-3473.383325));
expected = vec_3d().fill(-3473.383325);
result = trap.total_force_particles(0);
ASSERT(arma_vector_close_to(expected, result),
ASSERT(close_to(expected, result),
"Testing the total force of all particles on particle 0 "
"with 3 other particles.");
}

4
src/utils.cpp
View File

@@ -10,8 +10,6 @@
* @bug No known bugs
* */
#include <sys/stat.h>
#include "utils.hpp"
std::string scientific_format(double d, int width, int prec)
@@ -59,7 +57,7 @@ void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
}
}
bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol)
bool close_to(arma::vec &a, arma::vec &b, double tol)
{
if (a.n_elem != b.n_elem) {
return false;