added paralelization eval_rho and to diagnostics snapshot

nufi/grids.h

@@ -127,7 +127,7 @@ GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
 
       if (PRINT_GAUGE_DOF_POSITION)
         std::cout << " gauge_dof = " << gauge_dof
-                  << " gauge_point = " << point[0] << std::endl;
+                  << " gauge_point = " << point[0] << "\n";
       break;
     }
   }

src/main.cc

@@ -16,6 +16,7 @@ void clear_results_directory(const std::string &dir) {
 }
 
 int main() {
+  omp_set_max_active_levels(1);
   std::cout << "Threads: " << omp_get_max_threads() << "\n";
   try {
     clear_results_directory("results");

src/nufi_solver.cc

@@ -155,19 +155,24 @@ NuFISolver::eval_rho(unsigned int n, std::vector<double> &X,
 
   const double dv =
       (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
-
   const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv;
 
-  std::vector<double> integral(x_size, 0.0);
-  std::vector<double> tmp_int(x_size);
+  std::vector<double> partial(static_cast<size_t>(Nv) * x_size);
 
+#pragma omp parallel for
   for (unsigned int i = 0; i < Nv; ++i) {
-    tmp_int = eval_ftilda(n, X, v_min + i * dv, grid_struct,
-                          phi_history); // used eval_ftilda once per i
-    for (size_t ii = 0; ii < x_size; ++ii)
-      integral[ii] += tmp_int[ii];
+    std::vector<double> tmp_int =
+        eval_ftilda(n, X, v_min + i * dv, grid_struct,
+                    phi_history); // used eval_ftilda once per i
+    std::copy(tmp_int.begin(), tmp_int.end(),
+              partial.begin() + static_cast<size_t>(i) * x_size);
   }
 
+  std::vector<double> integral(x_size, 0.0);
+  for (unsigned int i = 0; i < Nv; ++i)
+    for (size_t ii = 0; ii < x_size; ++ii)
+      integral[ii] += partial[static_cast<size_t>(i) * x_size + ii];
+
   for (size_t i = 0; i < x_size; ++i)
     integral[i] = 1 - integral[i] * dv;
   return integral;
@@ -235,7 +240,7 @@ void NuFISolver::run() {
 
     std::cout << "Timestep " << it << " / " << Nt
               << " (simulation time = " << it * Parameters::DT << ")"
-              << std::endl;
+              << "\n";
 
     // START: diagnostics
     std::cout << "cells = " << poisson.get_triangulation().n_active_cells()

src/save_results.cc

@@ -4,7 +4,6 @@
 #include "nufi/grids.h"
 #include "nufi/nufi_solver.h"
 #include "nufi/parameters.h"
-#include "nufi/poisson_problem.h"
 #include <cstddef>
 #include <deal.II/numerics/solution_transfer.h>
 #include <filesystem>
@@ -57,6 +56,7 @@ compute_diagnostics(const NuFISolver &solver, unsigned int n,
     snap.v_eval[j] = vmin + (j + 0.5) * dv;
 
   snap.f.resize(static_cast<size_t>(Nx_out) * Nv_out);
+#pragma omp parallel for
   for (unsigned int j = 0; j < Nv_out; ++j) {
     std::vector<double> val =
         solver.eval_f(n, snap.x_eval, snap.v_eval[j], grid_struct, phi_history);