diff --git a/include/samurai/algorithm/graduation.hpp b/include/samurai/algorithm/graduation.hpp index d749f3201..47a74e50b 100644 --- a/include/samurai/algorithm/graduation.hpp +++ b/include/samurai/algorithm/graduation.hpp @@ -172,12 +172,13 @@ namespace samurai } template - void list_intervals_to_remove(const size_t grad_width, - const CellArray& ca, - [[maybe_unused]] const std::vector>& mpi_neighbourhood, - const std::array& is_periodic, - const std::array& nb_cells_finest_level, - std::array, CellArray::max_size>& out) + void list_interval_to_refine_for_graduation( + const size_t grad_width, + const CellArray& ca, + [[maybe_unused]] const std::vector>& mpi_neighbourhood, + const std::array& is_periodic, + const std::array& nb_cells_finest_level, + std::array, CellArray::max_size>& out) { const size_t max_level = ca.max_level(); const size_t min_level = ca.min_level(); @@ -218,7 +219,7 @@ namespace samurai } } } - xt::xtensor_fixed> translation = xt::xscalar(0); + DirectionVector translation = xt::xscalar(0); for (size_t d = 0; d != dim; ++d) { if (is_periodic[d]) @@ -260,6 +261,87 @@ namespace samurai #endif // SAMURAI_WITH_MPI } + template + void list_interval_to_refine_for_contiguous_boundary_cells( + const size_t half_stencil_width, + const CellArray& ca, + const LevelCellArray& domain, + const std::array& is_periodic, + std::array, CellArray::max_size>& out) + { + const size_t max_level = ca.max_level(); + const size_t min_level = ca.min_level(); + + // We want to avoid that a flux is computed with ghosts outside of the domain if the cell doesn't touch the boundary, + // because we only want to apply the B.C. on the cells that touch the boundary. + // 1. Case where the boundary is at level L and the jump is going down to L-1: + // We want to have enough contiguous boundary cells to ensure that the stencil at the lower level + // won't go outside the domain. + // To ensure half_stencil_width at L-1, we need 2*half_stencil_width at level L. + // However, since we project the B.C. in the first outside ghost at level L-1, we can reduce the number of contiguous + // cells by 1 at level L-1. This makes, at level L, 2*(half_stencil_width - 1) contiguous cells. + const int n_contiguous_boundary_cells = std::max(int(half_stencil_width), 2 * (int(half_stencil_width) - 2)); + + // 2. Case where the boundary is at level L and jump is going up: + // Then ensuring half_stencil_width contiguous cells at level L automatically ensures half_stencil_width + // at level L+1. + + for_each_cartesian_direction( + [&](const auto direction_idx, const auto& translation) + { + if (not is_periodic[direction_idx]) + { + // Jump level --> level-1 + for (size_t level = max_level; level != min_level; --level) + { + auto boundaryCells = difference(ca[level], translate(self(domain).on(level), -translation)).on(level); + + for (int i = 2; i <= n_contiguous_boundary_cells; i += 2) + { + // auto refine_subset = intersection(translate(boundaryCells, -i*translation), ca[level-1]).on(level-1); + LevelCellArray translated_boundary(translate(boundaryCells, -i * translation)); + auto refine_subset = intersection(translated_boundary, ca[level - 1]).on(level - 1); + refine_subset( + [&](const auto& x_interval, const auto& yz) + { + out[level - 1].push_back(x_interval, yz); + }); + } + } + // Jump level --> level+1 + for (size_t level = max_level - 1; level != min_level - 1; --level) + { + auto boundaryCells = difference(ca[level], translate(self(domain).on(level), -translation)); + for (size_t i = 1; i != half_stencil_width; ++i) + { + auto refine_subset = translate(intersection(translate(boundaryCells, -i * translation), ca[level + 1]).on(level), + i * translation) + .on(level); + refine_subset( + [&](const auto& x_interval, const auto& yz) + { + out[level].push_back(x_interval, yz); + }); + } + } + } + }); + } + + template + void list_intervals_to_refine(const size_t grad_width, + const size_t half_stencil_width, + const CellArray& ca, + const LevelCellArray& domain, + [[maybe_unused]] const std::vector>& mpi_neighbourhood, + const std::array& is_periodic, + const std::array& nb_cells_finest_level, + std::array, CellArray::max_size>& out) + { + list_interval_to_refine_for_graduation(grad_width, ca, mpi_neighbourhood, is_periodic, nb_cells_finest_level, out); + list_interval_to_refine_for_contiguous_boundary_cells(half_stencil_width, ca, domain, is_periodic, out); + } + // if add the intervals in add_m_interval // if dim = 2 then add_m_interval stores the y coord // if dim > 2 then add_intercal contains the 'inner_stencil' i.e. the coordinates y+s_x, z+s_z, etc. @@ -317,10 +399,12 @@ namespace samurai template size_t make_graduation(CellArray& ca, + const LevelCellArray& domain, [[maybe_unused]] const std::vector>& mpi_neighbourhood, const std::array& is_periodic, - const std::array& nb_cells_finest_level, - const size_t grad_width = 1) + const size_t grad_width = 1, + const size_t half_stencil_width = 1 // helf of width of the numerical scheme's stencil. + ) { using ca_type = CellArray; using coord_type = typename ca_type::lca_type::coord_type; @@ -328,6 +412,16 @@ namespace samurai const size_t max_level = ca.max_level(); const size_t min_level = ca.min_level(); + const auto& min_indices = domain.min_indices(); + const auto& max_indices = domain.max_indices(); + + std::array nb_cells_finest_level; + + for (size_t d = 0; d != max_indices.size(); ++d) + { + nb_cells_finest_level[d] = max_indices[d] - min_indices[d]; + } + std::vector add_p_interval; std::vector add_p_inner_stencil; std::vector add_p_idx; @@ -351,18 +445,18 @@ namespace samurai // Then, if the non-graduated is not taged as keep, we coarsen it ca_add_p.clear(); ca_remove_p.clear(); - list_intervals_to_remove(grad_width, ca, mpi_neighbourhood, is_periodic, nb_cells_finest_level, remove_m_all); + list_intervals_to_refine(grad_width, half_stencil_width, ca, domain, mpi_neighbourhood, is_periodic, nb_cells_finest_level, remove_m_all); add_p_interval.clear(); add_p_inner_stencil.clear(); add_p_idx.clear(); for (size_t level = min_level; level != max_level + 1; ++level) { -#ifdef SAMURAI_WITH_MPI + // #ifdef SAMURAI_WITH_MPI remove_m_all[level].remove_overlapping_intervals(); -#else - remove_m_all[level].sort_intervals(); -#endif // SAMURAI_WITH_MPI + // #else + // remove_m_all[level].sort_intervals(); + // #endif // SAMURAI_WITH_MPI const size_t imax = remove_m_all[level].size(); for (size_t i = 0; i != imax; ++i) { @@ -424,10 +518,10 @@ namespace samurai std::vector> mpi_neighbourhood; std::array is_periodic; - std::array nb_cells_finest_level; + LevelCellArray domain; is_periodic.fill(false); - return make_graduation(ca, mpi_neighbourhood, is_periodic, nb_cells_finest_level, grad_width); + return make_graduation(ca, domain, mpi_neighbourhood, is_periodic, grad_width); } template diff --git a/include/samurai/mr/adapt.hpp b/include/samurai/mr/adapt.hpp index 215a661b7..aa88660cc 100644 --- a/include/samurai/mr/adapt.hpp +++ b/include/samurai/mr/adapt.hpp @@ -311,18 +311,15 @@ namespace samurai // on test adapt_test/2.mutliple_fields with: // linux-mamba (clang-18, ubuntu-24.04, clang, clang-18, clang-18, clang++-18) // while the code bellow do not. - const auto& min_indices = mesh.domain().min_indices(); - const auto& max_indices = mesh.domain().max_indices(); - - std::array nb_cells_finest_level; - - for (size_t d = 0; d != max_indices.size(); ++d) - { - nb_cells_finest_level[d] = max_indices[d] - min_indices[d]; - } ca_type new_ca = update_cell_array_from_tag(mesh[mesh_id_t::cells], m_tag); - make_graduation(new_ca, mesh.mpi_neighbourhood(), mesh.periodicity(), nb_cells_finest_level, mesh_t::config::graduation_width); + make_graduation(new_ca, + mesh.domain(), + mesh.mpi_neighbourhood(), + mesh.periodicity(), + mesh_t::config::graduation_width, + mesh_t::config::max_stencil_width); + mesh_t new_mesh{new_ca, mesh}; #ifdef SAMURAI_WITH_MPI mpi::communicator world; diff --git a/include/samurai/stencil.hpp b/include/samurai/stencil.hpp index d7a58331a..0abaa3a9f 100644 --- a/include/samurai/stencil.hpp +++ b/include/samurai/stencil.hpp @@ -130,6 +130,21 @@ namespace samurai return 0; } + template + void for_each_cartesian_direction(Func&& f) + { + DirectionVector direction; + direction.fill(0); + for (std::size_t d = 0; d < dim; ++d) + { + direction[d] = 1; + f(d, direction); + direction[d] = -1; + f(d, direction); + direction[d] = 0; + } + } + //-----------------------// // Useful stencils // //-----------------------//