Merge pull request #14 from jvachier/jv/improvement

Improvement.

src/abp_3D_confine.cpp

@@ -13,8 +13,8 @@
 #include <stdio.h>
 #include <cmath>
 #include <time.h>
-#include <omp.h> //import library to use pragma
-#include <tuple> //to output multiple components of a function
+#include <omp.h>
+#include <tuple>
 
 #include "headers/print_file.h"
 #include "headers/cylindrical_reflective_boundary_conditions.h"
@@ -53,55 +53,50 @@ int main(int argc, char *argv[])
 	fscanf(parameter, "%lf\t%lf\t%d\t%lf\t%lf\t%lf\t%lf\t%lf\t%d\n", &epsilon, &delta, &Particles, &Dt, &De, &vs, &Wall, &height, &N);
 	printf("%lf\t%lf\t%d\t%lf\t%lf\t%lf\t%lf\t%lf\t%d\n", epsilon, delta, Particles, Dt, De, vs, Wall, height, N);
 
-    // Position
+	// Position
 	double *x = (double *)malloc(Particles * sizeof(double)); // x-position
 	double *y = (double *)malloc(Particles * sizeof(double)); // y-position
-    double *z = (double *)malloc(Particles * sizeof(double)); // z-position  
+	double *z = (double *)malloc(Particles * sizeof(double)); // z-position
 
-    // Orientation
-    double *ex = (double *)malloc(Particles * sizeof(double)); // ex-orientation
+	// Orientation
+	double *ex = (double *)malloc(Particles * sizeof(double)); // ex-orientation
 	double *ey = (double *)malloc(Particles * sizeof(double)); // ey-orientation
-    double *ez = (double *)malloc(Particles * sizeof(double)); // ez-orientation 
+	double *ez = (double *)malloc(Particles * sizeof(double)); // ez-orientation
 
 	// parameters
 	const int L = 1.0; // particle size
 
 	// initialization of the random generator
 	random_device rdev;
 	default_random_engine generator(rdev()); // random seed -> rdev
-	// default_random_engine generator(1); //same seed
 
 	// Distributions Gaussian
 	normal_distribution<double> Gaussdistribution(0.0, 1.0);
 	// Distribution Uniform for initialization
 	uniform_real_distribution<double> distribution(-Wall, Wall);
-    // Uniform distribution for the orientation - later on maybe take it from the unit sphere but normalized in update position
+	// Uniform distribution for the orientation
 	uniform_real_distribution<double> distribution_e(0.0, 1.0);
 
 	double xi_px = 0.0; // noise for x-position
-	double xi_py = 0.0; // noise for y-position 
-    double xi_pz = 0.0; // noise for z-position
-	double xi_ex = 0.0;  // noise ex ortientation
-    double xi_ey = 0.0;  // noise ey ortientation
-    double xi_ez = 0.0;  // noise ez ortientation
+	double xi_py = 0.0; // noise for y-position
+	double xi_pz = 0.0; // noise for z-position
+	double xi_ex = 0.0; // noise ex ortientation
+	double xi_ey = 0.0; // noise ey ortientation
+	double xi_ez = 0.0; // noise ez ortientation
 
 	// double phi = 0.0;
 	double prefactor_e = sqrt(2.0 * delta * De);
 	double prefactor_xi_px = sqrt(2.0 * delta * Dt);
 	double prefactor_xi_py = sqrt(2.0 * delta * Dt);
-    double prefactor_xi_pz = sqrt(2.0 * delta * Dt);
+	double prefactor_xi_pz = sqrt(2.0 * delta * Dt);
 	double prefactor_interaction = epsilon * 48.0;
 	double r = 5.0 * L;
 
-	/* does not work when using openmp
-	clock_t tStart = clock(); // check time for one trajectory 
-	*/
-
 	// Open MP to get execution time
 	double itime, ftime, exec_time;
