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 +102 -109
8 changed files with 102 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 79da936956 - Show all commits

8
include/Particle.hpp
View File

@ -14,12 +14,13 @@
#include <armadillo>
#include "typedefs.hpp"
#include "constants.hpp"
#include "typedefs.hpp"
/** @brief A class that holds attributes of a particle
* */
class Particle {
class Particle
{
private:
vec_3d r_vec; ///< position
vec_3d v_vec; ///< velocity
@ -37,7 +38,8 @@ public:
* @param r_vec The initial position of the particle
* @param v_vec The initial velocity of the particle
* */
Particle(vec_3d r_vec, vec_3d v_vec, double q = CA_CHARGE, double m = CA_MASS);
Particle(vec_3d r_vec, vec_3d v_vec, double q = CA_CHARGE,
double m = CA_MASS);
/** @brief Make private attributes available for PenningTrap.
* */

20
include/constants.hpp
View File

@ -12,14 +12,24 @@
#ifndef __CONST__
#define __CONST__
#define K_E 1.38935333*1e5 ///< Coulomb constant. unit: \f$\frac{u(\mu m)^3}{(\mu s)^2 e^2}\f$
/** @brief Coulomb constant. unit: \f$\frac{u(\mu m)^3}{(\mu s)^2 e^2}\f$
* */
#define K_E 1.38935333 * 1e5
#define T 9.64852558*1e1 ///< 1 Tesla. unit: \f$ \frac{u}{(\mu s) e} \f$
/** @brief 1 Tesla. unit: \f$ \frac{u}{(\mu s) e} \f$
* */
#define T 9.64852558 * 1e1
#define V 9.64852558*1e7 ///< 1 Volt. unit: \f$ \frac{u (\mu m)^2}{(\mu s)^2 e} \f$
/** @brief 1 Volt. unit: \f$ \frac{u (\mu m)^2}{(\mu s)^2 e} \f$
* */
#define V 9.64852558 * 1e7
#define CA_MASS 40.078 ///< Mass of a single calcium ion. unit: amu
/** @brief Mass of a single calcium ion. unit: amu
* */
#define CA_MASS 40.078
#define CA_CHARGE 1. ///< Charge of a singly charged calcium ion. unit: e
/** @brief Charge of a singly charged calcium ion. unit: e
* */
#define CA_CHARGE 1.
#endif

42
include/utils.hpp
View File

@ -31,10 +31,10 @@
* messages by adding the -DDBG flag when compiling.
* */
#ifdef DBG
#define DEBUG(msg) std::cout << __FILE__ << " " << __LINE__ << ": " \
<< msg << std::endl
#define DEBUG(msg) \
std::cout << __FILE__ << " " << __LINE__ << ": " << msg << std::endl
#else
#define DEBUG(msg)
#define DEBUG(msg)
#endif
/** @def ASSERT(expr)
@ -43,15 +43,14 @@
* This macro calls the m_assert function which is a more informative
* assertion function than the regular assert function from cassert.
* */
#define ASSERT(expr, msg) m_assert(expr, #expr, __METHOD_NAME__, __FILE__, \
__LINE__, msg)
#define ASSERT(expr, msg) \
m_assert(expr, #expr, __METHOD_NAME__, __FILE__, __LINE__, msg)
/** @def __METHOD_NAME__
* @brief Get the name of the current method/function without the return type.
* */
#define __METHOD_NAME__ methodName(__PRETTY_FUNCTION__)
/** @brief Turns a double into a string written in scientific format.
*
* @details The code is stolen from https://github.com/anderkve/FYS3150.
@ -62,10 +61,10 @@
*
* @return std::string
* */
std::string scientific_format(double d, int width=20, int prec=10);
std::string scientific_format(double d, int width = 20, int prec = 10);
/** @brief Turns a vector of doubles into a string written in scientific format.
/** @brief Turns a vector of doubles into a string written in scientific
* format.
*
* @details The code is stolen from https://github.com/anderkve/FYS3150.
*
@ -75,10 +74,8 @@ std::string scientific_format(double d, int width=20, int prec=10);
*
* @return std::string
* */
std::string scientific_format(const std::vector<double>& v,
int width=20,
int prec=10);
std::string scientific_format(const std::vector<double> &v, int width = 20,
int prec = 10);
/** @brief Test an expression, confirm that test is ok, or abort execution.
*
@ -92,13 +89,8 @@ std::string scientific_format(const std::vector<double>& v,
* @param line The line number where this function is called from
* @param msg The message to be displayed
* */
void m_assert(bool expr,
std::string expr_str,
std::string func,
std::string file,
int line,
std::string msg);
void m_assert(bool expr, std::string expr_str, std::string func,
std::string file, int line, std::string msg);
/** @brief Test if two armadillo vectors are close to each other.
*
@ -111,8 +103,7 @@ void m_assert(bool expr,
*
* @return bool
* */
bool close_to(arma::vec &a, arma::vec &b, double tol=1e-8);
bool close_to(arma::vec &a, arma::vec &b, double tol = 1e-8);
/** @brief Takes in the __PRETTY_FUNCTION__ string and removes the return type.
*
@ -124,16 +115,15 @@ bool close_to(arma::vec &a, arma::vec &b, double tol=1e-8);
*
* @return std::string
* */
static inline std::string methodName(const std::string& pretty_function)
static inline std::string methodName(const std::string &pretty_function)
{
size_t colons = pretty_function.find("::");
size_t begin = pretty_function.substr(0,colons).rfind(" ") + 1;
size_t begin = pretty_function.substr(0, colons).rfind(" ") + 1;
size_t end = pretty_function.rfind("(") - begin;
return pretty_function.substr(begin,end) + "()";
return pretty_function.substr(begin, end) + "()";
}
/** @brief Make path given.
*
* @details This tries to be the equivalent to "mkdir -p" and creates a new

44
src/PenningTrap.cpp
View File

@ -45,8 +45,9 @@ PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0,
void PenningTrap::set_pertubation(double f, double omega_V)
{
this->perturbation
= [f, omega_V](double t) { return f * std::cos(omega_V * t); };
this->perturbation = [f, omega_V](double t) {
return f * std::cos(omega_V * t);
};
}
void PenningTrap::reinitialize(double f, double omega_V, double t)
@ -167,8 +168,8 @@ 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)(uint)
= particle_interaction ? &PenningTrap::total_force
vec_3d (PenningTrap::*force)(uint) =
particle_interaction ? &PenningTrap::total_force
: &PenningTrap::total_force_external;
size_t size = this->particles.size();
@ -188,10 +189,10 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
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);
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 = tmp_particles;
}
@ -205,8 +206,8 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
vec_3d force_res[size];
Particle *p;
vec_3d (PenningTrap::*force)(uint)
= particle_interaction ? &PenningTrap::total_force
vec_3d (PenningTrap::*force)(uint) =
particle_interaction ? &PenningTrap::total_force
: &PenningTrap::total_force_external;
// Calculating the force for each particle is independent and therefore
@ -220,9 +221,8 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
// this as well
#pragma omp parallel for private(p)
for (size_t i = 0; i < size; i++) {
p = &this->particles[i];
p->r_vec += dt * p->v_vec;
p->v_vec += dt * force_res[i] / p->m;
p->r_vec += dt * this->particles[i].v_vec;
p->v_vec += dt * force_res[i] / this->particles[i].m;
}
this->t += dt;
@ -241,11 +241,9 @@ simulation_t PenningTrap::simulate(double time, uint steps, std::string method,
void (PenningTrap::*func)(double, bool);
if (method == "rk4") {
func = &PenningTrap::evolve_RK4;
}
else if (method == "euler") {
} else if (method == "euler") {
func = &PenningTrap::evolve_forward_euler;
}
else {
} else {
std::cout << "Not a valid method!" << std::endl;
abort();
}
@ -273,8 +271,8 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
abort();
}
simulation_t res
= this->simulate(time, steps, method, particle_interaction);
simulation_t res =
this->simulate(time, steps, method, particle_interaction);
std::ofstream ofile;
@ -308,11 +306,9 @@ double PenningTrap::fraction_of_particles_left(double time, uint steps,
void (PenningTrap::*func)(double, bool);
if (method == "rk4") {
func = &PenningTrap::evolve_RK4;
}
else if (method == "euler") {
} else if (method == "euler") {
func = &PenningTrap::evolve_forward_euler;
}
else {
} else {
std::cout << "Not a valid method!" << std::endl;
abort();
}
@ -323,7 +319,7 @@ double PenningTrap::fraction_of_particles_left(double time, uint steps,
int particles_left = 0;
// A reduction is perfect here
// 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) {

30
src/main.cpp
View File

@ -42,8 +42,8 @@ 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))
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)};
@ -156,8 +156,8 @@ void potential_resonance_wide_sweep()
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) + 1;
size_t freq_iterations =
(size_t)((freq_end - freq_start) / freq_increment) + 1;
double res[4][freq_iterations];
@ -176,14 +176,14 @@ void potential_resonance_wide_sweep()
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((unsigned int)100);
PenningTrap trap((uint)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(amplitudes[i], res[0][j]);
res[i + 1][j]
= trap.fraction_of_particles_left(time, N, "rk4", false);
res[i + 1][j] =
trap.fraction_of_particles_left(time, N, "rk4", false);
}
}
}
@ -213,8 +213,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);
double res[4][freq_iterations];
@ -233,14 +233,14 @@ void potential_resonance_narrow_sweep()
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((unsigned int)100);
PenningTrap trap((uint)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(amplitudes[i], res[0][j]);
res[i + 1][j]
= trap.fraction_of_particles_left(time, N, "rk4", false);
res[i + 1][j] =
trap.fraction_of_particles_left(time, N, "rk4", false);
}
}
}
@ -270,8 +270,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);
double res[4][freq_iterations];
@ -289,7 +289,7 @@ void potential_resonance_narrow_sweep_interaction()
{
// Each thread creates a PenningTrap instance and reuses it throughout
// the sweep.
PenningTrap trap((unsigned int)100);
PenningTrap trap((uint)100);
#pragma omp for collapse(2)
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < freq_iterations; j++) {

34
src/test_suite.cpp
View File

@ -20,7 +20,8 @@
/** @brief Test class for the Penning trap.
* */
class PenningTrapTest {
class PenningTrapTest
{
public:
/** @brief Test that the external E field gives correct values.
* */
@ -34,17 +35,17 @@ public:
tests.push_back(
std::make_pair(vec_3d{0., 0., 0.}, vec_3d{0., 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(vec_3d{10., 0., 0.}, vec_3d{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(vec_3d{10., 0., 0.}, vec_3d{96.4852558, 0., 0.}));
tests.push_back(std::make_pair(vec_3d{0., 10., 0.},
vec_3d{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(vec_3d{0., 0., 10.},
vec_3d{0., 0., -192.9705116}));
tests.push_back(
std::make_pair(vec_3d{0., 0., 10.}, vec_3d{0., 0., -192.9705116}));
vec_3d result;
vec_3d v;
@ -106,8 +107,7 @@ public:
static void test_total_force_external()
{
PenningTrap trap;
trap.add_particle(
Particle(vec_3d{1., 2., 3.}, vec_3d{3., 4., 5.}));
trap.add_particle(Particle(vec_3d{1., 2., 3.}, vec_3d{3., 4., 5.}));
vec_3d expected{395.58954878, -270.15871624, -57.89115348};
vec_3d result = trap.total_force_external(0);
@ -122,8 +122,7 @@ public:
static void test_total_force_particles()
{
PenningTrap trap;
trap.add_particle(
Particle(vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.}));
trap.add_particle(Particle(vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.}));
vec_3d expected{0., 0., 0.};
vec_3d result = trap.total_force_particles(0);
@ -131,12 +130,9 @@ public:
"Testing the total force of all particles on particle 0 "
"with only a single particle");
trap.add_particle(
Particle(vec_3d{1., 0., 0.}, vec_3d{0., 0., 0.}));
trap.add_particle(
Particle(vec_3d{0., 1., 0.}, vec_3d{0., 0., 0.}));
trap.add_particle(
Particle(vec_3d{0., 0., 1.}, vec_3d{0., 0., 0.}));
trap.add_particle(Particle(vec_3d{1., 0., 0.}, vec_3d{0., 0., 0.}));
trap.add_particle(Particle(vec_3d{0., 1., 0.}, vec_3d{0., 0., 0.}));
trap.add_particle(Particle(vec_3d{0., 0., 1.}, vec_3d{0., 0., 0.}));
expected = vec_3d().fill(-138935.333);
result = trap.total_force_particles(0);

15
src/utils.cpp
View File

@ -19,7 +19,8 @@ std::string scientific_format(double d, int width, int prec)
return ss.str();
}
std::string scientific_format(const std::vector<double> &v, int width, int prec)
std::string scientific_format(const std::vector<double> &v, int width,
int prec)
{
std::stringstream ss;
for (double elem : v) {
@ -32,8 +33,7 @@ static void print_message(std::string msg)
{
if (msg.size() > 0) {
std::cout << "message: " << msg << "\n\n";
}
else {
} else {
std::cout << "\n";
}
}
@ -47,8 +47,7 @@ void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
if (expr) {
std::cout << "\x1B[32mOK\033[0m\n";
print_message(msg);
}
else {
} else {
std::cout << "\x1B[31mFAIL\033[0m\n";
print_message(msg);
std::cout << file << " " << line << ": Assertion \"" << expr_str
@ -85,11 +84,11 @@ bool mkpath(std::string path, int mode)
pos = path.find('/', pos);
if (pos != std::string::npos) {
cur_dir = path.substr(0, pos);
if (mkdir(cur_dir.c_str(), mode) != 0 && stat(cur_dir.c_str(), &buf) != 0) {
if (mkdir(cur_dir.c_str(), mode) != 0
&& stat(cur_dir.c_str(), &buf) != 0) {
return -1;
}
}
else {
} else {
break;
}
}