diff --git a/.gitignore b/.gitignore index 59e2741..05e4185 100644 --- a/.gitignore +++ b/.gitignore @@ -3,5 +3,6 @@ CMakeCache.txt compile_commands.json CMakeFiles/ results/ +saved_sims/ nufi_poisson *.ipynb diff --git a/libnufi_lib.a b/libnufi_lib.a index adc4573..7b25e7c 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/fields.h b/nufi/fields.h index 2087f43..8d4e260 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -47,19 +47,19 @@ double eval(double x, const double *coeffs) noexcept x -= Parameters::X_DOMAIN_LEFT; // Get "periodic position" in box at origin. - x = x - Parameters::LX * floor( x/Parameters::LX ); + x = x - Parameters::LX * floor( x*Parameters::LX_INV ); // Knot number - double x_knot = floor( x/Parameters::SPLINE_DX); + 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 - x_knot; + 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; + for ( size_t i = 0; i < dx; ++i ) factor *= 1*Parameters::SPLINE_DX_INV; return factor*splines1d::eval(x, coeffs + ii); } @@ -142,7 +142,7 @@ void interpolate( real *coeffs, const real *values) coeffs[ i ] = tmp[ i % Parameters::SPLINE_NX ]; } -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 double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX) { double integral = 0.0; double x = Parameters::X_DOMAIN_LEFT; @@ -153,7 +153,7 @@ double integral_space_vector(const double *current_coeffs, double dx = Parameter return integral; }; -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 double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX) { double integral = 0.0; double x = Parameters::X_DOMAIN_LEFT; diff --git a/nufi/parameters.h b/nufi/parameters.h index 4dc598c..599c953 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -11,11 +11,13 @@ namespace Parameters 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.0; constexpr double V_DOMAIN_RIGHT = 10.0; - constexpr unsigned int NV = 562; + 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; @@ -25,15 +27,16 @@ namespace Parameters constexpr double EPS = 0.01; constexpr double WAVE_NR = 0.5; - constexpr double F0_FACTOR = 0.39894228040143267793994; + constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) // NUFI options - constexpr double DT=1./16.; - constexpr unsigned int TMAX = 20; + constexpr double DT=1./10.; + constexpr unsigned int TMAX = 50; //spline options - constexpr int SPLINE_NX = 562; - constexpr double SPLINE_DX = LX/(SPLINE_NX-1); + 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; //Plotting options diff --git a/nufi/splines.h b/nufi/splines.h index 8b1bd4d..73e8fd0 100644 --- a/nufi/splines.h +++ b/nufi/splines.h @@ -1,5 +1,5 @@ -#ifndef SPLINES_H -#define SPLINES_H +#ifndef SPLINES_HP +#define SPLINES_HP #include diff --git a/rho_E.mp4 b/rho_E.mp4 index 0f4af77..c907db0 100644 Binary files a/rho_E.mp4 and b/rho_E.mp4 differ diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index fd020a2..4d076e3 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -26,9 +26,10 @@ double NuFISolver::eval_ftilda(unsigned int n, 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 + Parameters::SPLINE_ORDER - 1); + const size_t stride_t = stride_x*(Nx + order - 1); double Ex; const double *c; @@ -107,26 +108,32 @@ void NuFISolver::run() std::vector sampled_potential = poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of FE - save_space_vector(sampled_potential, "potential", it); - // 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); + // 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); // //////////////////////////////////////////////// - double* current_coeffs = coeffs.get() + it*stride_t; - interpolate(current_coeffs, E_vals.data()); + // interpolate and save current field + double* current_coeffs = coeffs.get() + it*stride_t; + interpolate(current_coeffs, sampled_potential.data()); + + std::vector E_x(Nx,0.0) ; + for(size_t ix=0; ix(Parameters::X_DOMAIN_LEFT+ix*dx, current_coeffs); + } if (it % Parameters::PLOT_FREQUENCY == 0) { std::cout << "Saving results... \n\n"; 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_Efield(it, coeffs.get(), 128, "results/field_" + 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); double int_val = 0.5 * integral_space_vector_squared(current_coeffs); int_E_squared.push_back(int_val);