-    itime = omp_get_wtime(); 
-    
-    fprintf(datacsv, "Particles,x-position,y-position,z-position,ex-orientation,ey-orientation,ez-orientation,time\n");
+	itime = omp_get_wtime();
+
+	fprintf(datacsv, "Particles,x-position,y-position,z-position,ex-orientation,ey-orientation,ez-orientation,time\n");
 
 	// initialization position and activity
 	initialization(
@@ -119,10 +114,10 @@ int main(int argc, char *argv[])
 		update_position(
 			x, y, z, ex, ey, ez, prefactor_e, Particles,
 			delta, De, Dt, xi_ex, xi_ey, xi_ez, xi_px,
-			xi_py, xi_pz, vs, prefactor_xi_px, prefactor_xi_py, prefactor_xi_pz, 
+			xi_py, xi_pz, vs, prefactor_xi_px, prefactor_xi_py, prefactor_xi_pz,
 			r, prefactor_interaction,
 			generator, Gaussdistribution, distribution_e);
-		
+
 		cylindrical_reflective_boundary_conditions(
 			x, y, z, Particles,
 			Wall, height, L);
@@ -136,14 +131,9 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	/* Does not work when using openmp
-	printf("Time taken: %.2fs\n", ((double)(clock() - tStart) / CLOCKS_PER_SEC)); 
-	// printf("Time taken: %.2fs\n", ((double)(clock() - tStart) / CLOCKS_PER_SEC/N_thread));
-	*/
-
 	ftime = omp_get_wtime();
-    exec_time = ftime - itime;
-    printf("Time taken is %f", exec_time);
+	exec_time = ftime - itime;
+	printf("Time taken is %f", exec_time);
 
 	free(x);
 	free(y);

src/check_nooverlap.cpp

@@ -9,7 +9,7 @@ void check_nooverlap(
 {
 	int count = 0;
 	double R = 0.0;
-#pragma omp parallel for simd 
+#pragma omp parallel for simd
 	for (int k = 0; k < Particles; k++)
 	{
 		for (int j = 0; j < Particles; j++)

src/cylindrical_reflective_boundary_conditions.cpp

@@ -6,10 +6,10 @@ void cylindrical_reflective_boundary_conditions(
 	double *x, double *y, double *z, int Particles,
 	double Wall, double height, int L)
 {
-	double distance_squared = 0.0, Wall_squared = Wall * Wall; 
-	double height_L = height - L/2.0;
+	double distance_squared = 0.0, Wall_squared = Wall * Wall;
+	double height_L = height - L / 2.0;
 	double D_AW_z = 0.0;
-#pragma omp parallel for simd 
+#pragma omp parallel for simd
 	for (int k = 0; k < Particles; k++)
 	{
 		// x-y coordidnate circle
@@ -20,7 +20,7 @@ void cylindrical_reflective_boundary_conditions(
 			y[k] = (sqrt(Wall_squared) / sqrt(distance_squared)) * y[k];
 		}
 
-		// z coordinate 
+		// z coordinate
 		D_AW_z = 0.0;
 		if (abs(z[k]) > height_L)
 		{
@@ -50,5 +50,4 @@ void cylindrical_reflective_boundary_conditions(
 			}
 		}
 	}
-
 }

src/initialization.cpp

@@ -16,7 +16,7 @@ void initialization(
 		ez[k] = distribution_e(generator);
 
 		// Need to normalize the orientaional vector
-		norm_e = sqrt(ex[k]*ex[k] + ey[k]*ey[k] + ez[k]*ez[k]);
+		norm_e = sqrt(ex[k] * ex[k] + ey[k] * ey[k] + ez[k] * ez[k]);
 		invers_norm_e = 1.0 / norm_e;
 
 		ex[k] = ex[k] * invers_norm_e;

src/update_position.cpp

@@ -14,20 +14,20 @@ void update_position(
 	double F = 0.0, R = 0.0;
 
 // First orientation
-#pragma omp parallel for simd 
+#pragma omp parallel for simd
 	for (int k = 0; k < Particles; k++)
 	{
 		xi_ex = distribution_e(generator);
 		xi_ey = distribution_e(generator);
 		xi_ez = distribution_e(generator);
 
-// Ito formulation
-		ex[k] = prefactor_e * ( ey[k]*xi_ez - xi_ez*ez[k])  - ex[k]; 
-		ey[k] = prefactor_e * ( ex[k]*xi_ez - xi_ex*ez[k])  - ey[k];
-		ez[k] = prefactor_e * ( ex[k]*xi_ey - xi_ex*ey[k])  - ez[k];
+		// Ito formulation
+		ex[k] = prefactor_e * (ey[k] * xi_ez - xi_ez * ez[k]) - ex[k];
+		ey[k] = prefactor_e * (ex[k] * xi_ez - xi_ex * ez[k]) - ey[k];
+		ez[k] = prefactor_e * (ex[k] * xi_ey - xi_ex * ey[k]) - ez[k];
 
 		// Need to normalize the orientaional vector
-		norm_e = sqrt(ex[k]*ex[k] + ey[k]*ey[k] + ez[k]*ez[k]);
+		norm_e = sqrt(ex[k] * ex[k] + ey[k] * ey[k] + ez[k] * ez[k]);
 		invers_norm_e = 1.0 / norm_e;
 
 		ex[k] = ex[k] * invers_norm_e;
@@ -36,7 +36,7 @@ void update_position(
 	}
 
 // Second position
-#pragma omp parallel for simd 
+#pragma omp parallel for simd
 	for (int k = 0; k < Particles; k++)
 	{
 
@@ -63,6 +63,6 @@ void update_position(
 		}
 		x[k] = x[k] + vs * ex[k] * delta + F * x[k] * delta + xi_px * prefactor_xi_px;
 		y[k] = y[k] + vs * ey[k] * delta + F * y[k] * delta + xi_py * prefactor_xi_py;
-		z[k] = z[k] + vs * ez[k] * delta + F * z[k] * delta + xi_pz * prefactor_xi_pz; 
+		z[k] = z[k] + vs * ez[k] * delta + F * z[k] * delta + xi_pz * prefactor_xi_pz;
 	}
 }