From eeb26d31c14599b9ea90b3636709aab004195a66 Mon Sep 17 00:00:00 2001 From: VCB Ferreira Date: Mon, 3 Aug 2026 00:26:15 +0200 Subject: [PATCH] wip: Branch init Main changes done, eval_rho + eval_ftilda_batch need to be checked to work correcly together --- nufi/nufi_solver.h | 9 ++- src/main.cc | 25 ++++---- src/nufi_solver.cc | 150 +++++++++++++++++++++++---------------------- 3 files changed, 95 insertions(+), 89 deletions(-) diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h index c773f15..1e5eeb4 100644 --- a/nufi/nufi_solver.h +++ b/nufi/nufi_solver.h @@ -23,10 +23,13 @@ public: const std::vector> &grid_struct, const std::vector> &phi_history, const unsigned int Nv = Parameters::NV) const; + std::vector - eval_ftilda(unsigned int n, std::vector x, double u, - const std::vector> &grid_struct, - const std::vector> &phi_history) const; + eval_ftilda_batch(unsigned int n, std::vector X, + std::vector U, + const std::vector> &grid_struct, + const std::vector> &phi_history) const; + std::vector eval_f(unsigned int n, std::vector x, double u, const std::vector> &grid_struct, diff --git a/src/main.cc b/src/main.cc index 8b60171..ffec53a 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,3 +1,4 @@ +#include "nufi/stopwatch.h" #include #include #include @@ -5,18 +6,6 @@ #include #include -void clear_results_directory(const std::string &dir) { - if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir)) - return; - - for (const auto &entry : std::filesystem::directory_iterator(dir)) { - if (std::filesystem::is_regular_file(entry)) { - std::filesystem::remove(entry.path()); - std::cout << "Deleted: " << entry.path() << '\n'; - } - } -}; - template void run() { unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT); @@ -134,6 +123,18 @@ template void run() { std::cout << "NuFI simulation finished in " << total_time << " seconds.\n"; }; +void clear_results_directory(const std::string &dir) { + if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir)) + return; + + for (const auto &entry : std::filesystem::directory_iterator(dir)) { + if (std::filesystem::is_regular_file(entry)) { + std::filesystem::remove(entry.path()); + std::cout << "Deleted: " << entry.path() << '\n'; + } + } +}; + int main() { // omp_set_max_active_levels(1); std::cout << "Threads: " << omp_get_max_threads() << "\n"; diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 480f328..8d0511d 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -20,74 +20,10 @@ #include "nufi/fields.h" #include "nufi/grids.h" #include "nufi/parameters.h" -#include "nufi/poisson_problem.h" #include "nufi/save_results.h" -#include "nufi/stopwatch.h" using namespace dealii; -std::vector NuFISolver::eval_ftilda( - unsigned int n, std::vector X, double u, - const std::vector> &grid_struct, - const std::vector> &phi_history) const { - - size_t x_size = X.size(); - - std::vector U(x_size, u); - std::vector results(x_size); - if (n == 0) { - for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U[i]); - return results; - } - - std::vector Ex(x_size); - std::vector tmp(x_size); - - // We omit the initial half-step. - while (--n) { - for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U[i]; - - AssertThrow( - phi_history[n].solution.size() == - grid_struct[phi_history[n].grid_version].dof_handler->n_dofs(), - ExcMessage("In eval_ftilda: Solution size = " + - std::to_string(phi_history[n].solution.size()) + - ", expected by grid_struct = " + - std::to_string(grid_struct[phi_history[n].grid_version] - .dof_handler->n_dofs()))); - AssertThrow( - grid_struct[phi_history[n].grid_version].grid_version == - phi_history[n].grid_version, - ExcMessage( - "grid.grid_version not equal to phi_history[n].grid_version")); - - tmp = eval(X, grid_struct[phi_history[n].grid_version], - phi_history[n].solution); // call eval only once - - for (size_t i = 0; i < x_size; ++i) { - Ex[i] = -tmp[i]; - U[i] = U[i] + Parameters::DT * Ex[i]; - } - } - - // The final half-step. - for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U[i]; - - tmp = eval(X, grid_struct[phi_history[n].grid_version], - phi_history[n].solution); // call eval only once - - for (size_t i = 0; i < x_size; ++i) { - Ex[i] = -tmp[i]; - U[i] = U[i] + 0.5 * Parameters::DT * Ex[i]; - } - for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U[i]); - return results; -} - std::vector NuFISolver::eval_f(unsigned int n, std::vector X, double u, const std::vector> &grid_struct, @@ -142,32 +78,98 @@ NuFISolver::eval_f(unsigned int n, std::vector X, double u, return results; } +std::vector NuFISolver::eval_ftilda_batch( + unsigned int n, std::vector X, std::vector U, + const std::vector> &grid_struct, + const std::vector> &phi_history) const { + + const size_t x_size = X.size(); + std::vector results(x_size); + + if (n == 0) { + for (size_t k = 0; k < x_size; ++k) + results[k] = f0(X[k], U[k]); + return results; + } + + std::vector Ex(x_size); + std::vector tmp(x_size); + + // We omit the initial half-step. + while (--n) { + for (size_t k = 0; k < x_size; ++k) + X[k] = X[k] - Parameters::DT * U[k]; + + AssertThrow( + phi_history[n].solution.size() == + grid_struct[phi_history[n].grid_version].dof_handler->n_dofs(), + ExcMessage("In eval_ftilda_batch: Solution size = " + + std::to_string(phi_history[n].solution.size()) + + ", expected by grid_struct = " + + std::to_string(grid_struct[phi_history[n].grid_version] + .dof_handler->n_dofs()))); + AssertThrow( + grid_struct[phi_history[n].grid_version].grid_version == + phi_history[n].grid_version, + ExcMessage( + "grid.grid_version not equal to phi_history[n].grid_version")); + + tmp = eval(X, grid_struct[phi_history[n].grid_version], + phi_history[n].solution); // one call for the whole batch + + for (size_t k = 0; k < x_size; ++k) { + Ex[k] = -tmp[k]; + U[k] = U[k] + Parameters::DT * Ex[k]; + } + } + + // The final half-step. + for (size_t k = 0; k < x_size; ++k) + X[k] = X[k] - Parameters::DT * U[k]; + + tmp = eval(X, grid_struct[phi_history[n].grid_version], + phi_history[n].solution); + + for (size_t k = 0; k < x_size; ++k) { + Ex[k] = -tmp[k]; + U[k] = U[k] + 0.5 * Parameters::DT * Ex[k]; + } + for (size_t k = 0; k < x_size; ++k) + results[k] = f0(X[k], U[k]); + return results; +} + std::vector NuFISolver::eval_rho(unsigned int n, std::vector &X, const std::vector> &grid_struct, const std::vector> &phi_history, const unsigned int Nv) const { - size_t x_size = X.size(); + const size_t x_size = X.size(); + const size_t total = static_cast(Nv) * x_size; + + std::vector X_all(total); + std::vector U_all(total); const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; - const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv; - - std::vector partial(static_cast(Nv) * x_size); + const double v_min = Parameters::V_DOMAIN_LEFT; #pragma omp parallel for for (unsigned int i = 0; i < Nv; ++i) { - std::vector tmp_int = - eval_ftilda(n, X, v_min + i * dv, grid_struct, - phi_history); // used eval_ftilda once per i - std::copy(tmp_int.begin(), tmp_int.end(), - partial.begin() + static_cast(i) * x_size); + const double v_i = v_min + i * dv; + std::copy(X.begin(), X.end(), + X_all.begin() + static_cast(i) * x_size); + std::fill(U_all.begin() + static_cast(i) * x_size, + U_all.begin() + static_cast(i + 1) * x_size, v_i); } + std::vector ftilda_all = eval_ftilda_batch( + n, std::move(X_all), std::move(U_all), grid_struct, phi_history); + std::vector integral(x_size, 0.0); for (unsigned int i = 0; i < Nv; ++i) for (size_t ii = 0; ii < x_size; ++ii) - integral[ii] += partial[static_cast(i) * x_size + ii]; + integral[ii] += ftilda_all[static_cast(i) * x_size + ii]; for (size_t i = 0; i < x_size; ++i) integral[i] = 1 - integral[i] * dv;