Run clang-format with new rules

Coryab/create tests

.clang-format

@@ -0,0 +1,10 @@
+UseTab: Never
+IndentWidth: 4
+TabWidth: 4
+AccessModifierOffset: -4
+IndentAccessModifiers: false
+AllowShortFunctionsOnASingleLine: false
+BreakBeforeBraces: Custom
+BraceWrapping:
+  AfterFunction: true
+  BeforeElse: true

include/PenningTrap.hpp

@@ -13,9 +13,13 @@
 #define __PENNING_TRAP__
 
 #include <armadillo>
+#include <omp.h>
 
-#include "constants.hpp"
 #include "Particle.hpp"
+#include "constants.hpp"
+
+#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \
+  initializer( omp_priv = omp_orig )
 
 /** @brief A class that simulates a Penning trap.
  *

include/utils.hpp

@@ -15,15 +15,16 @@
 #ifndef __UTILS__
 #define __UTILS__
 
-#include <string>
-#include <vector>
+#include <armadillo>
 #include <iomanip>
 #include <sstream>
+#include <string>
+#include <vector>
 
 /** @def DEBUG(msg)
  *  @brief Writes a debug message
  *
- *  This function writes a debug message that includes the filename,
+ *  This macro writes a debug message that includes the filename,
  *  line number, and a custom message. The function is wrapped in an ifdef
  *  that checks if DBG is defined, so one can choose to display the debug 
  *  messages by adding the -DDBG flag when compiling.
@@ -35,6 +36,18 @@
     #define DEBUG(msg)
 #endif
 
+/** @def ASSERT(expr)
+ *  @brief A prettier assertion function.
+ *
+ *  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 __METHOD_NAME__ methodName(__PRETTY_FUNCTION__)
+
+
 /** Code stolen from https://github.com/anderkve/FYS3150
  *  Header: https://github.com/anderkve/FYS3150/blob/master/code_examples/compilation_linking/example_1/include/utils.hpp
  *  Source: https://github.com/anderkve/FYS3150/blob/master/code_examples/compilation_linking/example_1/src/utils.cpp
@@ -62,4 +75,47 @@ 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.
+ *
+ *  This function takes in an expression and prints an OK message if it's
+ *  true, or it prints a fail message and aborts execution if it fails.
+ *
+ *  @param expr The expression to be evaluated
+ *  @param expr_str The stringified version of the expression
+ *  @param func The function name of the caller
+ *  @param file The file of the caller
+ *  @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);
+
+
+/** @brief Test if two armadillo vectors are close to each other.
+ *
+ *  This function takes in 2 vectors and checks if they are approximately
+ *  equal to each other given a tolerance.
+ *
+ *  @param a Vector a
+ *  @param b Vector b
+ *  @param tol The tolerance
+ *
+ *  @return Boolean
+ * */
+bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8);
+
+
+static inline std::string methodName(const std::string& prettyFunction)
+{
+    size_t colons = prettyFunction.find("::");
+    size_t begin = prettyFunction.substr(0,colons).rfind(" ") + 1;
+    size_t end = prettyFunction.rfind("(") - begin;
+
+    return prettyFunction.substr(begin,end) + "()";
+}
+
 #endif

src/Makefile

@@ -27,7 +27,7 @@ all: test_suite main
 main: main.o $(LIBOBJS) $(CLASSOBJS)
 	$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) -I$(INCLUDE) $(OPENMP)
 
-test_suite: test_suite.o $(LIBOBJS)
+test_suite: test_suite.o $(LIBOBJS) $(CLASSOBJS)
 	$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) -I$(INCLUDE) $(OPENMP)
 
 # Rule for object files

src/PenningTrap.cpp

@@ -13,15 +13,9 @@
  *  @todo Implement evolve_forward_euler
  * */
 
-#include "utils.hpp"
 #include "PenningTrap.hpp"
 #include "constants.hpp"
-#include <algorithm>
-#include <stdexcept>
-#include <omp.h>
-
-#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \
-  initializer( omp_priv = omp_orig )
+#include "utils.hpp"
 
 PenningTrap::PenningTrap(double B_0, double V_0, double d)
 {
@@ -38,20 +32,17 @@ void PenningTrap::add_particle(Particle particle)
 arma::vec PenningTrap::external_E_field(arma::vec r)
 {
     arma::vec::fixed<3> res;
-    double f = this->V_0/(this->d*this->d);
+    double f = this->V_0 / (this->d * this->d);
     res(0) = r(0);
     res(1) = r(1);
-    res(2) = -2.*r(2);
+    res(2) = -2. * r(2);
 
-    return f*res;
+    return f * res;
 }
 
 arma::vec PenningTrap::external_B_field(arma::vec r)
 {
-    arma::vec::fixed<3> res;
-    res(0) = 0.;
-    res(1) = 0.;
-    res(2) = this->B_0;
+    arma::vec::fixed<3> res{0., 0., this->B_0};
 
     return res;
 }
@@ -59,14 +50,14 @@ arma::vec PenningTrap::external_B_field(arma::vec r)
 arma::vec PenningTrap::force_on_particle(int i, int j)
 {
     // Calculate the difference between the particles' position
-    arma::vec::fixed<3> res = this->particles.at(i).r_vec 
-        - this->particles.at(j).r_vec;
+    arma::vec::fixed<3> res =
+        this->particles.at(i).r_vec - this->particles.at(j).r_vec;
 
     // Get the distance between the particles
     double norm = arma::norm(res);
 
     // Multiply res with p_j's charge divided by the norm cubed
-    res *= this->particles.at(j).q/(norm*norm*norm);
+    res *= this->particles.at(j).q / (norm * norm * norm);
 
     return res;
 }
@@ -75,16 +66,13 @@ arma::vec PenningTrap::total_force_external(int i)
 {
     Particle p = this->particles.at(i);
 
-    arma::vec::fixed<3> v_cross_B;
-
     arma::vec::fixed<3> B = this->external_B_field(p.r_vec);
 
-    v_cross_B(0) = p.v_vec(1)*B(2) - p.v_vec(2)*B(1);
-    v_cross_B(1) = p.v_vec(2)*B(0) - p.v_vec(0)*B(2);
-    v_cross_B(2) = p.v_vec(0)*B(1) - p.v_vec(1)*B(0);
+    arma::vec::fixed<3> v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1),
+                                  p.v_vec(2) * B(0) - p.v_vec(0) * B(2),
+                                  p.v_vec(0) * B(1) - p.v_vec(1) * B(0)};
 
