diff --git a/Makefile b/Makefile index adbe8ea..73bbda3 100644 --- a/Makefile +++ b/Makefile @@ -67,10 +67,10 @@ RM = /usr/local/bin/cmake -E rm -f EQUALS = = # The top-level source directory on which CMake was run. -CMAKE_SOURCE_DIR = /home/vf/Thesis/dealII_NuFI/1x1v +CMAKE_SOURCE_DIR = /home/vf/Thesis/dealII_NuFI # The top-level build directory on which CMake was run. -CMAKE_BINARY_DIR = /home/vf/Thesis/dealII_NuFI/1x1v +CMAKE_BINARY_DIR = /home/vf/Thesis/dealII_NuFI #============================================================================= # Targets provided globally by CMake. @@ -89,7 +89,7 @@ rebuild_cache/fast: rebuild_cache # Special rule for the target edit_cache edit_cache: @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..." - /usr/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) + /usr/local/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) .PHONY : edit_cache # Special rule for the target edit_cache @@ -99,9 +99,9 @@ edit_cache/fast: edit_cache # The main all target all: cmake_check_build_system - $(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/1x1v/CMakeFiles /home/vf/Thesis/dealII_NuFI/1x1v//CMakeFiles/progress.marks + $(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/CMakeFiles /home/vf/Thesis/dealII_NuFI//CMakeFiles/progress.marks $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 all - $(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/1x1v/CMakeFiles 0 + $(CMAKE_COMMAND) -E cmake_progress_start /home/vf/Thesis/dealII_NuFI/CMakeFiles 0 .PHONY : all # The main clean target diff --git a/README.md b/README.md index 3e635be..6740fa2 100644 --- a/README.md +++ b/README.md @@ -5,9 +5,9 @@ This simulation of the Vlasov-Poisson system dimensions uses - [deal.ii](https://dealii.org/) FEM package --- -dimensions: 1x1v +dimensions: 1x2v -status: +status: working, needs testing, not optimized (at all!) Refinement working: - grid versions saved on a vector diff --git a/cmake_install.cmake b/cmake_install.cmake index ad9d58c..d973a14 100644 --- a/cmake_install.cmake +++ b/cmake_install.cmake @@ -1,4 +1,4 @@ -# Install script for directory: /home/vf/Thesis/dealII_NuFI/1x1v +# Install script for directory: /home/vf/Thesis/dealII_NuFI # Set the install prefix if(NOT DEFINED CMAKE_INSTALL_PREFIX) @@ -50,5 +50,5 @@ endif() string(REPLACE ";" "\n" CMAKE_INSTALL_MANIFEST_CONTENT "${CMAKE_INSTALL_MANIFEST_FILES}") -file(WRITE "/home/vf/Thesis/dealII_NuFI/1x1v/${CMAKE_INSTALL_MANIFEST}" +file(WRITE "/home/vf/Thesis/dealII_NuFI/${CMAKE_INSTALL_MANIFEST}" "${CMAKE_INSTALL_MANIFEST_CONTENT}") diff --git a/libnufi_lib.a b/libnufi_lib.a index b9e417c..d716ba6 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/fields.h b/nufi/fields.h index 0d75a8d..c91ca6f 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -1,7 +1,7 @@ #ifndef FIELDS_H #define FIELDS_H -#include "grids.h" +#include "nufi/grids.h" #include "nufi/parameters.h" #include "nufi/poisson_problem.h" #include @@ -21,26 +21,15 @@ inline std::vector make_x_eval(size_t Nx) { for (unsigned int i = 0; i < Nx; ++i) { x_eval_E.push_back(Parameters::X_DOMAIN_LEFT + (i + 0.5) * dx); } - // std::vector x_eval(Nx); - // const double dx = Parameters::LX / Nx; - // for (size_t i = 0; i < Nx; ++i) - // x_eval[i] = Parameters::X_DOMAIN_LEFT + i * dx; return x_eval_E; } -inline void reset_x_eval(std::vector &x_vals) { - const size_t Nx = x_vals.size(); - const double dx = Parameters::LX / Nx; - for (size_t i = 0; i < Nx; ++i) - x_vals[i] = Parameters::X_DOMAIN_LEFT + i * dx; -}; - -inline double f0(const double x, const double v, +inline double f0(const double x, const double v1, const double v2, 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 gaussian = v1 * v1 * std::exp(-0.5 * (v1 * v1 + v2 * v2)); return prefactor * gaussian; } diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h index 8aa75fc..3e30da7 100644 --- a/nufi/nufi_solver.h +++ b/nufi/nufi_solver.h @@ -24,13 +24,14 @@ public: eval_rho(unsigned int n, std::vector &x, const std::vector> &grid_struct, const std::vector> &phi_history, - const unsigned int Nv = Parameters::NV) const; + const unsigned int Nv1 = Parameters::NV_1, + const unsigned int Nv2 = Parameters::NV_2) const; std::vector - eval_ftilda(unsigned int n, std::vector x, double u, + eval_ftilda(unsigned int n, std::vector x, double u1, double u2, const std::vector> &grid_struct, const std::vector> &phi_history) const; std::vector - eval_f(unsigned int n, std::vector x, double u, + eval_f(unsigned int n, std::vector x, double u1, double u2, const std::vector> &grid_struct, const std::vector> &phi_history) const; @@ -38,7 +39,7 @@ public: eval_rho_points(unsigned int n, const std::vector> &points, const std::vector> &grid_struct, const std::vector> &phi_history, - const unsigned int Nv) const; + const unsigned int Nv1, const unsigned int Nv2) const; private: unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT); diff --git a/nufi/parameters.h b/nufi/parameters.h index 0f8926a..ce72a76 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -16,12 +16,17 @@ constexpr double LX_INV = 1 / LX; constexpr size_t CALC_NX = 256; constexpr double CALC_DX = LX / CALC_NX; -constexpr double V_DOMAIN_LEFT = -10.; -constexpr double V_DOMAIN_RIGHT = 10.; +constexpr double V_DOMAIN_LEFT_1 = -10.; +constexpr double V_DOMAIN_RIGHT_1 = 10.; -constexpr unsigned int NV = 128; -constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV; +constexpr unsigned int NV_1 = 128; +constexpr double DV_1 = std::abs(V_DOMAIN_RIGHT_1 - V_DOMAIN_LEFT_1) / NV_1; +constexpr double V_DOMAIN_LEFT_2 = -10.; +constexpr double V_DOMAIN_RIGHT_2 = 10.; + +constexpr unsigned int NV_2 = 128; +constexpr double DV_2 = std::abs(V_DOMAIN_RIGHT_2 - V_DOMAIN_LEFT_2) / NV_2; // deal.ii options constexpr unsigned int GLOBAL_REFINEMENT = 6; constexpr unsigned int FE_DEGREE = 3; @@ -45,6 +50,7 @@ constexpr unsigned int REFINE_FREQUENCY = 30; constexpr int PLOT_FREQUENCY = 10; constexpr size_t PLOT_NX = CALC_NX; constexpr double PLOT_DX = LX / PLOT_NX; +constexpr double PLOT_FIXED_V2 = 0.; const std::string PLOT_DIR = "results/"; } // namespace Parameters diff --git a/nufi/save_results.h b/nufi/save_results.h index c5f18fa..caad066 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -2,14 +2,13 @@ #define SAVE_RESULTS_H #include "nufi/grids.h" -#include "nufi/poisson_problem.h" #include #include #include class NuFISolver; -void save_f(const NuFISolver &solver, unsigned int n, +void save_f(const NuFISolver &solver, unsigned int n, const double v2_0, std::vector> &grid_struct, std::vector> &phi_history, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename); diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index 88733fd..8fde451 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -12,7 +12,6 @@ #include #include #include -#include #include #include #include @@ -27,18 +26,18 @@ using namespace dealii; std::vector NuFISolver::eval_ftilda( - unsigned int n, std::vector X, double u, + unsigned int n, std::vector X, double u1, double u2, const std::vector> &grid_struct, const std::vector> &phi_history) const { size_t x_size = X.size(); - std::vector U(x_size, u); + std::vector U1(x_size, u1); + std::vector U2(x_size, u2); std::vector results(x_size); if (n == 0) { for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U[i]); - // reset_x_eval(X); + results[i] = f0(X[i], U1[i], U2[i]); return results; } @@ -48,61 +47,53 @@ std::vector NuFISolver::eval_ftilda( // We omit the initial half-step. while (--n) { for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U[i]; - - AssertThrow( - phi_history[n].solution.size() == - grid_struct[phi_history[n].grid_version].dof_handler->n_dofs(), - ExcMessage("In eval_ftilda: Solution size = " + - std::to_string(phi_history[n].solution.size()) + - ", expected by grid_struct = " + - std::to_string(grid_struct[phi_history[n].grid_version] - .dof_handler->n_dofs()))); - AssertThrow( - grid_struct[phi_history[n].grid_version].grid_version == - phi_history[n].grid_version, - ExcMessage( - "grid.grid_version not equal to phi_history[n].grid_version")); + X[i] = X[i] - Parameters::DT * U1[i]; tmp = eval(X, grid_struct[phi_history[n].grid_version], phi_history[n].solution); // call eval only once for (size_t i = 0; i < x_size; ++i) { Ex[i] = -tmp[i]; - U[i] = U[i] + Parameters::DT * Ex[i]; + U1[i] = U1[i] + Parameters::DT * Ex[i]; + + // IMPORTANT!! LOOK at this before 1x2v, and check with Paul + // + // No change to U2 because dv1/dt = E_x + // but dv2/dt = E_y = 0 (because were in 1x not 2x) + // U2[i] = U2[i] + Parameters::DT * Ey[i]; TO_BUILD } } // The final half-step. for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U[i]; + X[i] = X[i] - Parameters::DT * U1[i]; tmp = eval(X, grid_struct[phi_history[n].grid_version], phi_history[n].solution); // call eval only once for (size_t i = 0; i < x_size; ++i) { Ex[i] = -tmp[i]; - U[i] = U[i] + 0.5 * Parameters::DT * Ex[i]; + U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i]; } for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U[i]); - // reset_x_eval(X); + results[i] = f0(X[i], U1[i], U2[i]); return results; } std::vector -NuFISolver::eval_f(unsigned int n, std::vector X, double u, +NuFISolver::eval_f(unsigned int n, std::vector X, double u1, double u2, const std::vector> &grid_struct, const std::vector> &phi_history) const { size_t x_size = X.size(); - std::vector U(x_size, u); + std::vector U1(x_size, u1); + std::vector U2(x_size, u2); + std::vector results(x_size); if (n == 0) { for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U[i]); - // reset_x_eval(X); + results[i] = f0(X[i], U1[i], U2[i]); return results; } @@ -114,35 +105,34 @@ NuFISolver::eval_f(unsigned int n, std::vector X, double u, phi_history[n].solution); // call eval only once for (size_t i = 0; i < x_size; ++i) { Ex[i] = -tmp[i]; - U[i] = U[i] + 0.5 * Parameters::DT * Ex[i]; + U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i]; } while (--n) { for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U[i]; + X[i] = X[i] - Parameters::DT * U1[i]; tmp = eval(X, grid_struct[phi_history[n].grid_version], phi_history[n].solution); // call eval only once for (size_t i = 0; i < x_size; ++i) { Ex[i] = -tmp[i]; - U[i] = U[i] + Parameters::DT * Ex[i]; + U1[i] = U1[i] + Parameters::DT * Ex[i]; } } // The final half-step. for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U[i]; + X[i] = X[i] - Parameters::DT * U1[i]; tmp = eval(X, grid_struct[phi_history[n].grid_version], phi_history[n].solution); // call eval only once for (size_t i = 0; i < x_size; ++i) { Ex[i] = -tmp[i]; - U[i] = U[i] + 0.5 * Parameters::DT * Ex[i]; + U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i]; } for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U[i]); - // reset_x_eval(X); + results[i] = f0(X[i], U1[i], U2[i]); return results; } @@ -150,26 +140,30 @@ std::vector NuFISolver::eval_rho(unsigned int n, std::vector &X, const std::vector> &grid_struct, const std::vector> &phi_history, - const unsigned int Nv) const { + const unsigned int Nv1, const unsigned int Nv2) const { size_t x_size = X.size(); - const double dv = - (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; + const double dv_0 = Parameters::DV_1; + const double v_min_0 = Parameters::V_DOMAIN_LEFT_1 + 0.5 * dv_0; - const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv; + const double dv_1 = Parameters::DV_2; + const double v_min_1 = Parameters::V_DOMAIN_LEFT_2 + 0.5 * dv_1; std::vector integral(x_size, 0.0); std::vector tmp_int(x_size); - for (unsigned int i = 0; i < Nv; ++i) { - tmp_int = eval_ftilda(n, X, v_min + i * dv, grid_struct, - phi_history); // used eval_ftilda once per i - for (size_t ii = 0; ii < x_size; ++ii) - integral[ii] += tmp_int[ii]; + for (unsigned int i = 0; i < Nv1; ++i) { + for (unsigned int j = 0; j < Nv2; ++j) { + tmp_int = + eval_ftilda(n, X, v_min_0 + i * dv_0, v_min_1 + j * dv_1, grid_struct, + phi_history); // used eval_ftilda once per i + for (size_t ii = 0; ii < x_size; ++ii) + integral[ii] += tmp_int[ii]; + } } for (size_t i = 0; i < x_size; ++i) - integral[i] = 1 - integral[i] * dv; + integral[i] = 1 - integral[i] * dv_0; return integral; } @@ -177,9 +171,11 @@ std::vector NuFISolver::eval_rho_points(unsigned int n, const std::vector> &points, const std::vector> &grid_struct, const std::vector> &phi_history, - const unsigned int Nv) const { + const unsigned int Nv1, + const unsigned int Nv2) const { std::vector point_vector = Point_vector_to_double_vector(points); - return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv); + return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv1, + Nv2); } void NuFISolver::run() { @@ -280,7 +276,8 @@ void NuFISolver::run() { for (size_t i = 0; i < points.size(); ++i) x[i] = points[i][0]; - return eval_rho(it, x, grid_versions, phi_history, Parameters::NV); + return eval_rho(it, x, grid_versions, phi_history, Parameters::NV_1, + Parameters::NV_2); }); if (it % Parameters::REFINE_FREQUENCY == 0) { @@ -306,8 +303,9 @@ void NuFISolver::run() { double plot_start = timer.elapsed(); std::cout << "Saving results... "; - save_f(*this, it, grid_versions, phi_history, Parameters::PLOT_NX, - Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat"); + save_f(*this, it, Parameters::PLOT_FIXED_V2, grid_versions, phi_history, + Parameters::PLOT_NX, Parameters::NV_1, + "results/ftilda_" + std::to_string(it) + ".dat"); save_rho(*this, it, grid_versions, phi_history, Parameters::PLOT_NX, "results/rho_" + std::to_string(it) + ".dat"); diff --git a/src/save_results.cc b/src/save_results.cc index 09ca5ca..bcc7559 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -12,7 +12,7 @@ #include #include -void save_f(const NuFISolver &solver, unsigned int n, +void save_f(const NuFISolver &solver, unsigned int n, const double v2_0, std::vector> &grid_struct, std::vector> &phi_history, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename) { @@ -21,21 +21,21 @@ void save_f(const NuFISolver &solver, unsigned int n, double xmin = Parameters::X_DOMAIN_LEFT; double xmax = Parameters::X_DOMAIN_RIGHT; - double vmin = Parameters::V_DOMAIN_LEFT; - double vmax = Parameters::V_DOMAIN_RIGHT; + double vmin_1 = Parameters::V_DOMAIN_LEFT_1; + double vmax_1 = Parameters::V_DOMAIN_RIGHT_1; - double dv = (vmax - vmin) / Nv_out; + double dv_1 = (vmax_1 - vmin_1) / Nv_out; file << Nx_out << " " << Nv_out << "\n"; file << xmin << " " << xmax << "\n"; - file << vmin << " " << vmax << "\n"; + file << vmin_1 << " " << vmax_1 << "\n"; std::vector x_eval = make_x_eval(Nx_out); std::vector val(Nx_out); for (unsigned int j = 0; j < Nv_out; ++j) { - double v = vmin + (j + 0.5) * dv; - val = solver.eval_f(n, x_eval, v, grid_struct, phi_history); + double v = vmin_1 + (j + 0.5) * dv_1; + val = solver.eval_f(n, x_eval, v, v2_0, grid_struct, phi_history); for (unsigned int i = 0; i < Nx_out; ++i) { file << val[i];