diff --git a/CMakeLists.txt b/CMakeLists.txt index 54666d7..a565c27 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,7 +23,7 @@ find_package(OpenMP REQUIRED) add_library(nufi_lib src/nufi_solver.cc src/save_results.cc - src/blas.cc + # src/blas.cc ) target_include_directories(nufi_lib PUBLIC diff --git a/Makefile b/Makefile index 1b6e392..38f5543 100644 --- a/Makefile +++ b/Makefile @@ -142,30 +142,6 @@ nufi_poisson/fast: $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_poisson.dir/build.make CMakeFiles/nufi_poisson.dir/build .PHONY : nufi_poisson/fast -src/blas.o: src/blas.cc.o -.PHONY : src/blas.o - -# target to build an object file -src/blas.cc.o: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.o -.PHONY : src/blas.cc.o - -src/blas.i: src/blas.cc.i -.PHONY : src/blas.i - -# target to preprocess a source file -src/blas.cc.i: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.i -.PHONY : src/blas.cc.i - -src/blas.s: src/blas.cc.s -.PHONY : src/blas.s - -# target to generate assembly for a file -src/blas.cc.s: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.s -.PHONY : src/blas.cc.s - src/main.o: src/main.cc.o .PHONY : src/main.o @@ -248,9 +224,6 @@ help: @echo "... rebuild_cache" @echo "... nufi_lib" @echo "... nufi_poisson" - @echo "... src/blas.o" - @echo "... src/blas.i" - @echo "... src/blas.s" @echo "... src/main.o" @echo "... src/main.i" @echo "... src/main.s" diff --git a/README.md b/README.md index 14ea93a..74c5412 100644 --- a/README.md +++ b/README.md @@ -6,4 +6,5 @@ This simulation of the Vlasov-Poisson system in 1x1v dimensions uses --- Todo: -- Correct spline interpolator/evaluation for discontinuity +- fix eval and saving fields. + - something with periodicity or eval range in x diff --git a/libnufi_lib.a b/libnufi_lib.a index 14ea604..4dbcb62 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/blas.h b/nufi/blas.h deleted file mode 100644 index 2f90e84..0000000 --- a/nufi/blas.h +++ /dev/null @@ -1,54 +0,0 @@ -#ifndef NUFI_BLAS_H -#define NUFI_BLAS_H - -#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/nufi/fields.h b/nufi/fields.h index 18e5e58..0e5b1fe 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -1,229 +1,112 @@ #ifndef FIELDS_H #define FIELDS_H +#include "nufi/parameters.h" +#include "nufi/poisson_problem.h" #include #include -#include "nufi/parameters.h" -#include "nufi/splines.h" -#include "nufi/lsmr.h" +#include using namespace dealii; -inline double f0(const double x, - const double v, +// inline std::vector Indices_of_points(const std::vector &points, +// double x_min, double x_max, double dx, int grid_type=0) +// { +// // grid type: +// // 0 => uniform +// // 1 => non uniform (TODO) +// +// if (dx <= 0.0) { +// throw std::invalid_argument("dx must be positive"); +// } +// if (x_max <= x_min) { +// throw std::invalid_argument("x_max must be > x_min"); +// } +// +// std::vector indices; +// indices.reserve(points.size()); +// +// switch (grid_type) { +// case 0: +// { +// const double L = x_max - x_min; +// const int N = std::floor(L/dx); +// +// +// for (double x : points) //GPT loop, to check +// { +// x-= x_min; +// x = x - L * std::floor(x/L); +// +// int i = static_cast(std::floor(x / dx)); +// +// // safety: handle rare edge case due to floating precision +// if (i == N) i = 0; +// +// indices.push_back(i); +// } +// } +// case 1: +// { +// throw std::invalid_argument("Case for non uniform grid is not +// completed"); +// } +// default: +// throw std::invalid_argument("Invalid grid_type argument"); +// +// } +// return indices; +// } + +inline double f0(const double x, const double v, const double eps = Parameters::EPS, - const double k = Parameters::WAVE_NR) -{ - const double prefactor = Parameters::F0_FACTOR * (1.0 + eps * std::cos(k*x)); - const double gaussian = v*v * std::exp(-0.5 * v*v); + const double k = Parameters::WAVE_NR) { + const double prefactor = + Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x)); + const double gaussian = v * v * std::exp(-0.5 * v * v); return prefactor * gaussian; } +// wrapper for eval_point() { VectorTools::point_values() } +inline double eval(double x, const PoissonProblem<1> &poisson, + const Vector &solution) noexcept { -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; + x -= Parameters::X_DOMAIN_LEFT; - double integral = 0.0; + x = x - Parameters::LX * std::floor(x * Parameters::LX_INV); // in domain - 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; + return eval_point_grad<1>(poisson.get_mapping(), poisson.get_dof_handler(), + solution, Point<1>(x)); } - -template -double eval(double x, const double *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_INV ); - - // Knot number - double x_knot = floor( x*Parameters::SPLINE_DX_INV); - - size_t ii = static_cast(x_knot); - - // Convert x to reference coordinates. - x = x*Parameters::SPLINE_DX_INV - x_knot; - - // Scale according to derivative. - double factor = 1; - for ( size_t i = 0; i < dx; ++i ) factor *= 1*Parameters::SPLINE_DX_INV; - - return factor*splines1d::eval(x, coeffs + ii); -} - -template -void interpolate( real *coeffs, const real *values) -{ - 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 - { - real N[ order ]; - - 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; - 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 - { - real N[ order ]; - - transposed_mat_t() - { - 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; transposed_mat_t Mt; - lsmr_options opt; opt.silent = true; - 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"; - - for ( size_t i = 0; i < Parameters::SPLINE_NX + order - 1; ++i ) - coeffs[ i ] = tmp[ i % Parameters::SPLINE_NX ]; -} - -inline double integral_space_vector(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX) -{ +inline double integral_space_vector(const PoissonProblem<1> &poisson, + const Vector &solution, + double dx = Parameters::PLOT_DX, + size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; double xmin = Parameters::X_DOMAIN_LEFT; - #pragma omp parallel for reduction (+:integral) - for (size_t i=0; i(x, current_coeffs); + integral += eval(x, poisson, solution); } - return integral*dx; + return integral * dx; }; -inline double integral_space_vector_squared(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX) -{ +inline double integral_space_vector_squared(const PoissonProblem<1> &poisson, + const Vector &solution, + double dx = Parameters::PLOT_DX, + size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; double xmin = Parameters::X_DOMAIN_LEFT; - #pragma omp parallel for reduction (+:integral) - for (size_t i=0; i(x, current_coeffs); - integral += val*val; +#pragma omp parallel for reduction(+ : integral) + for (size_t i = 0; i < Nx; ++i) { + double x = xmin + i * dx; + double val = eval(x, poisson, solution); + integral += val * val; } - return integral*dx; + 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 compute(const std::vector& values) const { - size_t n = values.size(); - if (n < 2) { - throw std::invalid_argument("Need at least 2 points"); - } - - std::vector 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 -class ChargeDensity : public Function // only uses f0 -{ -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/lsmr.h b/nufi/lsmr.h deleted file mode 100644 index 144793f..0000000 --- a/nufi/lsmr.h +++ /dev/null @@ -1,252 +0,0 @@ -#ifndef LSMR_H -#define LSMR_H - -#include -#include -#include -#include -#include "nufi/blas.h" - -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/nufi_solver.h b/nufi/nufi_solver.h index a023895..95fc301 100644 --- a/nufi/nufi_solver.h +++ b/nufi/nufi_solver.h @@ -2,33 +2,37 @@ #define NUFI_SOLVER_H #include -#include #include #include #include -#include +#include +#include +#include "nufi/fields.h" //dont remove #include "nufi/parameters.h" #include "nufi/poisson_problem.h" -#include "nufi/fields.h" //dont remove using namespace dealii; -class NuFISolver -{ +class NuFISolver { public: NuFISolver(); void run(); - double eval_rho(unsigned int n, double x, const double *E_coeffs, unsigned int Nv = Parameters::NV) const; - double eval_ftilda(unsigned int n, double x, double u, const double *E_coeffs) const; - double eval_f(unsigned int n, double x, double u, const double *E_coeffs) const; + double eval_rho(unsigned int n, const double x, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history, + const unsigned int Nv = Parameters::NV) const; + double eval_ftilda(unsigned int n, double x, double u, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history) const; + double eval_f(unsigned int n, double x, double u, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history) const; private: - - - unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT); - unsigned int Nx = Parameters::SPLINE_NX; + unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT); + unsigned int Nx = Parameters::CALC_NX; double Lx = Parameters::LX; @@ -40,36 +44,5 @@ private: unsigned int order; PoissonProblem<1> poisson; - }; - -template -class ChargeDensity_NuFI : public Function -{ - public: - ChargeDensity_NuFI(const double *rho_values, unsigned int Nx) - : Function(), rho(rho_values), Nx(Nx) {} - - virtual double value(const Point &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(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 diff --git a/nufi/parameters.h b/nufi/parameters.h index 97618e6..c379718 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -4,43 +4,41 @@ #include #include -namespace Parameters -{ - constexpr unsigned int DIMENSION = 1; +namespace Parameters { +constexpr unsigned int DIMENSION = 1; - constexpr double X_DOMAIN_LEFT = 0.0; - constexpr double X_DOMAIN_RIGHT = 4*M_PI; - constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT); - constexpr double LX_INV = 1/LX; +constexpr double X_DOMAIN_LEFT = 0.0; +constexpr double X_DOMAIN_RIGHT = 4 * M_PI; +constexpr double LX = std::abs(X_DOMAIN_RIGHT - X_DOMAIN_LEFT); +constexpr double LX_INV = 1 / LX; - constexpr double V_DOMAIN_LEFT = -10.; - constexpr double V_DOMAIN_RIGHT = 10.; +constexpr size_t CALC_NX = 128; +constexpr double CALC_DX = LX / CALC_NX; - constexpr unsigned int NV = 512; - constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT)/NV; +constexpr double V_DOMAIN_LEFT = -10.; +constexpr double V_DOMAIN_RIGHT = 10.; - // deal.ii options - constexpr unsigned int GLOBAL_REFINEMENT = 8; - constexpr unsigned int FE_DEGREE = 4; - constexpr unsigned int CONVERGENCE_ITERATIONS = 10000; - constexpr double CONVERGENCE_LIMIT = 1e-12; +constexpr unsigned int NV = 128; +constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV; - constexpr double EPS = 0.01; - constexpr double WAVE_NR = 0.5; - constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) +// deal.ii options +constexpr unsigned int GLOBAL_REFINEMENT = 6; +constexpr unsigned int FE_DEGREE = 2; +constexpr unsigned int CONVERGENCE_ITERATIONS = 5000; +constexpr double CONVERGENCE_LIMIT = 1e-8; - // NUFI options - constexpr double DT=1./16.; - constexpr unsigned int TMAX = 500; +constexpr double EPS = 0.01; +constexpr double WAVE_NR = 0.5; +constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) - //spline options - constexpr int SPLINE_NX = 256; - constexpr double SPLINE_DX = LX/(SPLINE_NX); - constexpr double SPLINE_DX_INV = 1/SPLINE_DX; - constexpr size_t SPLINE_ORDER = 4; +// NUFI options +constexpr double DT = 1. / 8.; +constexpr unsigned int TMAX = 100; - //Plotting options - constexpr int PLOT_FREQUENCY = 10; -} +// Plotting options +constexpr int PLOT_FREQUENCY = 4; +constexpr size_t PLOT_NX = CALC_NX; +constexpr double PLOT_DX = LX / PLOT_NX; +} // namespace Parameters #endif diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index 7571603..354903e 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -3,41 +3,45 @@ #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 #include -#include #include +#include #include #include #include -#include -#include #include +#include +#include #include #include +#include +#include +#include #include #include #include @@ -49,9 +53,7 @@ using namespace dealii; // =-=-=-=-= Poisson Solver =-=-=-=-= -template -class PoissonProblem -{ +template class PoissonProblem { public: PoissonProblem(unsigned int degree); @@ -59,13 +61,16 @@ public: void solve_step(); void run(); - void set_rhs_function(std::unique_ptr> rhs_function); + void set_rhs_function(std::function &)> f); const Vector &get_solution() const { return solution; } + const MappingQ &get_mapping() const { return mapping; } const DoFHandler &get_dof_handler() const { return dof_handler; } - std::vector sample_electric_field(double x_min, double x_max, unsigned int Nx); - std::vector sample_electric_potential(double x_min, double x_max, unsigned int Nx); + std::vector sample_electric_field(double x_min, double x_max, + unsigned int Nx); + std::vector sample_electric_potential(double x_min, double x_max, + unsigned int Nx); private: void create_mesh(); @@ -74,18 +79,18 @@ private: void solve(); Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; + FE_Q fe; + DoFHandler dof_handler; AffineConstraints constraints; - SparsityPattern sparsity_pattern; + SparsityPattern sparsity_pattern; SparseMatrix system_matrix; - Vector solution; // phi + Vector solution; // phi Vector system_rhs; - std::unique_ptr> rhs_function; + std::function &)> rhs_function; MappingQ mapping; }; @@ -93,81 +98,73 @@ private: // Utilities template -void PoissonProblem::set_rhs_function(std::unique_ptr> rhs) -{ - rhs_function = std::move(rhs); +void PoissonProblem::set_rhs_function( + std::function &)> f) { + rhs_function = std::move(f); } template PoissonProblem::PoissonProblem(unsigned int degree) - : fe(degree) - , dof_handler(triangulation) - , mapping(degree) -{} + : fe(degree), dof_handler(triangulation), mapping(degree) {} template -std::vector PoissonProblem::sample_electric_field(double x_min,double x_max,unsigned int Nx) -{ +std::vector +PoissonProblem::sample_electric_field(double x_min, double x_max, + unsigned int Nx) { std::vector E_values(Nx); const double dx = (x_max - x_min) / (Nx - 1); - for (unsigned int i = 0; i < Nx; ++i) - { - const double x = x_min + i * dx; - const Point point(x); + for (unsigned int i = 0; i < Nx; ++i) { + const double x = x_min + i * dx; + const Point point(x); - // 1. Find the active cell containing x - const auto cell_point_pair = - GridTools::find_active_cell_around_point(mapping, - dof_handler, - point); + // 1. Find the active cell containing x + const auto cell_point_pair = + GridTools::find_active_cell_around_point(mapping, dof_handler, point); - const auto cell = cell_point_pair.first; - const Point &unit_point = cell_point_pair.second; + const auto cell = cell_point_pair.first; + const Point &unit_point = cell_point_pair.second; - // 2. FEPointEvaluation expects an ArrayView of points - std::vector> points(1, unit_point); - ArrayView> point_view(points); + // 2. FEPointEvaluation expects an ArrayView of points + std::vector> points(1, unit_point); + ArrayView> point_view(points); - FEPointEvaluation<1, dim> evaluator(mapping, - dof_handler.get_fe(), - update_gradients); + FEPointEvaluation<1, dim> evaluator(mapping, dof_handler.get_fe(), + update_gradients); - // reinit with ArrayView of points - evaluator.reinit(cell, point_view); + // reinit with ArrayView of points + evaluator.reinit(cell, point_view); - Vector local_dofs(dof_handler.get_fe().dofs_per_cell); - cell->get_dof_values(solution, local_dofs); + Vector local_dofs(dof_handler.get_fe().dofs_per_cell); + cell->get_dof_values(solution, local_dofs); - // 3. Evaluate gradient at this point - evaluator.evaluate(local_dofs, EvaluationFlags::gradients); + // 3. Evaluate gradient at this point + evaluator.evaluate(local_dofs, EvaluationFlags::gradients); - const Tensor<1, dim> grad_phi = evaluator.get_gradient(0); + const Tensor<1, dim> grad_phi = evaluator.get_gradient(0); - // 4. Compute E = -grad(phi) - E_values[i] = -grad_phi[0]; + // 4. Compute E = -grad(phi) + E_values[i] = -grad_phi[0]; } return E_values; } template -std::vector PoissonProblem::sample_electric_potential( - double x_min, - double x_max, - unsigned int Nx) -{ +std::vector +PoissonProblem::sample_electric_potential(double x_min, double x_max, + unsigned int Nx) { std::vector values(Nx); std::vector> eval_points(Nx); double Lx = x_max - x_min; double dx = Lx / Nx; - for(unsigned int i=0 ; i(x_min + i * dx); - Utilities::MPI::RemotePointEvaluation cache; + Utilities::MPI::RemotePointEvaluation cache; cache.reinit(eval_points, triangulation, mapping); values = VectorTools::point_values(cache, dof_handler, solution); @@ -175,45 +172,72 @@ std::vector PoissonProblem::sample_electric_potential( return values; } -// dealii Poisson - -template -void PoissonProblem::create_mesh() -{ +template +double +eval_point_grad(const Mapping &mapping, const DoFHandler &dof_handler, + const Vector &solution, const Point &point) { - GridGenerator::hyper_cube(triangulation, - Parameters::X_DOMAIN_LEFT, + Tensor<1, dim> grad = + VectorTools::point_gradient(mapping, dof_handler, solution, point); + double Ex = grad[0]; + + return Ex; +} + +template +double eval_point_value(const Mapping &mapping, + const DoFHandler &dof_handler, + const Vector &solution, + const Point &point) { + + return VectorTools::point_value(mapping, dof_handler, solution, point); +} +// dealii Poisson +template void PoissonProblem::create_mesh() { + + GridGenerator::hyper_cube(triangulation, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT); + std::vector< + GridTools::PeriodicFacePair::cell_iterator>> + periodic_faces; - std::vector::cell_iterator>> periodic_faces; - - GridTools::collect_periodic_faces(triangulation, - 0, 1, // boundary IDs - 0, - periodic_faces); + GridTools::collect_periodic_faces(triangulation, 0, 1, // boundary IDs + 0, periodic_faces); triangulation.add_periodicity(periodic_faces); triangulation.refine_global(Parameters::GLOBAL_REFINEMENT); } -template -void PoissonProblem::setup_system() -{ +template void PoissonProblem::setup_system() { dof_handler.distribute_dofs(fe); constraints.clear(); - + DoFTools::make_hanging_node_constraints(dof_handler, constraints); - DoFTools::make_periodicity_constraints(dof_handler, - 0, 1, - 0, - constraints); - + DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints); + + // Gauge fix for periodic Poisson: + // remove the constant nullspace by pinning one unconstrained DoF. + // (by Paul Wilhelm) + types::global_dof_index gauge_dof = numbers::invalid_dof_index; + + for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) { + if (!constraints.is_constrained(i)) { + gauge_dof = i; + break; + } + } + + Assert(gauge_dof != numbers::invalid_dof_index, + ExcMessage("No unconstrained DoF found for gauge fixing.")); + + constraints.add_line(gauge_dof); + constraints.set_inhomogeneity(gauge_dof, 0.0); + constraints.close(); DynamicSparsityPattern dsp(dof_handler.n_dofs()); @@ -226,112 +250,77 @@ void PoissonProblem::setup_system() system_rhs.reinit(dof_handler.n_dofs()); } -template -void PoissonProblem::assemble_system() -{ - QGauss quadrature_formula(fe.degree + 1); +// Paul +template void PoissonProblem::assemble_system() { + Assert(system_matrix.m() == dof_handler.n_dofs(), + ExcMessage("Matrix not initialized correctly")); + system_matrix = 0; + system_rhs = 0; + + QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, - update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); + update_values | update_gradients | + update_quadrature_points | 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); + Vector cell_rhs(dofs_per_cell); std::vector local_dof_indices(dofs_per_cell); - Assert(rhs_function != nullptr, ExcMessage("RHS function not set")); - - for (const auto &cell : dof_handler.active_cell_iterators()) - { + for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); cell_matrix = 0; - cell_rhs = 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 auto q : fe_values.quadrature_point_indices()) { + 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()) - 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 + cell_matrix(i, j) += fe_values.shape_grad(i, q) * + fe_values.shape_grad(j, q) * fe_values.JxW(q); 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_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); - 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); + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); } - 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 PoissonProblem::solve() { -template -void PoissonProblem::solve() -{ - - SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS, Parameters::CONVERGENCE_LIMIT); + SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS, + Parameters::CONVERGENCE_LIMIT); SolverCG> solver(solver_control); // PreconditionSSOR> preconditioner; // 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); + constraints.distribute(solution); } -template -void PoissonProblem::initialize() -{ - create_mesh(); // build grid - setup_system(); // distribute DoFs and matrices +template void PoissonProblem::initialize() { + create_mesh(); // build grid + setup_system(); // distribute DoFs and matrices } -template -void PoissonProblem::solve_step() -{ +template void PoissonProblem::solve_step() { assemble_system(); solve(); } - // NuFI doesnt use this, kept only for testing PoissonProblem -template -void PoissonProblem::run() -{ +template void PoissonProblem::run() { create_mesh(); setup_system(); assemble_system(); diff --git a/nufi/save_results.h b/nufi/save_results.h index 2874edc..7e47ae7 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -1,28 +1,26 @@ #ifndef SAVE_RESULTS_H #define SAVE_RESULTS_H -#include #include "nufi/nufi_solver.h" +#include "nufi/poisson_problem.h" +#include +#include +void save_f(const NuFISolver &solver, unsigned int n, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history, unsigned int Nx_out, + unsigned int Nv_out, const std::string &filename); -void save_f( const NuFISolver &solver, - unsigned int n, - const double *E_coeffs, - unsigned int Nx_out, - unsigned int Nv_out, - const std::string &filename); +void save_rho(const NuFISolver &solver, unsigned int n, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history, + unsigned int Nx_out, const std::string &filename); -void save_rho(const NuFISolver &solver, - unsigned int n, - const double *E_coeffs, - unsigned int Nx_out, - const std::string &filename); +void save_Efield(unsigned int n, const PoissonProblem<1> &poisson, + const std::vector> &phi_history, + 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); +void save_space_vector(const std::vector &vals, + const std::string &filename, size_t it); -void save_space_vector(const std::vector& vals, const std::string& filename, size_t it); - -#endif +#endif diff --git a/nufi/splines.h b/nufi/splines.h deleted file mode 100644 index 73e8fd0..0000000 --- a/nufi/splines.h +++ /dev/null @@ -1,90 +0,0 @@ -#ifndef SPLINES_HP -#define SPLINES_HP - -#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 - - diff --git a/nufi/stopwatch.h b/nufi/stopwatch.h index 1d2b5ca..2c6bbe8 100644 --- a/nufi/stopwatch.h +++ b/nufi/stopwatch.h @@ -3,35 +3,27 @@ #include -template -class stopwatch -{ +template class stopwatch { public: - void reset(); - real elapsed(); + void reset(); + real elapsed(); private: - using clock = std::chrono::high_resolution_clock; - clock::time_point t0 { clock::now() }; + using clock = std::chrono::high_resolution_clock; + clock::time_point t0{clock::now()}; }; - -template inline -void stopwatch::reset() -{ - t0 = clock::now(); +template inline void stopwatch::reset() { + t0 = clock::now(); } -template inline -real stopwatch::elapsed() -{ - using seconds = std::chrono::duration>; +template inline real stopwatch::elapsed() { + using seconds = std::chrono::duration>; - auto tnow = clock::now(); - auto duration = std::chrono::duration_cast( tnow - t0 ); + auto tnow = clock::now(); + auto duration = std::chrono::duration_cast(tnow - t0); - return duration.count(); + return duration.count(); } #endif // STOPWATCH_H - diff --git a/rho_E.mp4 b/rho_E.mp4 deleted file mode 100644 index 2564d57..0000000 Binary files a/rho_E.mp4 and /dev/null differ diff --git a/src/blas.cc b/src/blas.cc deleted file mode 100644 index 9107ecc..0000000 --- a/src/blas.cc +++ /dev/null @@ -1,95 +0,0 @@ -#include "nufi/blas.h" -#include - -namespace blas -{ - -double dot( const size_t n, const double *x, size_t incx, - const double *y, size_t incy ) -{ - return cblas_ddot(n,x,incx,y,incy); -} - -float dot( const size_t n, const float *x, size_t incx, - const float *y, size_t incy ) -{ - return cblas_sdot(n,x,incx,y,incy); -} - -void axpy( size_t n, double alpha, const double *x, size_t incx, - double *y, size_t incy ) -{ - cblas_daxpy(n,alpha,x,incx,y,incy); -} - -void axpy( size_t n, float alpha, const float *x, size_t incx, - float *y, size_t incy ) -{ - cblas_saxpy(n,alpha,x,incx,y,incy); -} - -void scal( size_t n, double alpha, double *x, size_t incx ) -{ - cblas_dscal(n,alpha,x,incx); -} - -void scal( size_t n, float alpha, float *x, size_t incx ) -{ - cblas_sscal(n,alpha,x,incx); -} - -void copy( size_t n, const double *x, size_t incx, double *y, size_t incy ) -{ - cblas_dcopy(n,x,incx,y,incy); -} - -void copy( size_t n, const float *x, size_t incx, float *y, size_t incy ) -{ - cblas_scopy(n,x,incx,y,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) -{ - cblas_dger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, 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) -{ - cblas_sger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, 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 ) -{ - if ( trans == 'T' || trans == 'Y' ) - { - cblas_dgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy ); - } - else - { - cblas_dgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, 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 ) -{ - if ( trans == 'T' || trans == 'Y' ) - { - cblas_sgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy ); - } - else - { - cblas_sgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, incy ); - } -} - -} diff --git a/src/main.cc b/src/main.cc index 7e18060..e5f3c67 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,44 +1,35 @@ -#include #include - +#include #include +#include -void clear_results_directory(const std::string &dir) -{ - if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir)) - return; +void clear_results_directory(const std::string &dir) { + if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir)) + return; - for (const auto &entry : std::filesystem::directory_iterator(dir)) - { - if (std::filesystem::is_regular_file(entry)) - { - std::filesystem::remove(entry.path()); - std::cout << "Deleted: " << entry.path() << '\n'; - } + for (const auto &entry : std::filesystem::directory_iterator(dir)) { + if (std::filesystem::is_regular_file(entry)) { + std::filesystem::remove(entry.path()); + std::cout << "Deleted: " << entry.path() << '\n'; } + } } -int main() -{ - #include -std::cout << "Threads: " << omp_get_max_threads() << "\n"; - try - { - clear_results_directory("results"); +int main() { + std::cout << "Threads: " << omp_get_max_threads() << "\n"; + try { + clear_results_directory("results"); - NuFISolver solver; - solver.run(); - } - catch (const std::exception &exc) - { - std::cerr << "\nException:\n" << exc.what() << "\n"; - return 1; - } - catch (...) - { - std::cerr << "\nUnknown exception!\n"; - return 1; - } + NuFISolver solver; + solver.run(); - return 0; + } catch (const std::exception &exc) { + std::cerr << "\nException:\n" << exc.what() << "\n"; + return 1; + } catch (...) { + std::cerr << "\nUnknown exception!\n"; + return 1; + } + + return 0; } diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 83ca6a1..d776d61 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -1,207 +1,192 @@ #include "nufi/nufi_solver.h" #include +#include #include #include #include #include -#include +#include +#include #include #include #include -#include #include -#include "nufi/parameters.h" -#include "nufi/save_results.h" -#include "nufi/poisson_problem.h" #include "nufi/fields.h" +#include "nufi/parameters.h" +#include "nufi/poisson_problem.h" +#include "nufi/save_results.h" #include "nufi/stopwatch.h" using namespace dealii; -double NuFISolver::eval_ftilda(unsigned int n, - double x, - double u, - const double *E_coeffs) const -{ - if ( n == 0 ) return f0(x,u); - - const size_t order = Parameters::SPLINE_ORDER; - const size_t stride_x = 1; - const size_t stride_t = stride_x*(Nx + order - 1); +double +NuFISolver::eval_ftilda(unsigned int n, double x, double u, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history) const { + if (n == 0) + return f0(x, u); double Ex; - const double *c; // We omit the initial half-step. - while ( --n ) - { - x = x - Parameters::DT *u; - c = E_coeffs + n*stride_t; - Ex = -eval<1>(x, c); - u = u + Parameters::DT *Ex; + while (--n) { + x = x - Parameters::DT * u; + Ex = -eval(x, poisson, phi_history[n]); + u = u + Parameters::DT * Ex; } // 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; + x = x - Parameters::DT * u; + Ex = -eval(x, poisson, phi_history[n]); + u += 0.5 * Parameters::DT * Ex; - return f0(x,u); + return f0(x, u); } -double NuFISolver::eval_f(unsigned int n, - double x, - double u, - const double *E_coeffs) const -{ - if ( n == 0 ) return f0(x,u); - - const size_t order = Parameters::SPLINE_ORDER; - const size_t stride_x = 1; - const size_t stride_t = stride_x*(Nx + order - 1); +double +NuFISolver::eval_f(unsigned int n, double x, double u, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history) const { + if (n == 0) + return f0(x, u); double Ex; - const double *c; // Initial half-step. - c = E_coeffs + n*stride_t; - Ex = -eval<1>(x, c); - u += 0.5*Parameters::DT * Ex; + Ex = -eval(x, poisson, phi_history[n]); + u += 0.5 * Parameters::DT * Ex; - while ( --n ) - { - x = x - Parameters::DT *u; - c = E_coeffs + n*stride_t; - Ex = -eval<1>(x, c); - u = u + Parameters::DT *Ex; + while (--n) { + x = x - Parameters::DT * u; + Ex = -eval(x, poisson, phi_history[n]); + u = u + Parameters::DT * Ex; } // 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; + x = x - Parameters::DT * u; + Ex = -eval(x, poisson, phi_history[n]); + u += 0.5 * Parameters::DT * Ex; - return f0(x,u); + return f0(x, u); } -double NuFISolver::eval_rho(const unsigned int n, - const double x, - const double *E_coeffs, - const unsigned int Nv) const -{ - const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; - const double v_min = Parameters::V_DOMAIN_LEFT; +double NuFISolver::eval_rho(unsigned int n, const double x, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history, + const unsigned int Nv) const { + const double dv = + (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; + const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv; double integral = 0.0; -#pragma omp parallel for reduction (+ : integral) +#pragma omp parallel for reduction(+ : integral) for (unsigned int i = 0; i < Nv; ++i) - integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs); - - return 1.0 - integral*dv; + integral += eval_ftilda(n, x, v_min + i * dv, poisson, phi_history); + return 1.0 - integral * dv; } -void NuFISolver::run() -{ +void NuFISolver::run() { std::cout << "Building E_sline\n\n"; using std::abs; using std::max; - const size_t stride_t = Nx + order - 1; - - std::unique_ptr coeffs { new double[ Nt*stride_t ] {} }; - std::unique_ptr rho { reinterpret_cast(std::aligned_alloc(64,sizeof(double)*Nx)), std::free }; + std::unique_ptr rho{ + reinterpret_cast(std::aligned_alloc(64, sizeof(double) * Nx)), + std::free}; std::vector int_E_squared; int_E_squared.reserve(Nt); - if ( rho == nullptr ) throw std::bad_alloc {}; + std::vector> phi_history; - Gradient grad(x_min, x_max, Nx); + if (rho == nullptr) + throw std::bad_alloc{}; double total_time = 0; + std::ofstream time_file("results/simulation_time.txt"); + time_file << "# it step_time total_time" << "\n"; - for (unsigned int it = 0; it < Nt; ++it) - { - stopwatch timer; - - double time_elapsed_before = timer.elapsed(); + const double x_min = Parameters::X_DOMAIN_LEFT; + double dx = Parameters::CALC_DX; - std::cout << "Timestep " << it << " / " << Nt << " (simulation time = "<< it*Parameters::DT << ")"<< std::endl; + for (unsigned int it = 0; it < Nt; ++it) { + stopwatch timer; - // compute rho + double time_elapsed_before = timer.elapsed(); - double dx = Parameters::SPLINE_DX; - - #pragma omp parallel for - for(size_t i = 0; i>(rho.get(), Nx)); - poisson.solve_step(); +#pragma omp parallel for + for (size_t i = 0; i < Nx; i++) { + double x = Parameters::X_DOMAIN_LEFT + i * dx; + double ith_rho = eval_rho(it, x, poisson, phi_history, Parameters::NV); - std::vector sampled_potential = poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of FE + AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho")); + rho.get()[i] = ith_rho; + } - // These have been tested to be equivalent - // //////////////////////////////////////////////// - // std::vector E_vals = grad.compute(sampled_potential); // vector grad of FE solution - // save_space_vector(E_vals, "electric", it); - // std::vector E_vals_deal = poisson.sample_electric_field(x_min, x_max, Nx); // FE grad of soution - // save_space_vector(E_vals_deal, "electricdeal", it); - // //////////////////////////////////////////////// + poisson.set_rhs_function([&rho, x_min, dx, Nx = Nx](const Point<1> &p) { + double x = p[0]; + int i = static_cast(std::floor((x - x_min) / dx)); + i = (i % Nx + Nx) % Nx; + return rho.get()[i]; + }); + poisson.solve_step(); - // interpolate and save current field - double* current_coeffs = coeffs.get() + it*stride_t; - interpolate(current_coeffs, sampled_potential.data()); + phi_history.push_back(poisson.get_solution()); + // std::vector sampled_potential = + // poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of + // FE - std::vector E_x(Nx,0.0) ; - #pragma omp parallel for - for(size_t ix=0; ix(Parameters::X_DOMAIN_LEFT+ix*dx, current_coeffs); + double timer_elapsed = timer.elapsed(); + double step_time = timer_elapsed - time_elapsed_before; + total_time += timer_elapsed; + + time_file << it << " " << step_time << " " << total_time << "\n"; + time_file.flush(); + + std::cout << "step made in " << step_time << " seconds\n\n"; + if (it % Parameters::PLOT_FREQUENCY == 0) { + std::cout << "Saving results... "; + save_f(*this, it, poisson, phi_history, Parameters::PLOT_NX, + Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat"); + save_rho(*this, it, poisson, phi_history, Parameters::PLOT_NX, + "results/rho_" + std::to_string(it) + ".dat"); + // save_Efield(it, coeffs.get(), 128, "results/field_" + + // std::to_string(it) + ".dat"); + + std::vector E_x(Nx, 0.0); +#pragma omp parallel for + for (size_t ix = 0; ix < Nx; ++ix) { + E_x[ix] = -eval(x_min + ix * dx, poisson, phi_history[it]); } - double timer_elapsed = timer.elapsed(); - total_time += timer_elapsed; - - std::cout << "step made in "<< timer_elapsed-time_elapsed_before <<" seconds\n\n"; - if (it % Parameters::PLOT_FREQUENCY == 0) - { - std::cout << "Saving results... "; - save_f(*this, it, coeffs.get(), Parameters::SPLINE_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat"); - save_rho(*this, it, coeffs.get(), Parameters::SPLINE_NX, "results/rho_" + std::to_string(it) + ".dat"); - // save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat"); - save_space_vector(E_x, "field", it); + save_space_vector(E_x, "field", it); - double int_val = 0.5 * integral_space_vector_squared(current_coeffs); - int_E_squared.push_back(int_val); - save_space_vector(int_E_squared, "electricint", it); - std::cout << "Time since start = "<< total_time<<"\n\n"; - } + double int_val = + 0.5 * integral_space_vector_squared(poisson, phi_history[it]); + int_E_squared.push_back(int_val); + save_space_vector(int_E_squared, "electricint", it); + std::cout << "Time since start = " << total_time << "\n\n"; + } } - std::cout << "NuFI simulation finished in "<< total_time <<" seconds.\n"; + std::cout << "NuFI simulation finished in " << total_time << " seconds.\n"; } -NuFISolver::NuFISolver() - : order(Parameters::FE_DEGREE), - poisson(order) -{ +NuFISolver::NuFISolver() : order(Parameters::FE_DEGREE), poisson(order) { std::cout << "Initializing dealii Poisson Solver\n"; poisson.initialize(); } diff --git a/src/save_results.cc b/src/save_results.cc index 8c6d063..3f32acd 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -1,20 +1,18 @@ #include "nufi/save_results.h" +#include "nufi/fields.h" +#include "nufi/nufi_solver.h" #include "nufi/parameters.h" +#include "nufi/poisson_problem.h" #include #include #include #include -#include "nufi/nufi_solver.h" - -void save_f( const NuFISolver &solver, - unsigned int n, - const double *E_coeffs, - unsigned int Nx_out, - unsigned int Nv_out, - const std::string &filename) -{ +void save_f(const NuFISolver &solver, unsigned int n, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history, unsigned int Nx_out, + unsigned int Nv_out, const std::string &filename) { std::ofstream file(filename); double xmin = Parameters::X_DOMAIN_LEFT; @@ -30,34 +28,30 @@ void save_f( const NuFISolver &solver, file << xmin << " " << xmax << "\n"; file << vmin << " " << vmax << "\n"; - for (unsigned int i = 0; i < Nx_out; ++i) - { - double x = xmin + (i + 0.5)*dx; + for (unsigned int i = 0; i < Nx_out; ++i) { + double x = xmin + (i + 0.5) * dx; - for (unsigned int j = 0; j < Nv_out; ++j) - { - double v = vmin + (j + 0.5)*dv; + for (unsigned int j = 0; j < Nv_out; ++j) { + double v = vmin + (j + 0.5) * dv; - double val = solver.eval_f(n, x, v, E_coeffs); + double val = solver.eval_f(n, x, v, poisson, phi_history); - file << val; + file << val; - if (j < Nv_out - 1) - file << " "; - } + if (j < Nv_out - 1) + file << " "; + } - file << "\n"; + file << "\n"; } file.close(); } -void save_rho(const NuFISolver &solver, - unsigned int n, - const double *E_coeffs, - unsigned int Nx_out, - const std::string &filename) -{ +void save_rho(const NuFISolver &solver, unsigned int n, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history, + unsigned int Nx_out, const std::string &filename) { std::ofstream file(filename); double xmin = Parameters::X_DOMAIN_LEFT; @@ -67,20 +61,18 @@ void save_rho(const NuFISolver &solver, file << Nx_out << "\n"; file << xmin << " " << xmax << "\n"; - for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) - { - double val = solver.eval_rho(n, xmin, E_coeffs); - file << val; - file << "\n"; + for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) { + double val = solver.eval_rho(n, xmin, poisson, phi_history); + file << val; + file << "\n"; } file.close(); } -void save_Efield(unsigned int n, - const double *E_coeffs, - unsigned int Nx_out, - const std::string &filename) -{ +void save_Efield([[maybe_unused]] unsigned int n, + const PoissonProblem<1> &poisson, + const std::vector> &phi_history, + unsigned int Nx_out, const std::string &filename) { std::ofstream file(filename); double xmin = Parameters::X_DOMAIN_LEFT; @@ -88,31 +80,29 @@ void save_Efield(unsigned int n, 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"; + for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) { + double val = -eval(xmin, poisson, phi_history[n]); + file << val; + file << "\n"; } file.close(); } -void save_space_vector(const std::vector& vals, const std::string& filename, size_t it) -{ - std::ofstream file("results/" + filename + "_" + std::to_string(it) + ".dat"); +void save_space_vector(const std::vector &vals, + const std::string &filename, size_t it) { + std::ofstream file("results/" + filename + "_" + std::to_string(it) + ".dat"); - if (!file) throw std::runtime_error("failed to start file in results/"); + if (!file) + throw std::runtime_error("failed to start file in results/"); - file << vals.size() << "\n"; - file << Parameters::X_DOMAIN_LEFT << " " << Parameters::X_DOMAIN_RIGHT << "\n"; - file << std::fixed << std::setprecision(8); - for (double val : vals) file << val << "\n"; + file << vals.size() << "\n"; + file << Parameters::X_DOMAIN_LEFT << " " << Parameters::X_DOMAIN_RIGHT + << "\n"; + file << std::fixed << std::setprecision(8); + for (double val : vals) + file << val << "\n"; }