diff --git a/.clangd b/.clangd new file mode 100644 index 0000000..ec0b890 --- /dev/null +++ b/.clangd @@ -0,0 +1,2 @@ +CompileFlags: + CompilationDatabase: build diff --git a/.gitignore b/.gitignore index 6d7ee76..bf17354 100644 --- a/.gitignore +++ b/.gitignore @@ -1,10 +1,11 @@ -.ipynb_checkpoints/ -CMakeCache.txt -compile_commands.json -CMakeFiles/ +build/ + results/ + saved_sims/ -nufi_poisson + +.ipynb_checkpoints/ *.ipynb + .todo-txt/ Notes/ diff --git a/Makefile b/Makefile deleted file mode 100644 index 73bbda3..0000000 --- a/Makefile +++ /dev/null @@ -1,271 +0,0 @@ -# CMAKE generated file: DO NOT EDIT! -# Generated by "Unix Makefiles" Generator, CMake Version 3.18 - -# Default target executed when no arguments are given to make. -default_target: all - -.PHONY : default_target - -# Allow only one "make -f Makefile2" at a time, but pass parallelism. -.NOTPARALLEL: - - -#============================================================================= -# Special targets provided by cmake. - -# Disable implicit rules so canonical targets will work. -.SUFFIXES: - - -# Disable VCS-based implicit rules. -% : %,v - - -# Disable VCS-based implicit rules. -% : RCS/% - - -# Disable VCS-based implicit rules. -% : RCS/%,v - - -# Disable VCS-based implicit rules. -% : SCCS/s.% - - -# Disable VCS-based implicit rules. -% : s.% - - -.SUFFIXES: .hpux_make_needs_suffix_list - - -# Command-line flag to silence nested $(MAKE). -$(VERBOSE)MAKESILENT = -s - -#Suppress display of executed commands. -$(VERBOSE).SILENT: - -# A target that is always out of date. -cmake_force: - -.PHONY : cmake_force - -#============================================================================= -# Set environment variables for the build. - -# The shell in which to execute make rules. -SHELL = /bin/sh - -# The CMake executable. -CMAKE_COMMAND = /usr/local/bin/cmake - -# The command to remove a file. -RM = /usr/local/bin/cmake -E rm -f - -# Escaping for special characters. -EQUALS = = - -# The top-level source directory on which CMake was run. -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 - -#============================================================================= -# Targets provided globally by CMake. - -# Special rule for the target rebuild_cache -rebuild_cache: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..." - /usr/local/bin/cmake --regenerate-during-build -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) -.PHONY : rebuild_cache - -# Special rule for the target rebuild_cache -rebuild_cache/fast: rebuild_cache - -.PHONY : rebuild_cache/fast - -# Special rule for the target edit_cache -edit_cache: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..." - /usr/local/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) -.PHONY : edit_cache - -# Special rule for the target edit_cache -edit_cache/fast: edit_cache - -.PHONY : edit_cache/fast - -# The main all target -all: cmake_check_build_system - $(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/CMakeFiles 0 -.PHONY : all - -# The main clean target -clean: - $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 clean -.PHONY : clean - -# The main clean target -clean/fast: clean - -.PHONY : clean/fast - -# Prepare targets for installation. -preinstall: all - $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 preinstall -.PHONY : preinstall - -# Prepare targets for installation. -preinstall/fast: - $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 preinstall -.PHONY : preinstall/fast - -# clear depends -depend: - $(CMAKE_COMMAND) -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 1 -.PHONY : depend - -#============================================================================= -# Target rules for targets named nufi_poisson - -# Build rule for target. -nufi_poisson: cmake_check_build_system - $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 nufi_poisson -.PHONY : nufi_poisson - -# fast build rule for target. -nufi_poisson/fast: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_poisson.dir/build.make CMakeFiles/nufi_poisson.dir/build -.PHONY : nufi_poisson/fast - -#============================================================================= -# Target rules for targets named nufi_lib - -# Build rule for target. -nufi_lib: cmake_check_build_system - $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 nufi_lib -.PHONY : nufi_lib - -# fast build rule for target. -nufi_lib/fast: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/build -.PHONY : nufi_lib/fast - -src/main.o: src/main.cc.o - -.PHONY : src/main.o - -# target to build an object file -src/main.cc.o: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_poisson.dir/build.make CMakeFiles/nufi_poisson.dir/src/main.cc.o -.PHONY : src/main.cc.o - -src/main.i: src/main.cc.i - -.PHONY : src/main.i - -# target to preprocess a source file -src/main.cc.i: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_poisson.dir/build.make CMakeFiles/nufi_poisson.dir/src/main.cc.i -.PHONY : src/main.cc.i - -src/main.s: src/main.cc.s - -.PHONY : src/main.s - -# target to generate assembly for a file -src/main.cc.s: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_poisson.dir/build.make CMakeFiles/nufi_poisson.dir/src/main.cc.s -.PHONY : src/main.cc.s - -src/nufi_solver.o: src/nufi_solver.cc.o - -.PHONY : src/nufi_solver.o - -# target to build an object file -src/nufi_solver.cc.o: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/nufi_solver.cc.o -.PHONY : src/nufi_solver.cc.o - -src/nufi_solver.i: src/nufi_solver.cc.i - -.PHONY : src/nufi_solver.i - -# target to preprocess a source file -src/nufi_solver.cc.i: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/nufi_solver.cc.i -.PHONY : src/nufi_solver.cc.i - -src/nufi_solver.s: src/nufi_solver.cc.s - -.PHONY : src/nufi_solver.s - -# target to generate assembly for a file -src/nufi_solver.cc.s: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/nufi_solver.cc.s -.PHONY : src/nufi_solver.cc.s - -src/save_results.o: src/save_results.cc.o - -.PHONY : src/save_results.o - -# target to build an object file -src/save_results.cc.o: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/save_results.cc.o -.PHONY : src/save_results.cc.o - -src/save_results.i: src/save_results.cc.i - -.PHONY : src/save_results.i - -# target to preprocess a source file -src/save_results.cc.i: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/save_results.cc.i -.PHONY : src/save_results.cc.i - -src/save_results.s: src/save_results.cc.s - -.PHONY : src/save_results.s - -# target to generate assembly for a file -src/save_results.cc.s: - $(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/save_results.cc.s -.PHONY : src/save_results.cc.s - -# Help Target -help: - @echo "The following are some of the valid targets for this Makefile:" - @echo "... all (the default if no target is provided)" - @echo "... clean" - @echo "... depend" - @echo "... edit_cache" - @echo "... rebuild_cache" - @echo "... nufi_lib" - @echo "... nufi_poisson" - @echo "... src/main.o" - @echo "... src/main.i" - @echo "... src/main.s" - @echo "... src/nufi_solver.o" - @echo "... src/nufi_solver.i" - @echo "... src/nufi_solver.s" - @echo "... src/save_results.o" - @echo "... src/save_results.i" - @echo "... src/save_results.s" -.PHONY : help - - - -#============================================================================= -# Special targets to cleanup operation of make. - -# Special rule to run CMake to check the build system integrity. -# No rule that depends on this can have commands that come from listfiles -# because they might be regenerated. -cmake_check_build_system: - $(CMAKE_COMMAND) -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 0 -.PHONY : cmake_check_build_system - diff --git a/README.md b/README.md index 6740fa2..a164c93 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,28 @@ -# Vlasov-Poisson model solver +# Vlasov-Poisson model solver This simulation of the Vlasov-Poisson system dimensions uses + - [NuFI algorithm](https://doi.org/10.1002/pamm.202300162) - [deal.ii](https://dealii.org/) FEM package ---- -dimensions: 1x2v +______________________________________________________________________ -status: working, needs testing, not optimized (at all!) +dimensions: 1x1v + +notes: + +- Locator not optimized for 1d. +- works for higher dimensions +- locator up to O(dim * Max_depth) +- + +status: Working Refinement working: + - grid versions saved on a vector - solutions point to a version of the grid + +todo: + +- add ions diff --git a/cmake_install.cmake b/cmake_install.cmake deleted file mode 100644 index d973a14..0000000 --- a/cmake_install.cmake +++ /dev/null @@ -1,54 +0,0 @@ -# Install script for directory: /home/vf/Thesis/dealII_NuFI - -# Set the install prefix -if(NOT DEFINED CMAKE_INSTALL_PREFIX) - set(CMAKE_INSTALL_PREFIX "/usr/local") -endif() -string(REGEX REPLACE "/$" "" CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}") - -# Set the install configuration name. -if(NOT DEFINED CMAKE_INSTALL_CONFIG_NAME) - if(BUILD_TYPE) - string(REGEX REPLACE "^[^A-Za-z0-9_]+" "" - CMAKE_INSTALL_CONFIG_NAME "${BUILD_TYPE}") - else() - set(CMAKE_INSTALL_CONFIG_NAME "Debug") - endif() - message(STATUS "Install configuration: \"${CMAKE_INSTALL_CONFIG_NAME}\"") -endif() - -# Set the component getting installed. -if(NOT CMAKE_INSTALL_COMPONENT) - if(COMPONENT) - message(STATUS "Install component: \"${COMPONENT}\"") - set(CMAKE_INSTALL_COMPONENT "${COMPONENT}") - else() - set(CMAKE_INSTALL_COMPONENT) - endif() -endif() - -# Install shared libraries without execute permission? -if(NOT DEFINED CMAKE_INSTALL_SO_NO_EXE) - set(CMAKE_INSTALL_SO_NO_EXE "1") -endif() - -# Is this installation the result of a crosscompile? -if(NOT DEFINED CMAKE_CROSSCOMPILING) - set(CMAKE_CROSSCOMPILING "FALSE") -endif() - -# Set default install directory permissions. -if(NOT DEFINED CMAKE_OBJDUMP) - set(CMAKE_OBJDUMP "/usr/bin/objdump") -endif() - -if(CMAKE_INSTALL_COMPONENT) - set(CMAKE_INSTALL_MANIFEST "install_manifest_${CMAKE_INSTALL_COMPONENT}.txt") -else() - set(CMAKE_INSTALL_MANIFEST "install_manifest.txt") -endif() - -string(REPLACE ";" "\n" CMAKE_INSTALL_MANIFEST_CONTENT - "${CMAKE_INSTALL_MANIFEST_FILES}") -file(WRITE "/home/vf/Thesis/dealII_NuFI/${CMAKE_INSTALL_MANIFEST}" - "${CMAKE_INSTALL_MANIFEST_CONTENT}") diff --git a/libnufi_lib.a b/libnufi_lib.a deleted file mode 100644 index d716ba6..0000000 Binary files a/libnufi_lib.a and /dev/null differ diff --git a/nufi/cells.h b/nufi/cells.h index 19ee0c9..9f9d548 100644 --- a/nufi/cells.h +++ b/nufi/cells.h @@ -1,117 +1,146 @@ #ifndef CELLS_H #define CELLS_H -#include -#include +#include +#include #include #include #include -#include #include #include +#include "nufi/parameters.h" + using namespace dealii; -template struct CellInfo { - // what needs to be given to evaluator - typename DoFHandler::active_cell_iterator cell; - // usefull for locator - Point lower; - Point upper; - double h; -}; - template struct CellLocation { - const CellInfo *info; + typename DoFHandler::active_cell_iterator cell; Point reference_point; }; template class CellLocator { public: - using CellIterator = typename DoFHandler::active_cell_iterator; - void rebuild(const DoFHandler &dof_handler, const Triangulation &triangulation); + CellLocation locate(const Point &p) const; - const std::vector> &get_cell_centers() const; - private: - std::vector> cells; - std::vector> cell_centers; + const DoFHandler *dof_handler_ptr = nullptr; + Point lower; + Point upper; + + // Cached base level = Parameters::GLOBAL_REFINEMENT. + unsigned int base_level = 0; + unsigned int base_n_per_axis = 1; // 2^base_level + + std::vector::cell_iterator> base_cells; }; template void CellLocator::rebuild(const DoFHandler &dof_handler, const Triangulation &triangulation) { + dof_handler_ptr = &dof_handler; - cells.clear(); - cells.reserve(triangulation.n_active_cells()); + AssertThrow(triangulation.n_cells(0) == 1, + ExcMessage("CellLocator assumes exactly one coarse/root cell.")); - for (const auto &cell : dof_handler.active_cell_iterators()) { - CellInfo info; + typename Triangulation::cell_iterator root = triangulation.begin(0); + lower = root->vertex(0); + upper = root->vertex(GeometryInfo::vertices_per_cell - 1); - info.cell = cell; - info.lower = cell->vertex(0); - info.upper = cell->vertex(GeometryInfo::vertices_per_cell - 1); + base_level = Parameters::GLOBAL_REFINEMENT; + base_n_per_axis = 1u << base_level; // = 2^base_level - info.h = info.upper[0] - info.lower[0]; + const unsigned int n_base_cells = + 1u << (dim * base_level); // = 2^(dim*base_level) + base_cells.assign(n_base_cells, typename Triangulation::cell_iterator()); - cells.push_back(info); - } + std::vector::cell_iterator> stack; + std::vector index_stack; + std::vector depth_stack; + stack.push_back(root); + index_stack.push_back(0); + depth_stack.push_back(0); - std::sort(cells.begin(), cells.end(), - [](const CellInfo &a, const CellInfo &b) { - return a.lower[0] < b.lower[0]; - }); + while (!stack.empty()) { + auto cell = stack.back(); + unsigned int idx = index_stack.back(); + unsigned int depth = depth_stack.back(); + stack.pop_back(); + index_stack.pop_back(); + depth_stack.pop_back(); - cell_centers.clear(); - cell_centers.reserve(cells.size()); + if (depth == base_level) { + base_cells[idx] = cell; + continue; + } - for (const auto &cell : cells) { - Point center; - for (unsigned int d = 0; d < dim; ++d) - center[d] = 0.5 * (cell.lower[d] + cell.upper[d]); + AssertThrow(cell->has_children(), + ExcMessage("CellLocator: mesh is not uniformly refined to " + "Parameters::GLOBAL_REFINEMENT; base-level cache " + "cannot be built. Did you coarsen below the " + "global refinement level?")); - cell_centers.push_back(center); + const unsigned int n_children = GeometryInfo::max_children_per_cell; + for (unsigned int c = 0; c < n_children; ++c) { + stack.push_back(cell->child(c)); + index_stack.push_back(idx * n_children + c); + depth_stack.push_back(depth + 1); + } } } template CellLocation CellLocator::locate(const Point &p) const { - static_assert(dim == 1, - "Current CellLocator implementation only supports 1D."); - - AssertThrow(!cells.empty(), + AssertThrow(dof_handler_ptr != nullptr, ExcMessage("CellLocator::rebuild() has not been called.")); - const double x = p[0]; + Point p_wrapped; + for (unsigned int d = 0; d < dim; ++d) { + const double L = upper[d] - lower[d]; + double x = p[d] - lower[d]; + x = x - L * std::floor(x / L); + p_wrapped[d] = lower[d] + x; + } - auto it = std::upper_bound(cells.begin(), cells.end(), x, - [](double value, const CellInfo &cell) { - return value < cell.lower[0]; - }); // returns cell to the right of cell with x + Point xi; + for (unsigned int d = 0; d < dim; ++d) { + xi[d] = (p_wrapped[d] - lower[d]) / (upper[d] - lower[d]); + xi[d] = std::min(std::max(xi[d], 0.0), 1.0); // of cell at base level + } - if (it == cells.begin()) - it = cells.begin(); - else - --it; + unsigned int idx = 0; + Point xi_local = xi; + for (unsigned int l = 0; l < base_level; ++l) { + const unsigned int child_index = + GeometryInfo::child_cell_from_point(xi_local); + xi_local = + GeometryInfo::cell_to_child_coordinates(xi_local, child_index); + idx = idx * GeometryInfo::max_children_per_cell + child_index; + } - // Safety check: make sure the point is really inside this cell - AssertThrow(x >= it->lower[0] - 1e-12 && x <= it->upper[0] + 1e-12, - ExcMessage("CellLocator failed to find containing cell.")); + typename Triangulation::cell_iterator cell = base_cells[idx]; + xi = xi_local; + + // Step 4: continue descending only through ADAPTIVE refinement beyond + // the base level -- this loop now only runs `depth - base_level` times + // instead of `depth` times. + while (cell->has_children()) { + const unsigned int child_index = + GeometryInfo::child_cell_from_point(xi); + xi = GeometryInfo::cell_to_child_coordinates(xi, child_index); + cell = cell->child(child_index); + } + + typename DoFHandler::active_cell_iterator dof_cell( + &cell->get_triangulation(), cell->level(), cell->index(), + dof_handler_ptr); CellLocation location; - - location.info = &(*it); - location.reference_point[0] = (p[0] - it->lower[0]) / it->h; - + location.cell = dof_cell; + location.reference_point = xi; return location; } -template -const std::vector> &CellLocator::get_cell_centers() const { - return cell_centers; -} - #endif // !CELLS_H diff --git a/nufi/fields.h b/nufi/fields.h index c91ca6f..2693c5d 100644 --- a/nufi/fields.h +++ b/nufi/fields.h @@ -1,9 +1,9 @@ -#ifndef FIELDS_H -#define FIELDS_H +#ifndef NUFI_FIELDS_H_ +#define NUFI_FIELDS_H_ #include "nufi/grids.h" #include "nufi/parameters.h" -#include "nufi/poisson_problem.h" + #include #include #include @@ -13,25 +13,76 @@ using namespace dealii; inline std::vector make_x_eval(size_t Nx) { - std::vector x_eval_E; double dx = (Parameters::X_DOMAIN_RIGHT - Parameters::X_DOMAIN_LEFT) / Nx; - for (unsigned int i = 0; i < Nx; ++i) { + for (unsigned int i = 0; i < Nx; ++i) x_eval_E.push_back(Parameters::X_DOMAIN_LEFT + (i + 0.5) * dx); - } + return x_eval_E; } -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 = v1 * v1 * std::exp(-0.5 * (v1 * v1 + v2 * v2)); +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; +}; - return prefactor * gaussian; +inline double f0(const double x, const double v, + const size_t f0_type = Parameters::f0_TYPE) { + const double eps = Parameters::EPS; + const double k = Parameters::WAVE_NR; + + const double factor = Parameters::F0_FACTOR; + + auto maxwell = [factor](double u, double v = 0, double v_th = 1) { + return 1 / v_th * factor * + std::exp(-0.5 * (u - v) * (u - v) / (v_th * v_th)); + }; + + double prefactor; + double computed_max; + double result; + + switch (f0_type) { + case 0: // two-stream + { + computed_max = maxwell(v); + prefactor = (1.0 + eps * std::cos(k * x)) * v * v; + result = prefactor * computed_max; + break; + } + case 1: // Landau-damping + { + computed_max = maxwell(v); + prefactor = (1.0 + eps * std::cos(k * x)); + result = prefactor * computed_max; + break; + } + case 2: // Maxwellian + { + result = maxwell(v); + break; + } + case 3: // Bump-on tail + { + const double beam_density = 0.05; + const double beam_v_th = 0.2; + const double beam_v = 3; + const double alpha = beam_density / (1 - beam_density); + + const double beam_max = maxwell(v, beam_v, beam_v_th); + computed_max = maxwell(v); + result = (1 - alpha) * computed_max + alpha * beam_max; + break; + } + default: + throw std::invalid_argument("Invalid f0_type"); + } + + return result; } // wrapper for eval_point() { VectorTools::point_values() } @@ -42,50 +93,41 @@ inline std::vector eval(std::vector &X, AssertThrow(grid.dof_handler->n_dofs() == solution.size(), ExcMessage("@ eval(...) grid's number of DoFs doesn't correspond " "to solution's size")); - size_t x_size = X.size(); - std::vector> Points(x_size); - for (size_t i = 0; i < x_size; ++i) { - X[i] = X[i] - Parameters::X_DOMAIN_LEFT; - X[i] = X[i] - Parameters::LX * std::floor(X[i] * Parameters::LX_INV); + const size_t x_size = X.size(); + std::vector> points(x_size); - Points[i][0] = X[i]; - } + for (size_t i = 0; i < x_size; ++i) + points[i][0] = X[i]; - return grid.eval_vector_grad(solution, Points); + return grid.eval_vector_grad(solution, points); } inline double integral_space_vector(const GridStructure<1> &grid, const Vector &solution, - double dx = Parameters::PLOT_DX, size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; - double xmin = Parameters::X_DOMAIN_LEFT; - std::vector x_eval(Nx); - for (size_t i = 0; i < Nx; ++i) - x_eval[i] = xmin + i * dx; + std::vector x_eval = make_x_eval(Nx); + const double dx = Parameters::LX / Nx; std::vector tmp = eval(x_eval, grid, solution); for (size_t i = 0; i < Nx; ++i) integral += tmp[i]; return integral * dx; -}; +} inline double integral_space_vector_squared(const GridStructure<1> &grid, const Vector &solution, - double dx = Parameters::PLOT_DX, size_t Nx = Parameters::PLOT_NX) { double integral = 0.0; - double xmin = Parameters::X_DOMAIN_LEFT; - std::vector x_eval(Nx); - for (size_t i = 0; i < Nx; ++i) - x_eval[i] = xmin + i * dx; + std::vector x_eval = make_x_eval(Nx); + const double dx = Parameters::LX / Nx; std::vector tmp = eval(x_eval, grid, solution); for (size_t i = 0; i < Nx; ++i) integral += tmp[i] * tmp[i]; return integral * dx; -}; +} inline std::vector Point_vector_to_double_vector(const std::vector> &Points) { @@ -97,4 +139,4 @@ Point_vector_to_double_vector(const std::vector> &Points) { return vector; } -#endif +#endif // NUFI_FIELDS_H_ diff --git a/nufi/grids.h b/nufi/grids.h index d65fc90..65ce753 100644 --- a/nufi/grids.h +++ b/nufi/grids.h @@ -72,12 +72,12 @@ template struct GridStructure { const auto cell_location = locator->locate(points[p]); - cell_location.info->cell->get_dof_values(solution, - local_solution_buffer.begin(), - local_solution_buffer.end()); + cell_location.cell->get_dof_values(solution, + local_solution_buffer.begin(), + local_solution_buffer.end()); evaluator.reinit( - cell_location.info->cell, + cell_location.cell, ArrayView>(&cell_location.reference_point, 1)); evaluator.evaluate(local_solution_buffer, EvaluationFlags::gradients); @@ -127,7 +127,7 @@ GridStructure make_grid_snapshot(PoissonProblem &poisson) { if (PRINT_GAUGE_DOF_POSITION) std::cout << " gauge_dof = " << gauge_dof - << " gauge_point = " << point[0] << std::endl; + << " gauge_point = " << point[0] << "\n"; break; } } diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h index 3e30da7..c55aa12 100644 --- a/nufi/nufi_solver.h +++ b/nufi/nufi_solver.h @@ -8,7 +8,7 @@ #include #include -#include "nufi/fields.h" //dont remove +#include "nufi/fields.h" // dont remove #include "nufi/grids.h" #include "nufi/parameters.h" #include "nufi/poisson_problem.h" diff --git a/nufi/parameters.h b/nufi/parameters.h index ce72a76..cd8eeea 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -2,6 +2,7 @@ #define PARAMETERS_H #include +#include #include #include @@ -16,23 +17,30 @@ constexpr double LX_INV = 1 / LX; constexpr size_t CALC_NX = 256; constexpr double CALC_DX = LX / CALC_NX; -constexpr double V_DOMAIN_LEFT_1 = -10.; -constexpr double V_DOMAIN_RIGHT_1 = 10.; +constexpr double V_DOMAIN_LEFT = -10.; +constexpr double V_DOMAIN_RIGHT = 10.; -constexpr unsigned int NV_1 = 128; -constexpr double DV_1 = std::abs(V_DOMAIN_RIGHT_1 - V_DOMAIN_LEFT_1) / NV_1; +constexpr unsigned int NV = 128; +constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV; -constexpr double V_DOMAIN_LEFT_2 = -10.; -constexpr double V_DOMAIN_RIGHT_2 = 10.; +// f0_TYPE: +// 0 -> twos-stream +// 1 -> landau-damping +// 2 -> maxwellian +// 3 -> bump-on-tail +constexpr size_t f0_TYPE = 0; -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 GLOBAL_REFINEMENT = 7; constexpr unsigned int FE_DEGREE = 3; constexpr unsigned int CONVERGENCE_ITERATIONS = 5000; constexpr double CONVERGENCE_LIMIT = 1e-8; +// Adaptive refinement options +constexpr unsigned int REFINE_FREQUENCY = 30; +constexpr double REFINEMENT_TOP_FRACTION = 0.8; +constexpr double REFINEMENT_BOTTOM_FRACTION = 0.1; + // Gauge options constexpr double GAUGE_DOMAIN_LEFT = 3.2; constexpr double GAUGE_DOMAIN_RIGHT = 3.8; @@ -44,13 +52,11 @@ constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) // NUFI options constexpr double DT = 1. / 10.; constexpr unsigned int TMAX = 100; -constexpr unsigned int REFINE_FREQUENCY = 30; // Plotting options 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/poisson_problem.h b/nufi/poisson_problem.h index ac1e8b2..f3e0a9c 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -1,9 +1,7 @@ #ifndef POISSON_PROBLEM_H #define POISSON_PROBLEM_H -#include #include - #include #include #include @@ -42,6 +40,7 @@ #include #include +#include #include #include #include @@ -49,6 +48,7 @@ #include #include #include +#include #include #include #include @@ -56,6 +56,7 @@ #include "nufi/cells.h" #include "nufi/grids.h" #include "nufi/parameters.h" +#include "nufi/stopwatch.h" #include "omp.h" void save_space_vector(const std::vector &vals, @@ -70,8 +71,8 @@ public: PoissonProblem(unsigned int degree); void initialize(); - void solve_step(size_t it, std::vector> &grid_versions, - bool refining = false); + double solve_step(size_t it, std::vector> &grid_versions, + bool refining = false); void coarse_and_refine_grid(size_t it); void setup_constraints(AffineConstraints &constraints); void run(); @@ -91,6 +92,7 @@ public: return constraints; } const CellLocator &get_locator() const { return cell_locator; } + double get_error_estimate() const { return error_estimate; } std::vector sample_electric_field(double x_min, double x_max, unsigned int Nx); @@ -108,6 +110,7 @@ private: void setup_system(); void assemble_system(); void solve(size_t it); + void estimate_error(); std::function(const std::vector> &)> rhs_function; @@ -124,6 +127,7 @@ private: Vector system_rhs; const bool PRINT_GAUGE_DOF_POSITION = true; + double error_estimate = 0.0; }; //====//====// @@ -336,10 +340,6 @@ template void PoissonProblem::setup_system() { // used for evaluator to avoid running it anytime there is an eval cell_locator.rebuild(dof_handler, triangulation); - - // local_solution_buffer.resize(fe.n_dofs_per_cell()); - // evaluator = std::make_unique>(mapping, fe, - // update_gradients); } template void PoissonProblem::assemble_system() { @@ -375,7 +375,9 @@ template void PoissonProblem::assemble_system() { Assert(rhs_function, ExcMessage("Poisson RHS function has not been initialized.")); + std::cout << "Start of full rho eval..." << "\n"; std::vector all_rho = rhs_function(all_q_points); + std::cout << "End of full rho eval..." << "\n"; Assert(all_rho.size() == all_q_points.size(), ExcMessage("rhs_function returned wrong size")); @@ -417,21 +419,20 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { dof_handler, QGauss(fe.degree + 1), std::map *>(), solution, error_per_cell); - GridRefinement::refine_and_coarsen_fixed_number(triangulation, error_per_cell, - 0.3, 0.03); - // START: remove refinment flags from edges of domain to alow safe gauge - // fixing - // for (const auto &cell : triangulation.active_cell_iterators()) { - // const double x = cell->center()[0]; - // if (x >= Parameters::X_DOMAIN_RIGHT - .5) { - // cell->clear_refine_flag(); - // cell->clear_coarsen_flag(); - // } - // } - // END + // GridRefinement::refine_and_coarsen_fixed_number(triangulation, + // error_per_cell, + // 0.3, 0.03); + GridRefinement::refine_and_coarsen_fixed_fraction( + triangulation, error_per_cell, Parameters::REFINEMENT_TOP_FRACTION, + Parameters::REFINEMENT_BOTTOM_FRACTION, + std::numeric_limits::max(), VectorTools::L2_norm); + + // Avoid coarsing below GLOBAL_REFINEMENT level for CellLocator + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->level() <= static_cast(Parameters::GLOBAL_REFINEMENT)) + cell->clear_coarsen_flag(); - // triangulation.prepare_coarsening_and_refinement(); triangulation.execute_coarsening_and_refinement(); std::cout << "Refinement Finished..." << "\n"; @@ -439,13 +440,17 @@ template void PoissonProblem::coarse_and_refine_grid(size_t it) { std::string grid_file_name = Parameters::PLOT_DIR + "grid_" + std::to_string(it); save_grid_to_file(grid_file_name); +} - // std::vector Ex = - // sample_electric_potential(Parameters::X_DOMAIN_LEFT, - // Parameters::X_DOMAIN_RIGHT, - // Parameters::PLOT_NX); - // - // save_space_vector(Ex, "Ex_after_coarsed", it); +template void PoissonProblem::estimate_error() { + Vector error_per_cell(triangulation.n_active_cells()); + + KellyErrorEstimator::estimate( + dof_handler, QGauss(fe.degree + 1), + std::map *>(), solution, + error_per_cell); + + error_estimate = error_per_cell.l2_norm(); } template void PoissonProblem::solve(size_t it) { @@ -488,15 +493,21 @@ template void PoissonProblem::initialize() { } template -void PoissonProblem::solve_step( +double PoissonProblem::solve_step( size_t it, std::vector> &grid_versions, bool refining) { + double refining_time = 0.0; if (refining) { + stopwatch refining_timer; coarse_and_refine_grid(it); setup_system(); update_grid_versions(grid_versions, *this); + refining_time = refining_timer.elapsed(); } assemble_system(); solve(it); + estimate_error(); + + return refining_time; } // NuFI doesnt use this, kept only for testing PoissonProblem diff --git a/nufi/save_results.h b/nufi/save_results.h index caad066..40bbb49 100644 --- a/nufi/save_results.h +++ b/nufi/save_results.h @@ -8,21 +8,29 @@ class NuFISolver; -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); +struct DiagnosticsSnapshot { + unsigned int Nx = 0; + unsigned int Nv = 0; + std::vector x_eval; + std::vector v_eval; + std::vector f; + std::vector rho; + std::vector E; +}; -void save_rho(const NuFISolver &solver, unsigned int n, - std::vector> &grid_struct, - std::vector> &phi_history, - unsigned int Nx_out, const std::string &filename); +DiagnosticsSnapshot +compute_diagnostics(const NuFISolver &solver, unsigned int n, + std::vector> &grid_struct, + std::vector> &phi_history, + unsigned int Nx_out, unsigned int Nv_out); -void save_Efield(unsigned int it, std::vector> &grid_versions, - std::vector> &phi_history, - unsigned int Nx_out = Parameters::PLOT_NX); +void save_f(const DiagnosticsSnapshot &snap, const std::string &filepath); +void save_rho(const DiagnosticsSnapshot &snap, const std::string &filepath); +void save_Efield(const DiagnosticsSnapshot &snap, const std::string &filepath); +double compute_int_E_squared(const DiagnosticsSnapshot &snap); -void save_space_vector(const std::vector &vals, - const std::string &filename, size_t it); +void save_time_series(const std::vector &t, + const std::vector &values, + const std::string &filepath); #endif diff --git a/src/main.cc b/src/main.cc index e5f3c67..b6f9e42 100644 --- a/src/main.cc +++ b/src/main.cc @@ -16,6 +16,7 @@ void clear_results_directory(const std::string &dir) { } int main() { + // omp_set_max_active_levels(1); 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 8fde451..299161c 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -26,18 +27,18 @@ using namespace dealii; std::vector NuFISolver::eval_ftilda( - unsigned int n, std::vector X, double u1, double u2, + unsigned int n, std::vector X, double u, const std::vector> &grid_struct, const std::vector> &phi_history) const { size_t x_size = X.size(); - std::vector U1(x_size, u1); - std::vector U2(x_size, u2); + std::vector U(x_size, u); std::vector results(x_size); if (n == 0) { for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U1[i], U2[i]); + results[i] = f0(X[i], U[i]); + // reset_x_eval(X); return results; } @@ -47,53 +48,61 @@ 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 * U1[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")); 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]; - 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 + U[i] = U[i] + Parameters::DT * Ex[i]; } } // The final half-step. for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U1[i]; + X[i] = X[i] - Parameters::DT * U[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]; - U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i]; + U[i] = U[i] + 0.5 * Parameters::DT * Ex[i]; } for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U1[i], U2[i]); + results[i] = f0(X[i], U[i]); + // reset_x_eval(X); return results; } std::vector -NuFISolver::eval_f(unsigned int n, std::vector X, double u1, double u2, +NuFISolver::eval_f(unsigned int n, std::vector X, double u, const std::vector> &grid_struct, const std::vector> &phi_history) const { size_t x_size = X.size(); - std::vector U1(x_size, u1); - std::vector U2(x_size, u2); - + std::vector U(x_size, u); std::vector results(x_size); if (n == 0) { for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U1[i], U2[i]); + results[i] = f0(X[i], U[i]); + // reset_x_eval(X); return results; } @@ -105,34 +114,35 @@ NuFISolver::eval_f(unsigned int n, std::vector X, double u1, double u2, phi_history[n].solution); // call eval only once for (size_t i = 0; i < x_size; ++i) { Ex[i] = -tmp[i]; - U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i]; + U[i] = U[i] + 0.5 * Parameters::DT * Ex[i]; } while (--n) { for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U1[i]; + X[i] = X[i] - Parameters::DT * U[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]; - U1[i] = U1[i] + Parameters::DT * Ex[i]; + U[i] = U[i] + Parameters::DT * Ex[i]; } } // The final half-step. for (size_t i = 0; i < x_size; ++i) - X[i] = X[i] - Parameters::DT * U1[i]; + X[i] = X[i] - Parameters::DT * U[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]; - U1[i] = U1[i] + 0.5 * Parameters::DT * Ex[i]; + U[i] = U[i] + 0.5 * Parameters::DT * Ex[i]; } for (size_t i = 0; i < x_size; ++i) - results[i] = f0(X[i], U1[i], U2[i]); + results[i] = f0(X[i], U[i]); + // reset_x_eval(X); return results; } @@ -140,30 +150,31 @@ std::vector NuFISolver::eval_rho(unsigned int n, std::vector &X, const std::vector> &grid_struct, const std::vector> &phi_history, - const unsigned int Nv1, const unsigned int Nv2) const { + const unsigned int Nv) const { size_t x_size = X.size(); - const double dv_0 = Parameters::DV_1; - const double v_min_0 = Parameters::V_DOMAIN_LEFT_1 + 0.5 * dv_0; + const double dv = + (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv; + 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 partial(static_cast(Nv) * x_size); - std::vector integral(x_size, 0.0); - std::vector tmp_int(x_size); - - 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]; - } +#pragma omp parallel for + for (unsigned int i = 0; i < Nv; ++i) { + std::vector tmp_int = + eval_ftilda(n, X, v_min + i * dv, grid_struct, + phi_history); // used eval_ftilda once per i + std::copy(tmp_int.begin(), tmp_int.end(), + partial.begin() + static_cast(i) * x_size); } + std::vector integral(x_size, 0.0); + for (unsigned int i = 0; i < Nv; ++i) + for (size_t ii = 0; ii < x_size; ++ii) + integral[ii] += partial[static_cast(i) * x_size + ii]; + for (size_t i = 0; i < x_size; ++i) - integral[i] = 1 - integral[i] * dv_0; + integral[i] = 1 - integral[i] * dv; return integral; } @@ -171,11 +182,9 @@ 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 Nv1, - const unsigned int Nv2) const { + const unsigned int Nv) const { std::vector point_vector = Point_vector_to_double_vector(points); - return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv1, - Nv2); + return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv); } void NuFISolver::run() { @@ -187,29 +196,20 @@ void NuFISolver::run() { using std::abs; using std::max; - // std::unique_ptr rho{ - // reinterpret_cast(std::aligned_alloc(64, sizeof(double) * - // Nx)), - // std::free}; - // - // if (rho == nullptr) - // throw std::bad_alloc{}; - std::vector int_E_squared; int_E_squared.reserve(Nt); + std::vector int_E_squared_times; + int_E_squared_times.reserve(Nt); std::vector> grid_versions; std::vector> phi_history; - // update_grid_versions(grid_versions, poisson); - // update_solution_history(phi_history, poisson, - // grid_versions.back().grid_version); - std::vector x_eval(Parameters::CALC_NX); - std::ofstream time_file("results/simulation_time.dat"); + std::ofstream time_file(Parameters::PLOT_DIR + "simulation_time.dat"); double total_time = 0; + stopwatch total_timer; time_file << "it " << "step_time " @@ -219,6 +219,10 @@ void NuFISolver::run() { << "plot_time" << "\n"; + std::ofstream error_file(Parameters::PLOT_DIR + "error_estimate.dat"); + error_file << "# nufi Kelly l2 error estimate\n"; + error_file << "# it l2_error_estimate\n"; + [[maybe_unused]] const double x_min = Parameters::X_DOMAIN_LEFT; [[maybe_unused]] double dx = Parameters::CALC_DX; @@ -236,11 +240,12 @@ void NuFISolver::run() { std::cout << "Timestep " << it << " / " << Nt << " (simulation time = " << it * Parameters::DT << ")" - << std::endl; + << "\n"; // START: diagnostics - // std::cout << "cells = " << poisson.triangulation.n_active_cells() << "\n" - // << " dofs = " << poisson.dof_handler.n_dofs() << "\n"; + std::cout << "cells = " << poisson.get_triangulation().n_active_cells() + << "\n" + << " dofs = " << poisson.get_dof_handler().n_dofs() << "\n"; // double min_h = 1e100; // double max_h = 0; // @@ -276,14 +281,13 @@ 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_1, - Parameters::NV_2); + return eval_rho(it, x, grid_versions, phi_history, Parameters::NV); }); if (it % Parameters::REFINE_FREQUENCY == 0) { // if (it == 0) { - poisson.solve_step(it, grid_versions, true); - compute_time = timer.elapsed() - compute_start; + refine_time = poisson.solve_step(it, grid_versions, true); + compute_time = timer.elapsed() - compute_start - refine_time; } else { poisson.solve_step(it, grid_versions, false); compute_time = timer.elapsed() - compute_start; @@ -291,6 +295,9 @@ void NuFISolver::run() { update_solution_history(phi_history, poisson, grid_versions.back().grid_version); + error_file << it << " " << poisson.get_error_estimate() << "\n"; + error_file.flush(); + double timer_elapsed = timer.elapsed(); double step_time = timer_elapsed - time_elapsed_before; @@ -299,28 +306,31 @@ void NuFISolver::run() { //====//====// // Plotting // //====//====// + if (it % Parameters::PLOT_FREQUENCY == 0) { double plot_start = timer.elapsed(); std::cout << "Saving results... "; - 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"); + DiagnosticsSnapshot snap = + compute_diagnostics(*this, it, grid_versions, phi_history, + Parameters::PLOT_NX, Parameters::NV); - save_Efield(it, grid_versions, phi_history); + save_f(snap, Parameters::PLOT_DIR + "f_" + std::to_string(it) + ".dat"); + save_rho(snap, + Parameters::PLOT_DIR + "rho_" + std::to_string(it) + ".dat"); + save_Efield(snap, + Parameters::PLOT_DIR + "E_" + std::to_string(it) + ".dat"); - double int_val = 0.5 * integral_space_vector_squared( - grid_versions[phi_history[it].grid_version], - phi_history[it].solution); - int_E_squared.push_back(int_val); - save_space_vector(int_E_squared, "electricint", it); + int_E_squared.push_back(compute_int_E_squared(snap)); + int_E_squared_times.push_back(it * Parameters::DT); + save_time_series(int_E_squared_times, int_E_squared, + Parameters::PLOT_DIR + "int_E_sqr.dat"); plot_time = timer.elapsed() - plot_start; std::cout << "Results saved in " << plot_start << "[s]" << "\n"; } - total_time = timer.elapsed(); + + total_time = total_timer.elapsed(); std::cout << "Time since start = " << total_time << "\n\n"; time_file << it << " " << step_time << " " << total_time << " " diff --git a/src/save_results.cc b/src/save_results.cc index bcc7559..2bc3850 100644 --- a/src/save_results.cc +++ b/src/save_results.cc @@ -4,113 +4,128 @@ #include "nufi/grids.h" #include "nufi/nufi_solver.h" #include "nufi/parameters.h" -#include "nufi/poisson_problem.h" #include #include +#include #include #include #include #include -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) { - std::ofstream file(filename); +namespace { +void ensure_results_dir() { std::filesystem::create_directories("results"); } +} // namespace - double xmin = Parameters::X_DOMAIN_LEFT; - double xmax = Parameters::X_DOMAIN_RIGHT; - - double vmin_1 = Parameters::V_DOMAIN_LEFT_1; - double vmax_1 = Parameters::V_DOMAIN_RIGHT_1; - - double dv_1 = (vmax_1 - vmin_1) / Nv_out; - - file << Nx_out << " " << Nv_out << "\n"; - file << xmin << " " << xmax << "\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_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]; - - if (i < Nx_out - 1) - file << " "; - } - - file << "\n"; - } - - file.close(); -} - -void save_rho(const NuFISolver &solver, unsigned int n, - std::vector> &grid_struct, - std::vector> &phi_history, - unsigned int Nx_out, const std::string &filename) { - std::ofstream file(filename); - - double xmin = Parameters::X_DOMAIN_LEFT; - double xmax = Parameters::X_DOMAIN_RIGHT; - - std::vector x_eval = make_x_eval(Nx_out); - file << Nx_out << "\n"; - file << xmin << " " << xmax << "\n"; - - std::vector tmp = - solver.eval_rho(n, x_eval, grid_struct, phi_history); - for (size_t i = 0; i < Nx_out; ++i) { - file << tmp[i]; - file << "\n"; - } - file.close(); -} - -void save_Efield(unsigned int it, std::vector> &grid_versions, - std::vector> &phi_history, - unsigned int Nx_out) { - std::vector x_eval_E = make_x_eval(Nx_out); - - auto grad_phi = eval(x_eval_E, grid_versions[phi_history[it].grid_version], - phi_history[it].solution); - - std::vector> ordered_E; - - for (size_t i = 0; i < x_eval_E.size(); ++i) { - double E = -grad_phi[i]; - ordered_E.push_back({x_eval_E[i], E}); - } - - std::sort(ordered_E.begin(), ordered_E.end()); - - std::ofstream E_file("results/E_" + std::to_string(it) + ".dat"); - - for (auto [x, E] : ordered_E) - E_file << x << " " << E << "\n"; - - std::cout << "Saving iteration " << it << " using grid version " - << phi_history[it].grid_version << "\n\n"; - - E_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_field_1d(const std::vector &x, + const std::vector &values, + const std::string &filepath) { + if (x.size() != values.size()) + throw std::runtime_error( + "save_field_1d: size mismatch (x=" + std::to_string(x.size()) + + ", values=" + std::to_string(values.size()) + ") writing " + filepath); + ensure_results_dir(); + std::ofstream file(filepath); if (!file) - throw std::runtime_error("failed to start file in results/"); + throw std::runtime_error("save_field_1d: failed to open " + filepath); - 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 << "# nufi 1d field\n"; + file << "# n = " << x.size() << "\n"; + file << "# columns: x value\n"; + file << std::setprecision(10); + for (size_t i = 0; i < x.size(); ++i) + file << x[i] << " " << values[i] << "\n"; +} + +DiagnosticsSnapshot +compute_diagnostics(const NuFISolver &solver, unsigned int n, + std::vector> &grid_struct, + std::vector> &phi_history, + unsigned int Nx_out, unsigned int Nv_out) { + DiagnosticsSnapshot snap; + snap.Nx = Nx_out; + snap.Nv = Nv_out; + + const double vmin = Parameters::V_DOMAIN_LEFT; + const double vmax = Parameters::V_DOMAIN_RIGHT; + const double dv = (vmax - vmin) / Nv_out; + + snap.x_eval = make_x_eval(Nx_out); + snap.v_eval.resize(Nv_out); + for (unsigned int j = 0; j < Nv_out; ++j) + snap.v_eval[j] = vmin + (j + 0.5) * dv; + + snap.f.resize(static_cast(Nx_out) * Nv_out); +#pragma omp parallel for + for (unsigned int j = 0; j < Nv_out; ++j) { + std::vector val = + solver.eval_f(n, snap.x_eval, snap.v_eval[j], grid_struct, phi_history); + std::copy(val.begin(), val.end(), snap.f.begin() + j * Nx_out); + } + + snap.rho.assign(Nx_out, 1.0); + for (unsigned int j = 0; j < Nv_out; ++j) + for (unsigned int i = 0; i < Nx_out; ++i) + snap.rho[i] -= snap.f[j * Nx_out + i] * dv; + + std::vector x_copy = snap.x_eval; + auto grad_phi = eval(x_copy, grid_struct[phi_history[n].grid_version], + phi_history[n].solution); + snap.E.resize(Nx_out); + for (unsigned int i = 0; i < Nx_out; ++i) + snap.E[i] = -grad_phi[i]; + + return snap; +} + +void save_f(const DiagnosticsSnapshot &snap, const std::string &filepath) { + ensure_results_dir(); + std::ofstream file(filepath); + if (!file) + throw std::runtime_error("save_f: failed to open " + filepath); + + file << "# nufi f(x,v) ascii dump\n"; + file << "# Nx = " << snap.Nx << " Nv = " << snap.Nv << "\n"; + file << "# columns: x v f\n"; + file << std::setprecision(10); + for (unsigned int j = 0; j < snap.Nv; ++j) + for (unsigned int i = 0; i < snap.Nx; ++i) + file << snap.x_eval[i] << " " << snap.v_eval[j] << " " + << snap.f[j * snap.Nx + i] << "\n"; +} + +void save_rho(const DiagnosticsSnapshot &snap, const std::string &filepath) { + save_field_1d(snap.x_eval, snap.rho, filepath); +} + +void save_Efield(const DiagnosticsSnapshot &snap, const std::string &filepath) { + save_field_1d(snap.x_eval, snap.E, filepath); +} + +double compute_int_E_squared(const DiagnosticsSnapshot &snap) { + const double dx = Parameters::LX / snap.Nx; + double integral = 0.0; + for (double e : snap.E) + integral += e * e; + return 0.5 * integral * dx; +} + +void save_time_series(const std::vector &t, + const std::vector &values, + const std::string &filepath) { + if (t.size() != values.size()) + throw std::runtime_error( + "save_time_series: size mismatch (t=" + std::to_string(t.size()) + + ", values=" + std::to_string(values.size()) + ") writing " + filepath); + + ensure_results_dir(); + std::ofstream file(filepath); + if (!file) + throw std::runtime_error("save_time_series: failed to open " + filepath); + + file << "# nufi time series\n"; + file << "# columns: t value\n"; + file << std::setprecision(10); + for (size_t i = 0; i < t.size(); ++i) + file << t[i] << " " << values[i] << "\n"; }