-    arma::vec force = p.q
-        *(this->external_E_field(p.r_vec) + v_cross_B);
+    arma::vec force = p.q * (this->external_E_field(p.r_vec) + v_cross_B);
 
     return force;
 }
@@ -95,7 +83,7 @@ arma::vec PenningTrap::total_force_particles(int i)
 
     arma::vec res(3);
 
-    for (int j=0; j < this->particles.size(); j++) {
+    for (int j = 0; j < this->particles.size(); j++) {
         if (i == j) {
             continue;
         }
@@ -103,7 +91,7 @@ arma::vec PenningTrap::total_force_particles(int i)
         res += this->force_on_particle(i, j);
     }
 
-    res *= K_E*(p.q/p.m);
+    res *= K_E * (p.q / p.m);
 
     return res;
 }
@@ -115,7 +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)
@@ -124,11 +170,11 @@ void PenningTrap::evolve_forward_euler(double dt)
 
     Particle *p;
 
-    #pragma omp parallel for private(p)
-    for (int i=0; i < this->particles.size(); i++) {
+#pragma omp parallel for private(p)
+    for (int i = 0; i < this->particles.size(); i++) {
         p = &new_state.at(i);
-        p->v_vec += dt*this->total_force(i)/new_state.at(i).m;
-        p->r_vec += dt*this->particles.at(i).v_vec;
+        p->v_vec += dt * this->total_force(i) / new_state.at(i).m;
+        p->r_vec += dt * this->particles.at(i).v_vec;
     }
 
     this->particles = new_state;

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

@@ -10,9 +10,7 @@
  *  @bug No known bugs
  * */
 
-#include <filesystem>
 #include <fstream>
-#include <string>
 #include <omp.h>
 #include <sys/stat.h>
 
@@ -23,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();
     }
@@ -78,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

@@ -9,7 +9,135 @@
  *
  *  @bug No known bugs
  * */
-int main() 
+
+#include "PenningTrap.hpp"
+#include "utils.hpp"
+
+#include <iomanip>
+#include <sstream>
+#include <string>
+
+class PenningTrapTest {
+public:
+    static void test_external_E_field()
+    {
+        PenningTrap trap;
+
+        // 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.}, 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{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++) {
+            v = tests.at(i).first;
+            result = trap.external_E_field(v);
+
+            msg.str("");
+            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.
+        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),
+               "Testing the external B field at (0,0,0)");
+    }
+
+    static void test_force_on_particle()
+    {
+        PenningTrap trap;
+        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));
+
+        // 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),
+               "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),
+               "Testing the force on a particle at (0,0,0) from a "
+               "particle at (0,3,4).");
+    }
+
+    static void test_total_force_external()
+    {
+        PenningTrap trap;
+        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 result = trap.total_force_external(0);
+        ASSERT(arma_vector_close_to(expected, result),
+               "Testing the total external force on a particle at "
+               "(1,2,3) with velocity (3,4,5)");
+    }
+
+    static void test_total_force_particles()
+    {
+        PenningTrap trap;
+        trap.add_particle(
+            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);
+        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.}));
+
+        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 "
+               "with 3 other particles.");
+    }
+};
+
+int main()
 {
+    PenningTrapTest::test_external_E_field();
+    PenningTrapTest::test_external_B_field();
+    PenningTrapTest::test_force_on_particle();
+    PenningTrapTest::test_total_force_external();
+    PenningTrapTest::test_total_force_particles();
     return 0;
 }

src/utils.cpp

@@ -18,11 +18,54 @@ 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) {
+    for (double elem : v) {
         ss << scientific_format(elem, width, prec);
     }
     return ss.str();
 }
+
+static void print_message(std::string msg)
+{
+    if (msg.size() > 0) {
+        std::cout << "message: " << msg << "\n\n";
+    }
+    else {
+        std::cout << "\n";
+    }
+}
+
+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";
+    std::cout << f << ": ";
+    if (expr) {
+        std::cout << "\x1B[32mOK\033[0m\n";
+        print_message(msg);
+    }
+    else {
+        std::cout << "\x1B[31mFAIL\033[0m\n";
+        print_message(msg);
+        std::cout << file << " " << line << ": Assertion \"" << expr_str
+                  << "\" Failed\n\n";
+        abort();
+    }
+}
+
+bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol)
+{
+    if (a.n_elem != b.n_elem) {
+        return false;
+    }
+
+    for (int i = 0; i < a.n_elem; i++) {
+        if (std::abs(a(i) - b(i)) >= tol) {
+            return false;
+        }
+    }
+    return true;
+}