wip: Branch init

Main changes done,
eval_rho + eval_ftilda_batch need to be checked to work correcly together

nufi/nufi_solver.h

@@ -23,10 +23,13 @@ public:
            const std::vector<GridStructure<1>> &grid_struct,
            const std::vector<SolutionSnapshot<1>> &phi_history,
            const unsigned int Nv = Parameters::NV) const;
+
   std::vector<double>
-  eval_ftilda(unsigned int n, std::vector<double> x, double u,
+  eval_ftilda_batch(unsigned int n, std::vector<double> X,
+                    std::vector<double> U,
                     const std::vector<GridStructure<1>> &grid_struct,
                     const std::vector<SolutionSnapshot<1>> &phi_history) const;
+
   std::vector<double>
   eval_f(unsigned int n, std::vector<double> x, double u,
          const std::vector<GridStructure<1>> &grid_struct,

src/main.cc

@@ -1,3 +1,4 @@
+#include "nufi/stopwatch.h"
 #include <filesystem>
 #include <iostream>
 #include <nufi/nufi_solver.h>
@@ -5,18 +6,6 @@
 #include <nufi/save_results.h>
 #include <omp.h>
 
-void clear_results_directory(const std::string &dir) {
-  if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir))
-    return;
-
-  for (const auto &entry : std::filesystem::directory_iterator(dir)) {
-    if (std::filesystem::is_regular_file(entry)) {
-      std::filesystem::remove(entry.path());
-      std::cout << "Deleted: " << entry.path() << '\n';
-    }
-  }
-};
-
 template <int dim> void run() {
   unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT);
 
@@ -134,6 +123,18 @@ template <int dim> void run() {
   std::cout << "NuFI simulation finished in " << total_time << " seconds.\n";
 };
 
