better save results method, to test

nufi/fields.h

@@ -17,13 +17,9 @@ inline std::vector<double> make_x_eval(size_t Nx) {
 
   double dx = (Parameters::X_DOMAIN_RIGHT - Parameters::X_DOMAIN_LEFT) / Nx;
 
-  for (unsigned int i = 0; i < Nx; ++i) {
+  for (unsigned int i = 0; i < Nx; ++i)
     x_eval_E.push_back(Parameters::X_DOMAIN_LEFT + (i + 0.5) * dx);
-  }
-  // std::vector<double> x_eval(Nx);
-  // const double dx = Parameters::LX / Nx;
-  // for (size_t i = 0; i < Nx; ++i)
-  //   x_eval[i] = Parameters::X_DOMAIN_LEFT + i * dx;
+
   return x_eval_E;
 }
 
@@ -75,13 +71,10 @@ inline std::vector<double> eval(std::vector<double> &X,
 
 inline double integral_space_vector(const GridStructure<1> &grid,
                                     const Vector<double> &solution,
-                                    double dx = Parameters::PLOT_DX,
                                     size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
-  double xmin = Parameters::X_DOMAIN_LEFT;
-  std::vector<double> x_eval(Nx);
-  for (size_t i = 0; i < Nx; ++i)
-    x_eval[i] = xmin + i * dx;
+  std::vector<double> x_eval = make_x_eval(Nx);
+  const double dx = Parameters::LX / Nx;
 
   std::vector<double> tmp = eval(x_eval, grid, solution);
   for (size_t i = 0; i < Nx; ++i)
@@ -91,13 +84,10 @@ inline double integral_space_vector(const GridStructure<1> &grid,
 
 inline double integral_space_vector_squared(const GridStructure<1> &grid,
                                             const Vector<double> &solution,
-                                            double dx = Parameters::PLOT_DX,
                                             size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
-  double xmin = Parameters::X_DOMAIN_LEFT;
-  std::vector<double> x_eval(Nx);
-  for (size_t i = 0; i < Nx; ++i)
-    x_eval[i] = xmin + i * dx;
+  std::vector<double> x_eval = make_x_eval(Nx);
+  const double dx = Parameters::LX / Nx;
 
   std::vector<double> tmp = eval(x_eval, grid, solution);
   for (size_t i = 0; i < Nx; ++i)

nufi/save_results.h

@@ -2,17 +2,25 @@
 #define SAVE_RESULTS_H
 
 #include "nufi/grids.h"
-#include "nufi/poisson_problem.h"
 #include <deal.II/lac/vector.h>
 #include <string>
 #include <vector>
 
 class NuFISolver;
 
-void save_f(const NuFISolver &solver, unsigned int n,
-            std::vector<GridStructure<1>> &grid_struct,
-            std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,
-            unsigned int Nv_out, const std::string &filename);
+void save_field_1d(const std::vector<double> &x,
+                   const std::vector<double> &values,
+                   const std::string &filepath);
+
+void save_time_series(const std::vector<double> &t,
+                      const std::vector<double> &values,
+                      const std::string &filepath);
+
+void save_f_binary(const NuFISolver &solver, unsigned int n,
+                   std::vector<GridStructure<1>> &grid_struct,
+                   std::vector<SolutionSnapshot<1>> &phi_history,
+                   unsigned int Nx_out, unsigned int Nv_out,
+                   const std::string &filepath);
 
 void save_rho(const NuFISolver &solver, unsigned int n,
               std::vector<GridStructure<1>> &grid_struct,
@@ -23,7 +31,4 @@ void save_Efield(unsigned int it, std::vector<GridStructure<1>> &grid_versions,
                  std::vector<SolutionSnapshot<1>> &phi_history,
                  unsigned int Nx_out = Parameters::PLOT_NX);
 
-void save_space_vector(const std::vector<double> &vals,
-                       const std::string &filename, size_t it);
-
 #endif

src/nufi_solver.cc

@@ -193,6 +193,8 @@ void NuFISolver::run() {
 
   std::vector<double> int_E_squared;
   int_E_squared.reserve(Nt);
+  std::vector<double> int_E_squared_times; // <-- add this
+  int_E_squared_times.reserve(Nt);
 
   std::vector<GridStructure<1>> grid_versions;
   std::vector<SolutionSnapshot<1>> phi_history;
@@ -301,10 +303,11 @@ void NuFISolver::run() {
       double plot_start = timer.elapsed();
 
       std::cout << "Saving results...   ";
-      save_f(*this, it, grid_versions, phi_history, Parameters::PLOT_NX,
-             Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
+      save_f_binary(*this, it, grid_versions, phi_history, Parameters::PLOT_NX,
+                    Parameters::NV,
+                    Parameters::PLOT_DIR + "f_" + std::to_string(it) + ".bin");
       save_rho(*this, it, grid_versions, phi_history, Parameters::PLOT_NX,
-               "results/rho_" + std::to_string(it) + ".dat");
+               Parameters::PLOT_DIR + "rho_" + std::to_string(it) + ".dat");
 
       save_Efield(it, grid_versions, phi_history);
 
@@ -312,7 +315,9 @@ void NuFISolver::run() {
                                  grid_versions[phi_history[it].grid_version],
                                  phi_history[it].solution);
       int_E_squared.push_back(int_val);
-      save_space_vector(int_E_squared, "electricint", it);
+      int_E_squared_times.push_back(it * Parameters::DT);
+      save_time_series(int_E_squared_times, int_E_squared,
+                       Parameters::PLOT_DIR + "int_E_sqr.dat");
 
       plot_time = timer.elapsed() - plot_start;
       std::cout << "Results saved in " << plot_start << "[s]" << "\n";

src/save_results.cc

@@ -7,69 +7,104 @@
 #include "nufi/poisson_problem.h"
 #include <cstddef>
 #include <deal.II/numerics/solution_transfer.h>
+#include <filesystem>
 #include <fstream>
 #include <stdexcept>
 #include <string>
 #include <vector>
 
-void save_f(const NuFISolver &solver, unsigned int n,
-            std::vector<GridStructure<1>> &grid_struct,
-            std::vector<SolutionSnapshot<1>> &phi_history, unsigned int Nx_out,
-            unsigned int Nv_out, const std::string &filename) {
-  std::ofstream file(filename);
+namespace {
+void ensure_results_dir() { std::filesystem::create_directories("results"); }
+} // namespace
 
-  double xmin = Parameters::X_DOMAIN_LEFT;
-  double xmax = Parameters::X_DOMAIN_RIGHT;
+void save_field_1d(const std::vector<double> &x,
+                   const std::vector<double> &values,
+                   const std::string &filepath) {
+  if (x.size() != values.size())
+    throw std::runtime_error(
+        "save_field_1d: size mismatch (x=" + std::to_string(x.size()) +
+        ", values=" + std::to_string(values.size()) + ") writing " + filepath);
 
-  double vmin = Parameters::V_DOMAIN_LEFT;
-  double vmax = Parameters::V_DOMAIN_RIGHT;
+  ensure_results_dir();
+  std::ofstream file(filepath);
+  if (!file)
+    throw std::runtime_error("save_field_1d: failed to open " + filepath);
 
-  double dv = (vmax - vmin) / Nv_out;
+  file << "# nufi 1d field\n";
+  file << "# n = " << x.size() << "\n";
+  file << "# columns: x value\n";
+  file << std::setprecision(10);
+  for (size_t i = 0; i < x.size(); ++i)
+    file << x[i] << " " << values[i] << "\n";
+}
 
-  file << Nx_out << " " << Nv_out << "\n";
-  file << xmin << " " << xmax << "\n";
-  file << vmin << " " << vmax << "\n";
+void save_time_series(const std::vector<double> &t,
+                      const std::vector<double> &values,
+                      const std::string &filepath) {
+  if (t.size() != values.size())
+    throw std::runtime_error(
+        "save_time_series: size mismatch (t=" + std::to_string(t.size()) +
+        ", values=" + std::to_string(values.size()) + ") writing " + filepath);
+
+  ensure_results_dir();
+  std::ofstream file(filepath);
+  if (!file)
+    throw std::runtime_error("save_time_series: failed to open " + filepath);
+
+  file << "# nufi time series\n";
+  file << "# columns: t value\n";
+  file << std::setprecision(10);
+  for (size_t i = 0; i < t.size(); ++i)
+    file << t[i] << " " << values[i] << "\n";
+}
+
+void save_f_binary(const NuFISolver &solver, unsigned int n,
+                   std::vector<GridStructure<1>> &grid_struct,
+                   std::vector<SolutionSnapshot<1>> &phi_history,
+                   unsigned int Nx_out, unsigned int Nv_out,
+                   const std::string &filepath) {
+  ensure_results_dir();
+  std::ofstream file(filepath, std::ios::binary);
+  if (!file)
+    throw std::runtime_error("save_f_binary: failed to open " + filepath);
+
+  const double vmin = Parameters::V_DOMAIN_LEFT;
+  const double vmax = Parameters::V_DOMAIN_RIGHT;
+  const double dv = (vmax - vmin) / Nv_out;
 
   std::vector<double> x_eval = make_x_eval(Nx_out);
+  std::vector<double> v_eval(Nv_out);
+  for (unsigned int j = 0; j < Nv_out; ++j)
+    v_eval[j] = vmin + (j + 0.5) * dv;
+
+  // simple fixed binary layout: magic, Nx, Nv, x[], v[], f[Nv*Nx]
+  const uint32_t magic = 0x4E554649; // "NUFI"
+  const uint64_t nx64 = Nx_out, nv64 = Nv_out;
+
+  file.write(reinterpret_cast<const char *>(&magic), sizeof(magic));
+  file.write(reinterpret_cast<const char *>(&nx64), sizeof(nx64));
+  file.write(reinterpret_cast<const char *>(&nv64), sizeof(nv64));
+  file.write(reinterpret_cast<const char *>(x_eval.data()),
+             x_eval.size() * sizeof(double));
+  file.write(reinterpret_cast<const char *>(v_eval.data()),
+             v_eval.size() * sizeof(double));
+
   std::vector<double> val(Nx_out);
-
   for (unsigned int j = 0; j < Nv_out; ++j) {
-    double v = vmin + (j + 0.5) * dv;
-    val = solver.eval_f(n, x_eval, v, grid_struct, phi_history);
-
-    for (unsigned int i = 0; i < Nx_out; ++i) {
-      file << val[i];
-
-      if (i < Nx_out - 1)
-        file << " ";
-    }
-
-    file << "\n";
+    val = solver.eval_f(n, x_eval, v_eval[j], grid_struct, phi_history);
+    file.write(reinterpret_cast<const char *>(val.data()),
+               val.size() * sizeof(double));
   }
-
-  file.close();
 }
 
 void save_rho(const NuFISolver &solver, unsigned int n,
               std::vector<GridStructure<1>> &grid_struct,
               std::vector<SolutionSnapshot<1>> &phi_history,
               unsigned int Nx_out, const std::string &filename) {
-  std::ofstream file(filename);
-
-  double xmin = Parameters::X_DOMAIN_LEFT;
-  double xmax = Parameters::X_DOMAIN_RIGHT;
-
   std::vector<double> x_eval = make_x_eval(Nx_out);
-  file << Nx_out << "\n";
-  file << xmin << " " << xmax << "\n";
-
-  std::vector<double> tmp =
+  std::vector<double> rho =
       solver.eval_rho(n, x_eval, grid_struct, phi_history);
-  for (size_t i = 0; i < Nx_out; ++i) {
-    file << tmp[i];
-    file << "\n";
-  }
-  file.close();
+  save_field_1d(x_eval, rho, filename);
 }
 
 void save_Efield(unsigned int it, std::vector<GridStructure<1>> &grid_versions,
@@ -80,37 +115,12 @@ void save_Efield(unsigned int it, std::vector<GridStructure<1>> &grid_versions,
   auto grad_phi = eval(x_eval_E, grid_versions[phi_history[it].grid_version],
                        phi_history[it].solution);
 
-  std::vector<std::pair<double, double>> ordered_E;
+  std::vector<double> E(x_eval_E.size());
+  for (size_t i = 0; i < x_eval_E.size(); ++i)
+    E[i] = -grad_phi[i];
 
-  for (size_t i = 0; i < x_eval_E.size(); ++i) {
-    double E = -grad_phi[i];
-    ordered_E.push_back({x_eval_E[i], E});
-  }
-
-  std::sort(ordered_E.begin(), ordered_E.end());
-
-  std::ofstream E_file("results/E_" + std::to_string(it) + ".dat");
-
-  for (auto [x, E] : ordered_E)
-    E_file << x << " " << E << "\n";
+  save_field_1d(x_eval_E, E, "results/E_" + std::to_string(it) + ".dat");
 
   std::cout << "Saving iteration " << it << " using grid version "
             << phi_history[it].grid_version << "\n\n";
-
-  E_file.close();
-}
-
-void save_space_vector(const std::vector<double> &vals,
-                       const std::string &filename, size_t it) {
-  std::ofstream file("results/" + filename + "_" + std::to_string(it) + ".dat");
-
-  if (!file)
-    throw std::runtime_error("failed to start file in results/");
-
-  file << vals.size() << "\n";
-  file << Parameters::X_DOMAIN_LEFT << " " << Parameters::X_DOMAIN_RIGHT
-       << "\n";
-  file << std::fixed << std::setprecision(8);
-  for (double val : vals)
-    file << val << "\n";
 }