diff --git a/libnufi_lib.a b/libnufi_lib.a index 21d0353..e581d10 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/fields.h b/nufi/fields.h index a91341a..c8951be 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -97,8 +97,7 @@ inline std::vector eval(std::vector &X, Points[i][0] = X[i]; } - return eval_vector_grad(poisson.get_mapping(), poisson.get_dof_handler(), - solution, Points); + return poisson.eval_vector_grad(solution, Points); } inline double integral_space_vector(const PoissonProblem<1> &poisson, diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index fecbe53..fb26183 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -22,6 +22,7 @@ #include #include +#include #include #include @@ -41,7 +42,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -72,6 +75,12 @@ public: std::vector sample_electric_potential(double x_min, double x_max, unsigned int Nx); + std::vector + eval_vector_grad(const Vector &solution, + const std::vector> &points) const; + + void save_grid_to_file(const std::string &filename) const; + private: void create_mesh(); void setup_system(); @@ -93,6 +102,14 @@ private: std::function &)> rhs_function; MappingQ mapping; + + // typename DoFHandler::active_cell_iterator + // find_cell(const DoFHandler &dof_handler, const Point x); + + std::vector::active_cell_iterator> active_cells; + + mutable std::vector local_solution_buffer; + mutable std::unique_ptr> evaluator; }; // Utilities @@ -173,30 +190,56 @@ PoissonProblem::sample_electric_potential(double x_min, double x_max, } template -double -eval_point_grad(const Mapping &mapping, const DoFHandler &dof_handler, - const Vector &solution, const Point &point) { +std::vector PoissonProblem::eval_vector_grad( + const Vector &solution, + const std::vector> &points) const { - Tensor<1, dim> grad = - VectorTools::point_gradient(mapping, dof_handler, solution, point); - double Ex = grad[0]; + const unsigned int n_cells = active_cells.size(); + const double xmin = Parameters::X_DOMAIN_LEFT; + const double h = Parameters::LX / static_cast(n_cells); - return Ex; + std::vector values(points.size()); + + for (unsigned int p = 0; p < points.size(); ++p) { + const double x = points[p][0]; + + const unsigned int cell_index = + std::min(static_cast(std::floor((x - xmin) / h)), + n_cells - 1); // index is at most n_cells - 1 + + const auto cell = active_cells[cell_index]; + + const double xi = (x - cell->vertex(0)[0]) / h; + + Point unit_point; + unit_point[0] = xi; + + cell->get_dof_values(solution, local_solution_buffer.begin(), + local_solution_buffer.end()); + + evaluator->reinit(cell, ArrayView>(&unit_point, 1)); + + evaluator->evaluate(local_solution_buffer, EvaluationFlags::gradients); + + values[p] = evaluator->get_gradient(0)[0]; + } + + return values; } -template -std::vector eval_vector_grad(const Mapping &mapping, - const DoFHandler &dof_handler, - const Vector &solution, - const std::vector> &points) { - size_t p_size = points.size(); - - std::vector Ex(p_size); - for (size_t i = 0; i < p_size; ++i) - Ex[i] = eval_point_grad(mapping, dof_handler, solution, points[i]); - - return Ex; -} +// template +// std::vector eval_vector_grad(const Mapping &mapping, +// const DoFHandler &dof_handler, +// const Vector &solution, +// const std::vector> &points) { +// size_t p_size = points.size(); +// +// std::vector Ex(p_size); +// for (size_t i = 0; i < p_size; ++i) +// Ex[i] = eval_point_grad(mapping, dof_handler, solution, points[i]); +// +// return Ex; +// } template double eval_point_value(const Mapping &mapping, @@ -206,6 +249,16 @@ double eval_point_value(const Mapping &mapping, return VectorTools::point_value(mapping, dof_handler, solution, point); } + +template +void PoissonProblem::save_grid_to_file(const std::string &filename) const { + std::ofstream out(filename); + GridOut grid_out; + grid_out.write_svg(triangulation, out); + + std::cout << "Grid written to " << filename << "\n"; +} + // dealii Poisson template void PoissonProblem::create_mesh() { @@ -262,6 +315,17 @@ template void PoissonProblem::setup_system() { solution.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); + + // used for evaluator to avoid running it anytime there is an eval + active_cells.clear(); + + for (auto cell = dof_handler.begin_active(); cell != dof_handler.end(); + ++cell) + active_cells.push_back(cell); + + local_solution_buffer.resize(fe.n_dofs_per_cell()); + evaluator = std::make_unique>(mapping, fe, + update_gradients); } // Paul diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index ec4edd9..d5aa1dc 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -9,6 +9,7 @@ #include #include +#include #include #include #include @@ -34,6 +35,7 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector &X, double u, if (n == 0) { for (size_t i = 0; i < x_size; ++i) results[i] = f0(X[i], U[i]); + reset_x_eval(X); return results; } @@ -65,6 +67,7 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector &X, double u, } for (size_t i = 0; i < x_size; ++i) results[i] = f0(X[i], U[i]); + reset_x_eval(X); return results; } @@ -80,6 +83,7 @@ NuFISolver::eval_f(unsigned int n, std::vector &X, double u, if (n == 0) { for (size_t i = 0; i < x_size; ++i) results[i] = f0(X[i], U[i]); + reset_x_eval(X); return results; } @@ -116,6 +120,7 @@ NuFISolver::eval_f(unsigned int n, std::vector &X, double u, for (size_t i = 0; i < x_size; ++i) results[i] = f0(X[i], U[i]); + reset_x_eval(X); return results; } @@ -134,14 +139,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); + tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson, + phi_history); // used eval_ftilda once per i 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,20 +160,20 @@ 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; - size_t CALC_NX = Parameters::CALC_NX; - std::vector x_eval = make_x_eval(CALC_NX); + 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; - std::ofstream time_file("results/simulation_time.txt"); - time_file << "# it step_time total_time" << "\n"; + time_file << "it step_time total_time" << "\n"; const double x_min = Parameters::X_DOMAIN_LEFT; double dx = Parameters::CALC_DX; @@ -183,10 +188,11 @@ void NuFISolver::run() { << std::endl; // compute rho + std::vector x_eval = make_x_eval(Nx); std::vector tmp_rho = eval_rho(it, x_eval, poisson, phi_history, Parameters::NV); - for (size_t i = 0; i < CALC_NX; ++i) { + for (size_t i = 0; i < Nx; i++) { AssertThrow(std::isfinite(tmp_rho[i]), ExcMessage("NaN detected in rho")); rho.get()[i] = tmp_rho[i]; } @@ -223,12 +229,12 @@ void NuFISolver::run() { // save_Efield(it, coeffs.get(), 128, "results/field_" + // std::to_string(it) + ".dat"); - 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]; + tmp_rho = eval(x_eval_Ex, poisson, phi_history[it]); + + 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 = diff --git a/src/save_results.cc b/src/save_results.cc index 7f248f5..bb6fbf5 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -38,7 +38,7 @@ void save_f(const NuFISolver &solver, unsigned int n, for (unsigned int i = 0; i < Nx_out; ++i) { file << val[i]; - if (j < Nx_out - 1) + if (i < Nx_out - 1) file << " "; }