+void clear_results_directory(const std::string &dir) {
+  if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir))
+    return;
+
+  for (const auto &entry : std::filesystem::directory_iterator(dir)) {
+    if (std::filesystem::is_regular_file(entry)) {
+      std::filesystem::remove(entry.path());
+      std::cout << "Deleted: " << entry.path() << '\n';
+    }
+  }
+};
+
 int main() {
   // omp_set_max_active_levels(1);
   std::cout << "Threads: " << omp_get_max_threads() << "\n";

src/nufi_solver.cc

@@ -20,74 +20,10 @@
 #include "nufi/fields.h"
 #include "nufi/grids.h"
 #include "nufi/parameters.h"
-#include "nufi/poisson_problem.h"
 #include "nufi/save_results.h"
-#include "nufi/stopwatch.h"
 
 using namespace dealii;
 
-std::vector<double> NuFISolver::eval_ftilda(
-    unsigned int n, std::vector<double> X, double u,
-    const std::vector<GridStructure<1>> &grid_struct,
-    const std::vector<SolutionSnapshot<1>> &phi_history) const {
-
-  size_t x_size = X.size();
-
-  std::vector<double> U(x_size, u);
-  std::vector<double> results(x_size);
-  if (n == 0) {
-    for (size_t i = 0; i < x_size; ++i)
-      results[i] = f0(X[i], U[i]);
-    return results;
-  }
-
-  std::vector<double> Ex(x_size);
-  std::vector<double> tmp(x_size);
-
-  // We omit the initial half-step.
-  while (--n) {
-    for (size_t i = 0; i < x_size; ++i)
-      X[i] = X[i] - Parameters::DT * U[i];
-
-    AssertThrow(
-        phi_history[n].solution.size() ==
-            grid_struct[phi_history[n].grid_version].dof_handler->n_dofs(),
-        ExcMessage("In eval_ftilda: Solution size = " +
-                   std::to_string(phi_history[n].solution.size()) +
-                   ", expected by grid_struct = " +
-                   std::to_string(grid_struct[phi_history[n].grid_version]
-                                      .dof_handler->n_dofs())));
-    AssertThrow(
-        grid_struct[phi_history[n].grid_version].grid_version ==
-            phi_history[n].grid_version,
-        ExcMessage(
-            "grid.grid_version not equal to phi_history[n].grid_version"));
-
-    tmp = eval(X, grid_struct[phi_history[n].grid_version],
-               phi_history[n].solution); // call eval only once
-
-    for (size_t i = 0; i < x_size; ++i) {
-      Ex[i] = -tmp[i];
-      U[i] = U[i] + Parameters::DT * Ex[i];
-    }
-  }
-
-  // The final half-step.
-  for (size_t i = 0; i < x_size; ++i)
-    X[i] = X[i] - Parameters::DT * U[i];
-
-  tmp = eval(X, grid_struct[phi_history[n].grid_version],
-             phi_history[n].solution); // call eval only once
-
-  for (size_t i = 0; i < x_size; ++i) {
-    Ex[i] = -tmp[i];
-    U[i] = U[i] + 0.5 * Parameters::DT * Ex[i];
-  }
-  for (size_t i = 0; i < x_size; ++i)
-    results[i] = f0(X[i], U[i]);
-  return results;
-}
-
 std::vector<double>
 NuFISolver::eval_f(unsigned int n, std::vector<double> X, double u,
                    const std::vector<GridStructure<1>> &grid_struct,
@@ -142,32 +78,98 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> X, double u,
   return results;
 }
 
+std::vector<double> NuFISolver::eval_ftilda_batch(
+    unsigned int n, std::vector<double> X, std::vector<double> U,
+    const std::vector<GridStructure<1>> &grid_struct,
+    const std::vector<SolutionSnapshot<1>> &phi_history) const {
+
+  const size_t x_size = X.size();
+  std::vector<double> results(x_size);
+
+  if (n == 0) {
+    for (size_t k = 0; k < x_size; ++k)
+      results[k] = f0(X[k], U[k]);
+    return results;
+  }
+
+  std::vector<double> Ex(x_size);
+  std::vector<double> tmp(x_size);
+
+  // We omit the initial half-step.
+  while (--n) {
+    for (size_t k = 0; k < x_size; ++k)
+      X[k] = X[k] - Parameters::DT * U[k];
+
+    AssertThrow(
+        phi_history[n].solution.size() ==
+            grid_struct[phi_history[n].grid_version].dof_handler->n_dofs(),
+        ExcMessage("In eval_ftilda_batch: Solution size = " +
+                   std::to_string(phi_history[n].solution.size()) +
+                   ", expected by grid_struct = " +
+                   std::to_string(grid_struct[phi_history[n].grid_version]
+                                      .dof_handler->n_dofs())));
+    AssertThrow(
+        grid_struct[phi_history[n].grid_version].grid_version ==
+            phi_history[n].grid_version,
+        ExcMessage(
+            "grid.grid_version not equal to phi_history[n].grid_version"));
+
+    tmp = eval(X, grid_struct[phi_history[n].grid_version],
+               phi_history[n].solution); // one call for the whole batch
+
+    for (size_t k = 0; k < x_size; ++k) {
+      Ex[k] = -tmp[k];
+      U[k] = U[k] + Parameters::DT * Ex[k];
+    }
+  }
+
+  // The final half-step.
+  for (size_t k = 0; k < x_size; ++k)
+    X[k] = X[k] - Parameters::DT * U[k];
+
+  tmp = eval(X, grid_struct[phi_history[n].grid_version],
+             phi_history[n].solution);
+
+  for (size_t k = 0; k < x_size; ++k) {
+    Ex[k] = -tmp[k];
+    U[k] = U[k] + 0.5 * Parameters::DT * Ex[k];
+  }
+  for (size_t k = 0; k < x_size; ++k)
+    results[k] = f0(X[k], U[k]);
+  return results;
+}
+
 std::vector<double>
 NuFISolver::eval_rho(unsigned int n, std::vector<double> &X,
                      const std::vector<GridStructure<1>> &grid_struct,
                      const std::vector<SolutionSnapshot<1>> &phi_history,
                      const unsigned int Nv) const {
-  size_t x_size = X.size();
+  const size_t x_size = X.size();
+  const size_t total = static_cast<size_t>(Nv) * x_size;
+
+  std::vector<double> X_all(total);
+  std::vector<double> U_all(total);
 
   const double dv =
       (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
-  const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv;
+  const double v_min = Parameters::V_DOMAIN_LEFT;
-
-  std::vector<double> partial(static_cast<size_t>(Nv) * x_size);
 
 #pragma omp parallel for
   for (unsigned int i = 0; i < Nv; ++i) {
-    std::vector<double> tmp_int =
+    const double v_i = v_min + i * dv;
-        eval_ftilda(n, X, v_min + i * dv, grid_struct,
+    std::copy(X.begin(), X.end(),
-                    phi_history); // used eval_ftilda once per i
+              X_all.begin() + static_cast<size_t>(i) * x_size);
-    std::copy(tmp_int.begin(), tmp_int.end(),
+    std::fill(U_all.begin() + static_cast<size_t>(i) * x_size,
-              partial.begin() + static_cast<size_t>(i) * x_size);
+              U_all.begin() + static_cast<size_t>(i + 1) * x_size, v_i);
   }
 
+  std::vector<double> ftilda_all = eval_ftilda_batch(
+      n, std::move(X_all), std::move(U_all), grid_struct, phi_history);
+
   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];
+      integral[ii] += ftilda_all[static_cast<size_t>(i) * x_size + ii];
 
   for (size_t i = 0; i < x_size; ++i)
     integral[i] = 1 - integral[i] * dv;