From a99eeab8b51fd7d1d957e10eac97ecb0f0494fb8 Mon Sep 17 00:00:00 2001
From: Cory <cory@balaton.dev>
Date: Sun, 22 Oct 2023 03:50:48 +0200
Subject: [PATCH] Optimize

- Use the new reinitialize method to reuse PenningTrap instances
---
 src/main.cpp | 187 +++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 138 insertions(+), 49 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index b6760cf..bdc967f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -107,7 +107,7 @@ void simulate_single_particle_with_different_steps()
         for (int i = 0; i < steps; i++) {
             ofile << arma::norm(res.r_vecs[0][i]
                                 - analytical_solution_particle_1(dt * i))
-                  << "\n";
+                  << '\n';
         }
         ofile.close();
     }
@@ -125,7 +125,7 @@ void simulate_single_particle_with_different_steps()
         for (int i = 0; i < steps; i++) {
             ofile << arma::norm(res.r_vecs[0][i]
                                 - analytical_solution_particle_1(dt * i))
-                  << "\n";
+                  << '\n';
         }
         ofile.close();
     }
@@ -150,7 +150,7 @@ void simulate_100_particles()
  * MHz.
  *
  * */
-void simulate_100_particles_with_time_potential_wide_sweep()
+void potential_resonance_wide_sweep()
 {
     double time = 500.;
 
@@ -160,7 +160,7 @@ void simulate_100_particles_with_time_potential_wide_sweep()
     double freq_end = 2.5;
     double freq_increment = .02;
     size_t freq_iterations
-        = (size_t)((freq_end - freq_start) / freq_increment);
+        = (size_t)((freq_end - freq_start) / freq_increment) + 1;
 
     double res[4][freq_iterations];
 
@@ -174,29 +174,27 @@ void simulate_100_particles_with_time_potential_wide_sweep()
         res[0][i] = freq_start + freq_increment * i;
     }
 
