Merge pull request #5 from FYS3150-G2-2023/coryab/implement-RK4

Coryab/implement rk4

.clang-format

@@ -1,4 +1,8 @@
+UseTab: Never
 IndentWidth: 4
+TabWidth: 4
+AccessModifierOffset: -4
+IndentAccessModifiers: false
 AllowShortFunctionsOnASingleLine: false
 BreakBeforeBraces: Custom
 BraceWrapping:

src/PenningTrap.cpp

@@ -103,6 +103,65 @@ arma::vec PenningTrap::total_force(int i)
 
 void PenningTrap::evolve_RK4(double dt)
 {
+    std::vector<Particle> 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; i<size; i++) {
+        k_v[i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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; i<size; i++) {
+        k_v[1*size + i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[1*size + i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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; i<size; i++) {
+        k_v[2*size + i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[2*size + i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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; i<size; i++) {
+        k_v[3*size + i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[3*size + i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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)

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()
 
 

src/main.cpp

@@ -21,14 +21,16 @@
 #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 r = arma::vec(3).randn() * 0.1 *
-        arma::vec v = arma::vec(3).randn() * 0.1 * trap.get_d();  // random initial velocity
+                      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));
     }
 
@@ -45,7 +47,7 @@ void euler_100_particles()
         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;
@@ -54,6 +56,7 @@ void euler_100_particles()
     arma::vec::fixed<3> cur_elem;
 
     mkdir("output", 0777);
+    mkdir("output/simulate_100_particles", 0777);
 
     std::ofstream ofile;
 
@@ -61,12 +64,11 @@ void euler_100_particles()
 #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");
+        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) << "," 
+            ofile << cur_elem(0) << "," << cur_elem(1) << "," << cur_elem(2)
-                << cur_elem(1) << "," 
+                  << "\n";
-                << 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();
 

src/test_suite.cpp

@@ -26,8 +26,8 @@ public:
         // Vector containing inputs and expected results
         std::vector<std::pair<arma::vec, arma::vec>> tests;
 
-        tests.push_back(std::make_pair(arma::vec{0.,0.,0.}, 
+        tests.push_back(
-                                       arma::vec{0.,0.,0.}));
+            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.}));
@@ -41,7 +41,6 @@ public:
         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;
@@ -57,11 +56,10 @@ public:
         }
     }
 
-
     static void test_external_B_field()
     {
-        // No point in testing at different points since it's not dependent on 
+        // No point in testing at different points since it's not dependent
-        // position.
+        // on position.
         PenningTrap trap;
         arma::vec expected{0., 0., T};
         arma::vec result = trap.external_B_field(arma::vec{0., 0., 0.});
@@ -69,7 +67,6 @@ public:
                "Testing the external B field at (0,0,0)");
     }
 
-
     static void test_force_on_particle()
     {
         PenningTrap trap;
@@ -98,8 +95,8 @@ public:
     static void test_total_force_external()
     {
         PenningTrap trap;
-        trap.add_particle(Particle(1.,40.,arma::vec{1.,2.,3.},
+        trap.add_particle(
-                                   arma::vec{3.,4.,5.}));
+            Particle(1., 40., arma::vec{1., 2., 3.}, arma::vec{3., 4., 5.}));
 
         arma::vec expected{395.58954878, -270.15871624, -57.89115348};
         arma::vec result = trap.total_force_external(0);
@@ -111,8 +108,8 @@ public:
     static void test_total_force_particles()
     {
         PenningTrap trap;
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,0.,0.},
+        trap.add_particle(
-                                   arma::vec{0.,0.,0.}));
+            Particle(1., 40., arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.}));
 
         arma::vec expected{0., 0., 0.};
         arma::vec result = trap.total_force_particles(0);
@@ -120,12 +117,12 @@ public:
                "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.},
+        trap.add_particle(
-                                   arma::vec{0.,0.,0.}));
+            Particle(1., 40., arma::vec{1., 0., 0.}, arma::vec{0., 0., 0.}));
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,1.,0.},
+        trap.add_particle(
-                                   arma::vec{0.,0.,0.}));
+            Particle(1., 40., arma::vec{0., 1., 0.}, arma::vec{0., 0., 0.}));
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,0.,1.},
+        trap.add_particle(
-                                   arma::vec{0.,0.,0.}));
+            Particle(1., 40., arma::vec{0., 0., 1.}, arma::vec{0., 0., 0.}));
 
         expected = arma::vec(3, arma::fill::value(-3473.383325));
         result = trap.total_force_particles(0);
@@ -135,7 +132,6 @@ public:
     }
 };
 
-
 int main()
 {
     PenningTrapTest::test_external_E_field();

src/utils.cpp

@@ -37,12 +37,8 @@ static void print_message(std::string msg)
     }
 }
 
-void m_assert(bool expr, 
+void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
-              std::string expr_str, 
+              int line, std::string msg)
-              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 \"" 
+        std::cout << file << " " << line << ": Assertion \"" << expr_str
-            << expr_str << "\" Failed\n\n";
+                  << "\" Failed\n\n";
         abort();
     }
 }