diff --git a/libnufi_lib.a b/libnufi_lib.a index a9045c0..5703b05 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/cells.h b/nufi/cells.h index 95a8ac2..f30c239 100644 --- a/nufi/cells.h +++ b/nufi/cells.h @@ -6,7 +6,9 @@ #include #include #include +#include #include +#include #include using namespace dealii; diff --git a/nufi/fields.h b/nufi/fields.h index 46e58f6..7cb1f83 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -1,6 +1,7 @@ #ifndef FIELDS_H #define FIELDS_H +#include "grids.h" #include "nufi/parameters.h" #include "nufi/poisson_problem.h" #include @@ -11,54 +12,6 @@ using namespace dealii; -// inline std::vector Indices_of_points(const std::vector &points, -// double x_min, double x_max, double dx, int grid_type=0) -// { -// // grid type: -// // 0 => uniform -// // 1 => non uniform (TODO) -// -// if (dx <= 0.0) { -// throw std::invalid_argument("dx must be positive"); -// } -// if (x_max <= x_min) { -// throw std::invalid_argument("x_max must be > x_min"); -// } -// -// std::vector indices; -// indices.reserve(points.size()); -// -// switch (grid_type) { -// case 0: -// { -// const double L = x_max - x_min; -// const int N = std::floor(L/dx); -// -// -// for (double x : points) //GPT loop, to check -// { -// x-= x_min; -// x = x - L * std::floor(x/L); -// -// int i = static_cast(std::floor(x / dx)); -// -// // safety: handle rare edge case due to floating precision -// if (i == N) i = 0; -// -// indices.push_back(i); -// } -// } -// case 1: -// { -// throw std::invalid_argument("Case for non uniform grid is not -// completed"); -// } -// default: -// throw std::invalid_argument("Invalid grid_type argument"); -// -// } -// return indices; -// } inline std::vector make_x_eval(size_t Nx) { std::vector x_eval(Nx); const double dx = Parameters::LX / Nx; @@ -86,7 +39,7 @@ inline double f0(const double x, const double v, // wrapper for eval_point() { VectorTools::point_values() } inline std::vector eval(std::vector &X, - const PoissonProblem<1> &poisson, + const GridStructure<1> &grid, const Vector &solution) noexcept { size_t x_size = X.size(); std::vector evals(x_size); @@ -99,10 +52,10 @@ inline std::vector eval(std::vector &X, Points[i][0] = X[i]; } - return poisson.eval_vector_grad(solution, Points); + return grid.eval_vector_grad(solution, Points); } -inline double integral_space_vector(const PoissonProblem<1> &poisson, +inline double integral_space_vector(const GridStructure<1> &grid, const Vector &solution, double dx = Parameters::PLOT_DX, size_t Nx = Parameters::PLOT_NX) { @@ -112,13 +65,13 @@ inline double integral_space_vector(const PoissonProblem<1> &poisson, for (size_t i = 0; i < Nx; ++i) x_eval[i] = xmin + i * dx; - std::vector tmp = eval(x_eval, poisson, solution); + std::vector tmp = eval(x_eval, grid, solution); for (size_t i = 0; i < Nx; ++i) integral += tmp[i]; return integral * dx; }; -inline double integral_space_vector_squared(const PoissonProblem<1> &poisson, +inline double integral_space_vector_squared(const GridStructure<1> &grid, const Vector &solution, double dx = Parameters::PLOT_DX, size_t Nx = Parameters::PLOT_NX) { @@ -128,7 +81,7 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson, for (size_t i = 0; i < Nx; ++i) x_eval[i] = xmin + i * dx; - std::vector tmp = eval(x_eval, poisson, solution); + std::vector tmp = eval(x_eval, grid, solution); for (size_t i = 0; i < Nx; ++i) integral += tmp[i] * tmp[i]; return integral * dx; diff --git a/nufi/grids.h b/nufi/grids.h new file mode 100644 index 0000000..66a5e0b --- /dev/null +++ b/nufi/grids.h @@ -0,0 +1,116 @@ +#ifndef GRIDS_H +#define GRIDS_H + +#include "nufi/cells.h" +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +template class PoissonProblem; + +template struct GridStructure { + //==//==// + // Vars // + //==//==// + std::unique_ptr> triangulation; + std::unique_ptr> dof_handler; + std::unique_ptr> mapping; + std::unique_ptr> fe; + + CellLocator locator; + + unsigned int grid_version; + + // === // === // + // Evaluator // + // === // === // + std::vector + eval_vector_grad(const Vector &solution, + const std::vector> &points) const { + + std::vector values(points.size()); +#pragma omp parallel + { + std::vector local_solution_buffer(fe->n_dofs_per_cell()); + FEPointEvaluation evaluator(*mapping, *fe, update_gradients); +#pragma omp for + for (unsigned int p = 0; p < points.size(); ++p) { + + const auto cell_location = locator.locate(points[p]); + + cell_location.info->cell->get_dof_values(solution, + local_solution_buffer.begin(), + local_solution_buffer.end()); + + evaluator.reinit( + cell_location.info->cell, + ArrayView>(&cell_location.reference_point, 1)); + + evaluator.evaluate(local_solution_buffer, EvaluationFlags::gradients); + + values[p] = evaluator.get_gradient(0)[0]; + } + } + return values; + } +}; + +template +GridStructure make_grid_snapshot(const PoissonProblem &poisson) { + GridStructure grid; + grid.grid_version = 0; + + grid.triangulation = std::make_unique>(); + grid.triangulation->copy_triangulation(poisson.get_triangulation()); + + grid.mapping = std::make_unique>(poisson.get_mapping()); + + grid.dof_handler = std::make_unique>(*grid.triangulation); + grid.dof_handler->distribute_dofs(poisson.get_dof_handler().get_fe()); + + grid.locator.rebuild(*grid.dof_handler, *grid.triangulation); + + grid.fe = std::make_unique>(poisson.get_fe()); + + return grid; +} + +template struct SolutionSnapshot { + unsigned int grid_version; + Vector solution; +}; + +template +inline void update_grid_versions(std::vector> &grid_versions, + PoissonProblem &poisson) { + auto grid = make_grid_snapshot(poisson); + + if (!grid_versions.empty()) + grid.grid_version = grid_versions.back().grid_version + 1; + + grid_versions.push_back(std::move(grid)); +} + +template +inline void +update_solution_history(std::vector> &solution_history, + PoissonProblem &poisson, + unsigned int current_grid_version) { + + SolutionSnapshot snapshot; + + snapshot.grid_version = current_grid_version; + + Vector solution = poisson.get_solution(); + snapshot.solution = solution; + + solution_history.push_back(std::move(snapshot)); +} + +#endif // !GRIDS_H diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h index 5dc46fa..f9a9d18 100644 --- a/nufi/nufi_solver.h +++ b/nufi/nufi_solver.h @@ -9,6 +9,7 @@ #include #include "nufi/fields.h" //dont remove +#include "nufi/grids.h" #include "nufi/parameters.h" #include "nufi/poisson_problem.h" @@ -19,18 +20,19 @@ public: NuFISolver(); void run(); - std::vector eval_rho(unsigned int n, std::vector &x, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history, - const unsigned int Nv = Parameters::NV) const; + std::vector + eval_rho(unsigned int n, std::vector &x, + const std::vector> &grid_struct, + const std::vector> &phi_history, + const unsigned int Nv = Parameters::NV) const; std::vector eval_ftilda(unsigned int, std::vector &x, double u, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history) const; + const std::vector> &grid_struct, + const std::vector> &phi_history) const; std::vector eval_f(unsigned int n, std::vector &x, double u, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history) const; + const std::vector> &grid_struct, + const std::vector> &phi_history) const; private: unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT); diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index 641cdb2..4a2517f 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -54,7 +54,9 @@ #include #include "nufi/cells.h" +#include "nufi/grids.h" #include "nufi/parameters.h" +#include "omp.h" using namespace dealii; @@ -66,8 +68,7 @@ public: void initialize(); void solve_step(); - void coarse_and_refine_grid(size_t it, - std::vector> &solution_history); + void coarse_and_refine_grid(size_t it); void run(); unsigned int get_rhs_size(); @@ -78,6 +79,9 @@ public: const Vector &get_solution() const { return solution; } const MappingQ &get_mapping() const { return mapping; } const DoFHandler &get_dof_handler() const { return dof_handler; } + const Triangulation &get_triangulation() const { return triangulation; } + const FE_Q &get_fe() const { return fe; } + const CellLocator &get_locator() const { return cell_locator; } std::vector sample_electric_field(double x_min, double x_max, unsigned int Nx); @@ -117,8 +121,8 @@ private: MappingQ mapping; - mutable std::vector local_solution_buffer; - mutable std::unique_ptr> evaluator; + // mutable std::vector local_solution_buffer; + // mutable std::unique_ptr> evaluator; }; //====//====// @@ -211,31 +215,6 @@ PoissonProblem::sample_electric_potential(double x_min, double x_max, return values; } -template -std::vector PoissonProblem::eval_vector_grad( - const Vector &solution, - const std::vector> &points) const { - - std::vector values(points.size()); - - for (unsigned int p = 0; p < points.size(); ++p) { - - const auto cell_location = cell_locator.locate(points[p]); - - cell_location.info->cell->get_dof_values( - solution, local_solution_buffer.begin(), local_solution_buffer.end()); - - evaluator->reinit( - cell_location.info->cell, - ArrayView>(&cell_location.reference_point, 1)); - evaluator->evaluate(local_solution_buffer, EvaluationFlags::gradients); - - values[p] = evaluator->get_gradient(0)[0]; - } - - return values; -} - template std::vector eval_point_grad(const Mapping &mapping, const DoFHandler &dof_handler, @@ -347,9 +326,9 @@ template void PoissonProblem::setup_system() { // used for evaluator to avoid running it anytime there is an eval cell_locator.rebuild(dof_handler, triangulation); - local_solution_buffer.resize(fe.n_dofs_per_cell()); - evaluator = std::make_unique>(mapping, fe, - update_gradients); + // local_solution_buffer.resize(fe.n_dofs_per_cell()); + // evaluator = std::make_unique>(mapping, fe, + // update_gradients); } // Paul @@ -402,9 +381,7 @@ template void PoissonProblem::assemble_system() { } } -template -void PoissonProblem::coarse_and_refine_grid( - size_t it, std::vector> &solution_history) { +template void PoissonProblem::coarse_and_refine_grid(size_t it) { std::cout << "Refinement Started" << "\n"; Vector error_per_cell(triangulation.n_active_cells()); @@ -417,29 +394,20 @@ void PoissonProblem::coarse_and_refine_grid( 0.3, 0.03); triangulation.prepare_coarsening_and_refinement(); + SolutionTransfer> transfer(dof_handler); - transfer.prepare_for_coarsening_and_refinement(solution_history); + const Vector refined_solution = solution; + + transfer.prepare_for_coarsening_and_refinement(refined_solution); triangulation.execute_coarsening_and_refinement(); setup_system(); - std::vector> new_solution_history(solution_history.size()); - - for (auto &vec : new_solution_history) - vec.reinit(dof_handler.n_dofs()); - - transfer.interpolate(solution_history, new_solution_history); - - solution_history.swap(new_solution_history); - - solution = solution_history.back(); + transfer.interpolate(refined_solution, solution); constraints.distribute(solution); - // cell_locator.rebuild(dof_handler, triangulation); // No need to be called - // again because its in setup_system(); - std::cout << "Refinement Finished" << "\n"; std::string grid_file_name = diff --git a/nufi/save_results.h b/nufi/save_results.h index 7e47ae7..8bb7d53 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -1,23 +1,25 @@ #ifndef SAVE_RESULTS_H #define SAVE_RESULTS_H +#include "nufi/grids.h" #include "nufi/nufi_solver.h" #include "nufi/poisson_problem.h" #include #include +#include void save_f(const NuFISolver &solver, unsigned int n, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history, unsigned int Nx_out, + std::vector> &grid_struct, + std::vector> &phi_history, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename); void save_rho(const NuFISolver &solver, unsigned int n, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history, + std::vector> &grid_struct, + std::vector> &phi_history, unsigned int Nx_out, const std::string &filename); -void save_Efield(unsigned int n, const PoissonProblem<1> &poisson, - const std::vector> &phi_history, +void save_Efield(unsigned int n, GridStructure<1> &grid_struct, + std::vector> &phi_history, unsigned int Nx_out, const std::string &filename); void save_space_vector(const std::vector &vals, diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 0d6131f..84f68fe 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -17,6 +17,7 @@ #include #include "nufi/fields.h" +#include "nufi/grids.h" #include "nufi/parameters.h" #include "nufi/poisson_problem.h" #include "nufi/save_results.h" @@ -24,10 +25,10 @@ using namespace dealii; -std::vector -NuFISolver::eval_ftilda(unsigned int n, std::vector &X, double u, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history) const { +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(); @@ -48,7 +49,8 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector &X, double u, for (size_t i = 0; i < x_size; ++i) X[i] = X[i] - Parameters::DT * U[i]; - tmp = eval(X, poisson, phi_history[n]); // call eval only once + 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]; @@ -60,7 +62,8 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector &X, double u, for (size_t i = 0; i < x_size; ++i) X[i] = X[i] - Parameters::DT * U[i]; - tmp = eval(X, poisson, phi_history[n]); // call eval only once + 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]; @@ -74,8 +77,8 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector &X, double u, std::vector NuFISolver::eval_f(unsigned int n, std::vector &X, double u, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history) const { + const std::vector> &grid_struct, + const std::vector> &phi_history) const { size_t x_size = X.size(); @@ -92,7 +95,8 @@ NuFISolver::eval_f(unsigned int n, std::vector &X, double u, std::vector tmp(x_size); // Initial half-step. - tmp = eval(X, poisson, phi_history[n]); // call eval only once + 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]; @@ -102,7 +106,8 @@ NuFISolver::eval_f(unsigned int n, std::vector &X, double u, for (size_t i = 0; i < x_size; ++i) X[i] = X[i] - Parameters::DT * U[i]; - tmp = eval(X, poisson, phi_history[n]); // call eval only once + 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]; @@ -113,7 +118,8 @@ NuFISolver::eval_f(unsigned int n, std::vector &X, double u, for (size_t i = 0; i < x_size; ++i) X[i] = X[i] - Parameters::DT * U[i]; - tmp = eval(X, poisson, phi_history[n]); // call eval only once + 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]; @@ -127,20 +133,21 @@ NuFISolver::eval_f(unsigned int n, std::vector &X, double u, std::vector NuFISolver::eval_rho(unsigned int n, std::vector &X, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history, + const std::vector> &grid_struct, + const std::vector> &phi_history, const unsigned int Nv) const { size_t x_size = X.size(); + 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 integral(x_size, 0.0); - 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, + tmp_int = eval_ftilda(n, X, v_min + i * dv, grid_struct, phi_history); // used eval_ftilda once per i for (size_t ii = 0; ii < x_size; ++ii) integral[ii] += tmp_int[ii]; @@ -164,16 +171,21 @@ 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; + std::vector> grid_versions; + std::vector> phi_history; + + update_grid_versions(grid_versions, poisson); + update_solution_history(phi_history, poisson, + grid_versions.back().grid_version); std::vector x_eval(Parameters::CALC_NX); - if (rho == nullptr) - throw std::bad_alloc{}; - std::ofstream time_file("results/simulation_time.dat"); double total_time = 0; @@ -224,7 +236,7 @@ void NuFISolver::run() { // compute rho std::vector x_eval = make_x_eval(poisson.get_dof_size()); std::vector rho_values = - eval_rho(it, x_eval, poisson, phi_history, Parameters::NV); + eval_rho(it, x_eval, grid_versions, phi_history, Parameters::NV); Vector rhs(x_eval.size()); for (unsigned int i = 0; i < rhs.size(); ++i) @@ -233,18 +245,23 @@ void NuFISolver::run() { poisson.set_rhs(rhs); poisson.solve_step(); - phi_history.push_back(poisson.get_solution()); compute_time = timer.elapsed() - compute_start; if (it % Parameters::REFINE_FREQUENCY == 0 && it != 0) { double refine_start = timer.elapsed(); - poisson.coarse_and_refine_grid(it, phi_history); + + poisson.coarse_and_refine_grid(it); + + update_grid_versions(grid_versions, poisson); + refine_time = timer.elapsed() - refine_start; std::cout << "Refinement step done in " << std::to_string(std::round(std::floor(refine_time))) << "[s]" << "\n"; } + update_solution_history(phi_history, poisson, + grid_versions.back().grid_version); double timer_elapsed = timer.elapsed(); double step_time = timer_elapsed - time_elapsed_before; @@ -253,30 +270,33 @@ void NuFISolver::run() { if (it % Parameters::PLOT_FREQUENCY == 0) { double plot_start = timer.elapsed(); std::cout << "Saving results... "; - save_f(*this, it, poisson, phi_history, Parameters::PLOT_NX, + save_f(*this, it, grid_versions, phi_history, Parameters::PLOT_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat"); - save_rho(*this, it, poisson, phi_history, Parameters::PLOT_NX, + save_rho(*this, it, grid_versions, phi_history, Parameters::PLOT_NX, "results/rho_" + std::to_string(it) + ".dat"); // save_Efield(it, coeffs.get(), 128, "results/field_" + // std::to_string(it) + ".dat"); std::vector x_eval_Ex = make_x_eval(Parameters::PLOT_NX); std::vector tmp_rho(x_eval_Ex.size()); - tmp_rho = eval(x_eval_Ex, poisson, phi_history[it]); + tmp_rho = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version], + phi_history[it].solution); std::vector E_x(Parameters::PLOT_NX); for (size_t i = 0; i < Parameters::PLOT_NX; ++i) E_x[i] = -tmp_rho[i]; save_space_vector(E_x, "field", it); - double int_val = - 0.5 * integral_space_vector_squared(poisson, phi_history[it]); + 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); save_space_vector(int_E_squared, "electricint", it); std::cout << "Time since start = " << total_time << "\n\n"; plot_time = timer.elapsed() - plot_start; } + total_time = timer.elapsed(); time_file << it << " " << step_time << " " << total_time << " " << compute_time << " " << refine_time << " " << plot_time << "\n"; diff --git a/src/save_results.cc b/src/save_results.cc index bb6fbf5..bd787d3 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -1,18 +1,20 @@ #include "nufi/save_results.h" #include "nufi/fields.h" +#include "nufi/grids.h" #include "nufi/nufi_solver.h" #include "nufi/parameters.h" #include "nufi/poisson_problem.h" #include +#include #include #include #include #include void save_f(const NuFISolver &solver, unsigned int n, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history, unsigned int Nx_out, + std::vector> &grid_struct, + std::vector> &phi_history, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename) { std::ofstream file(filename); @@ -33,7 +35,7 @@ void save_f(const NuFISolver &solver, unsigned int n, 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); + val = solver.eval_f(n, x_eval, v, grid_struct, phi_history); for (unsigned int i = 0; i < Nx_out; ++i) { file << val[i]; @@ -49,8 +51,8 @@ void save_f(const NuFISolver &solver, unsigned int n, } void save_rho(const NuFISolver &solver, unsigned int n, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history, + std::vector> &grid_struct, + std::vector> &phi_history, unsigned int Nx_out, const std::string &filename) { std::ofstream file(filename); @@ -61,7 +63,8 @@ void save_rho(const NuFISolver &solver, unsigned int n, file << Nx_out << "\n"; file << xmin << " " << xmax << "\n"; - std::vector tmp = solver.eval_rho(n, x_eval, poisson, phi_history); + std::vector tmp = + solver.eval_rho(n, x_eval, grid_struct, phi_history); for (size_t i = 0; i < Nx_out; ++i) { file << tmp[i]; file << "\n"; @@ -69,9 +72,8 @@ void save_rho(const NuFISolver &solver, unsigned int n, file.close(); } -void save_Efield([[maybe_unused]] unsigned int n, - const PoissonProblem<1> &poisson, - const std::vector> &phi_history, +void save_Efield([[maybe_unused]] unsigned int n, GridStructure<1> &grid_struct, + std::vector> &phi_history, unsigned int Nx_out, const std::string &filename) { std::ofstream file(filename); @@ -85,7 +87,7 @@ void save_Efield([[maybe_unused]] unsigned int n, file << Nx_out << "\n"; file << xmin << " " << xmax << "\n"; - std::vector tmp = eval(x_eval, poisson, phi_history[n]); + std::vector tmp = eval(x_eval, grid_struct, phi_history[n].solution); for (size_t i = 0; i < Nx_out; ++i) { file << -tmp[i]; file << "\n";