From ed6b4bf373f40ceee84624e0803a7d46e0615160 Mon Sep 17 00:00:00 2001 From: VCB Ferreira Date: Tue, 4 Aug 2026 19:05:11 +0200 Subject: [PATCH] Branch init use of already calulated values needed --- nufi/save_results.h | 20 +++++++++++++- src/main.cc | 20 +++++--------- src/save_results.cc | 65 ++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 87 insertions(+), 18 deletions(-) diff --git a/nufi/save_results.h b/nufi/save_results.h index 40bbb49..ef94fb9 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -8,6 +8,24 @@ class NuFISolver; +struct RhoESnapshot { + unsigned int it = 0; + double time = 0.0; + std::vector x_eval; + std::vector rho; + std::vector E; + double int_E_sqr = 0; +}; + +RhoESnapshot +compute_rho_E_snapshot(const NuFISolver &solver, unsigned int n, + std::vector> &grid_struct, + std::vector> &phi_history, + unsigned int Nx_out, unsigned int Nv = Parameters::NV); + +void flush_rho_E_history(std::vector &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 &E, double Lx); void save_time_series(const std::vector &t, const std::vector &values, diff --git a/src/main.cc b/src/main.cc index 8b60171..670a743 100644 --- a/src/main.cc +++ b/src/main.cc @@ -28,13 +28,10 @@ template void run() { std::cout << "Initializing dealii Poisson Solver\n"; poisson.initialize(); - std::vector int_E_squared; - int_E_squared.reserve(Nt); - std::vector int_E_squared_times; - int_E_squared_times.reserve(Nt); + std::vector> grid_versions; + std::vector> phi_history; - std::vector> grid_versions; - std::vector> phi_history; + std::vector rho_e_history; std::ofstream time_file(Parameters::PLOT_DIR + "simulation_time.dat"); @@ -91,6 +88,8 @@ template 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 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"; diff --git a/src/save_results.cc b/src/save_results.cc index 2bc3850..a55aaa1 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -37,6 +37,65 @@ void save_field_1d(const std::vector &x, file << x[i] << " " << values[i] << "\n"; } +RhoESnapshot +compute_rho_E_snapshot(const NuFISolver &solver, unsigned int n, + std::vector> &grid_struct, + std::vector> &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 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 &history, + const std::string &dir) { + ensure_results_dir(); + + std::vector 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> &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 &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; }