diff --git a/blas.hpp b/blas.hpp new file mode 100644 index 0000000..d61537f --- /dev/null +++ b/blas.hpp @@ -0,0 +1,54 @@ +#ifndef NUFI_BLAS_HPP +#define NUFI_BLAS_HPP + +#include + +/*! + * \brief Convenience wrappers for BLAS, with overloads for single and double + * precision. + */ +namespace blas +{ + +double dot( const size_t n, const double *x, size_t incx, + const double *y, size_t incy ); + +float dot( const size_t n, const float *x, size_t incx, + const float *y, size_t incy ); + +void axpy( size_t n, double alpha, const double *x, size_t incx, + double *y, size_t incy ); + +void axpy( size_t n, float alpha, const float *x, size_t incx, + float *y, size_t incy ); + + +void scal( size_t n, double alpha, double *x, size_t incx ); +void scal( size_t n, float alpha, float *x, size_t incx ); + +void copy( size_t n, const double *x, size_t incx, double *y, size_t incy ); +void copy( size_t n, const float *x, size_t incx, float *y, size_t incy ); + +void ger( const size_t M, const size_t N, const double alpha, + const double *X, const size_t incX, const double *Y, const size_t incY, + double *A, const size_t lda); + +void ger( const size_t M, const size_t N, const float alpha, + const float *X, const size_t incX, const float *Y, const size_t incY, + float *A, const size_t lda); + + +void gemv( const char trans, size_t m, size_t n, + double alpha, const double *a, size_t lda, + const double *x, size_t incx, double beta, + double *y, size_t incy ); + +void gemv( const char trans, size_t m, size_t n, + float alpha, const float *a, size_t lda, + const float *x, size_t incx, float beta, + float *y, size_t incy ); +} + + +#endif + diff --git a/fields.hpp b/fields.hpp index d7df45a..1dc9b28 100644 --- a/fields.hpp +++ b/fields.hpp @@ -4,6 +4,8 @@ #include #include #include "parameters.hpp" +#include "splines.hpp" +#include "lsmr.hpp" using namespace dealii; @@ -36,8 +38,8 @@ inline double compute_rho(const double x, } -template -real eval( real x, const real *coeffs) noexcept +template +double eval(double x, const double *coeffs) noexcept { using std::floor; @@ -48,7 +50,7 @@ real eval( real x, const real *coeffs) noexcept x = x - Parameters::LX * floor( x/Parameters::LX ); // Knot number - real x_knot = floor( x/Parameters::SPLINE_DX); + double x_knot = floor( x/Parameters::SPLINE_DX); size_t ii = static_cast(x_knot); @@ -56,10 +58,10 @@ real eval( real x, const real *coeffs) noexcept x = x/Parameters::SPLINE_DX - x_knot; // Scale according to derivative. - real factor = 1; + double factor = 1; for ( size_t i = 0; i < dx; ++i ) factor *= 1/Parameters::SPLINE_DX; - return factor*splines1d::eval( x, coeffs + ii ); + return factor*splines1d::eval(x, coeffs + ii); } template @@ -72,17 +74,15 @@ void interpolate( real *coeffs, const real *values) // Least Squares needs to be 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 } + mat_t() { 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; @@ -103,10 +103,9 @@ void interpolate( real *coeffs, const real *values) // Least Squares needs to be 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 } + transposed_mat_t() { splines1d::N(0,N); } @@ -132,9 +131,9 @@ void interpolate( real *coeffs, const real *values) // Least Squares needs to be } }; - mat_t M { config }; transposed_mat_t Mt { config }; + mat_t M; transposed_mat_t Mt; lsmr_options opt; opt.silent = true; - lsmr( config.Nx, config.Nx, M, Mt, values, tmp.get(), opt ); + lsmr( Parameters::SPLINE_NX, Parameters::SPLINE_NX , M, Mt, values, tmp.get(), opt ); if ( opt.iter == opt.max_iter ) std::cerr << "Warning. LSMR did not converge.\n"; diff --git a/lsmr.hpp b/lsmr.hpp new file mode 100644 index 0000000..e3703e5 --- /dev/null +++ b/lsmr.hpp @@ -0,0 +1,252 @@ +#ifndef LSMR_HPP +#define LSMR_HPP + +#include +#include +#include +#include +#include "blas.hpp" + +template +struct lsmr_options +{ + /////////// + // INPUT // + /////////// + + // Whether to print messages to std::cout. + bool silent = true; + + // Residual of normal equations AᵀAx = Aᵀb + bool relative_residual = true; + real target_residual = std::numeric_limits::epsilon(); + size_t max_iter = 1000; + + // How many Lánczos vectors to keep for local reorthogonalisation. + // Choose zero for no reorthogonalisation, pure LSMR. + // Choose a large value for complete reorthognalisation. + // + // In an ideal world without roundoff errors, this would have no effect + // at all, as the Lánczos vectors would be perfectly orthogonal. In practice + // this property is lost rather quickly. One may choose to store some of + // the most recent Lánczos vectors to enforce this property manually. This + // increase convergence speed at the cost of additional memory requirements. + size_t reorthogonalise_u = 50; + size_t reorthogonalise_v = 50; + + //////////// + // OUTPUT // + //////////// + + // Iteration count and reached residual. + // Estimates of ‖A‖ and cond(A) + size_t iter; real residual; + real norm_A_estimate, cond_estimate; +}; + +template +void lsmr( size_t m, size_t n, const mat& A, const transposed_mat& At, + const real *b, real *x, lsmr_options &S ); + +namespace lsmr_impl +{ + +template +real norm( size_t n, const real *x ) +{ + using std::hypot; + + real result = 0; + for ( size_t i = 0; i < n; ++i ) + result = hypot(result,x[i]); + + return result; +} + +// Reorthognalise u with respect to the previous vectors in buffer, +// using the modified Gram–Schmidt process. Overwrite the oldest vector +// in buffer when full. +template +void reorthogonalise( real *buf, size_t n, size_t buffer_max, + real *u, size_t iter ) +{ + using std::min; + using blas::dot; + using blas::axpy; + using blas::scal; + using blas::copy; + + size_t n_buffered = min( iter+1, buffer_max ); + for ( size_t i = 0; i < n_buffered; ++i ) + { + real fac = -dot( n, u, 1, buf + i*n, 1 ); + axpy( n, fac, buf + i*n, 1, u, 1 ); + } + + scal( n, 1/norm(n,u), u, 1 ); + copy( n, u, 1, buf + ((iter+1)%buffer_max)*n, 1 ); +} + +} + +template +void lsmr( size_t m, size_t n, const mat& A, const transposed_mat& At, + const real *b, real *x, lsmr_options &S ) +{ + using std::min; + using std::max; + using std::abs; + using std::swap; + using std::hypot; + using blas::axpy; + using blas::scal; + using blas::copy; + using lsmr_impl::norm; + using lsmr_impl::reorthogonalise; + + + // Allocation of buffers. + size_t max_buf = min(n,m)-1; + S.reorthogonalise_u = min(S.reorthogonalise_u,max_buf); + S.reorthogonalise_v = min(S.reorthogonalise_v,max_buf); + size_t u_buffer_size = max( S.reorthogonalise_u, size_t(1) ); + size_t v_buffer_size = max( S.reorthogonalise_v, size_t(1) ); + + std::unique_ptr data { new real[ n*( 4 + v_buffer_size ) + + m*( 2 + u_buffer_size ) ] {} }; + + real *u = data.get(); + real *utmp = u + m; + real *ubuf = utmp + m; + real *v = ubuf + m*u_buffer_size; + real *vtmp = v + n; + real *h = vtmp + n; + real *h_bar = h + n; + real *vbuf = h_bar + n; + + At(b,v); + const real norm_ATb = norm(n,v); + + + A(x,u); axpy(m,real(-1),b,1,u,1); + scal(m, real(-1), u, 1 ); // u = b - Ax; + + real alpha = 0; + real beta = norm(m,u); + + if ( beta > real(0) ) + { + scal(m, real(1)/beta, u, 1 ); // u = b - Ax / norm(b-Ax) + At(u,v); // v = At*u + alpha = norm(n,v); + } + + if ( alpha > real(0) ) + scal(n, real(1)/alpha, v, 1 ); // v = At*u/norm(At*u) + + copy(n,u,1,ubuf,1); // u_buf.col(0) = u_buf + copy(n,v,1,vbuf,1); // v_buf.col(0) = v + copy(n,v,1,h,1); // h = v + + if ( alpha * beta == real(0) ) return; + + + real alpha_bar = alpha, zeta_bar = alpha*beta; + real rho = 1, rho_bar = 1, c_bar = 1, s_bar = 0; + real c, s, theta, zeta, theta_bar, rho_prev, rho_bar_prev; + + // For estimating the condition number. + real sigma_max = 0, sigma_min = std::numeric_limits::max(); + real rho_bar_max = 0, rho_bar_min = std::numeric_limits::max(); + + S.norm_A_estimate = 0; + for ( S.iter = 0; S.iter < S.max_iter; ++S.iter ) + { + // Continue the bidiagonalisation. + A(v,utmp); axpy(m,-alpha,u,1,utmp,1); swap(u,utmp); // u = A*v - alpha*u + beta = norm(m,u); + + if ( beta > 0 ) + { + scal(m, real(1)/beta, u, 1 ); + if ( S.reorthogonalise_u ) + reorthogonalise( ubuf, m, u_buffer_size, u, S.iter ); + + S.norm_A_estimate = hypot( alpha, S.norm_A_estimate ); + S.norm_A_estimate = hypot( beta , S.norm_A_estimate ); + + At(u,vtmp); axpy(n,-beta,v,1,vtmp,1); swap(v,vtmp); // v = At*u - beta*v + alpha = norm(n,v); + + if ( alpha > 0 ) + { + scal(n,real(1)/alpha, v, 1 ); + if ( S.reorthogonalise_v ) + reorthogonalise( vbuf, n, v_buffer_size, v, S.iter ); + } + } + + // Construct and apply rotation P_k + rho_prev = rho; + rho = hypot(alpha_bar,beta); + c = alpha_bar/rho; + s = beta/rho; + theta = s*alpha; + alpha_bar = c*alpha; + + // Construct and apply rotation \bar{P}_k + rho_bar_prev = rho_bar; + if ( S.iter ) + { + rho_bar_max = max( rho_bar, rho_bar_max ); + rho_bar_min = min( rho_bar, rho_bar_min ); + } + theta_bar = s_bar*rho; + rho_bar = hypot( c_bar*rho, theta ); + if ( S.iter ) + { + sigma_max = max( rho_bar_max, c_bar*rho ); + sigma_min = min( rho_bar_min, c_bar*rho ); + } + c_bar = c_bar * rho/rho_bar; + s_bar = theta/rho_bar; + zeta = c_bar * zeta_bar; + zeta_bar = -s_bar*zeta_bar; + + + // Update h, h_bar, x + scal(n, -(theta_bar*rho)/(rho_prev*rho_bar_prev), h_bar, 1 ) ; + axpy(n, real(1), h, 1, h_bar, 1 ); // h_bar = h - factor*h_bar + + axpy( n, zeta/(rho*rho_bar), h_bar, 1, x, 1 ); // x += factor * h_bar + + scal(n, -theta/rho, h, 1 ); + axpy(n, real(1), v, 1, h, 1 ); // h = v - factor*h; + + // Estimate quantities. + if ( S.relative_residual ) S.residual = abs(zeta_bar)/norm_ATb; + else S.residual = abs(zeta_bar); + S.cond_estimate = sigma_max / sigma_min; + + if ( S.residual <= S.target_residual ) + { + if ( S.silent == false ) + { + std::cout << "LSMR: Iteration: " << std::setw(4) << S.iter << ", " + << "Residual: " << std::setw(12) << std::scientific << S.residual << ", " + << "cond estimate: " << std::setw(12) << std::scientific << S.cond_estimate << ".\n"; + } + return; + } + + if ( S.silent == false && (S.iter%10) == 0 ) + { + std::cout << "LSMR: Iteration: " << std::setw(4) << S.iter << ", " + << "Residual: " << std::setw(12) << std::scientific << S.residual << ", " + << "cond estimate: " << std::setw(12) << std::scientific << S.cond_estimate << ".\n"; + } + } +} + +#endif + diff --git a/nufi_solver.hpp b/nufi_solver.hpp index adb5078..e3502b5 100644 --- a/nufi_solver.hpp +++ b/nufi_solver.hpp @@ -8,15 +8,14 @@ #include #include +#include #include #include #include #include "parameters.hpp" #include "poisson_problem.hpp" -#include "fields.hpp" // holds f0(x,v), and compute_rho(x) -#include "spline_field.hpp" // old GPT splines -#include "splines.hpp" //new splines +#include "fields.hpp" using namespace dealii; @@ -26,15 +25,15 @@ public: NuFISolver(); void run(); - 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); + 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); private: unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT); - [[maybe_unused]] unsigned int Nx = Parameters::SPLINE_NX; + unsigned int Nx = Parameters::SPLINE_NX; double Lx = Parameters::LX; @@ -42,8 +41,6 @@ private: unsigned int order; - double dt = Parameters::DT; - PoissonProblem<1> poisson; }; @@ -51,72 +48,82 @@ private: inline double NuFISolver::eval_ftilda(unsigned int n, double x, double u, - const std::vector E_coeffs) + const double *E_coeffs) + { - double Lu = std::abs(Parameters::V_DOMAIN_LEFT - Parameters::V_DOMAIN_RIGHT); + if ( n == 0 ) return f0(x,u); - if (n == 0) - return f0(x, u); + const size_t stride_x = 1; + const size_t stride_t = stride_x*(Parameters::SPLINE_NX + Parameters::SPLINE_ORDER - 1); - // Initial half-step. - u += 0.5*dt*E_spline.eval(x); + double Ex; + const double *c; + + // We omit the initial half-step. while ( --n ) { - x -= dt*u; - u += dt*E_spline.eval(x); + x = x - Parameters::DT *u; + c = E_coeffs + n*stride_t; + Ex = -eval<1>(x, c); + u = u + Parameters::DT *Ex; } - // Final half-step. - x -= dt*u; - u += 0.5*dt*E_spline.eval(x); + // The final half-step. + x -= Parameters::DT*u; + c = E_coeffs + n*stride_t; + Ex = -eval<1>(x, c); + u += 0.5*Parameters::DT*Ex; - double x_periodic = x - Lx * std::floor(x / Lx); - double u_periodic = u - Lu * std::floor(u / Lu); - - return f0(x_periodic, u_periodic); + return f0(x,u); } inline double NuFISolver::eval_rho(const unsigned int n, const double x, - const std::vector E_coeffs, + const double *E_coeffs, const unsigned int Nv) { const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; + const double v_min = Parameters::V_DOMAIN_LEFT; double integral = 0.0; for (unsigned int i = 0; i < Nv; ++i) - { - const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv; - AssertThrow(std::isfinite(E_spline.eval(x)), ExcMessage("NaN detected in E_spline.eval(x) inside NuFISolver::eval_rho integral loop")); - integral += eval_ftilda(n, x, v, E_spline) * dv; - } + integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs) * dv; return 1.0 - integral; } -class ChargeDensity_NuFI : public Function<1> +template +class ChargeDensity_NuFI : public Function { -public: - ChargeDensity_NuFI(NuFISolver &solver, size_t n, const std::vector E_coeffs) - : solver(solver), n(n), E_coeffs(E_coeffs) {} + public: + ChargeDensity_NuFI(const double *rho_values, unsigned int Nx) + : Function(), rho(rho_values), Nx(Nx) {} - virtual double value(const Point<1> &p, - [[maybe_unused]] const unsigned int component = 0) const override - { - double x = p[0]; + virtual double value(const Point &p, + [[maybe_unused]] const unsigned int component = 0) const override + { + const double x = p[0]; - return solver.eval_rho(n, x, E_coeffs); - } + // Map x -> grid index + const double L = Parameters::LX; + const double dx = L / (Nx-1); -private: - NuFISolver &solver; - size_t n; - const std::vector E_coeffs; + int i = static_cast(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; }; inline void NuFISolver::save_ftilda(unsigned int n, - const std::vector E_coeffs, + const double *E_coeffs, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename) @@ -144,7 +151,7 @@ inline void NuFISolver::save_ftilda(unsigned int n, { double v = vmin + (j + 0.5)*dv; - double val = eval_ftilda(n, x, v, E_spline); + double val = eval_ftilda(n, x, v, E_coeffs); file << val; @@ -160,60 +167,39 @@ inline void NuFISolver::save_ftilda(unsigned int n, inline void NuFISolver::run() { - std::cout << "Start of NuFISolver::run()\n\n"; - - // init E_spline + std::cout << "Building E_sline\n\n"; - unsigned int Nx = Parameters::SPLINE_NX; - // Nx grid points - double dx = Lx / (Nx-1); + using std::abs; + using std::max; - std::vector E_grid(Nx); + const size_t stride_t = Nx + order - 1; - //set initial E points - for (unsigned int i=0; i coeffs { new double[ Nt*stride_t ] {} }; + std::unique_ptr rho { reinterpret_cast(std::aligned_alloc(64,sizeof(double)*Nx)), std::free }; - std::vector E_coeffs(E_grid, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT); // Needs correction + if ( rho == nullptr ) throw std::bad_alloc {}; for (unsigned int it = 0; it < Nt; ++it) { std::cout << "Timestep " << it << " / " << Nt << std::endl << std::endl; - // Step 1: Evaluate rho^n(x) using current E_spline - - rho.resize(Nx); - - for (unsigned int i = 0; i < (Nx); ++i) - { - double x = (i + 0.5) * dx; - rho[i] = eval_rho(it, x, E_spline, Parameters::NV); - } - - ChargeDensity_NuFI rho_function(*this, it, E_spline); - - poisson.set_rhs_function(rho_function); - - for (unsigned int i=0; i< rho.size(); ++i) // check for bad rho[i] - { - AssertThrow(std::isfinite(rho[i]), ExcMessage("NaN detected in rho")); - } + // compute rho + for(size_t i = 0; i>(rho.get(), Parameters::SPLINE_NX)); poisson.solve_step(); if (it % Parameters::PLOT_FREQUENCY == 0) { std::cout << "Saving results... \n\n"; - save_ftilda(it, E_spline, 128, 128, "results/ftilda_" + std::to_string(it) + ".dat"); + save_ftilda(it, coeffs.get(), 128, 128, "results/ftilda_" + std::to_string(it) + ".dat"); poisson.output_results(it); } - - E_grid = poisson.sample_electric_field(poisson, Nx, 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 1860973..a74a274 100644 --- a/parameters.hpp +++ b/parameters.hpp @@ -3,6 +3,7 @@ #include #include + namespace Parameters { constexpr unsigned int DIMENSION = 1; @@ -27,9 +28,11 @@ namespace Parameters constexpr double DT=1./16.; constexpr unsigned int TMAX = 10; + //spline options constexpr int SPLINE_NX = 512; constexpr double SPLINE_DX = LX/SPLINE_NX; + constexpr size_t SPLINE_ORDER = 4; //Plotting options constexpr int PLOT_FREQUENCY = 2; diff --git a/poisson_problem.hpp b/poisson_problem.hpp index eb5a768..81c494e 100644 --- a/poisson_problem.hpp +++ b/poisson_problem.hpp @@ -33,7 +33,9 @@ #include #include +#include #include +#include #include #include "parameters.hpp" @@ -53,7 +55,7 @@ public: void solve_step(); void run(); - void set_rhs_function(const Function &rhs); + 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; } @@ -83,7 +85,7 @@ private: Vector solution; // phi Vector system_rhs; - const Function *rhs_function; + std::unique_ptr> rhs_function; MappingQ mapping; }; @@ -91,9 +93,9 @@ private: // Utilities template -void PoissonProblem::set_rhs_function(const Function &rhs) +void PoissonProblem::set_rhs_function(std::unique_ptr> rhs) { - rhs_function = &rhs; + rhs_function = std::move(rhs); } template @@ -214,9 +216,6 @@ public: } }; - -// =-=-=-=-= Poisson equation solver =-=-=-=-= - template void PoissonProblem::assemble_system() {