Bunch of stuff

- Add reinitialize to reuse PenningTrap instance and its particles
- Reformatting
- Refactor fraction_of_particles_left to not use simulate

include/PenningTrap.hpp

@@ -74,8 +74,8 @@ public:
      * */
     PenningTrap(
         double B_0 = T,
-        std::function<double(double)> V_0 =
+        std::function<double(double)> V_0
-            [](double t) { return 25. * V / 1000.; },
+        = [](double t) { return 25. * V / 1000.; },
         double d = 500., double t = 0.);
 
     /** @brief Constructor for the PenningTrap class
@@ -88,8 +88,8 @@ public:
      * */
     PenningTrap(
         unsigned int i, double B_0 = T,
-        std::function<double(double)> V_0 =
+        std::function<double(double)> V_0
-            [](double t) { return 25. * V / 1000.; },
+        = [](double t) { return 25. * V / 1000.; },
         double d = 500., double t = 0.);
 
     /** @brief Constructor for the PenningTrap class
@@ -102,10 +102,21 @@ public:
      * */
     PenningTrap(
         std::vector<Particle> particles, double B_0 = T,
-        std::function<double(double)> V_0 =
+        std::function<double(double)> V_0
-            [](double t) { return 25. * V / 1000.; },
+        = [](double t) { return 25. * V / 1000.; },
         double d = 500., double t = 0.);
 
+    /** @brief Give all particles new positions and velocities, and change t 
+     *  and V_0.
+     *
+     *  @param V_0 The tiome dependent applied potential
+     *  @param t The starting time
+     * */
+    void reinitialize(
+        std::function<double(double)> V_0
+        = [](double t) { return 25. * V / 1000.; },
+        double t = 0.);
+
     /** @brief Add a particle to the system
      *
      *  @param particle The particle to add to the Penning trap
@@ -151,7 +162,8 @@ public:
      * */
     vec_3d total_force_external(unsigned int i);
 
-    /** @brief Calculate the total force on a particle p_i from other particles.
+    /** @brief Calculate the total force on a particle p_i from other
+     * particles.
      *
      *  @param i The index of particle p_i
      *
@@ -204,7 +216,8 @@ public:
      *  @param particle_interaction Turn particle interactions on/off
      * */
     void write_simulation_to_dir(std::string path, double time,
-                                 unsigned int steps, std::string method = "rk4",
+                                 unsigned int steps,
+                                 std::string method = "rk4",
                                  bool particle_interaction = true);
 
     /** @brief Simulate and calculate what fraction of particles are still

src/PenningTrap.cpp

@@ -11,7 +11,10 @@
  * */
 
 #include "PenningTrap.hpp"
+#include "typedefs.hpp"
+#include <algorithm>
 #include <functional>
+#include <vector>
 
 PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0,
                          double d, double t)
@@ -26,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;
+        this->particles.push_back(Particle(1., 40., vec_3d(vec_3d().randn() * .1 * this->d),
-        v = vec_3d().randn() * .1 * this->d;
+                                    vec_3d(vec_3d().randn() * .1 * this->d)));
-        this->add_particle(Particle(1., 40., r, v));
     }
 }
 
@@ -41,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) {
@@ -109,15 +120,16 @@ vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
 
 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.};
     }
 
-    return vec_3d(p.q
+    return vec_3d(
-                  * (this->external_E_field(p.r_vec)
+        p->q
-                     + 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,9 +146,10 @@ vec_3d PenningTrap::total_force_particles(unsigned int i)
 
 vec_3d PenningTrap::total_force(unsigned int i)
 {
-    return arma::norm(this->particles[i].r_vec) > this->d
+    if (arma::norm(this->particles[i].r_vec) > this->d) {
-               ? vec_3d{0., 0., 0.}
+        return vec_3d{0., 0., 0.};
-               : vec_3d(this->total_force_external(i)
+    }
+    return vec_3d(this->total_force_external(i)
                   - this->total_force_particles(i));
 }
 
@@ -146,13 +159,9 @@ 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);
+    vec_3d (PenningTrap::*force)(unsigned int)
-    if (particle_interaction) {
+        = particle_interaction ? &PenningTrap::total_force
-        force = &PenningTrap::total_force;
+                               : &PenningTrap::total_force_external;
-    }
-    else {
-        force = &PenningTrap::total_force_external;
-    }
 
     size_t size = this->particles.size();
 
@@ -163,7 +172,9 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
         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 < size; j++) {
             this->k_v[i][j] = (this->*force)(j) / this->particles[j].m;
@@ -176,6 +187,7 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
         }
         this->particles = tmp_particles;
     }
+
     this->t += dt;
 }
 
@@ -185,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);
+    vec_3d (PenningTrap::*force)(unsigned int)
-    if (particle_interaction) {
+        = particle_interaction ? &PenningTrap::total_force
-        force = &PenningTrap::total_force_external;
+                               : &PenningTrap::total_force_external;
-    }
-    else {
-        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];
@@ -251,8 +263,8 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
         path += '/';
     }
     if (mkpath(path, 0777) != 0) {
-        std::cout << "Hello" << std::endl;
+        std::cout << "Failed to make path" << std::endl;
-        return;
+        abort();
     }
 
     simulation_t res
@@ -260,18 +272,22 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
 
     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) << ","
+            ofile << scientific_format(vec(0), 10, 8) << ','
-                  << scientific_format(vec(1), 8, 10) << ","
+                  << scientific_format(vec(1), 10, 8) << ','
-                  << scientific_format(vec(2), 8, 10) << "\n";
+                  << scientific_format(vec(2), 10, 8) << '\n';
         }
         ofile.close();
     }
@@ -281,26 +297,33 @@ double PenningTrap::fraction_of_particles_left(double time, unsigned int steps,
                                                std::string method,
                                                bool particle_interaction)
 {
-    simulation_t res
+    double dt = time / (double)steps;
-        = this->simulate(time, steps, method, particle_interaction);
+
+    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) {
+    // A reduction is perfect here
-        if (arma::norm(p.r_vec) < this->d) {
+    #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;
-}