Small fix

src/PenningTrap.cpp

@@ -178,15 +178,15 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
             p->v_vec = original_particles[j].v_vec + this->v_func(i, j, dt);
             p->r_vec = original_particles[j].r_vec + this->r_func(i, j, dt);
         }
-        this->particles = tmp_particles;
+        this->particles.swap(tmp_particles);
     }
     this->t += dt;
 }
 
 void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
 {
-    std::vector<Particle> new_state = this->particles;
-
+    size_t size = this->particles.size();
+    vec_3d force_res[size];
     Particle *p;
 
     vec_3d (PenningTrap::*force)(unsigned int);
@@ -197,14 +197,18 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
         force = &PenningTrap::total_force_external;
     }
 
-#pragma omp parallel for private(p)
-    for (size_t i = 0; i < this->particles.size(); i++) {
-        p = &new_state[i];
-        p->v_vec += dt * (this->*force)(i) / p->m;
-        p->r_vec += dt * this->particles[i].v_vec;
+#pragma omp parallel for
+    for (size_t i = 0; i < size; i++) {
+        force_res[i] = (this->*force)(i);
+    }
+
+#pragma omp parallel for private(p)
+    for (size_t i = 0; i < size; i++) {
+        p = &this->particles[i];
+        p->r_vec += dt * p->v_vec;
+        p->v_vec += dt * force_res[i] / p->m;
     }
 
-    this->particles = new_state;
     this->t += dt;
 }
 

src/main.cpp

@@ -85,7 +85,7 @@ void simulate_100_particles()
 
     double time = 500.; // microseconds
 
-    trap.write_simulation_to_dir("output/simulate_100_particles", time, N*5);
+    trap.write_simulation_to_dir("output/simulate_100_particles", time, N*4);
 }
 
 void simulate_100_particles_with_time_potential()
@@ -132,13 +132,15 @@ int main()
 
     double start = omp_get_wtime();
 
-    //simulate_100_particles();
+    for (int i=0; i<5; i++) {
+        simulate_100_particles();
+    }
 
-    simulate_100_particles_with_time_potential();
+    //simulate_100_particles_with_time_potential();
 
     double end = omp_get_wtime();
 
-    std::cout << "Time: " << end - start << " seconds" << std::endl;
+    std::cout << "Average time: " << (end - start) / 5. << " seconds" << std::endl;
 
     return 0;
 }