From c91d39d9e11cc2188abf85f1e5e1010119db12dd Mon Sep 17 00:00:00 2001 From: "Vasco C. B. Ferreira" Date: Thu, 30 Jul 2026 23:17:20 +0200 Subject: [PATCH] seems to be faster, to test --- README.md | 14 +++- nufi/cells.h | 170 ++++++++++++++++++++++------------------- nufi/grids.h | 3 +- nufi/poisson_problem.h | 3 +- 4 files changed, 105 insertions(+), 85 deletions(-) diff --git a/README.md b/README.md index ef9832b..b5151a7 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,25 @@ -# 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: 1x1v -status: Working, to be re-re-viewed +notes: + +- about 3 times slower than previous locator + +status: Working, to be re-reviewed Refinement working: + - grid versions saved on a vector - solutions point to a version of the grid todo: + - add ions diff --git a/nufi/cells.h b/nufi/cells.h index 8bcb9b7..91d1312 100644 --- a/nufi/cells.h +++ b/nufi/cells.h @@ -2,8 +2,10 @@ #define CELLS_H #include +#include #include #include +#include #include #include @@ -13,19 +15,18 @@ using namespace dealii; -// Current Locator needs Grid to be a hypercube !! - +// Current Locator needs Grid to be a single hypercube coarse cell, +// uniformly refined to `base_level`, with isotropic local refinement on top. template struct CellLocation { - using CellIterator = typename DoFHandler::cell_iterator; - - CellIterator cell; + typename DoFHandler::active_cell_iterator cell; Point reference_point; }; template class CellLocator { public: void rebuild(const DoFHandler &dof_handler, - const Triangulation &triangulation); + const Triangulation &triangulation, + unsigned int base_level); CellLocation locate(const Point &point) const; @@ -34,125 +35,134 @@ private: Point domain_lower_; Point domain_upper_; + std::array base_n_{}; + std::array base_dx_{}; + + // Flat, O(1)-indexable array of the base_level cells. + std::vector::cell_iterator> base_cells_; }; template void CellLocator::rebuild(const DoFHandler &dof_handler, - const Triangulation &triangulation) { + const Triangulation &triangulation, + unsigned int base_level) { AssertThrow(&dof_handler.get_triangulation() == &triangulation, ExcMessage("CellLocator::rebuild(): DoFHandler and Triangulation " "do not refer to the same mesh.")); - auto root = dof_handler.begin(0); - - AssertThrow( - root != dof_handler.end(0), - ExcMessage("CellLocator::rebuild(): triangulation has no level-zero " - "cell.")); - - auto after_root = root; - ++after_root; - - AssertThrow( - after_root == dof_handler.end(0), - ExcMessage("CellLocator currently requires exactly one coarse cell.")); - dof_handler_ = &dof_handler; + // --- bounding box of the (single) coarse cell --- + auto root = dof_handler.begin(0); + AssertThrow(root != dof_handler.end(0), + ExcMessage("CellLocator::rebuild(): no level-zero cell.")); + { + auto after_root = root; + ++after_root; + AssertThrow(after_root == dof_handler.end(0), + ExcMessage("CellLocator currently requires exactly one " + "coarse cell.")); + } + for (unsigned int d = 0; d < dim; ++d) { domain_lower_[d] = std::numeric_limits::max(); domain_upper_[d] = std::numeric_limits::lowest(); } - for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) for (unsigned int d = 0; d < dim; ++d) { domain_lower_[d] = std::min(domain_lower_[d], root->vertex(v)[d]); - domain_upper_[d] = std::max(domain_upper_[d], root->vertex(v)[d]); } + // --- build the O(1)-indexable array of base_level cells --- + unsigned int total = 1; for (unsigned int d = 0; d < dim; ++d) { const double length = domain_upper_[d] - domain_lower_[d]; - AssertThrow(length > 0.0, - ExcMessage("CellLocator::rebuild(): coarse-cell domain has " - "non-positive extent.")); - - const double scale = - std::max({1.0, std::abs(domain_lower_[d]), std::abs(domain_upper_[d])}); - - const double tolerance = - 100.0 * std::numeric_limits::epsilon() * scale; - - for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) { - const double coordinate = root->vertex(v)[d]; - - const bool lies_on_lower = - std::abs(coordinate - domain_lower_[d]) <= tolerance; - - const bool lies_on_upper = - std::abs(coordinate - domain_upper_[d]) <= tolerance; - - AssertThrow( - lies_on_lower || lies_on_upper, - ExcMessage( - "CellLocator requires an axis-aligned hypercube coarse cell.")); - } + ExcMessage("CellLocator::rebuild(): non-positive domain " + "extent.")); + base_n_[d] = + 1u << base_level; // refine_global(base_level) -> 2^level per axis + base_dx_[d] = length / base_n_[d]; + total *= base_n_[d]; } + + base_cells_.assign(total, typename DoFHandler::cell_iterator()); + + for (auto cell = dof_handler.begin(base_level); + cell != dof_handler.end(base_level); ++cell) { + const auto bb = cell->bounding_box(); + const auto c_lower = bb.get_boundary_points().first; + + unsigned int idx = 0, stride = 1; + for (unsigned int d = 0; d < dim; ++d) { + unsigned int i = static_cast( + std::round((c_lower[d] - domain_lower_[d]) / base_dx_[d])); + i = std::min(i, base_n_[d] - 1); + idx += i * stride; + stride *= base_n_[d]; + } + base_cells_[idx] = cell; + } + + for (const auto &c : base_cells_) + AssertThrow(c.state() == IteratorState::valid, + ExcMessage("CellLocator::rebuild(): failed to fill the base " + "grid — mesh isn't uniformly refined to " + "'base_level' as expected.")); } template CellLocation CellLocator::locate(const Point &point) const { - AssertThrow( - dof_handler_ != nullptr, - ExcMessage("CellLocator::locate(): rebuild() has not been called.")); + AssertThrow(dof_handler_ != nullptr, + ExcMessage("CellLocator::locate(): rebuild() has not been " + "called.")); - Point reference_point; + std::array base_index; + Point reference_point; // local coords, updated at each descent step + // 1) periodic wrap + O(1) base-cell index per axis for (unsigned int d = 0; d < dim; ++d) { const double length = domain_upper_[d] - domain_lower_[d]; + double shifted = point[d] - domain_lower_[d]; + shifted -= length * std::floor(shifted / length); + if (shifted >= length) + shifted = 0.0; - const double shifted = point[d] - domain_lower_[d]; - - double periodic_offset = shifted - length * std::floor(shifted / length); - - if (periodic_offset < 0.0) - periodic_offset += length; - - if (periodic_offset >= length) - periodic_offset = 0.0; - - reference_point[d] = periodic_offset / length; + const double xi = shifted / base_dx_[d]; + const unsigned int i = + std::min(base_n_[d] - 1, static_cast(std::floor(xi))); + base_index[d] = i; + reference_point[d] = xi - i; // fraction within the base cell } - auto cell = dof_handler_->begin(0); + // 2) O(1) lookup of the base-level cell + unsigned int idx = 0, stride = 1; + for (unsigned int d = 0; d < dim; ++d) { + idx += base_index[d] * stride; + stride *= base_n_[d]; + } + auto cell = base_cells_[idx]; + // 3) O(R_level_max) descent — only costs anything for cells that are + // actually locally refined beyond base_level while (cell->has_children()) { - AssertThrow( - cell->n_children() == GeometryInfo::max_children_per_cell, - ExcMessage( - "CellLocator currently supports isotropic refinement only.")); - + AssertThrow(cell->n_children() == GeometryInfo::max_children_per_cell, + ExcMessage("CellLocator currently supports isotropic " + "refinement only.")); const unsigned int child_index = GeometryInfo::child_cell_from_point(reference_point); - reference_point = GeometryInfo::cell_to_child_coordinates( reference_point, child_index); - cell = cell->child(child_index); } - AssertThrow( - cell->is_active(), - ExcMessage("CellLocator tree traversal did not finish on an active " - "cell.")); + AssertThrow(cell->is_active(), + ExcMessage("CellLocator tree traversal did not finish on an " + "active cell.")); - AssertThrow( - GeometryInfo::is_inside_unit_cell(reference_point, 1e-12), - ExcMessage( - "CellLocator produced a reference point outside the unit cell.")); - - return CellLocation{cell, reference_point}; + return CellLocation{typename DoFHandler::active_cell_iterator(cell), + reference_point}; } #endif // CELLS_H diff --git a/nufi/grids.h b/nufi/grids.h index 65ce753..5540eac 100644 --- a/nufi/grids.h +++ b/nufi/grids.h @@ -140,7 +140,8 @@ GridStructure make_grid_snapshot(PoissonProblem &poisson) { grid.constraints->close(); grid.locator = std::make_unique>(); - grid.locator->rebuild(*grid.dof_handler, *grid.triangulation); + grid.locator->rebuild(*grid.dof_handler, *grid.triangulation, + Parameters::GLOBAL_REFINEMENT); // START: diagnostics AssertThrow( diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h index 0d15e69..332c309 100644 --- a/nufi/poisson_problem.h +++ b/nufi/poisson_problem.h @@ -338,7 +338,8 @@ template void PoissonProblem::setup_system() { system_rhs.reinit(dof_handler.n_dofs()); // used for evaluator to avoid running it anytime there is an eval - cell_locator.rebuild(dof_handler, triangulation); + cell_locator.rebuild(dof_handler, triangulation, + Parameters::GLOBAL_REFINEMENT); } template void PoissonProblem::assemble_system() {