From eadf717dfca66e239f2bf4ce54201162b2a955ce Mon Sep 17 00:00:00 2001 From: "Vasco C. B. Ferreira" Date: Sat, 14 Mar 2026 18:34:48 +0100 Subject: [PATCH] corrected results output, E is 0 all the time everywhere?? --- nufi_solver.hpp | 60 ++++++++- poisson_non_periodic.hpp | 257 +++++++++++++++++++++++++++++++++++++++ poisson_problem.hpp | 57 --------- 3 files changed, 316 insertions(+), 58 deletions(-) create mode 100644 poisson_non_periodic.hpp diff --git a/nufi_solver.hpp b/nufi_solver.hpp index 036bfe2..cb817a6 100644 --- a/nufi_solver.hpp +++ b/nufi_solver.hpp @@ -27,7 +27,10 @@ public: void run(); double eval_rho(unsigned int n, double x, const double *E_coeffs, unsigned int Nv = Parameters::NV); double eval_ftilda(unsigned int n, double x, double u, const double *E_coeffs); + void save_ftilda(unsigned int n, const double *E_coeffs, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename); + void save_rho(unsigned int n, const double *E_coeffs, unsigned int Nx_out, const std::string &filename); + void save_Efield(unsigned int n, const double *E_coeffs, unsigned int Nx_out, const std::string &filename); private: @@ -165,6 +168,60 @@ inline void NuFISolver::save_ftilda(unsigned int n, file.close(); } +inline void NuFISolver::save_rho(unsigned int n, + const double *E_coeffs, + unsigned int Nx_out, + const std::string &filename) +{ + std::ofstream file(filename); + + double xmin = Parameters::X_DOMAIN_LEFT; + double xmax = Parameters::X_DOMAIN_RIGHT; + double dx = (xmax - xmin) / Nx_out; + + file << Nx_out << "\n"; + file << xmin << " " << xmax << "\n"; + + for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) + { + double val = eval_rho(n, xmin, E_coeffs); + file << val; + file << "\n"; + } + file.close(); +} + +inline void NuFISolver::save_Efield(unsigned int n, + const double *E_coeffs, + unsigned int Nx_out, + const std::string &filename) +{ + std::ofstream file(filename); + + double xmin = Parameters::X_DOMAIN_LEFT; + double xmax = Parameters::X_DOMAIN_RIGHT; + double dx = (xmax - xmin) / Nx_out; + + // select from E_coeffs + const size_t stride_x = 1; + const size_t stride_t = stride_x*(Parameters::SPLINE_NX + Parameters::SPLINE_ORDER - 1); + const double *c; + c = E_coeffs + n*stride_t; + + file << Nx_out << "\n"; + file << xmin << " " << xmax << "\n"; + + for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) + { + double val = -eval<1>(xmin, c); + file << val; + file << "\n"; + } + file.close(); +} + + + inline void NuFISolver::run() { std::cout << "Building E_sline\n\n"; @@ -220,7 +277,8 @@ inline void NuFISolver::run() { std::cout << "Saving results... \n\n"; save_ftilda(it, coeffs.get(), 128, 128, "results/ftilda_" + std::to_string(it) + ".dat"); - poisson.output_results(it); + save_rho(it, coeffs.get(), 128, "results/rho_" + std::to_string(it) + ".dat"); + save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat"); } } diff --git a/poisson_non_periodic.hpp b/poisson_non_periodic.hpp new file mode 100644 index 0000000..cb30a8a --- /dev/null +++ b/poisson_non_periodic.hpp @@ -0,0 +1,257 @@ +#ifndef POISSON_NON_PERIODIC_HPP +#define POISSON_NON_PERIODIC_HPP + +#include "parameters.hpp" +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace dealii; + + +template +class Poisson_non_periodic +{ +public: + Poisson_non_periodic (); + + void run(); + + void initialize(); + void solve_step(); + + void set_rhs_function(std::unique_ptr> rhs_function); + + const Vector &get_solution() const { return solution; } + const DoFHandler &get_dof_handler() const { return dof_handler; } + + std::vector sample_electric_field(const Poisson_non_periodic &problem, // sampling to save as spline + unsigned int Nx, + double x_min, + double x_max); + void output_results(unsigned int n); +private: + + void make_grid(); + void setup_system(); + void assemble_system(); + void solve(); + void output_results() const; + + Triangulation<1> triangulation; + const FE_Q<1> fe; + DoFHandler<1> dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + std::unique_ptr> rhs_function; + +}; + +template +Poisson_non_periodic::Poisson_non_periodic() + : fe(/* polynomial degree = */ 1) + , dof_handler(triangulation) +{} + +template +void Poisson_non_periodic::set_rhs_function(std::unique_ptr> rhs) +{ + rhs_function = std::move(rhs); +} + + +template +void Poisson_non_periodic::make_grid() +{ + Point x0 = Parameters::X_DOMAIN_RIGHT; + Point x1 = Parameters::X_DOMAIN_RIGHT; + GridGenerator::hyper_rectangle(triangulation, x0, x1); + triangulation.refine_global(Parameters::GLOBAL_REFINEMENT); + + std::cout << "Number of active cells: " << triangulation.n_active_cells() + << std::endl; +} + + + +template +void Poisson_non_periodic::setup_system() +{ + dof_handler.distribute_dofs(fe); + std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); +} + + +template +void Poisson_non_periodic::assemble_system() +{ + const QGauss<1> quadrature_formula(fe.degree + 1); + FEValues<1> fe_values(fe, + quadrature_formula, + update_values | update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + + cell_matrix = 0; + cell_rhs = 0; + + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + + const double rho = rhs_function->value(fe_values.quadrature_point(q_index)); + + for (const unsigned int i : fe_values.dof_indices()) + for (const unsigned int j : fe_values.dof_indices()) + cell_matrix(i, j) += + (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q) + fe_values.shape_grad(j, q_index) * // grad phi_j(x_q) + fe_values.JxW(q_index)); // dx + + for (const unsigned int i : fe_values.dof_indices()) + cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q) + rho * // f(x_q) + fe_values.JxW(q_index)); // dx + } + cell->get_dof_indices(local_dof_indices); + + for (const unsigned int i : fe_values.dof_indices()) + for (const unsigned int j : fe_values.dof_indices()) + system_matrix.add(local_dof_indices[i], + local_dof_indices[j], + cell_matrix(i, j)); + + for (const unsigned int i : fe_values.dof_indices()) + system_rhs(local_dof_indices[i]) += cell_rhs(i); + } + + + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + types::boundary_id(0), + Functions::ZeroFunction<1>(), + boundary_values); + MatrixTools::apply_boundary_values(boundary_values, + system_matrix, + solution, + system_rhs); +} + + +template +void Poisson_non_periodic::solve() +{ + SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm()); + SolverCG> solver(solver_control); + solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); + + std::cout << solver_control.last_step() + << " CG iterations needed to obtain convergence." << std::endl; +} + +template +void Poisson_non_periodic::output_results(unsigned int n) +{ + + // --- extract DoF coordinates --- + std::vector> support_points(dof_handler.n_dofs()); + Vector x_coordinate(dof_handler.n_dofs()); + + for (unsigned int i = 0; i < support_points.size(); ++i) + x_coordinate[i] = support_points[i][0]; // x-component in 1D + + //---- Output density ---- + ChargeDensity rho(Parameters::EPS, Parameters::WAVE_NR, Parameters::NV); + + DataOut data_out_rho; + data_out_rho.attach_dof_handler(dof_handler); + + Vector density(solution.size()); + VectorTools::interpolate(dof_handler, rho, density); + + data_out_rho.add_data_vector(density, "density"); + data_out_rho.add_data_vector(x_coordinate, "x_coordinate"); + + data_out_rho.build_patches(); + + std::ofstream out1("results/density_" + std::to_string(n) + ".vtk"); + data_out_rho.write_vtk(out1); + + //---- Output electric field & potential ---- + DataOut data_out_E; + data_out_E.attach_dof_handler(dof_handler); + + ElectricFieldPostprocessor electric_field; + Vector dummy(solution.size() * dim); + + data_out_E.add_data_vector(solution, "potential"); + data_out_E.add_data_vector(solution, electric_field); + data_out_E.add_data_vector(x_coordinate, "x_coordinate"); + + data_out_E.build_patches(); + + std::ofstream out2("results/electric_field_"+ std::to_string(n)+".vtk"); + data_out_E.write_vtk(out2); +} + +template +void Poisson_non_periodic::initialize() +{ + make_mesh(); // build grid + setup_system(); // distribute DoFs and matrices +} + +template +void Poisson_non_periodic::solve_step() +{ + assemble_system(); + solve(); +} + +#endif diff --git a/poisson_problem.hpp b/poisson_problem.hpp index 81c494e..2f4c12b 100644 --- a/poisson_problem.hpp +++ b/poisson_problem.hpp @@ -39,7 +39,6 @@ #include #include "parameters.hpp" -#include "fields.hpp" using namespace dealii; @@ -64,8 +63,6 @@ public: unsigned int Nx, double x_min, double x_max); - void output_results(unsigned int n); - private: void create_mesh(); @@ -284,57 +281,6 @@ void PoissonProblem::solve() constraints.distribute(solution); } - -template -void PoissonProblem::output_results(unsigned int n) -{ - - // --- extract DoF coordinates --- - std::vector> support_points(dof_handler.n_dofs()); - - DoFTools::map_dofs_to_support_points(mapping, - dof_handler, - support_points); - - Vector x_coordinate(dof_handler.n_dofs()); - - for (unsigned int i = 0; i < support_points.size(); ++i) - x_coordinate[i] = support_points[i][0]; // x-component in 1D - - //---- Output density ---- - ChargeDensity rho(Parameters::EPS, Parameters::WAVE_NR, Parameters::NV); - - DataOut data_out_rho; - data_out_rho.attach_dof_handler(dof_handler); - - Vector density(solution.size()); - VectorTools::interpolate(dof_handler, rho, density); - - data_out_rho.add_data_vector(density, "density"); - data_out_rho.add_data_vector(x_coordinate, "x_coordinate"); - - data_out_rho.build_patches(); - - std::ofstream out1("results/density_" + std::to_string(n) + ".vtk"); - data_out_rho.write_vtk(out1); - - //---- Output electric field & potential ---- - DataOut data_out_E; - data_out_E.attach_dof_handler(dof_handler); - - ElectricFieldPostprocessor electric_field; - Vector dummy(solution.size() * dim); - - data_out_E.add_data_vector(solution, "potential"); - data_out_E.add_data_vector(solution, electric_field); - data_out_E.add_data_vector(x_coordinate, "x_coordinate"); - - data_out_E.build_patches(); - - std::ofstream out2("results/electric_field_"+ std::to_string(n)+".vtk"); - data_out_E.write_vtk(out2); -} - template void PoissonProblem::initialize() { @@ -345,8 +291,6 @@ void PoissonProblem::initialize() template void PoissonProblem::solve_step() { - system_matrix = 0; - system_rhs = 0; assemble_system(); solve(); } @@ -360,7 +304,6 @@ void PoissonProblem::run() setup_system(); assemble_system(); solve(); - output_results(); } #endif