Branch init use of already calulated values needed

nufi/save_results.h

@@ -8,6 +8,24 @@
 
 class NuFISolver;
 
+struct RhoESnapshot {
+  unsigned int it = 0;
+  double time = 0.0;
+  std::vector<double> x_eval;
+  std::vector<double> rho;
+  std::vector<double> E;
+  double int_E_sqr = 0;
+};
+
+RhoESnapshot
+compute_rho_E_snapshot(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 = Parameters::NV);
+
+void flush_rho_E_history(std::vector<RhoESnapshot> &history,
+                         const std::string &dir);
+
 struct DiagnosticsSnapshot {
   unsigned int Nx = 0;
   unsigned int Nv = 0;
@@ -27,7 +45,7 @@ compute_diagnostics(const NuFISolver &solver, unsigned int n,
 void save_f(const DiagnosticsSnapshot &snap, const std::string &filepath);
 void save_rho(const DiagnosticsSnapshot &snap, const std::string &filepath);
 void save_Efield(const DiagnosticsSnapshot &snap, const std::string &filepath);
-double compute_int_E_squared(const DiagnosticsSnapshot &snap);
+double compute_int_E_squared(const std::vector<double> &E, double Lx);
 
 void save_time_series(const std::vector<double> &t,
                       const std::vector<double> &values,

src/main.cc

@@ -28,13 +28,10 @@ template <int dim> void run() {
   std::cout << "Initializing dealii Poisson Solver\n";
   poisson.initialize();
 
-  std::vector<double> int_E_squared;
-  int_E_squared.reserve(Nt);
-  std::vector<double> int_E_squared_times;
-  int_E_squared_times.reserve(Nt);
+  std::vector<GridStructure<dim>> grid_versions;
+  std::vector<SolutionSnapshot<dim>> phi_history;
 
-  std::vector<GridStructure<1>> grid_versions;
-  std::vector<SolutionSnapshot<1>> phi_history;
+  std::vector<RhoESnapshot> rho_e_history;
 
   std::ofstream time_file(Parameters::PLOT_DIR + "simulation_time.dat");
 
@@ -91,6 +88,8 @@ template <int dim> void run() {
     }
     update_solution_history(phi_history, poisson,
                             grid_versions.back().grid_version);
+    rho_e_history.push_back(compute_rho_E_snapshot(
+        solver, it, grid_versions, phi_history, Parameters::PLOT_NX));
 
     error_file << it << " " << poisson.get_error_estimate() << "\n";
     error_file.flush();
@@ -109,15 +108,8 @@ template <int dim> void run() {
                               Parameters::PLOT_NX, Parameters::NV);
 
       save_f(snap, Parameters::PLOT_DIR + "f_" + std::to_string(it) + ".dat");
-      save_rho(snap,
-               Parameters::PLOT_DIR + "rho_" + std::to_string(it) + ".dat");
-      save_Efield(snap,
-                  Parameters::PLOT_DIR + "E_" + std::to_string(it) + ".dat");
 
-      int_E_squared.push_back(compute_int_E_squared(snap));
-      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");
+      flush_rho_E_history(rho_e_history, Parameters::PLOT_DIR);
 
       plot_time = timer.elapsed() - plot_start;
       std::cout << "Results saved in " << plot_start << "[s]" << "\n";

src/save_results.cc

@@ -37,6 +37,65 @@ void save_field_1d(const std::vector<double> &x,
     file << x[i] << " " << values[i] << "\n";
 }
 
+RhoESnapshot
+compute_rho_E_snapshot(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) {
+  RhoESnapshot snap;
+  snap.it = n;
+  snap.time = n * Parameters::DT;
+  snap.x_eval = make_x_eval(Nx_out);
+
+  // rho at the plot grid (the expensive part — full characteristic trace)
+  snap.rho = solver.eval_rho(n, snap.x_eval, grid_struct, phi_history, Nv);
+
+  // E is just -grad(phi), cheap — reuse the already-solved phi for this step
+  std::vector<double> x_copy = snap.x_eval;
+  auto grad_phi = eval(x_copy, grid_struct[phi_history[n].grid_version],
+                       phi_history[n].solution);
+  snap.E.resize(Nx_out);
+  for (unsigned int i = 0; i < Nx_out; ++i)
+    snap.E[i] = -grad_phi[i];
+
+  snap.int_E_sqr = compute_int_E_squared(snap.E, Parameters::LX);
+
+  return snap;
+}
+
+void flush_rho_E_history(std::vector<RhoESnapshot> &history,
+                         const std::string &dir) {
+  ensure_results_dir();
+
+  std::vector<double> t, int_E_sqr;
+  t.reserve(history.size());
+  int_E_sqr.reserve(history.size());
+
+  for (const auto &snap : history) {
+    save_field_1d(snap.x_eval, snap.rho,
+                  dir + "rho_" + std::to_string(snap.it) + ".dat");
+    save_field_1d(snap.x_eval, snap.E,
+                  dir + "E_" + std::to_string(snap.it) + ".dat");
+    t.push_back(snap.time);
+    int_E_sqr.push_back(snap.int_E_sqr);
+  }
+
+  // Append rather than overwrite, since flush is called repeatedly.
+  static bool wrote_header = false;
+  std::ofstream file(dir + "int_E_sqr.dat",
+                     wrote_header ? std::ios::app : std::ios::trunc);
+  if (!wrote_header) {
+    file << "# nufi time series\n";
+    file << "# columns: t value\n";
+    wrote_header = true;
+  }
+  file << std::setprecision(10);
+  for (size_t i = 0; i < t.size(); ++i)
+    file << t[i] << " " << int_E_sqr[i] << "\n";
+
+  history.clear();
+}
+
 DiagnosticsSnapshot
 compute_diagnostics(const NuFISolver &solver, unsigned int n,
                     std::vector<GridStructure<1>> &grid_struct,
@@ -102,10 +161,10 @@ void save_Efield(const DiagnosticsSnapshot &snap, const std::string &filepath) {
   save_field_1d(snap.x_eval, snap.E, filepath);
 }
 
-double compute_int_E_squared(const DiagnosticsSnapshot &snap) {
-  const double dx = Parameters::LX / snap.Nx;
+double compute_int_E_squared(const std::vector<double> &E, double Lx) {
+  const double dx = Lx / E.size();
   double integral = 0.0;
-  for (double e : snap.E)
+  for (double e : E)
     integral += e * e;
   return 0.5 * integral * dx;
 }