vcbferreira/NuFI_deal.ii
1
0
Fork 0
mirror of https://codeberg.org/vcbferreira/NuFI_deal.ii synced 2026-08-12 14:33:18 +02:00
Code Issues Packages Projects Releases Wiki Activity

field history added

Browse Source
This commit is contained in:
Vasco C. B. Ferreira
2026-06-29 11:41:20 +02:00
parent 77392e4731
commit 9f6cca2aa2
7 changed files with 38 additions and 276 deletions
Download Patch File Download Diff File Expand all files Collapse all files
libnufi_lib.a
BIN
View File
Binary file not shown.
nufi/fields.h
+4 -71
View File
@@ -2,7 +2,7 @@
#define FIELDS_H
#include "nufi/parameters.h"
#include "poisson_problem.h"
#include "nufi/poisson_problem.h"
#include <cmath>
#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
@@ -68,27 +68,14 @@ inline double f0(const double x, const double 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;
}
// wrapper for eval_point() { VectorTools::point_values() }
inline double eval(double x, const PoissonProblem<1> &poisson,
const Vector<double> &solution) noexcept {
x -= Parameters::X_DOMAIN_LEFT;
x = x - Parameters::LX * std::floor(x * Parameters::LX_INV);
x = x - Parameters::LX * std::floor(x * Parameters::LX_INV); // in domain
return eval_point<1>(poisson.get_mapping(), poisson.get_dof_handler(),
solution, Point<1>(x));
}
@@ -122,58 +109,4 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
return integral * dx;
};
// class Gradient {
// public:
// Gradient(double xmin, double xmax, unsigned int Nx)
// : xmin_(xmin), xmax_(xmax), Nx_(Nx) {
// if (xmax_ <= xmin_) {
// throw std::invalid_argument("xmax must be greater than xmin");
// }
// }
//
// std::vector<double> compute(const std::vector<double> &values) const {
// size_t n = values.size();
// if (n < 2) {
// throw std::invalid_argument("Need at least 2 points");
// }
//
// std::vector<double> grad(n);
//
// double dx = (xmax_ - xmin_) / (n - 1);
// // periodic boundaries
// grad[0] = -(values[1] - values[n - 1]) / (2.0 * dx);
// grad[n - 1] = -(values[0] - values[n - 2]) / (2.0 * dx);
//
// for (size_t i = 1; i < n - 1; ++i) {
// grad[i] = -(values[i + 1] - values[i - 1]) / (2.0 * dx);
// }
//
// return grad;
// }
//
// private:
// double xmin_;
// double xmax_;
// [[maybe_unused]] unsigned int Nx_;
// };
template <int dim>
class ChargeDensity : public Function<dim> // only uses f0
{
public:
ChargeDensity(double eps, double k, unsigned int Nv)
: Function<dim>(1), eps(eps), k(k), Nv(Nv) {}
virtual double
value(const Point<dim> &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
nufi/nufi_solver.h
+4 -31
View File
@@ -5,7 +5,7 @@
#include <cmath>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/vector_tools.h>
#include <vector>
#include "nufi/fields.h" //dont remove
@@ -19,9 +19,10 @@ public:
NuFISolver();
void run();
double eval_rho(unsigned int n, double x, const PoissonProblem<1> &poisson,
double eval_rho(unsigned int n, const double x,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
unsigned int Nv = Parameters::NV) const;
const unsigned int Nv = Parameters::NV) const;
double eval_ftilda(unsigned int n, double x, double u,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history) const;
@@ -44,32 +45,4 @@ private:
PoissonProblem<1> poisson;
};
template <unsigned int dim> class ChargeDensity_NuFI : public Function<dim> {
public:
ChargeDensity_NuFI(const double *rho_values, unsigned int Nx)
: Function<dim>(), rho(rho_values), Nx(Nx) {}
virtual double
value(const Point<dim> &p,
[[maybe_unused]] const unsigned int component = 0) const override {
const double x = p[0];
// Map x -> grid index
const double L = Parameters::LX;
const double dx = L / (Nx - 1);
int i = static_cast<int>(std::floor((x - Parameters::X_DOMAIN_LEFT) / dx));
// periodic wrap
i = (i % Nx + Nx) % Nx;
return rho[i];
}
private:
const double *rho;
const unsigned int Nx;
};
#endif
nufi/parameters.h
+1 -1
View File
@@ -32,7 +32,7 @@ constexpr double WAVE_NR = 0.5;
constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
// NUFI options
constexpr double DT = 1. / 10.;
constexpr double DT = 1. / 8.;
constexpr unsigned int TMAX = 100;
// Plotting options
nufi/poisson_problem.h
+11 -158
View File
@@ -33,12 +33,12 @@
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/vector_tools_evaluate.h>
#include <deal.II/numerics/vector_tools_interpolate.h>
#include <functional>
#include <memory>
#include <string>
#include <utility>
@@ -58,11 +58,11 @@ public:
void solve_step();
void run();
void set_rhs_function(std::unique_ptr<Function<dim>> rhs_function);
void set_rhs_function(std::function<double(const Point<dim> &)> f);
const Vector<double> &get_solution() const { return solution; }
const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
const MappingQ<dim> &get_mapping() const { return mapping; }
const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
std::vector<double> sample_electric_field(double x_min, double x_max,
unsigned int Nx);
@@ -87,18 +87,17 @@ private:
Vector<double> solution; // phi
Vector<double> system_rhs;
std::unique_ptr<const Function<dim>> rhs_function;
std::function<double(const Point<dim> &)> rhs_function;
MappingQ<dim> mapping;
// std::unique_ptr<Functions::FEFieldFunction<dim>> fe_field_function;
};
// Utilities
template <int dim>
void PoissonProblem<dim>::set_rhs_function(std::unique_ptr<Function<dim>> rhs) {
rhs_function = std::move(rhs);
void PoissonProblem<dim>::set_rhs_function(
std::function<double(const Point<dim> &)> f) {
rhs_function = std::move(f);
}
template <int dim>
@@ -170,71 +169,6 @@ PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
return values;
}
// // by GPT to re-re-re-check
// template <int dim> std::vector<double> eval_solution_on_points(
// const std::vector<Vector<double>> &solutions,
// const unsigned int n,
// const std::vector<Point<dim>> &points, // need to be in [x_min,
// x_max]. I think.... const std::vector<unsigned int> &cell_indices,
// const DoFHandler<dim> &dof_handler,
// const MappingQ<dim> &mapping)
// {
// AssertIndexRange(n, solutions.size());
// Assert(points.size() == cell_indices.size(),
// ExcMessage("points and cell_indices must have same size"));
//
// const Vector<double> &solution = solutions[n];
//
// std::vector<double> result(points.size());
//
// // Group points by cell (required for FEPointEvaluation efficiency)
// std::map<unsigned int, std::vector<unsigned int>> cell_to_point_ids;
//
// for (unsigned int i = 0; i < points.size(); ++i)
// cell_to_point_ids[cell_indices[i]].push_back(i);
//
// FEPointEvaluation<1, dim> evaluator(mapping,
// dof_handler.get_fe(),
// update_values);
//
// std::vector<Point<dim>> cell_points;
// Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
//
// for (const auto &entry : cell_to_point_ids)
// {
// const unsigned int cell_id = entry.first;
// const auto &point_ids = entry.second;
//
// // these two lines bellow assume some order not sure how or why
// auto cell = dof_handler.begin_active();
// std::advance(cell, cell_id);
//
// // extract points belonging to this cell
// cell_points.clear();
// cell_points.reserve(point_ids.size());
//
// for (unsigned int id : point_ids)
// cell_points.push_back(points[id]);
//
// std::vector<types::global_dof_index>
// indices(dof_handler.get_fe().n_dofs_per_cell());
// cell->get_dof_indices(indices);
//
// for (unsigned int i=0;i<indices.size();++i)
// local_dofs[i] = solution[indices[i]];
//
// // initialize evaluator on this cell
// evaluator.reinit(cell, cell_points);
//
// evaluator.evaluate(local_dofs, EvaluationFlags::values);
//
// for (unsigned int k = 0; k < point_ids.size(); ++k)
// result[point_ids[k]] = evaluator.get_value(k);
// }
//
// return result;
// }
template <int dim>
double eval_point(const Mapping<dim> &mapping,
const DoFHandler<dim> &dof_handler,
@@ -299,90 +233,9 @@ template <int dim> void PoissonProblem<dim>::setup_system() {
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
// fe_field_function =
// std::make_unique<Functions::FEFieldFunction<dim>>(
// dof_handler, solution, mapping);
}
/* (Mine)
template <int dim>
void PoissonProblem<dim>::assemble_system()
{
system_matrix = 0;
system_rhs = 0;
QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> 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();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
Assert(rhs_function != nullptr, ExcMessage("RHS function not set"));
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
for (const auto q : fe_values.quadrature_point_indices())
{
const double rho = rhs_function->value(fe_values.quadrature_point(q));
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) * // grad phi_i(x_q)
fe_values.shape_grad(j, q) * // grad phi_j(x_q)
fe_values.JxW(q)); // dx
for (const unsigned int i : fe_values.dof_indices())
cell_rhs(i) += (fe_values.shape_value(i, q) * // phi_i(x_q)
rho * // f(x_q)
fe_values.JxW(q)); // dx
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
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<types::global_dof_index, double> 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);
}
*/
// Paul's, mine's above
// Paul
template <int dim> void PoissonProblem<dim>::assemble_system() {
Assert(system_matrix.m() == dof_handler.n_dofs(),
ExcMessage("Matrix not initialized correctly"));
@@ -400,8 +253,6 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
Assert(rhs_function != nullptr, ExcMessage("RHS function not set"));
for (const auto &cell : dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
@@ -409,7 +260,8 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
cell_rhs = 0;
for (const auto q : fe_values.quadrature_point_indices()) {
const double rho = rhs_function->value(fe_values.quadrature_point(q));
const double rho = rhs_function(
fe_values.quadrature_point(q)); // Eval rhs_function at q points
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
@@ -437,6 +289,7 @@ template <int dim> void PoissonProblem<dim>::solve() {
// preconditioner.initialize(system_matrix, 1.2);
// solver.solve(system_matrix, solution, system_rhs, preconditioner);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
constraints.distribute(solution);
}
rho_E.mp4
BIN
View File
Binary file not shown.
src/nufi_solver.cc
+16 -13
View File
@@ -6,7 +6,7 @@
#include <cstdlib>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/vector_tools.h>
#include <cstddef>
#include <iostream>
@@ -74,7 +74,7 @@ NuFISolver::eval_f(unsigned int n, double x, double u,
return f0(x, u);
}
double NuFISolver::eval_rho(const unsigned int n, const double x,
double NuFISolver::eval_rho(unsigned int n, const double x,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
const unsigned int Nv) const {
@@ -84,14 +84,9 @@ double NuFISolver::eval_rho(const unsigned int n, const double x,
double integral = 0.0;
if (n == 0) {
for (unsigned int i = 0; i < Nv; ++i)
integral += f0(x, v_min + i * dv);
} else {
#pragma omp parallel for reduction(+ : integral)
for (unsigned int i = 0; i < Nv; ++i)
integral += eval_ftilda(n, x, v_min + i * dv, poisson, phi_history);
}
return 1.0 - integral * dv;
}
@@ -115,6 +110,10 @@ void NuFISolver::run() {
double total_time = 0;
std::ofstream time_file("results/simulation_time.txt");
time_file << "# it step_time total_time" << "\n";
const double x_min = Parameters::X_DOMAIN_LEFT;
double dx = Parameters::CALC_DX;
for (unsigned int it = 0; it < Nt; ++it) {
stopwatch<double> timer;
@@ -127,8 +126,6 @@ void NuFISolver::run() {
// compute rho
double dx = Parameters::CALC_DX;
#pragma omp parallel for
for (size_t i = 0; i < Nx; i++) {
double x = Parameters::X_DOMAIN_LEFT + i * dx;
@@ -138,8 +135,13 @@ void NuFISolver::run() {
rho.get()[i] = ith_rho;
}
poisson.set_rhs_function(
std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Nx));
poisson.set_rhs_function([&rho, x_min, dx, Nx = Nx](const Point<1> &p) {
double x = p[0];
int i = static_cast<int>(std::floor((x - x_min) / dx));
i = (i % Nx + Nx) % Nx;
return rho.get()[i];
});
poisson.solve_step();
phi_history.push_back(poisson.get_solution());
@@ -151,6 +153,8 @@ void NuFISolver::run() {
double step_time = timer_elapsed - time_elapsed_before;
total_time += timer_elapsed;
time_file << it << " " << step_time << " " << total_time << "\n";
std::cout << "step made in " << step_time << " seconds\n\n";
if (it % Parameters::PLOT_FREQUENCY == 0) {
std::cout << "Saving results... ";
@@ -164,8 +168,7 @@ void NuFISolver::run() {
std::vector<double> E_x(Nx, 0.0);
#pragma omp parallel for
for (size_t ix = 0; ix < Nx; ++ix) {
E_x[ix] = -eval(Parameters::X_DOMAIN_LEFT + ix * dx, poisson,
phi_history[it]);
E_x[ix] = -eval(x_min + ix * dx, poisson, phi_history[it]);
}
save_space_vector(E_x, "field", it);