diff --git a/.clang-format b/.clang-format index 2298bfb..d4d13ad 100644 --- a/.clang-format +++ b/.clang-format @@ -1,4 +1,8 @@ +UseTab: Never IndentWidth: 4 +TabWidth: 4 +AccessModifierOffset: -4 +IndentAccessModifiers: false AllowShortFunctionsOnASingleLine: false BreakBeforeBraces: Custom BraceWrapping: diff --git a/images/100_particles.gif b/images/100_particles.gif new file mode 100644 index 0000000..3dde053 Binary files /dev/null and b/images/100_particles.gif differ diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index b4b2dd8..5d9a16b 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -103,6 +103,65 @@ arma::vec PenningTrap::total_force(int i) void PenningTrap::evolve_RK4(double dt) { + std::vector tmp_particles = this->particles; + + arma::vec::fixed<3> *k_v = new arma::vec::fixed<3>[this->particles.size()*4]; + arma::vec::fixed<3> *k_r = new arma::vec::fixed<3>[this->particles.size()*4]; + + int size = this->particles.size(); + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[i]; + p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[i]; + } + + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[1*size + i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[1*size + i]; + p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[1*size + i]; + } + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[2*size + i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + dt*k_v[2*size + i]; + p->r_vec = tmp_particles.at(i).r_vec + dt*k_r[2*size + i]; + } + + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[3*size + i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + dt*(k_v[i] + k_v[size + i] + k_v[2*size + i] + k_v[3*size + i])/6; + p->r_vec = tmp_particles.at(i).r_vec + dt*(k_r[i] + k_r[size + i] + k_r[2*size + i] + k_r[3*size + i])/6; + } + + delete [] k_v; + delete [] k_r; } void PenningTrap::evolve_forward_euler(double dt) diff --git a/src/test.py b/src/animate_100_particles.py similarity index 80% rename from src/test.py rename to src/animate_100_particles.py index a215cad..a62510f 100644 --- a/src/test.py +++ b/src/animate_100_particles.py @@ -7,11 +7,11 @@ def get_data(files): res = [] for file in files: arr = [[], [], []] - with open(file) as f: + with open(file, encoding="utf8") as f: lines = f.readlines() for line in lines: - xi,yi,zi = map(lambda x: float(x), line.strip().split(",")) + xi,yi,zi = map(float, line.strip().split(",")) arr[0].append(xi) arr[1].append(yi) arr[2].append(zi) @@ -27,11 +27,12 @@ def update(num, lines, arr): def animate(): + plt.style.use("dark_background") fig = plt.figure() ax = fig.add_subplot(projection="3d") - arr = get_data([f"output/p{i}.txt" for i in range(100)]) + arr = get_data([f"output/p{i}_RK4.txt" for i in range(100)]) arr = arr[:,:,::10] @@ -53,7 +54,7 @@ def animate(): blit=False) - # ani.save("100_particles.gif", writer=animation.FFMpegFileWriter(fps=30)) + ani.save("100_particles.gif", writer=animation.FFMpegFileWriter(fps=30)) plt.show() diff --git a/src/main.cpp b/src/main.cpp index 14e8fa7..1c55279 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -21,52 +21,54 @@ #define CHARGE 1. #define MASS 40. // unit: amu -void euler_100_particles() +void simulate_100_particles() { PenningTrap trap; // Add particles inside trap - for (int i=0; i < PARTICLES; i++) { - arma::vec r = arma::vec(3).randn() * 0.1 * trap.get_d(); // random initial position - arma::vec v = arma::vec(3).randn() * 0.1 * trap.get_d(); // random initial velocity + for (int i = 0; i < PARTICLES; i++) { + arma::vec r = arma::vec(3).randn() * 0.1 * + trap.get_d(); // random initial position + arma::vec v = arma::vec(3).randn() * 0.1 * + trap.get_d(); // random initial velocity trap.add_particle(Particle(CHARGE, MASS, r, v)); } double time = 50.; // microseconds - double dt = time / (double) N; + double dt = time / (double)N; auto res = new arma::vec::fixed<3>[PARTICLES][N]; int counter = 0; // Get the path of all particles - for (int j=0; j < N; j++) { - #pragma omp parallel for - for (int i=0; i < PARTICLES; i++) { + for (int j = 0; j < N; j++) { +#pragma omp parallel for + for (int i = 0; i < PARTICLES; i++) { res[i][j] = trap.get_particle(i); } - trap.evolve_forward_euler(dt); + trap.evolve_RK4(dt); } std::cout << counter << std::endl; - arma::vec::fixed<3>* cur_row; + arma::vec::fixed<3> *cur_row; arma::vec::fixed<3> cur_elem; mkdir("output", 0777); + mkdir("output/simulate_100_particles", 0777); std::ofstream ofile; - // Write particle paths to file - #pragma omp parallel for private(cur_row, cur_elem, ofile) - for (int i=0; i < PARTICLES; i++) { +// Write particle paths to file +#pragma omp parallel for private(cur_row, cur_elem, ofile) + for (int i = 0; i < PARTICLES; i++) { cur_row = res[i]; - ofile.open("output/p" + std::to_string(i) + ".txt"); - for (int j=0; j < N; j++) { + ofile.open("output/simulate_100_particles/p" + std::to_string(i) + ".txt"); + for (int j = 0; j < N; j++) { cur_elem = cur_row[j]; - ofile << cur_elem(0) << "," - << cur_elem(1) << "," - << cur_elem(2) << "\n"; + ofile << cur_elem(0) << "," << cur_elem(1) << "," << cur_elem(2) + << "\n"; } ofile.close(); } @@ -76,7 +78,7 @@ int main() { double start = omp_get_wtime(); - euler_100_particles(); + simulate_100_particles(); double end = omp_get_wtime(); diff --git a/src/test_suite.cpp b/src/test_suite.cpp index 17f9c0e..81ea5c1 100644 --- a/src/test_suite.cpp +++ b/src/test_suite.cpp @@ -26,71 +26,68 @@ public: // Vector containing inputs and expected results std::vector> tests; - tests.push_back(std::make_pair(arma::vec{0.,0.,0.}, - arma::vec{0.,0.,0.})); + tests.push_back( + std::make_pair(arma::vec{0., 0., 0.}, arma::vec{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(arma::vec{10., 0., 0.}, + arma::vec{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(arma::vec{10., 0., 0.}, + arma::vec{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(arma::vec{0.,0.,10.}, - arma::vec{0.,0.,-192.9705116})); + tests.push_back(std::make_pair(arma::vec{0., 10., 0.}, + arma::vec{0., 96.4852558, 0.})); + tests.push_back(std::make_pair(arma::vec{0., 0., 10.}, + arma::vec{0., 0., -192.9705116})); arma::vec result; arma::vec v; std::stringstream msg; - for (int i=0; i < tests.size(); i++) { + for (int i = 0; i < tests.size(); i++) { v = tests.at(i).first; result = trap.external_E_field(v); msg.str(""); - msg << "Testing the external E field at (" << std::setprecision(2) + 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()); } } - static void test_external_B_field() { - // No point in testing at different points since it's not dependent on - // position. + // 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.}); + 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), "Testing the external B field at (0,0,0)"); } - static void test_force_on_particle() { PenningTrap trap; - arma::vec v{0.,0.,0.}; + arma::vec 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., 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)); // Test p0 and p1 - arma::vec expected{-1.,0.,0.}; + arma::vec expected{-1., 0., 0.}; arma::vec result = trap.force_on_particle(0, 1); ASSERT(arma_vector_close_to(expected, result), - "Testing the force on a particle at (0,0,0) from a " - "particle at (1,0,0)."); + "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}; result = trap.force_on_particle(0, 2); - ASSERT(arma_vector_close_to(expected, result), + ASSERT(arma_vector_close_to(expected, result), "Testing the force on a particle at (0,0,0) from a " "particle at (0,3,4)."); } @@ -98,10 +95,10 @@ public: static void test_total_force_external() { PenningTrap trap; - trap.add_particle(Particle(1.,40.,arma::vec{1.,2.,3.}, - arma::vec{3.,4.,5.})); + trap.add_particle( + Particle(1., 40., arma::vec{1., 2., 3.}, arma::vec{3., 4., 5.})); - arma::vec expected{395.58954878,-270.15871624,-57.89115348}; + arma::vec expected{395.58954878, -270.15871624, -57.89115348}; arma::vec result = trap.total_force_external(0); ASSERT(arma_vector_close_to(expected, result), "Testing the total external force on a particle at " @@ -111,23 +108,23 @@ public: static void test_total_force_particles() { PenningTrap trap; - trap.add_particle(Particle(1.,40.,arma::vec{0.,0.,0.}, - arma::vec{0.,0.,0.})); + trap.add_particle( + Particle(1., 40., arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.})); - arma::vec expected{0.,0.,0.}; + arma::vec expected{0., 0., 0.}; arma::vec result = trap.total_force_particles(0); ASSERT(arma_vector_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.})); - trap.add_particle(Particle(1.,40.,arma::vec{0.,1.,0.}, - arma::vec{0.,0.,0.})); - trap.add_particle(Particle(1.,40.,arma::vec{0.,0.,1.}, - arma::vec{0.,0.,0.})); + trap.add_particle( + Particle(1., 40., arma::vec{1., 0., 0.}, arma::vec{0., 0., 0.})); + trap.add_particle( + Particle(1., 40., arma::vec{0., 1., 0.}, arma::vec{0., 0., 0.})); + trap.add_particle( + Particle(1., 40., arma::vec{0., 0., 1.}, arma::vec{0., 0., 0.})); - expected = arma::vec(3,arma::fill::value(-3473.383325)); + expected = arma::vec(3, arma::fill::value(-3473.383325)); result = trap.total_force_particles(0); ASSERT(arma_vector_close_to(expected, result), "Testing the total force of all particles on particle 0 " @@ -135,8 +132,7 @@ public: } }; - -int main() +int main() { PenningTrapTest::test_external_E_field(); PenningTrapTest::test_external_B_field(); diff --git a/src/utils.cpp b/src/utils.cpp index 9694610..b55fa34 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -18,10 +18,10 @@ std::string scientific_format(double d, int width, int prec) return ss.str(); } -std::string scientific_format(const std::vector& v, int width, int prec) +std::string scientific_format(const std::vector &v, int width, int prec) { std::stringstream ss; - for(double elem : v) { + for (double elem : v) { ss << scientific_format(elem, width, prec); } return ss.str(); @@ -37,12 +37,8 @@ static void print_message(std::string msg) } } -void m_assert(bool expr, - std::string expr_str, - std::string f, - std::string file, - int line, - std::string msg) +void m_assert(bool expr, std::string expr_str, std::string f, std::string file, + int line, std::string msg) { std::string new_assert(f.size() + (expr ? 4 : 6), '-'); std::cout << "\x1B[36m" << new_assert << "\033[0m\n"; @@ -54,8 +50,8 @@ void m_assert(bool expr, else { std::cout << "\x1B[31mFAIL\033[0m\n"; print_message(msg); - std::cout << file << " " << line << ": Assertion \"" - << expr_str << "\" Failed\n\n"; + std::cout << file << " " << line << ": Assertion \"" << expr_str + << "\" Failed\n\n"; abort(); } } @@ -66,7 +62,7 @@ bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol) return false; } - for (int i=0; i < a.n_elem; i++) { + for (int i = 0; i < a.n_elem; i++) { if (std::abs(a(i) - b(i)) >= tol) { return false; }