From 908c90053bc585245d4ee72229f00c2e3846a9b5 Mon Sep 17 00:00:00 2001
From: "Vasco C. B. Ferreira" <vascocbf@gmail.com>
Date: Wed, 29 Jul 2026 12:28:46 +0200
Subject: [PATCH] added paralelization eval_rho and to diagnostics snapshot

---
 nufi/grids.h        |  2 +-
 src/main.cc         |  1 +
 src/nufi_solver.cc  | 21 +++++++++++++--------
 src/save_results.cc |  2 +-
 4 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/nufi/grids.h b/nufi/grids.h
index d65fc90..2d90bdf 100644
--- a/nufi/grids.h
+++ b/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;
     }
   }
diff --git a/src/main.cc b/src/main.cc
index e5f3c67..11bf2a6 100644
--- a/src/main.cc
+++ b/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");
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index 418c061..df25d41 100644
--- a/src/nufi_solver.cc
+++ b/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()
diff --git a/src/save_results.cc b/src/save_results.cc
index cee8cf3..2bc3850 100644
--- a/src/save_results.cc
+++ b/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);