diff --git a/libnufi_lib.a b/libnufi_lib.a index 65a365c..e57c4b4 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/fields.h b/nufi/fields.h index 4e33a2d..54151e3 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -17,13 +17,9 @@ inline std::vector 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 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 eval(std::vector &X, inline double integral_space_vector(const GridStructure<1> &grid, const Vector &solution, - double dx = Parameters::PLOT_DX, size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; - double xmin = Parameters::X_DOMAIN_LEFT; - std::vector x_eval(Nx); - for (size_t i = 0; i < Nx; ++i) - x_eval[i] = xmin + i * dx; + std::vector x_eval = make_x_eval(Nx); + const double dx = Parameters::LX / Nx; std::vector 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 &solution, - double dx = Parameters::PLOT_DX, size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; - double xmin = Parameters::X_DOMAIN_LEFT; - std::vector x_eval(Nx); - for (size_t i = 0; i < Nx; ++i) - x_eval[i] = xmin + i * dx; + std::vector x_eval = make_x_eval(Nx); + const double dx = Parameters::LX / Nx; std::vector tmp = eval(x_eval, grid, solution); for (size_t i = 0; i < Nx; ++i) diff --git a/nufi/save_results.h b/nufi/save_results.h index c5f18fa..03d1cb8 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -2,17 +2,25 @@ #define SAVE_RESULTS_H #include "nufi/grids.h" -#include "nufi/poisson_problem.h" #include #include #include class NuFISolver; -void save_f(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 &filename); +void save_field_1d(const std::vector &x, + const std::vector &values, + const std::string &filepath); + +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, @@ -23,7 +31,4 @@ void save_Efield(unsigned int it, std::vector> &grid_versions, std::vector> &phi_history, unsigned int Nx_out = Parameters::PLOT_NX); -void save_space_vector(const std::vector &vals, - const std::string &filename, size_t it); - #endif diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index c876e28..51c9ad3 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -193,6 +193,8 @@ void NuFISolver::run() { std::vector int_E_squared; int_E_squared.reserve(Nt); + std::vector int_E_squared_times; // <-- add this + int_E_squared_times.reserve(Nt); std::vector> grid_versions; std::vector> 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"; diff --git a/src/save_results.cc b/src/save_results.cc index 09ca5ca..8b0ea0e 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -7,69 +7,104 @@ #include "nufi/poisson_problem.h" #include #include +#include #include #include #include #include -void save_f(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 &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 &x, + const std::vector &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 &t, + const std::vector &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> &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) { - 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(val.data()), + val.size() * sizeof(double)); } - - file.close(); } 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::ofstream file(filename); - - double xmin = Parameters::X_DOMAIN_LEFT; - double xmax = Parameters::X_DOMAIN_RIGHT; - std::vector x_eval = make_x_eval(Nx_out); - file << Nx_out << "\n"; - file << xmin << " " << xmax << "\n"; - - std::vector tmp = + std::vector 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> &grid_versions, @@ -80,37 +115,12 @@ void save_Efield(unsigned int it, std::vector> &grid_versions, auto grad_phi = eval(x_eval_E, grid_versions[phi_history[it].grid_version], phi_history[it].solution); - std::vector> ordered_E; + std::vector 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 &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"; }