From 2b42bb7d02c164ec1127353587b865e79cb88fe3 Mon Sep 17 00:00:00 2001 From: "Vasco C. B. Ferreira" Date: Wed, 1 Jul 2026 16:34:27 +0200 Subject: [PATCH] update save_results part 1 --- nufi/fields.h | 37 ++++++++++++++++++++++++++----------- src/nufi_solver.cc | 35 ++++++++++++++++++----------------- src/save_results.cc | 34 +++++++++++++++++----------------- 3 files changed, 61 insertions(+), 45 deletions(-) diff --git a/nufi/fields.h b/nufi/fields.h index 252e819..a91341a 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -59,6 +59,18 @@ using namespace dealii; // } // return indices; // } +inline std::vector make_x_eval(size_t Nx) { + std::vector x_eval(Nx); + for (size_t i = 0; i < Nx; ++i) + x_eval[i] = Parameters::X_DOMAIN_LEFT + i * Parameters::CALC_DX; + return x_eval; +} + +inline void reset_x_eval(std::vector &x_vals) { + const size_t Nx = x_vals.size(); + for (size_t i = 0; i < Nx; ++i) + x_vals[i] = Parameters::X_DOMAIN_LEFT + i * Parameters::CALC_DX; +}; inline double f0(const double x, const double v, const double eps = Parameters::EPS, @@ -95,11 +107,13 @@ inline double integral_space_vector(const PoissonProblem<1> &poisson, size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; double xmin = Parameters::X_DOMAIN_LEFT; -#pragma omp parallel for reduction(+ : integral) - for (size_t i = 0; i < Nx; ++i) { - double x = xmin + i * dx; - integral += eval(x, poisson, solution); - } + std::vector x_eval(Nx); + for (size_t i = 0; i < Nx; ++i) + x_eval[i] = xmin + i * dx; + + std::vector tmp = eval(x_eval, poisson, solution); + for (size_t i = 0; i < Nx; ++i) + integral += tmp[i]; return integral * dx; }; @@ -109,12 +123,13 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson, size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; double xmin = Parameters::X_DOMAIN_LEFT; -#pragma omp parallel for reduction(+ : integral) - for (size_t i = 0; i < Nx; ++i) { - double x = xmin + i * dx; - double val = eval(x, poisson, solution); - integral += val * val; - } + std::vector x_eval(Nx); + for (size_t i = 0; i < Nx; ++i) + x_eval[i] = xmin + i * dx; + + std::vector tmp = eval(x_eval, poisson, solution); + for (size_t i = 0; i < Nx; ++i) + integral += tmp[i] * tmp[i]; return integral * dx; }; diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 7a512be..ec4edd9 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -134,14 +134,14 @@ NuFISolver::eval_rho(unsigned int n, std::vector &X, std::vector tmp_int(x_size); for (unsigned int i = 0; i < Nv; ++i) { - tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson, - phi_history); // used eval_ftilda once per i + tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson, phi_history); for (size_t ii = 0; ii < x_size; ++ii) integral[ii] += tmp_int[ii]; } for (size_t i = 0; i < x_size; ++i) integral[i] = 1 - integral[i] * dv; + reset_x_eval(X); return integral; } @@ -155,13 +155,16 @@ void NuFISolver::run() { reinterpret_cast(std::aligned_alloc(64, sizeof(double) * Nx)), std::free}; + if (rho == nullptr) + throw std::bad_alloc{}; + std::vector int_E_squared; int_E_squared.reserve(Nt); std::vector> phi_history; - if (rho == nullptr) - throw std::bad_alloc{}; + size_t CALC_NX = Parameters::CALC_NX; + std::vector x_eval = make_x_eval(CALC_NX); double total_time = 0; std::ofstream time_file("results/simulation_time.txt"); @@ -180,14 +183,12 @@ void NuFISolver::run() { << std::endl; // compute rho + std::vector tmp_rho = + eval_rho(it, x_eval, poisson, phi_history, Parameters::NV); -#pragma omp parallel for - for (size_t i = 0; i < Nx; i++) { - double x = Parameters::X_DOMAIN_LEFT + i * dx; - double ith_rho = eval_rho(it, x, poisson, phi_history, Parameters::NV); - - AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho")); - rho.get()[i] = ith_rho; + for (size_t i = 0; i < CALC_NX; ++i) { + AssertThrow(std::isfinite(tmp_rho[i]), ExcMessage("NaN detected in rho")); + rho.get()[i] = tmp_rho[i]; } poisson.set_rhs_function([&rho, x_min, dx, Nx = Nx](const Point<1> &p) { @@ -222,12 +223,12 @@ void NuFISolver::run() { // save_Efield(it, coeffs.get(), 128, "results/field_" + // std::to_string(it) + ".dat"); - std::vector E_x(Nx, 0.0); -#pragma omp parallel for - for (size_t ix = 0; ix < Nx; ++ix) { - E_x[ix] = -eval(x_min + ix * dx, poisson, phi_history[it]); - } - + std::vector E_x(CALC_NX); + std::vector x_eval_Ex = make_x_eval(Parameters::PLOT_NX); + std::vector tmp_Ex = eval(x_eval_Ex, poisson, phi_history[it]); + reset_x_eval(x_eval_Ex); + for (size_t i = 0; i < CALC_NX; ++i) + E_x[i] = -tmp_Ex[i]; save_space_vector(E_x, "field", it); double int_val = diff --git a/src/save_results.cc b/src/save_results.cc index 3f32acd..1223501 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -4,6 +4,7 @@ #include "nufi/nufi_solver.h" #include "nufi/parameters.h" #include "nufi/poisson_problem.h" +#include #include #include #include @@ -21,22 +22,20 @@ void save_f(const NuFISolver &solver, unsigned int n, double vmin = Parameters::V_DOMAIN_LEFT; double vmax = Parameters::V_DOMAIN_RIGHT; - double dx = (xmax - xmin) / Nx_out; double dv = (vmax - vmin) / Nv_out; file << Nx_out << " " << Nv_out << "\n"; file << xmin << " " << xmax << "\n"; file << vmin << " " << vmax << "\n"; - for (unsigned int i = 0; i < Nx_out; ++i) { - double x = xmin + (i + 0.5) * dx; + std::vector x_eval = make_x_eval(Nx_out); + std::vector val(Nx_out); - for (unsigned int j = 0; j < Nv_out; ++j) { - double v = vmin + (j + 0.5) * dv; - - double val = solver.eval_f(n, x, v, poisson, phi_history); - - file << val; + for (unsigned int j = 0; j < Nv_out; ++j) { + double v = vmin + (j + 0.5) * dv; + val = solver.eval_f(n, x_eval, v, poisson, phi_history); + for (unsigned int i = 0; i < Nx_out; ++i) { + file << val[i]; if (j < Nv_out - 1) file << " "; @@ -56,14 +55,14 @@ void save_rho(const NuFISolver &solver, unsigned int n, double xmin = Parameters::X_DOMAIN_LEFT; double xmax = Parameters::X_DOMAIN_RIGHT; - double dx = (xmax - xmin) / Nx_out; + std::vector x_eval = make_x_eval(Nx_out); file << Nx_out << "\n"; file << xmin << " " << xmax << "\n"; - for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) { - double val = solver.eval_rho(n, xmin, poisson, phi_history); - file << val; + std::vector tmp = solver.eval_rho(n, x_eval, poisson, phi_history); + for (size_t i = 0; i < Nx_out; ++i) { + file << tmp[i]; file << "\n"; } file.close(); @@ -77,16 +76,17 @@ void save_Efield([[maybe_unused]] unsigned int n, double xmin = Parameters::X_DOMAIN_LEFT; double xmax = Parameters::X_DOMAIN_RIGHT; - double dx = (xmax - xmin) / Nx_out; + + std::vector x_eval = make_x_eval(Nx_out); // select from E_coeffs file << Nx_out << "\n"; file << xmin << " " << xmax << "\n"; - for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) { - double val = -eval(xmin, poisson, phi_history[n]); - file << val; + std::vector tmp = eval(x_eval, poisson, phi_history[n]); + for (size_t i = 0; i < Nx_out; ++i) { + file << -tmp[i]; file << "\n"; } file.close();