From 61599fbce01e0e661aa3bb00609c3e5c5fabb2ea Mon Sep 17 00:00:00 2001 From: "Vasco C. B. Ferreira" Date: Wed, 11 Mar 2026 21:01:43 +0100 Subject: [PATCH] advanced with new splines, updated readme todo --- README.md | 4 +- fields.hpp | 109 ++++++++++++++++++++++++++++++++++++++++++++++++ nufi_solver.hpp | 27 ++++++------ parameters.hpp | 5 ++- splines.hpp | 90 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 219 insertions(+), 16 deletions(-) create mode 100644 splines.hpp diff --git a/README.md b/README.md index fe5057c..0aa4092 100644 --- a/README.md +++ b/README.md @@ -7,4 +7,6 @@ This simulation of the Vlasov-Poisson system in 1x1v dimensions uses --- Todo: - Correct dealii solver, check ftilda results to see whats happening -- ... + - implement least squares thingy + - marry it to fields interpolation + - make nufi solver use new spline interpolation and evaluation diff --git a/fields.hpp b/fields.hpp index 04d8cf9..d7df45a 100644 --- a/fields.hpp +++ b/fields.hpp @@ -35,6 +35,115 @@ inline double compute_rho(const double x, return 1.0 - integral; } + +template +real eval( real x, const real *coeffs) noexcept +{ + using std::floor; + + // Shift to a box that starts at 0. + x -= Parameters::X_DOMAIN_LEFT; + + // Get "periodic position" in box at origin. + x = x - Parameters::LX * floor( x/Parameters::LX ); + + // Knot number + real x_knot = floor( x/Parameters::SPLINE_DX); + + size_t ii = static_cast(x_knot); + + // Convert x to reference coordinates. + x = x/Parameters::SPLINE_DX - x_knot; + + // Scale according to derivative. + real factor = 1; + for ( size_t i = 0; i < dx; ++i ) factor *= 1/Parameters::SPLINE_DX; + + return factor*splines1d::eval( x, coeffs + ii ); +} + +template +void interpolate( real *coeffs, const real *values) // Least Squares needs to be made +{ + std::unique_ptr tmp { new real[ Parameters::SPLINE_NX ] }; + + for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i ) + tmp[ i ] = coeffs[ i ]; + + struct mat_t // STRUCT AND CONFIG NEEDS TO BE REVIEWED + { + const config_t &config; + real N[ order ]; + + mat_t( const config_t &conf ): config { conf } + { + splines1d::N(0,N); + } + + void operator()( const real *in, real *out ) const + { + #pragma omp parallel for + for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i ) + { + real result = 0; + if ( i + order <= Parameters::SPLINE_NX ) + { + for ( size_t ii = 0; ii < order; ++ii ) + result += N[ii] * in[ i + ii ]; + } + else + { + for ( size_t ii = 0; ii < order; ++ii ) + result += N[ii]*in[ (i+ii) % Parameters::SPLINE_NX]; + } + out[ i ] = result; + } + } + }; + + struct transposed_mat_t // STRUCT AND CONFIG NEEDS TO BE REVIEWED + { + const config_t &config; + real N[ order ]; + + transposed_mat_t( const config_t &conf ): config { conf } + { + splines1d::N(0,N); + } + + void operator()( const real *in, real *out ) const + { + for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i ) + out[ i ] = 0; + + for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i ) + { + if ( i + order <= Parameters::SPLINE_NX ) + { + for ( size_t ii = 0; ii < order; ++ii ) + out[ i + ii ] += N[ii] * in[ i ]; + } + else + { + for ( size_t ii = 0; ii < order; ++ii ) + out[ (i+ii) % Parameters::SPLINE_NX ] += N[ii]*in[ i ]; + } + } + } + }; + + mat_t M { config }; transposed_mat_t Mt { config }; + lsmr_options opt; opt.silent = true; + lsmr( config.Nx, config.Nx, M, Mt, values, tmp.get(), opt ); + + if ( opt.iter == opt.max_iter ) + std::cerr << "Warning. LSMR did not converge.\n"; + + for ( size_t i = 0; i < Parameters::SPLINE_NX + order - 1; ++i ) + coeffs[ i ] = tmp[ i % Parameters::SPLINE_NX ]; +} + + template class ChargeDensity : public Function // only uses f0 { diff --git a/nufi_solver.hpp b/nufi_solver.hpp index c8b7ae2..adb5078 100644 --- a/nufi_solver.hpp +++ b/nufi_solver.hpp @@ -15,7 +15,8 @@ #include "parameters.hpp" #include "poisson_problem.hpp" #include "fields.hpp" // holds f0(x,v), and compute_rho(x) -#include "spline_field.hpp" +#include "spline_field.hpp" // old GPT splines +#include "splines.hpp" //new splines using namespace dealii; @@ -25,9 +26,9 @@ public: NuFISolver(); void run(); - double eval_rho(unsigned int n, double x, const UniformSpline1D& E_spline, unsigned int Nv = Parameters::NV); - double eval_ftilda(unsigned int n, double x, double u, const UniformSpline1D& E_spline); - void save_ftilda(unsigned int n, const UniformSpline1D& E_spline, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename); + double eval_rho(unsigned int n, double x, const std::vector E_coeffs, unsigned int Nv = Parameters::NV); + double eval_ftilda(unsigned int n, double x, double u, const std::vector E_coeffs); + void save_ftilda(unsigned int n, const std::vector E_coeffs, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename); private: @@ -50,7 +51,7 @@ private: inline double NuFISolver::eval_ftilda(unsigned int n, double x, double u, - const UniformSpline1D& E_spline) + const std::vector E_coeffs) { double Lu = std::abs(Parameters::V_DOMAIN_LEFT - Parameters::V_DOMAIN_RIGHT); @@ -78,7 +79,7 @@ inline double NuFISolver::eval_ftilda(unsigned int n, inline double NuFISolver::eval_rho(const unsigned int n, const double x, - const UniformSpline1D& E_spline, + const std::vector E_coeffs, const unsigned int Nv) { const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; @@ -97,25 +98,25 @@ inline double NuFISolver::eval_rho(const unsigned int n, class ChargeDensity_NuFI : public Function<1> { public: - ChargeDensity_NuFI(NuFISolver &solver, size_t n, const UniformSpline1D &E_spline) - : solver(solver), n(n), E_spline(E_spline) {} + ChargeDensity_NuFI(NuFISolver &solver, size_t n, const std::vector E_coeffs) + : solver(solver), n(n), E_coeffs(E_coeffs) {} virtual double value(const Point<1> &p, [[maybe_unused]] const unsigned int component = 0) const override { double x = p[0]; - return solver.eval_rho(n, x, E_spline); + return solver.eval_rho(n, x, E_coeffs); } private: NuFISolver &solver; size_t n; - const UniformSpline1D &E_spline; + const std::vector E_coeffs; }; inline void NuFISolver::save_ftilda(unsigned int n, - const UniformSpline1D& E_spline, + const std::vector E_coeffs, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename) @@ -176,7 +177,7 @@ inline void NuFISolver::run() E_grid[i] = 0; } - UniformSpline1D E_spline(E_grid, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT); + std::vector E_coeffs(E_grid, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT); // Needs correction for (unsigned int it = 0; it < Nt; ++it) { @@ -212,7 +213,7 @@ inline void NuFISolver::run() E_grid = poisson.sample_electric_field(poisson, Nx, 0.0, Lx); - E_spline = UniformSpline1D(E_grid, 0.0, Lx); + E_spline = std::vector E_coeffs; // needs correction } std::cout << "NuFI simulation finished.\n"; diff --git a/parameters.hpp b/parameters.hpp index 8fe7c69..1860973 100644 --- a/parameters.hpp +++ b/parameters.hpp @@ -28,10 +28,11 @@ namespace Parameters constexpr unsigned int TMAX = 10; //spline options - constexpr int SPLINE_NX = 256; + constexpr int SPLINE_NX = 512; + constexpr double SPLINE_DX = LX/SPLINE_NX; //Plotting options - constexpr int PLOT_FREQUENCY = 3; + constexpr int PLOT_FREQUENCY = 2; } #endif diff --git a/splines.hpp b/splines.hpp new file mode 100644 index 0000000..eb488d4 --- /dev/null +++ b/splines.hpp @@ -0,0 +1,90 @@ +#ifndef SPLINES_HPP +#define SPLINES_HPP + +#include + +namespace splines1d +{ + +template +constexpr real faculty( size_t n ) noexcept +{ + return (n > 1) ? real(n)*faculty(n-1) : real(1); +} + +template +void N( real x, real *result, size_t stride = 1 ) noexcept +{ + static_assert( order > 0, "Splines must have order greater than zero." ); + constexpr int n { order }; + constexpr int d { derivative }; + + if ( derivative >= order ) + for ( size_t i = 0; i < order; ++i ) + result[ i*stride ] = 0; + + if ( n == 1 ) + { + *result = 1; + return; + } + + real v[n]; v[n-1] = 1; + for ( int k = 1; k < n - d; ++k ) + { + v[n-k-1] = (1-x)*v[n-k]; + + for ( int i = 1-k; i < 0; ++i ) + v[n-1+i] = (x-i)*v[n-1+i] + (k+1+i-x)*v[n+i]; + + v[n-1] *= x; + } + + // Differentiate if necessary. + for ( size_t j = derivative; j-- > 0; ) + { + v[j] = -v[j+1]; + for ( size_t i = j + 1; i < order - 1; ++i ) + v[i] = v[i] - v[i+1]; + } + + constexpr real factor = real(1) / faculty(order-derivative-1); + for ( size_t i = 0; i < order; ++i ) + result[i*stride] = v[i]*factor; +} + +template +real eval( real x, const real *coefficients, size_t stride = 1 ) noexcept +{ + static_assert( order > 0, "Splines must have order greater than zero." ); + static_assert( order > derivative, "Too high derivative requested." ); + constexpr size_t n { order }; + constexpr size_t d { derivative }; + + if ( d >= n ) return 0; + if ( n == 1 ) return *coefficients; + + // Gather coefficients. + real c[ order ]; + for ( size_t j = 0; j < order; ++j ) + c[j] = coefficients[ stride * j ]; + + // Differentiate if necessary. + for ( size_t j = 1; j <= d; ++j ) + for ( size_t i = n; i-- > j; ) + c[i] = c[i] - c[i-1]; + + // Evaluate using de Boor’s algorithm. + for ( size_t j = 1; j < n-d; ++j ) + for ( size_t i = n-d; i-- > j; ) + c[d+i] = (x+n-d-1-i)*c[d+i] + (i-j+1-x)*c[d+i-1]; + + constexpr real factor = real(1) / faculty(order-derivative-1); + return factor*c[n-1]; +} + +} + +#endif + +