.

CMakeLists.txt

@@ -16,6 +16,8 @@ endif()
 
 deal_ii_initialize_cached_variables()
 
+find_package(OpenMP REQUIRED)
+
 # -------------------------
 
 add_library(nufi_lib
@@ -30,6 +32,10 @@ target_include_directories(nufi_lib PUBLIC
 
 deal_ii_setup_target(nufi_lib)
 
+target_link_libraries(nufi_lib
+  OpenMP::OpenMP_CXX
+) 
+
 # -------------------------
 
 add_executable(nufi_poisson
@@ -38,6 +44,8 @@ add_executable(nufi_poisson
 
 target_link_libraries(nufi_poisson
   nufi_lib
+  OpenMP::OpenMP_CXX
 )
 
+
 deal_ii_setup_target(nufi_poisson)

nufi/fields.h

@@ -83,6 +83,7 @@ void interpolate( real *coeffs, const real *values)
 
         void operator()( const real *in, real *out ) const
         {
+            #pragma omp parallel for
             for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
             {
                 real result = 0;
@@ -145,24 +146,26 @@ void interpolate( real *coeffs, const real *values)
 inline double integral_space_vector(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
 {
   double integral = 0.0;
-  double x = Parameters::X_DOMAIN_LEFT;
+  double xmin = Parameters::X_DOMAIN_LEFT;
+  #pragma omp parallel for reduction (+:integral)
   for (size_t i=0; i<Nx ; ++i) {
-    x += dx;
+    double x = xmin + i * dx;
-    integral += eval<1>(x, current_coeffs)*dx;
+    integral += eval<1>(x, current_coeffs);
   }
-  return integral;
+  return integral*dx;
 };
 
 inline double integral_space_vector_squared(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
 {
   double integral = 0.0;
-  double x = Parameters::X_DOMAIN_LEFT;
+  double xmin = Parameters::X_DOMAIN_LEFT;
+  #pragma omp parallel for reduction (+:integral)
   for (size_t i=0; i<Nx ; ++i) {
-    x += dx;
+    double x = xmin + i*dx;
     double val = eval<1>(x, current_coeffs); 
-    integral += val*val*dx;
+    integral += val*val;
   }
-  return integral;
+  return integral*dx;
 };
 
 class Gradient {

nufi/nufi_solver.h

@@ -22,6 +22,7 @@ public:
   void run();
   double eval_rho(unsigned int n, double x, const double *E_coeffs, unsigned int Nv = Parameters::NV) const;
   double eval_ftilda(unsigned int n, double x, double u, const double *E_coeffs) const;
+  double eval_f(unsigned int n, double x, double u, const double *E_coeffs) const;
 
 private:
 

nufi/parameters.h

@@ -13,16 +13,16 @@ namespace Parameters
   constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT);
   constexpr double LX_INV = 1/LX;
 
-  constexpr double V_DOMAIN_LEFT = -10.0;
+  constexpr double V_DOMAIN_LEFT = -10.;
-  constexpr double V_DOMAIN_RIGHT = 10.0;
+  constexpr double V_DOMAIN_RIGHT = 10.;
 
   constexpr unsigned int NV = 512;
   constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT)/NV;
 
   // deal.ii options
-  constexpr unsigned int GLOBAL_REFINEMENT = 7;
+  constexpr unsigned int GLOBAL_REFINEMENT = 8;
   constexpr unsigned int FE_DEGREE = 4;
-  constexpr unsigned int CONVERGENCE_ITERATIONS = 20000;
+  constexpr unsigned int CONVERGENCE_ITERATIONS = 10000;
   constexpr double CONVERGENCE_LIMIT = 1e-12;
 
   constexpr double EPS = 0.01;
@@ -30,8 +30,8 @@ namespace Parameters
   constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
 
   // NUFI options
-  constexpr double DT=1./10.;
+  constexpr double DT=1./16.;
-  constexpr unsigned int TMAX = 50;
+  constexpr unsigned int TMAX = 500;
 
   //spline options
   constexpr int SPLINE_NX = 256;
@@ -40,7 +40,7 @@ namespace Parameters
   constexpr size_t SPLINE_ORDER = 4;
 
   //Plotting options
-  constexpr int PLOT_FREQUENCY = 5;
+  constexpr int PLOT_FREQUENCY = 10;
 }
 
 #endif

nufi/save_results.h

@@ -5,7 +5,7 @@
 #include "nufi/nufi_solver.h"
 
 
-void save_ftilda( const NuFISolver &solver,
+void save_f( const NuFISolver &solver,
     unsigned int n,
                                     const double *E_coeffs,
                                     unsigned int Nx_out,

src/main.cc

@@ -20,6 +20,8 @@ void clear_results_directory(const std::string &dir)
 
 int main()
 {
+  #include <omp.h>
+std::cout << "Threads: " << omp_get_max_threads() << "\n";
     try
     {
         clear_results_directory("results");

src/nufi_solver.cc

@@ -54,6 +54,42 @@ double NuFISolver::eval_ftilda(unsigned int n,
   return f0(x,u);
 }
 
+double NuFISolver::eval_f(unsigned int n,
+                                      double x,
+                                      double u,
+                                      const double *E_coeffs) const
+{
+  if ( n == 0 ) return f0(x,u);
+  
+  const size_t order = Parameters::SPLINE_ORDER;
+  const size_t stride_x = 1;
+  const size_t stride_t = stride_x*(Nx + order - 1);
+
+  double Ex;
+  const double *c;
+
+  // Initial half-step.
+  c  = E_coeffs + n*stride_t;
+  Ex = -eval<1>(x, c);
+  u += 0.5*Parameters::DT * Ex;
+
+  while ( --n )
+  {
+      x  = x - Parameters::DT *u;
+      c  = E_coeffs + n*stride_t;
+      Ex = -eval<1>(x, c);
+      u  = u + Parameters::DT *Ex;
+  }
+
+  // The final half-step.
+  x -= Parameters::DT*u;
+  c  = E_coeffs + n*stride_t;
+  Ex = -eval<1>(x, c);
+  u += 0.5*Parameters::DT*Ex;
+
+  return f0(x,u);
+}
+
 double NuFISolver::eval_rho(const unsigned int n,
                           const double x,
                           const double *E_coeffs,
@@ -63,10 +99,12 @@ double NuFISolver::eval_rho(const unsigned int n,
   const double v_min = Parameters::V_DOMAIN_LEFT;
 
   double integral = 0.0;
-  for (unsigned int i = 0; i < Nv; ++i)
-    integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs) * dv;
 
-  return 1.0 - integral;
+#pragma omp parallel for reduction (+ : integral)
+  for (unsigned int i = 0; i < Nv; ++i)
+    integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs);
+  
+  return 1.0 - integral*dv;
 }
 
 void NuFISolver::run()
@@ -94,16 +132,20 @@ void NuFISolver::run()
   {
       stopwatch<double> timer;
       
-      std::cout << "Timestep " << it << " / " << Nt << std::endl << std::endl;
+      double time_elapsed_before = timer.elapsed();
+
+      std::cout << "Timestep " << it << " / " << Nt << " (simulation time = "<< it*Parameters::DT << ")"<< std::endl;
 
       // compute rho
 
       double dx = Parameters::SPLINE_DX;
-      double x = Parameters::X_DOMAIN_LEFT;
       
-    	for(size_t i = 0; i<Nx; i++, x+=dx)
+      #pragma omp parallel for
+    	for(size_t i = 0; i<Nx; i++)
     	{
+        double x = Parameters::X_DOMAIN_LEFT + i*dx;
         double ith_rho = eval_rho(it, x, coeffs.get(),Parameters::NV); 
+
         AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
     		rho.get()[i] =  ith_rho;
     	}
@@ -126,7 +168,9 @@ void NuFISolver::run()
       double* current_coeffs = coeffs.get() + it*stride_t;
       interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, sampled_potential.data());
 
+
       std::vector<double> E_x(Nx,0.0) ;
+      #pragma omp parallel for
       for(size_t ix=0; ix<Nx; ++ix)
       {
         E_x[ix] = -eval<1>(Parameters::X_DOMAIN_LEFT+ix*dx, current_coeffs);
@@ -135,11 +179,12 @@ void NuFISolver::run()
       double timer_elapsed = timer.elapsed();
       total_time += timer_elapsed;
       
+      std::cout << "step made in "<< timer_elapsed-time_elapsed_before <<" seconds\n\n";
       if (it % Parameters::PLOT_FREQUENCY == 0)
       {
         std::cout << "Saving results...   ";
-        save_ftilda(*this, it, coeffs.get(), 128, 128, "results/ftilda_" + std::to_string(it) + ".dat");
+        save_f(*this, it, coeffs.get(), Parameters::SPLINE_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
-        save_rho(*this, it, coeffs.get(), 128, "results/rho_" + std::to_string(it) + ".dat");
+        save_rho(*this, it, coeffs.get(), Parameters::SPLINE_NX, "results/rho_" + std::to_string(it) + ".dat");
         // save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat");
         save_space_vector(E_x, "field", it);
 

src/save_results.cc

@@ -8,7 +8,7 @@
 #include "nufi/nufi_solver.h"
 
 
-void save_ftilda( const NuFISolver &solver,
+void save_f( const NuFISolver &solver,
     unsigned int n,
                                     const double *E_coeffs,
                                     unsigned int Nx_out,
@@ -38,7 +38,7 @@ void save_ftilda( const NuFISolver &solver,
       {
           double v = vmin + (j + 0.5)*dv;
 
-          double val = solver.eval_ftilda(n, x, v, E_coeffs);
+          double val = solver.eval_f(n, x, v, E_coeffs);
 
           file << val;