update save_results part 1

nufi/fields.h

@@ -59,6 +59,18 @@ using namespace dealii;
 //     }
 //     return indices;
 // }
+inline std::vector<double> make_x_eval(size_t Nx) {
+  std::vector<double> x_eval(Nx);
+  for (size_t i = 0; i < Nx; ++i)
+    x_eval[i] = Parameters::X_DOMAIN_LEFT + i * Parameters::CALC_DX;
+  return x_eval;
+}
+
+inline void reset_x_eval(std::vector<double> &x_vals) {
+  const size_t Nx = x_vals.size();
+  for (size_t i = 0; i < Nx; ++i)
+    x_vals[i] = Parameters::X_DOMAIN_LEFT + i * Parameters::CALC_DX;
+};
 
 inline double f0(const double x, const double v,
                  const double eps = Parameters::EPS,
@@ -95,11 +107,13 @@ inline double integral_space_vector(const PoissonProblem<1> &poisson,
                                     size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
   double xmin = Parameters::X_DOMAIN_LEFT;
-#pragma omp parallel for reduction(+ : integral)
-  for (size_t i = 0; i < Nx; ++i) {
-    double x = xmin + i * dx;
-    integral += eval(x, poisson, solution);
-  }
+  std::vector<double> x_eval(Nx);
+  for (size_t i = 0; i < Nx; ++i)
+    x_eval[i] = xmin + i * dx;
+
+  std::vector<double> tmp = eval(x_eval, poisson, solution);
+  for (size_t i = 0; i < Nx; ++i)
+    integral += tmp[i];
   return integral * dx;
 };
 
@@ -109,12 +123,13 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
                                             size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
   double xmin = Parameters::X_DOMAIN_LEFT;
-#pragma omp parallel for reduction(+ : integral)
-  for (size_t i = 0; i < Nx; ++i) {
-    double x = xmin + i * dx;
-    double val = eval(x, poisson, solution);
-    integral += val * val;
-  }
+  std::vector<double> x_eval(Nx);
+  for (size_t i = 0; i < Nx; ++i)
+    x_eval[i] = xmin + i * dx;
+
+  std::vector<double> tmp = eval(x_eval, poisson, solution);
+  for (size_t i = 0; i < Nx; ++i)
+    integral += tmp[i] * tmp[i];
   return integral * dx;
 };
 

src/nufi_solver.cc

@@ -134,14 +134,14 @@ NuFISolver::eval_rho(unsigned int n, std::vector<double> &X,
   std::vector<double> tmp_int(x_size);
 
   for (unsigned int i = 0; i < Nv; ++i) {
-    tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson,
-                          phi_history); // used eval_ftilda once per i
+    tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson, phi_history);
     for (size_t ii = 0; ii < x_size; ++ii)
       integral[ii] += tmp_int[ii];
   }
 
   for (size_t i = 0; i < x_size; ++i)
     integral[i] = 1 - integral[i] * dv;
+  reset_x_eval(X);
   return integral;
 }
 
@@ -155,13 +155,16 @@ void NuFISolver::run() {
       reinterpret_cast<double *>(std::aligned_alloc(64, sizeof(double) * Nx)),
       std::free};
 
+  if (rho == nullptr)
+    throw std::bad_alloc{};
+
   std::vector<double> int_E_squared;
   int_E_squared.reserve(Nt);
 
   std::vector<Vector<double>> phi_history;
 
-  if (rho == nullptr)
-    throw std::bad_alloc{};
+  size_t CALC_NX = Parameters::CALC_NX;
+  std::vector<double> x_eval = make_x_eval(CALC_NX);
 
   double total_time = 0;
   std::ofstream time_file("results/simulation_time.txt");
