From dbb3b21a67ec439da6da6c400fafbd19aec189dd Mon Sep 17 00:00:00 2001 From: "V. Ferreira" Date: Thu, 26 Feb 2026 02:48:25 +0100 Subject: [PATCH] reorganized files --- density.vtk | 2 +- electric_field.vtk | 2 +- fields.hpp | 61 +++++++ nufi_poisson.cc | 401 +++----------------------------------------- parameters.hpp | 23 +++ poisson_problem.hpp | 298 ++++++++++++++++++++++++++++++++ 6 files changed, 403 insertions(+), 384 deletions(-) create mode 100644 fields.hpp create mode 100644 parameters.hpp create mode 100644 poisson_problem.hpp diff --git a/density.vtk b/density.vtk index 97aae28..8d4520a 100644 --- a/density.vtk +++ b/density.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 3.0 -#This file was generated by the deal.II library on 2026/2/26 at 1:54:49 +#This file was generated by the deal.II library on 2026/2/26 at 2:47:45 ASCII DATASET UNSTRUCTURED_GRID diff --git a/electric_field.vtk b/electric_field.vtk index f3ae802..7e853a5 100644 --- a/electric_field.vtk +++ b/electric_field.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 3.0 -#This file was generated by the deal.II library on 2026/2/26 at 1:54:49 +#This file was generated by the deal.II library on 2026/2/26 at 2:47:45 ASCII DATASET UNSTRUCTURED_GRID diff --git a/fields.hpp b/fields.hpp new file mode 100644 index 0000000..4d8ee0b --- /dev/null +++ b/fields.hpp @@ -0,0 +1,61 @@ +#ifndef FIELDS_HPP +#define FIELDS_HPP + +#include +#include +#include "parameters.hpp" + +using namespace dealii; + +inline double f0(const double x, + const double v, + const double eps = Parameters::EPS, + const double k = Parameters::WAVE_NR) +{ + const double prefactor = (1.0 + eps * std::cos(k*x)); + const double gaussian = + (v*v / std::sqrt(2.0 * M_PI)) * + std::exp(-0.5 * v*v); + + return prefactor * gaussian; +} + +inline double compute_rho(const double x, + const unsigned int Nv = Parameters::NV) +{ + const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; + + double integral = 0.0; + + for (unsigned int i = 0; i < Nv; ++i) + { + const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv; + integral += f0(x, v) * dv; + } + + return 1.0 - integral; +} + +template +class ChargeDensity : public Function +{ +public: + ChargeDensity(double eps, + double k, + unsigned int Nv) + : Function(1), eps(eps), k(k), Nv(Nv) {} + + virtual double value(const Point &p, + [[maybe_unused]] const unsigned int component = 0) const override + { + return compute_rho(p[0], Nv); + } + +private: + const double eps; + const double k; + const unsigned int Nv; +}; + + +#endif diff --git a/nufi_poisson.cc b/nufi_poisson.cc index e37e707..2d1dfae 100644 --- a/nufi_poisson.cc +++ b/nufi_poisson.cc @@ -1,397 +1,34 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include - -#include -#include - -#include -#include - -#include -#include #include -using namespace dealii; - -// =-=-=-=-=-= Parameter choice =-=-=-=-=-=-= - -// Domain dimension -constexpr unsigned int DIMENSION = 1; - -// Domain boundaries -constexpr double X_DOMAIN_LEFT = 0.0; -constexpr double X_DOMAIN_RIGHT = 12.0; - -constexpr double V_DOMAIN_LEFT = -6.0; -constexpr double V_DOMAIN_RIGHT = 6.0; -constexpr unsigned int NV = 1e3; // used only to evaluate rho(x), independent of deal.ii grid - -// Global refinement level -constexpr unsigned int GLOBAL_REFINEMENT = 7; - -// Polynomial degree -constexpr unsigned int FE_DEGREE = 3; - -// f0 parameters -constexpr double EPS = 0.01; -constexpr double WAVE_NR = 0.5; - - -// =-=-=-=-= f_0(x,v) =-=-=-=-= - -double f0(const double x, - const double v, - const double eps=EPS, - const double k=WAVE_NR) -{ - const double prefactor = (1.0 + eps * std::cos(k*x)); - const double gaussian = (v*v / std::sqrt((2.0 * M_PI))) - * std::exp(- 0.5 * v*v); - - return prefactor*gaussian; -} - -// =-=-=-=-= Compute rho(x) =-=-=-=-= - -double compute_rho(const double x, const unsigned int Nv=NV) -{ - const double v_min = V_DOMAIN_LEFT; - const double v_max = V_DOMAIN_RIGHT; - - const double dv = std::abs(v_min - v_max)/Nv; - - double integral = 0.0; - - for (unsigned int i=0; i -class ChargeDensity : public Function -{ -public: - ChargeDensity(double eps, - double k, - unsigned int Nv) - : Function(1), eps(eps), k(k), Nv(Nv) {} - - virtual double value(const Point &p, - [[maybe_unused]] const unsigned int component = 0) const override - { - return compute_rho(p[0], Nv); - } - -private: - const double eps; - const double k; - const unsigned int Nv; -}; - - -// =-=-=-=-= Poisson Solver =-=-=-=-= - -template -class PoissonProblem -{ -public: - PoissonProblem(unsigned int degree, unsigned int Nv); - void run(); - - void set_Nv(unsigned int new_Nv); - -private: - void create_mesh(); - void setup_system(); - void assemble_system(); - void solve(); - void output_results() const; - - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; - - AffineConstraints constraints; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; // phi - Vector system_rhs; - - MappingQ mapping; - - unsigned int Nv; -}; - - -template -void PoissonProblem::set_Nv(unsigned int new_Nv) -{ - Nv = new_Nv; -} - -template -PoissonProblem::PoissonProblem(unsigned int degree, unsigned int Nv) - : fe(degree) - , dof_handler(triangulation) - , mapping(degree) - , Nv(Nv) -{} - - -// =-=-=-=-= Make Grid =-=-=-=-= - -template -void PoissonProblem::create_mesh() -{ - GridGenerator::hyper_cube(triangulation, - X_DOMAIN_LEFT, - X_DOMAIN_RIGHT); - - // Make x-dim boundaries periodic - Tensor<1, dim> offset; - std::vector::cell_iterator>> periodicity_vector; - - GridTools::collect_periodic_faces(triangulation, - 0, - 1, - 0, - periodicity_vector, - offset); - - triangulation.add_periodicity(periodicity_vector); - - triangulation.refine_global(GLOBAL_REFINEMENT); -} - -template -void PoissonProblem::setup_system() -{ - dof_handler.distribute_dofs(fe); - - constraints.clear(); - DoFTools::make_hanging_node_constraints(dof_handler, constraints); - - // 'boundary' condition phi(x_0) = 0 - constraints.add_line(0); - constraints.set_inhomogeneity(0, 0.0); - - constraints.close(); - - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); - sparsity_pattern.copy_from(dsp); - - system_matrix.reinit(sparsity_pattern); - solution.reinit(dof_handler.n_dofs()); - system_rhs.reinit(dof_handler.n_dofs()); -} - -// =-=-=-=-= E_field = -dPhi/dx =-=-=-=-= - -template -class ElectricFieldPostprocessor : public DataPostprocessorVector -{ -public: - ElectricFieldPostprocessor() - : DataPostprocessorVector("electric_field", update_gradients) - {} - - virtual void evaluate_scalar_field( - const DataPostprocessorInputs::Scalar &input_data, - std::vector> &computed_quantities) const override - { - AssertDimension(input_data.solution_gradients.size(), - computed_quantities.size()); - - for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p) - { - AssertDimension(computed_quantities[p].size(), dim); - for (unsigned int d = 0; d < dim; ++d) - computed_quantities[p][d] = -input_data.solution_gradients[p][d]; - } - } -}; - - -// =-=-=-=-= Poisson equation solver =-=-=-=-= - -template -void PoissonProblem::assemble_system() -{ - QGauss quadrature_formula(fe.degree + 1); - FEValues fe_values(fe, quadrature_formula, - update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - Vector cell_rhs(dofs_per_cell); - std::vector local_dof_indices(dofs_per_cell); - - ChargeDensity rhs_function(EPS, WAVE_NR, Nv); - - for (const auto &cell : dof_handler.active_cell_iterators()) - { - fe_values.reinit(cell); - cell_matrix = 0; - cell_rhs = 0; - - for (unsigned int q = 0; q < n_q_points; ++q) - { - const double rho = rhs_function.value(fe_values.quadrature_point(q)); - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < dofs_per_cell; ++j) - cell_matrix(i, j) += - fe_values.shape_grad(i, q) * - fe_values.shape_grad(j, q) * - fe_values.JxW(q); - - cell_rhs(i) += - fe_values.shape_value(i, q) * - rho * - fe_values.JxW(q); - } - } - - cell->get_dof_indices(local_dof_indices); - constraints.distribute_local_to_global(cell_matrix, - cell_rhs, - local_dof_indices, - system_matrix, - system_rhs); - } -} - - -template -void PoissonProblem::solve() -{ - SolverControl solver_control(1000, 1e-12); - SolverCG> solver(solver_control); - - PreconditionSSOR> preconditioner; - preconditioner.initialize(system_matrix, 1.2); - - solver.solve(system_matrix, solution, system_rhs, preconditioner); - constraints.distribute(solution); -} - - -template -void PoissonProblem::output_results() const -{ - - // --- 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(EPS, WAVE_NR, Nv); - Vector density(triangulation.n_active_cells()); - - DataOut data_out_rho; - data_out_rho.attach_dof_handler(dof_handler); - - Vector density_nodal(solution.size()); - VectorTools::interpolate(dof_handler, rho, density_nodal); - - data_out_rho.add_data_vector(density_nodal, "density"); - data_out_rho.add_data_vector(x_coordinate, "x_coordinate"); - - data_out_rho.build_patches(); - - std::ofstream out1("density.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("electric_field.vtk"); - data_out_E.write_vtk(out2); -} - - -template -void PoissonProblem::run() -{ - create_mesh(); - setup_system(); - assemble_system(); - solve(); - output_results(); -} - +#include "parameters.hpp" +#include "poisson_problem.hpp" int main() { try { - PoissonProblem poisson_problem(FE_DEGREE, NV); - poisson_problem.set_Nv(NV); + PoissonProblem poisson_problem(Parameters::FE_DEGREE, + Parameters::NV); + poisson_problem.run(); } catch (std::exception &exc) { - std::cerr << exc.what() << std::endl; + std::cerr << std::endl + << "Exception: " + << std::endl + << exc.what() + << std::endl; return 1; } + catch (...) + { + std::cerr << std::endl + << "Unknown exception!" + << std::endl; + return 1; + } + return 0; } + diff --git a/parameters.hpp b/parameters.hpp new file mode 100644 index 0000000..b7c7add --- /dev/null +++ b/parameters.hpp @@ -0,0 +1,23 @@ +#ifndef PARAMETERS_HPP +#define PARAMETERS_HPP + +namespace Parameters +{ + constexpr unsigned int DIMENSION = 1; + + constexpr double X_DOMAIN_LEFT = 0.0; + constexpr double X_DOMAIN_RIGHT = 12.0; + + constexpr double V_DOMAIN_LEFT = -6.0; + constexpr double V_DOMAIN_RIGHT = 6.0; + + constexpr unsigned int NV = 1000; + + constexpr unsigned int GLOBAL_REFINEMENT = 7; + constexpr unsigned int FE_DEGREE = 3; + + constexpr double EPS = 0.01; + constexpr double WAVE_NR = 0.5; +} + +#endif diff --git a/poisson_problem.hpp b/poisson_problem.hpp new file mode 100644 index 0000000..9c05ab7 --- /dev/null +++ b/poisson_problem.hpp @@ -0,0 +1,298 @@ +#ifndef POISSON_PROBLEM_HPP +#define POISSON_PROBLEM_HPP + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "parameters.hpp" +#include "fields.hpp" + +using namespace dealii; + +// =-=-=-=-= Poisson Solver =-=-=-=-= + +template +class PoissonProblem +{ +public: + PoissonProblem(unsigned int degree, unsigned int Nv); + void run(); + + void set_Nv(unsigned int new_Nv); + +private: + void create_mesh(); + void setup_system(); + void assemble_system(); + void solve(); + void output_results() const; + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; // phi + Vector system_rhs; + + MappingQ mapping; + + unsigned int Nv; +}; + + +template +void PoissonProblem::set_Nv(unsigned int new_Nv) +{ + Nv = new_Nv; +} + +template +PoissonProblem::PoissonProblem(unsigned int degree, unsigned int Nv) + : fe(degree) + , dof_handler(triangulation) + , mapping(degree) + , Nv(Nv) +{} + + +// =-=-=-=-= Make Grid =-=-=-=-= + +template +void PoissonProblem::create_mesh() +{ + GridGenerator::hyper_cube(triangulation, + Parameters::X_DOMAIN_LEFT, + Parameters::X_DOMAIN_RIGHT); + + // Make x-dim boundaries periodic + Tensor<1, dim> offset; + std::vector::cell_iterator>> periodicity_vector; + + GridTools::collect_periodic_faces(triangulation, + 0, + 1, + 0, + periodicity_vector, + offset); + + triangulation.add_periodicity(periodicity_vector); + + triangulation.refine_global(Parameters::GLOBAL_REFINEMENT); +} + +template +void PoissonProblem::setup_system() +{ + dof_handler.distribute_dofs(fe); + + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + // 'boundary' condition phi(x_0) = 0 + constraints.add_line(0); + constraints.set_inhomogeneity(0, 0.0); + + constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); +} + +// =-=-=-=-= E_field = -dPhi/dx =-=-=-=-= + +template +class ElectricFieldPostprocessor : public DataPostprocessorVector +{ +public: + ElectricFieldPostprocessor() + : DataPostprocessorVector("electric_field", update_gradients) + {} + + virtual void evaluate_scalar_field( + const DataPostprocessorInputs::Scalar &input_data, + std::vector> &computed_quantities) const override + { + AssertDimension(input_data.solution_gradients.size(), + computed_quantities.size()); + + for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p) + { + AssertDimension(computed_quantities[p].size(), dim); + for (unsigned int d = 0; d < dim; ++d) + computed_quantities[p][d] = -input_data.solution_gradients[p][d]; + } + } +}; + + +// =-=-=-=-= Poisson equation solver =-=-=-=-= + +template +void PoissonProblem::assemble_system() +{ + QGauss quadrature_formula(fe.degree + 1); + FEValues fe_values(fe, quadrature_formula, + update_values | + update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + + ChargeDensity rhs_function(Parameters::EPS, Parameters::WAVE_NR, Nv); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double rho = rhs_function.value(fe_values.quadrature_point(q)); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += + fe_values.shape_grad(i, q) * + fe_values.shape_grad(j, q) * + fe_values.JxW(q); + + cell_rhs(i) += + fe_values.shape_value(i, q) * + rho * + fe_values.JxW(q); + } + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } +} + + +template +void PoissonProblem::solve() +{ + SolverControl solver_control(1000, 1e-12); + SolverCG> solver(solver_control); + + PreconditionSSOR> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + solver.solve(system_matrix, solution, system_rhs, preconditioner); + constraints.distribute(solution); +} + + +template +void PoissonProblem::output_results() const +{ + + // --- 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, 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("density.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("electric_field.vtk"); + data_out_E.write_vtk(out2); +} + + +template +void PoissonProblem::run() +{ + create_mesh(); + setup_system(); + assemble_system(); + solve(); + output_results(); +} + +#endif