diff --git a/CMakeLists.txt b/CMakeLists.txt index 60ddc93..54666d7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,6 +16,8 @@ endif() deal_ii_initialize_cached_variables() +find_package(OpenMP REQUIRED) + # ------------------------- add_library(nufi_lib @@ -30,6 +32,10 @@ target_include_directories(nufi_lib PUBLIC deal_ii_setup_target(nufi_lib) +target_link_libraries(nufi_lib + OpenMP::OpenMP_CXX +) + # ------------------------- add_executable(nufi_poisson @@ -38,6 +44,8 @@ add_executable(nufi_poisson target_link_libraries(nufi_poisson nufi_lib + OpenMP::OpenMP_CXX ) + deal_ii_setup_target(nufi_poisson) diff --git a/libnufi_lib.a b/libnufi_lib.a index 00e4ef9..14ea604 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/fields.h b/nufi/fields.h index 8d4e260..18e5e58 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -83,6 +83,7 @@ void interpolate( real *coeffs, const real *values) 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; @@ -145,24 +146,26 @@ void interpolate( real *coeffs, const real *values) inline double integral_space_vector(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX) { double integral = 0.0; - double x = Parameters::X_DOMAIN_LEFT; + double xmin = Parameters::X_DOMAIN_LEFT; + #pragma omp parallel for reduction (+:integral) for (size_t i=0; i(x, current_coeffs)*dx; + double x = xmin + i * dx; + integral += eval<1>(x, current_coeffs); } - return integral; + return integral*dx; }; inline double integral_space_vector_squared(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX) { double integral = 0.0; - double x = Parameters::X_DOMAIN_LEFT; + double xmin = Parameters::X_DOMAIN_LEFT; + #pragma omp parallel for reduction (+:integral) for (size_t i=0; i(x, current_coeffs); - integral += val*val*dx; + integral += val*val; } - return integral; + return integral*dx; }; class Gradient { diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h index 3d8669a..a023895 100644 --- a/nufi/nufi_solver.h +++ b/nufi/nufi_solver.h @@ -22,6 +22,7 @@ public: 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; private: diff --git a/nufi/parameters.h b/nufi/parameters.h index 599c953..97618e6 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -13,16 +13,16 @@ namespace Parameters constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT); constexpr double LX_INV = 1/LX; - constexpr double V_DOMAIN_LEFT = -10.0; - constexpr double V_DOMAIN_RIGHT = 10.0; + constexpr double V_DOMAIN_LEFT = -10.; + constexpr double V_DOMAIN_RIGHT = 10.; constexpr unsigned int NV = 512; constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT)/NV; // deal.ii options - constexpr unsigned int GLOBAL_REFINEMENT = 7; + constexpr unsigned int GLOBAL_REFINEMENT = 8; constexpr unsigned int FE_DEGREE = 4; - constexpr unsigned int CONVERGENCE_ITERATIONS = 20000; + constexpr unsigned int CONVERGENCE_ITERATIONS = 10000; constexpr double CONVERGENCE_LIMIT = 1e-12; constexpr double EPS = 0.01; @@ -30,8 +30,8 @@ namespace Parameters constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) // NUFI options - constexpr double DT=1./10.; - constexpr unsigned int TMAX = 50; + constexpr double DT=1./16.; + constexpr unsigned int TMAX = 500; //spline options constexpr int SPLINE_NX = 256; @@ -40,7 +40,7 @@ namespace Parameters constexpr size_t SPLINE_ORDER = 4; //Plotting options - constexpr int PLOT_FREQUENCY = 5; + constexpr int PLOT_FREQUENCY = 10; } #endif diff --git a/nufi/save_results.h b/nufi/save_results.h index 26b13af..2874edc 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -5,7 +5,7 @@ #include "nufi/nufi_solver.h" -void save_ftilda( const NuFISolver &solver, +void save_f( const NuFISolver &solver, unsigned int n, const double *E_coeffs, unsigned int Nx_out, diff --git a/rho_E.mp4 b/rho_E.mp4 index c907db0..2564d57 100644 Binary files a/rho_E.mp4 and b/rho_E.mp4 differ diff --git a/src/main.cc b/src/main.cc index 9ea40c4..7e18060 100644 --- a/src/main.cc +++ b/src/main.cc @@ -20,6 +20,8 @@ void clear_results_directory(const std::string &dir) int main() { + #include +std::cout << "Threads: " << omp_get_max_threads() << "\n"; try { clear_results_directory("results"); diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index d97c628..83ca6a1 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -54,6 +54,42 @@ double NuFISolver::eval_ftilda(unsigned int n, 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 Ex; + const double *c; + + // Initial half-step. + c = E_coeffs + n*stride_t; + Ex = -eval<1>(x, c); + 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; + } + + // 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; + + return f0(x,u); +} + double NuFISolver::eval_rho(const unsigned int n, const double x, const double *E_coeffs, @@ -63,10 +99,12 @@ double NuFISolver::eval_rho(const unsigned int n, const double v_min = Parameters::V_DOMAIN_LEFT; double integral = 0.0; + +#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) * dv; + integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs); - return 1.0 - integral; + return 1.0 - integral*dv; } void NuFISolver::run() @@ -93,17 +131,21 @@ void NuFISolver::run() for (unsigned int it = 0; it < Nt; ++it) { stopwatch timer; + + double time_elapsed_before = timer.elapsed(); - std::cout << "Timestep " << it << " / " << Nt << std::endl << std::endl; + std::cout << "Timestep " << it << " / " << Nt << " (simulation time = "<< it*Parameters::DT << ")"<< std::endl; // compute rho double dx = Parameters::SPLINE_DX; - double x = Parameters::X_DOMAIN_LEFT; - - for(size_t i = 0; i(current_coeffs, sampled_potential.data()); + 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); @@ -134,12 +178,13 @@ void NuFISolver::run() 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_ftilda(*this, it, coeffs.get(), 128, 128, "results/ftilda_" + std::to_string(it) + ".dat"); - save_rho(*this, it, coeffs.get(), 128, "results/rho_" + std::to_string(it) + ".dat"); + 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); diff --git a/src/save_results.cc b/src/save_results.cc index dfb8f7d..8c6d063 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -8,7 +8,7 @@ #include "nufi/nufi_solver.h" -void save_ftilda( const NuFISolver &solver, +void save_f( const NuFISolver &solver, unsigned int n, const double *E_coeffs, unsigned int Nx_out, @@ -38,7 +38,7 @@ void save_ftilda( const NuFISolver &solver, { double v = vmin + (j + 0.5)*dv; - double val = solver.eval_ftilda(n, x, v, E_coeffs); + double val = solver.eval_f(n, x, v, E_coeffs); file << val;