diff --git a/libnufi_lib.a b/libnufi_lib.a index e57c4b4..5f3f7a0 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/save_results.h b/nufi/save_results.h index 03d1cb8..40bbb49 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -8,27 +8,29 @@ class NuFISolver; -void save_field_1d(const std::vector &x, - const std::vector &values, - const std::string &filepath); +struct DiagnosticsSnapshot { + unsigned int Nx = 0; + unsigned int Nv = 0; + std::vector x_eval; + std::vector v_eval; + std::vector f; + std::vector rho; + std::vector E; +}; + +DiagnosticsSnapshot +compute_diagnostics(const NuFISolver &solver, unsigned int n, + std::vector> &grid_struct, + std::vector> &phi_history, + unsigned int Nx_out, unsigned int Nv_out); + +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); void save_time_series(const std::vector &t, const std::vector &values, const std::string &filepath); -void save_f_binary(const NuFISolver &solver, unsigned int n, - std::vector> &grid_struct, - std::vector> &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> &grid_struct, - std::vector> &phi_history, - unsigned int Nx_out, const std::string &filename); - -void save_Efield(unsigned int it, std::vector> &grid_versions, - std::vector> &phi_history, - unsigned int Nx_out = Parameters::PLOT_NX); - #endif diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 51c9ad3..95fc568 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -299,22 +299,22 @@ void NuFISolver::run() { //====//====// // Plotting // //====//====// + if (it % Parameters::PLOT_FREQUENCY == 0) { double plot_start = timer.elapsed(); std::cout << "Saving results... "; - 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, + DiagnosticsSnapshot snap = + compute_diagnostics(*this, it, grid_versions, phi_history, + 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"); - save_Efield(it, grid_versions, phi_history); - - double int_val = 0.5 * integral_space_vector_squared( - grid_versions[phi_history[it].grid_version], - phi_history[it].solution); - int_E_squared.push_back(int_val); + 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"); @@ -322,6 +322,7 @@ void NuFISolver::run() { plot_time = timer.elapsed() - plot_start; std::cout << "Results saved in " << plot_start << "[s]" << "\n"; } + total_time = total_timer.elapsed(); std::cout << "Time since start = " << total_time << "\n\n"; diff --git a/src/save_results.cc b/src/save_results.cc index 8b0ea0e..cee8cf3 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -38,6 +38,78 @@ void save_field_1d(const std::vector &x, file << x[i] << " " << values[i] << "\n"; } +DiagnosticsSnapshot +compute_diagnostics(const NuFISolver &solver, unsigned int n, + std::vector> &grid_struct, + std::vector> &phi_history, + unsigned int Nx_out, unsigned int Nv_out) { + DiagnosticsSnapshot snap; + snap.Nx = Nx_out; + snap.Nv = Nv_out; + + const double vmin = Parameters::V_DOMAIN_LEFT; + const double vmax = Parameters::V_DOMAIN_RIGHT; + const double dv = (vmax - vmin) / Nv_out; + + snap.x_eval = make_x_eval(Nx_out); + snap.v_eval.resize(Nv_out); + for (unsigned int j = 0; j < Nv_out; ++j) + snap.v_eval[j] = vmin + (j + 0.5) * dv; + + snap.f.resize(static_cast(Nx_out) * Nv_out); + for (unsigned int j = 0; j < Nv_out; ++j) { + std::vector val = + solver.eval_f(n, snap.x_eval, snap.v_eval[j], grid_struct, phi_history); + std::copy(val.begin(), val.end(), snap.f.begin() + j * Nx_out); + } + + snap.rho.assign(Nx_out, 1.0); + for (unsigned int j = 0; j < Nv_out; ++j) + for (unsigned int i = 0; i < Nx_out; ++i) + snap.rho[i] -= snap.f[j * Nx_out + i] * dv; + + 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]; + + return snap; +} + +void save_f(const DiagnosticsSnapshot &snap, const std::string &filepath) { + ensure_results_dir(); + std::ofstream file(filepath); + if (!file) + throw std::runtime_error("save_f: failed to open " + filepath); + + file << "# nufi f(x,v) ascii dump\n"; + file << "# Nx = " << snap.Nx << " Nv = " << snap.Nv << "\n"; + file << "# columns: x v f\n"; + file << std::setprecision(10); + for (unsigned int j = 0; j < snap.Nv; ++j) + for (unsigned int i = 0; i < snap.Nx; ++i) + file << snap.x_eval[i] << " " << snap.v_eval[j] << " " + << snap.f[j * snap.Nx + i] << "\n"; +} + +void save_rho(const DiagnosticsSnapshot &snap, const std::string &filepath) { + save_field_1d(snap.x_eval, snap.rho, filepath); +} + +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 integral = 0.0; + for (double e : snap.E) + integral += e * e; + return 0.5 * integral * dx; +} + void save_time_series(const std::vector &t, const std::vector &values, const std::string &filepath) { @@ -57,70 +129,3 @@ void save_time_series(const std::vector &t, 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> &grid_struct, - std::vector> &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 x_eval = make_x_eval(Nx_out); - std::vector 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(&magic), sizeof(magic)); - file.write(reinterpret_cast(&nx64), sizeof(nx64)); - file.write(reinterpret_cast(&nv64), sizeof(nv64)); - file.write(reinterpret_cast(x_eval.data()), - x_eval.size() * sizeof(double)); - file.write(reinterpret_cast(v_eval.data()), - v_eval.size() * sizeof(double)); - - std::vector val(Nx_out); - for (unsigned int j = 0; j < Nv_out; ++j) { - val = solver.eval_f(n, x_eval, v_eval[j], grid_struct, phi_history); - file.write(reinterpret_cast(val.data()), - val.size() * sizeof(double)); - } -} - -void save_rho(const NuFISolver &solver, unsigned int n, - std::vector> &grid_struct, - std::vector> &phi_history, - unsigned int Nx_out, const std::string &filename) { - std::vector x_eval = make_x_eval(Nx_out); - std::vector rho = - solver.eval_rho(n, x_eval, grid_struct, phi_history); - save_field_1d(x_eval, rho, filename); -} - -void save_Efield(unsigned int it, std::vector> &grid_versions, - std::vector> &phi_history, - unsigned int Nx_out) { - std::vector x_eval_E = make_x_eval(Nx_out); - - auto grad_phi = eval(x_eval_E, grid_versions[phi_history[it].grid_version], - phi_history[it].solution); - - std::vector E(x_eval_E.size()); - for (size_t i = 0; i < x_eval_E.size(); ++i) - E[i] = -grad_phi[i]; - - 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"; -}