@@ -180,14 +183,12 @@ void NuFISolver::run() {
               << std::endl;
 
     // compute rho
+    std::vector<double> tmp_rho =
+        eval_rho(it, x_eval, poisson, phi_history, Parameters::NV);
 
-#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, poisson, phi_history, Parameters::NV);
-
-      AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
-      rho.get()[i] = ith_rho;
+    for (size_t i = 0; i < CALC_NX; ++i) {
+      AssertThrow(std::isfinite(tmp_rho[i]), ExcMessage("NaN detected in rho"));
+      rho.get()[i] = tmp_rho[i];
     }
 
     poisson.set_rhs_function([&rho, x_min, dx, Nx = Nx](const Point<1> &p) {
@@ -222,12 +223,12 @@ void NuFISolver::run() {
       // save_Efield(it, coeffs.get(), 128, "results/field_" +
       // std::to_string(it) + ".dat");
 
-      std::vector<double> E_x(Nx, 0.0);
-#pragma omp parallel for
-      for (size_t ix = 0; ix < Nx; ++ix) {
-        E_x[ix] = -eval(x_min + ix * dx, poisson, phi_history[it]);
-      }
-
+      std::vector<double> E_x(CALC_NX);
+      std::vector<double> x_eval_Ex = make_x_eval(Parameters::PLOT_NX);
+      std::vector<double> tmp_Ex = eval(x_eval_Ex, poisson, phi_history[it]);
+      reset_x_eval(x_eval_Ex);
+      for (size_t i = 0; i < CALC_NX; ++i)
+        E_x[i] = -tmp_Ex[i];
       save_space_vector(E_x, "field", it);
 
       double int_val =

src/save_results.cc

@@ -4,6 +4,7 @@
 #include "nufi/nufi_solver.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
+#include <cstddef>
 #include <fstream>
 #include <stdexcept>
 #include <string>
@@ -21,22 +22,20 @@ void save_f(const NuFISolver &solver, unsigned int n,
   double vmin = Parameters::V_DOMAIN_LEFT;
   double vmax = Parameters::V_DOMAIN_RIGHT;
 
-  double dx = (xmax - xmin) / Nx_out;
   double dv = (vmax - vmin) / Nv_out;
 
   file << Nx_out << " " << Nv_out << "\n";
   file << xmin << " " << xmax << "\n";
   file << vmin << " " << vmax << "\n";
 
-  for (unsigned int i = 0; i < Nx_out; ++i) {
-    double x = xmin + (i + 0.5) * dx;
+  std::vector<double> x_eval = make_x_eval(Nx_out);
+  std::vector<double> val(Nx_out);
 
   for (unsigned int j = 0; j < Nv_out; ++j) {
     double v = vmin + (j + 0.5) * dv;
-
-      double val = solver.eval_f(n, x, v, poisson, phi_history);
-
-      file << val;
+    val = solver.eval_f(n, x_eval, v, poisson, phi_history);
+    for (unsigned int i = 0; i < Nx_out; ++i) {
+      file << val[i];
 
       if (j < Nv_out - 1)
         file << " ";
@@ -56,14 +55,14 @@ void save_rho(const NuFISolver &solver, unsigned int n,
 
   double xmin = Parameters::X_DOMAIN_LEFT;
   double xmax = Parameters::X_DOMAIN_RIGHT;
-  double dx = (xmax - xmin) / Nx_out;
 
+  std::vector<double> x_eval = make_x_eval(Nx_out);
   file << Nx_out << "\n";
   file << xmin << " " << xmax << "\n";
 
-  for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) {
-    double val = solver.eval_rho(n, xmin, poisson, phi_history);
-    file << val;
+  std::vector<double> tmp = solver.eval_rho(n, x_eval, poisson, phi_history);
+  for (size_t i = 0; i < Nx_out; ++i) {
+    file << tmp[i];
     file << "\n";
   }
   file.close();
@@ -77,16 +76,17 @@ void save_Efield([[maybe_unused]] unsigned int n,
 
   double xmin = Parameters::X_DOMAIN_LEFT;
   double xmax = Parameters::X_DOMAIN_RIGHT;
-  double dx = (xmax - xmin) / Nx_out;
+
+  std::vector<double> x_eval = make_x_eval(Nx_out);
 
   // select from E_coeffs
 
   file << Nx_out << "\n";
   file << xmin << " " << xmax << "\n";
 
-  for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) {
-    double val = -eval(xmin, poisson, phi_history[n]);
-    file << val;
+  std::vector<double> tmp = eval(x_eval, poisson, phi_history[n]);
+  for (size_t i = 0; i < Nx_out; ++i) {
+    file << -tmp[i];
     file << "\n";
   }
   file.close();