-// Using num_threads() is usually bad practice, but not having it
-// causes a SIGKILL.
-#pragma omp parallel for collapse(2) num_threads(4)
-    for (size_t i = 0; i < 3; i++) {
-        for (size_t j = 0; j < freq_iterations; j++) {
-            res[i + 1][j]
-                = PenningTrap((unsigned)100, T,
-                              std::bind(
-                                  [](double f, double r, double t) {
-                                      return (25. * V / 1000.)
-                                             * (1. + f * std::cos(r * t));
-                                  },
-                                  amplitudes[i], res[0][j],
-                                  std::placeholders::_1),
-                              500., 0.)
-                      .fraction_of_particles_left(time, N, "rk4", false);
+#pragma omp parallel
+    {
+        PenningTrap trap((unsigned int)100);
+#pragma omp for collapse(2)
+        for (size_t i = 0; i < 3; i++) {
+            for (size_t j = 0; j < freq_iterations; j++) {
+                trap.reinitialize(std::bind(
+                    [](double f, double r, double t) {
+                        return (25. * V / 1000.) * (1. + f * std::cos(r * t));
+                    },
+                    amplitudes[i], res[0][j], std::placeholders::_1));
+                res[i + 1][j]
+                    = trap.fraction_of_particles_left(time, N, "rk4", false);
+            }
         }
     }
 
-    ofile.open(path + "res.txt");
+    ofile.open(path + "wide_sweep.txt");
     for (size_t i = 0; i < freq_iterations; i++) {
-        ofile << res[0][i] << "," << res[1][i] << "," << res[2][i] << ","
-              << res[3][i] << "\n";
+        ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
+              << res[3][i] << '\n';
     }
     ofile.close();
 }
@@ -208,7 +206,7 @@ void simulate_100_particles_with_time_potential_wide_sweep()
  * MHz.
  *
  * */
-void simulate_100_particles_with_time_potential_narrow_sweep()
+void potential_resonance_no_interaction_narrow_sweep()
 {
     double time = 500.;
 
@@ -216,7 +214,7 @@ void simulate_100_particles_with_time_potential_narrow_sweep()
 
     double freq_start = 1.;
     double freq_end = 1.7;
-    double freq_increment = .001;
+    double freq_increment = .002;
     size_t freq_iterations
         = (size_t)((freq_end - freq_start) / freq_increment);
 
@@ -232,36 +230,90 @@ void simulate_100_particles_with_time_potential_narrow_sweep()
         res[0][i] = freq_start + freq_increment * i;
     }
 
-// Using num_threads() is usually bad practice, but not having it
-// causes a SIGKILL.
-#pragma omp parallel for collapse(2) num_threads(4)
-    for (size_t i = 0; i < 3; i++) {
-        for (size_t j = 0; j < freq_iterations; j++) {
-            res[i + 1][j]
-                = PenningTrap((unsigned)100, T,
-                              std::bind(
-                                  [](double f, double r, double t) {
-                                      return (25. * V / 1000.)
-                                             * (1. + f * std::cos(r * t));
-                                  },
-                                  amplitudes[i], res[0][j],
-                                  std::placeholders::_1),
-                              500., 0.)
-                      .fraction_of_particles_left(time, N, "rk4", false);
+#pragma omp parallel
+    {
+        PenningTrap trap((unsigned int)100);
+#pragma omp for collapse(2)
+        for (size_t i = 0; i < 3; i++) {
+            for (size_t j = 0; j < freq_iterations; j++) {
+                trap.reinitialize(std::bind(
+                    [](double f, double r, double t) {
+                        return (25. * V / 1000.) * (1. + f * std::cos(r * t));
+                    },
+                    amplitudes[i], res[0][j], std::placeholders::_1));
+                res[i + 1][j]
+                    = trap.fraction_of_particles_left(time, N, "rk4", false);
+            }
         }
     }
 
-    ofile.open(path + "res.txt");
+    ofile.open(path + "narrow_sweep.txt");
     for (size_t i = 0; i < freq_iterations; i++) {
-        ofile << res[0][i] << "," << res[1][i] << "," << res[2][i] << ","
-              << res[3][i] << "\n";
+        ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
+              << res[3][i] << '\n';
+    }
+    ofile.close();
+}
+
+/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
+ *  dependent potential.
+ *
+ *  @details The simulation sweeps over different frequencies in [1., 1.7]
+ * MHz.
+ *
+ * */
+void potential_resonance_with_interaction_narrow_sweep()
+{
+    double time = 500.;
+
+    double amplitudes[]{.1, .4, .7};
+
+    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);
+
+    double res[4][freq_iterations];
+
+    std::string path = "output/time_dependent_potential/";
+    mkpath(path);
+
+    std::ofstream ofile;
+
+#pragma omp parallel for
+    for (size_t i = 0; i < freq_iterations; i++) {
+        res[0][i] = freq_start + freq_increment * i;
+    }
+
+#pragma omp parallel
+    {
+        PenningTrap trap((unsigned int)100);
+#pragma omp for collapse(2)
+        for (size_t i = 0; i < 3; i++) {
+            for (size_t j = 0; j < freq_iterations; j++) {
+                trap.reinitialize(std::bind(
+                    [](double f, double r, double t) {
+                        return (25. * V / 1000.) * (1. + f * std::cos(r * t));
+                    },
+                    amplitudes[i], res[0][j], std::placeholders::_1));
+                res[i + 1][j] = trap.fraction_of_particles_left(time, N);
+            }
+        }
+    }
+
+    ofile.open(path + "narrow_sweep_interactions.txt");
+    for (size_t i = 0; i < freq_iterations; i++) {
+        ofile << res[0][i] << ',' << res[1][i] << ',' << res[2][i] << ','
+              << res[3][i] << '\n';
     }
     ofile.close();
 }
 
 int main()
 {
-    double t0 = omp_get_wtime();
+    double start,end,t1,t2;
+    start = omp_get_wtime();
 
     simulate_single_particle();
 
@@ -269,14 +321,51 @@ int main()
 
     simulate_single_particle_with_different_steps();
 
+    t2 = omp_get_wtime();
+
+    std::cout << "Time single and double            : " << (t2 - start)
+              << " seconds" << std::endl;
+
+    t1 = omp_get_wtime();
+
     simulate_100_particles();
 
-    simulate_100_particles_with_time_potential_wide_sweep();
-    simulate_100_particles_with_time_potential_narrow_sweep();
+    t2 = omp_get_wtime();
 
-    double end = omp_get_wtime();
+    std::cout << "Time 100 particles                : " << (t2 - t1)
+              << " seconds" << std::endl;
 
-    std::cout << "Time: " << (end - t0) << " seconds" << std::endl;
+    t1 = omp_get_wtime();
+
+    potential_resonance_wide_sweep();
+
+    t2 = omp_get_wtime();
+
+    std::cout << "Time wide sweep                   : " << (t2 - t1)
+              << " seconds" << std::endl;
+
+    t1 = omp_get_wtime();
+
+    potential_resonance_no_interaction_narrow_sweep();
+
+    t2 = omp_get_wtime();
+
+    std::cout << "Time narrow sweep no interaction  : " << (t2 - t1)
+              << " seconds" << std::endl;
+
+    t1 = omp_get_wtime();
+
+    potential_resonance_with_interaction_narrow_sweep();
+
+    t2 = omp_get_wtime();
+
+    std::cout << "Time narrow sweep with interaction: " << (t2 - t1)
+              << " seconds" << std::endl;
+    
+    end = omp_get_wtime();
+
+    std::cout << "Time                              : " << (end - start)
+              << " seconds" << std::endl;
 
     return 0;
 }