From ba3b1cf1c08d090fab21330ab57896b98dbb5ce4 Mon Sep 17 00:00:00 2001 From: Loic Gouarin Date: Mon, 23 Feb 2026 11:23:44 +0100 Subject: [PATCH 1/4] init --- include/samurai/mr/adapt.hpp | 39 +++++++++++++++++++++--------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/include/samurai/mr/adapt.hpp b/include/samurai/mr/adapt.hpp index 475a6ca89..4400e98c5 100644 --- a/include/samurai/mr/adapt.hpp +++ b/include/samurai/mr/adapt.hpp @@ -125,6 +125,7 @@ namespace samurai PredictionFn m_prediction_fn; fields_t m_fields; // NOLINT(cppcoreguidelines-avoid-const-or-ref-data-members) detail_t m_detail; + typename mesh_t::ca_type m_interest_cells; tag_t m_tag; }; @@ -152,6 +153,7 @@ namespace samurai } cfg.parse_args(); + m_interest_cells = mesh[mesh_id_t::cells]; for (std::size_t i = 0; i < max_level - min_level; ++i) { // std::cout << "MR mesh adaptation " << i << std::endl; @@ -275,23 +277,22 @@ namespace samurai // We compute the detail in the cells and ghosts below the cells, except near the (non-periodic) boundaries, where we compute // the detail only in the cells (justification in the comments below). - bool periodic_in_all_directions = true; - std::array contract_directions; - for (std::size_t d = 0; d < dim; ++d) - { - periodic_in_all_directions = periodic_in_all_directions && mesh.is_periodic(d); - contract_directions[d] = !mesh.is_periodic(d); - } + // bool periodic_in_all_directions = true; + // std::array contract_directions; + // for (std::size_t d = 0; d < dim; ++d) + // { + // periodic_in_all_directions = periodic_in_all_directions && mesh.is_periodic(d); + // contract_directions[d] = !mesh.is_periodic(d); + // } times::timers.start("detail computation"); - for (std::size_t level = ((min_level > 0) ? min_level - 1 : 0); level < max_level - ite; ++level) + for (std::size_t level = ((min_level > 0) ? min_level - 1 : 0); level < max_level; ++level) { // 1. detail computation in the cells (at level+1) - auto ghosts_below_cells = intersection(mesh[mesh_id_t::all_cells][level], mesh[mesh_id_t::cells][level + 1]).on(level); + auto ghosts_below_cells = intersection(mesh[mesh_id_t::all_cells][level], m_interest_cells[level + 1]).on(level); ghosts_below_cells.apply_op(compute_detail(m_detail, m_fields)); // 'compute_detail' applies 1 level above the set // it is applied to, i.e. level+1 - // 2. detail computation in the ghosts below cells (at level) if (level >= min_level) { if (periodic_in_all_directions) @@ -310,7 +311,7 @@ namespace samurai // contract the domain only in non-periodic directions auto domain_without_bdry = contract(self(mesh.domain()).on(level - 1), 1, contract_directions); auto cells_without_bdry = intersection(intersection(mesh[mesh_id_t::all_cells][level - 1], domain_without_bdry), - mesh[mesh_id_t::cells][level + 1]) + m_interest_cells[level + 1]) .on(level - 1); cells_without_bdry.apply_op(compute_detail(m_detail, m_fields)); } @@ -326,14 +327,14 @@ namespace samurai update_ghost_subdomains(m_detail); times::timers.start("tag computation"); - for (std::size_t level = min_level; level <= max_level - ite; ++level) + for (std::size_t level = min_level; level <= max_level; ++level) { std::size_t exponent = dim * (max_level - level); double eps_l = cfg.epsilon() / (1 << exponent); double regularity_to_use = cfg.regularity() + dim; - auto subset_1 = intersection(mesh[mesh_id_t::cells][level], mesh[mesh_id_t::all_cells][level - 1]).on(level - 1); + auto subset_1 = intersection(m_interest_cells[level], mesh[mesh_id_t::all_cells][level - 1]).on(level - 1); subset_1.apply_op(to_coarsen_mr(m_detail, m_tag, eps_l, min_level), to_refine_mr(m_detail, @@ -348,9 +349,9 @@ namespace samurai keep_boundary_refined(mesh, m_tag); } - for (std::size_t level = min_level; level <= max_level - ite; ++level) + for (std::size_t level = min_level; level <= max_level; ++level) { - auto subset_2 = intersection(mesh[mesh_id_t::cells][level], mesh[mesh_id_t::cells][level]); + auto subset_2 = intersection(m_interest_cells[level], mesh[mesh_id_t::cells][level]); subset_2.apply_op(keep_around_refine(m_tag)); @@ -367,7 +368,7 @@ namespace samurai for (std::size_t level = max_level; level > 0; --level) { - auto keep_subset = intersection(mesh[mesh_id_t::cells][level], mesh[mesh_id_t::all_cells][level - 1]).on(level - 1); + auto keep_subset = intersection(m_interest_cells[level], mesh[mesh_id_t::all_cells][level - 1]).on(level - 1); update_tag_periodic(level, m_tag); update_tag_subdomains(level, m_tag); @@ -392,6 +393,12 @@ namespace samurai times::timers.stop("mesh update"); return true; } + for (std::size_t level = min_level; level <= max_level; ++level) + { + auto subset = difference(new_mesh[mesh_id_t::cells][level], mesh[mesh_id_t::cells][level]); + m_interest_cells[level] = subset.to_lca(); + } + times::timers.stop("mesh update"); update_ghost_mr(other_fields...); From f1f8c9d791508f2540542faae063eba452e37340 Mon Sep 17 00:00:00 2001 From: Loic Gouarin Date: Tue, 24 Feb 2026 06:51:29 +0100 Subject: [PATCH 2/4] ... --- include/samurai/mr/adapt.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/include/samurai/mr/adapt.hpp b/include/samurai/mr/adapt.hpp index 4400e98c5..520aa4c12 100644 --- a/include/samurai/mr/adapt.hpp +++ b/include/samurai/mr/adapt.hpp @@ -293,6 +293,7 @@ namespace samurai ghosts_below_cells.apply_op(compute_detail(m_detail, m_fields)); // 'compute_detail' applies 1 level above the set // it is applied to, i.e. level+1 + // 2. detail computation in the ghosts below cells (at level) if (level >= min_level) { if (periodic_in_all_directions) @@ -302,8 +303,6 @@ namespace samurai } else { - // We don't want to compute the detail in the ghosts below the boundary cells. In those ghosts, we want to keep the - // detail to 0. We do that because that detail would use the outer ghost cells at level L-2, which holds the BC // projected 2 times, and this method actually does not work well. So we're removing a layer of 1 boundary cells // from the domain at level L-2. This action ensures that the outer ghost at level L-2 will not be used in the // prediction stencil of interior ghosts. Note: where we don't compute the detail, it stays at its initial value of 0. @@ -311,7 +310,7 @@ namespace samurai // contract the domain only in non-periodic directions auto domain_without_bdry = contract(self(mesh.domain()).on(level - 1), 1, contract_directions); auto cells_without_bdry = intersection(intersection(mesh[mesh_id_t::all_cells][level - 1], domain_without_bdry), - m_interest_cells[level + 1]) + mesh[mesh_id_t::cells][level + 1]) .on(level - 1); cells_without_bdry.apply_op(compute_detail(m_detail, m_fields)); } From d42527723eb0368ccc7947070c6aa5ad773cd2a0 Mon Sep 17 00:00:00 2001 From: Loic Gouarin Date: Wed, 25 Feb 2026 18:18:04 +0100 Subject: [PATCH 3/4] ... --- include/samurai/mr/adapt.hpp | 71 +++++++++++++++++++++++++++++------- 1 file changed, 57 insertions(+), 14 deletions(-) diff --git a/include/samurai/mr/adapt.hpp b/include/samurai/mr/adapt.hpp index 520aa4c12..f66c44ce7 100644 --- a/include/samurai/mr/adapt.hpp +++ b/include/samurai/mr/adapt.hpp @@ -277,16 +277,16 @@ namespace samurai // We compute the detail in the cells and ghosts below the cells, except near the (non-periodic) boundaries, where we compute // the detail only in the cells (justification in the comments below). - // bool periodic_in_all_directions = true; - // std::array contract_directions; - // for (std::size_t d = 0; d < dim; ++d) - // { - // periodic_in_all_directions = periodic_in_all_directions && mesh.is_periodic(d); - // contract_directions[d] = !mesh.is_periodic(d); - // } + bool periodic_in_all_directions = true; + std::array contract_directions; + for (std::size_t d = 0; d < dim; ++d) + { + periodic_in_all_directions = periodic_in_all_directions && mesh.is_periodic(d); + contract_directions[d] = !mesh.is_periodic(d); + } times::timers.start("detail computation"); - for (std::size_t level = ((min_level > 0) ? min_level - 1 : 0); level < max_level; ++level) + for (std::size_t level = ((min_level > 0) ? min_level - 1 : 0); level < max_level - ite; ++level) { // 1. detail computation in the cells (at level+1) auto ghosts_below_cells = intersection(mesh[mesh_id_t::all_cells][level], m_interest_cells[level + 1]).on(level); @@ -310,7 +310,7 @@ namespace samurai // contract the domain only in non-periodic directions auto domain_without_bdry = contract(self(mesh.domain()).on(level - 1), 1, contract_directions); auto cells_without_bdry = intersection(intersection(mesh[mesh_id_t::all_cells][level - 1], domain_without_bdry), - mesh[mesh_id_t::cells][level + 1]) + m_interest_cells[level + 1]) .on(level - 1); cells_without_bdry.apply_op(compute_detail(m_detail, m_fields)); } @@ -326,7 +326,7 @@ namespace samurai update_ghost_subdomains(m_detail); times::timers.start("tag computation"); - for (std::size_t level = min_level; level <= max_level; ++level) + for (std::size_t level = min_level; level <= max_level - ite; ++level) { std::size_t exponent = dim * (max_level - level); double eps_l = cfg.epsilon() / (1 << exponent); @@ -348,7 +348,7 @@ namespace samurai keep_boundary_refined(mesh, m_tag); } - for (std::size_t level = min_level; level <= max_level; ++level) + for (std::size_t level = min_level; level <= max_level - ite; ++level) { auto subset_2 = intersection(m_interest_cells[level], mesh[mesh_id_t::cells][level]); @@ -367,7 +367,7 @@ namespace samurai for (std::size_t level = max_level; level > 0; --level) { - auto keep_subset = intersection(m_interest_cells[level], mesh[mesh_id_t::all_cells][level - 1]).on(level - 1); + auto keep_subset = intersection(mesh[mesh_id_t::cells][level], mesh[mesh_id_t::all_cells][level - 1]).on(level - 1); update_tag_periodic(level, m_tag); update_tag_subdomains(level, m_tag); @@ -392,12 +392,55 @@ namespace samurai times::timers.stop("mesh update"); return true; } + + std::vector> directions; + if (mesh.is_periodic()) + { + std::array nb_cells_finest_level; + const auto& min_indices = mesh.domain().min_indices(); + const auto& max_indices = mesh.domain().max_indices(); + + for (size_t d = 0; d != max_indices.size(); ++d) + { + nb_cells_finest_level[d] = max_indices[d] - min_indices[d]; + } + directions = detail::get_periodic_directions(nb_cells_finest_level, 0, mesh.periodicity()); + } + + cl_type cell_list(mesh.origin_point(), mesh.scaling_factor()); + auto update_interest_cells = [&](std::size_t level, auto&& old_cells, auto&& new_cells) + { + auto& lcl = cell_list[level]; + auto diff = nestedExpand(difference(new_cells, old_cells).on(level - 1), mesh.cfg().prediction_stencil_radius); + auto subset1 = intersection(diff, new_mesh[mesh_id_t::cells][level]); + subset1( + [&](auto& i, auto& index) + { + lcl[index].add_interval(i); + }); + auto& lclm1 = cell_list[level - 1]; + auto subset2 = intersection(diff.on(level - 2), new_mesh[mesh_id_t::cells][level - 1]); + subset2( + [&](auto& i, auto& index) + { + lclm1[index].add_interval(i); + }); + }; + for (std::size_t level = min_level; level <= max_level; ++level) { - auto subset = difference(new_mesh[mesh_id_t::cells][level], mesh[mesh_id_t::cells][level]); - m_interest_cells[level] = subset.to_lca(); + update_interest_cells(level, mesh[mesh_id_t::cells][level], new_mesh[mesh_id_t::cells][level]); + const int delta_l = int(mesh.domain().level() - level); + for (auto& d : directions) + { + update_interest_cells(level, + translate(mesh[mesh_id_t::cells][level], d >> delta_l), + translate(new_mesh[mesh_id_t::cells][level], d >> delta_l)); + } } + m_interest_cells = {cell_list, true}; + times::timers.stop("mesh update"); update_ghost_mr(other_fields...); From 61cbba5bfd3ea117f08fc6b9dafa7745486fc4c8 Mon Sep 17 00:00:00 2001 From: Loic Gouarin Date: Thu, 26 Feb 2026 16:45:56 +0100 Subject: [PATCH 4/4] ... --- include/samurai/mr/adapt.hpp | 79 +++++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 32 deletions(-) diff --git a/include/samurai/mr/adapt.hpp b/include/samurai/mr/adapt.hpp index f66c44ce7..50b4be112 100644 --- a/include/samurai/mr/adapt.hpp +++ b/include/samurai/mr/adapt.hpp @@ -249,7 +249,8 @@ namespace samurai template template - bool Adapt::harten(std::size_t ite, const mra_config& cfg, Fields&... other_fields) + bool + Adapt::harten([[maybe_unused]] std::size_t ite, const mra_config& cfg, Fields&... other_fields) { // ScopedTimer timer(fmt::format("harten criterion {}", ite)); auto& mesh = m_fields.mesh(); @@ -296,24 +297,28 @@ namespace samurai // 2. detail computation in the ghosts below cells (at level) if (level >= min_level) { - if (periodic_in_all_directions) - { - auto ghosts_2_levels_below_cells = intersection(mesh[mesh_id_t::all_cells][level - 1], ghosts_below_cells).on(level - 1); - ghosts_2_levels_below_cells.apply_op(compute_detail(m_detail, m_fields)); - } - else - { - // projected 2 times, and this method actually does not work well. So we're removing a layer of 1 boundary cells - // from the domain at level L-2. This action ensures that the outer ghost at level L-2 will not be used in the - // prediction stencil of interior ghosts. Note: where we don't compute the detail, it stays at its initial value of 0. - - // contract the domain only in non-periodic directions - auto domain_without_bdry = contract(self(mesh.domain()).on(level - 1), 1, contract_directions); - auto cells_without_bdry = intersection(intersection(mesh[mesh_id_t::all_cells][level - 1], domain_without_bdry), - m_interest_cells[level + 1]) - .on(level - 1); - cells_without_bdry.apply_op(compute_detail(m_detail, m_fields)); - } + auto ghosts_2_levels_below_cells = intersection(mesh[mesh_id_t::all_cells][level - 1], ghosts_below_cells).on(level - 1); + ghosts_2_levels_below_cells.apply_op(compute_detail(m_detail, m_fields)); + + // if (periodic_in_all_directions) + // { + // auto ghosts_2_levels_below_cells = intersection(mesh[mesh_id_t::all_cells][level - 1], ghosts_below_cells).on(level - + // 1); ghosts_2_levels_below_cells.apply_op(compute_detail(m_detail, m_fields)); + // } + // else + // { + // // projected 2 times, and this method actually does not work well. So we're removing a layer of 1 boundary cells + // // from the domain at level L-2. This action ensures that the outer ghost at level L-2 will not be used in the + // // prediction stencil of interior ghosts. Note: where we don't compute the detail, it stays at its initial value of + // 0. + + // // contract the domain only in non-periodic directions + // auto domain_without_bdry = contract(self(mesh.domain()).on(level - 1), 1, contract_directions); + // auto cells_without_bdry = intersection(intersection(mesh[mesh_id_t::all_cells][level - 1], domain_without_bdry), + // m_interest_cells[level + 1]) + // .on(level - 1); + // cells_without_bdry.apply_op(compute_detail(m_detail, m_fields)); + // } } } @@ -348,22 +353,22 @@ namespace samurai keep_boundary_refined(mesh, m_tag); } - for (std::size_t level = min_level; level <= max_level - ite; ++level) - { - auto subset_2 = intersection(m_interest_cells[level], mesh[mesh_id_t::cells][level]); + // for (std::size_t level = min_level; level <= max_level - ite; ++level) + // { + // auto subset_2 = intersection(m_interest_cells[level], mesh[mesh_id_t::cells][level]); - subset_2.apply_op(keep_around_refine(m_tag)); + // subset_2.apply_op(keep_around_refine(m_tag)); - if constexpr (enlarge_v) - { - auto subset_3 = intersection(mesh[mesh_id_t::cells_and_ghosts][level], mesh[mesh_id_t::cells_and_ghosts][level]); - subset_2.apply_op(enlarge(m_tag)); - subset_3.apply_op(tag_to_keep<0>(m_tag, CellFlag::enlarge)); - } + // if constexpr (enlarge_v) + // { + // auto subset_3 = intersection(mesh[mesh_id_t::cells_and_ghosts][level], mesh[mesh_id_t::cells_and_ghosts][level]); + // subset_2.apply_op(enlarge(m_tag)); + // subset_3.apply_op(tag_to_keep<0>(m_tag, CellFlag::enlarge)); + // } - update_tag_periodic(level, m_tag); - update_tag_subdomains(level, m_tag); - } + // update_tag_periodic(level, m_tag); + // update_tag_subdomains(level, m_tag); + // } for (std::size_t level = max_level; level > 0; --level) { @@ -437,6 +442,16 @@ namespace samurai translate(mesh[mesh_id_t::cells][level], d >> delta_l), translate(new_mesh[mesh_id_t::cells][level], d >> delta_l)); } + for (auto& neighbour : new_mesh.mpi_neighbourhood()) + { + update_interest_cells(level, mesh[mesh_id_t::cells][level], neighbour.mesh[mesh_id_t::cells][level]); + for (auto& d : directions) + { + update_interest_cells(level, + translate(mesh[mesh_id_t::cells][level], d >> delta_l), + translate(neighbour.mesh[mesh_id_t::cells][level], d >> delta_l)); + } + } } m_interest_cells = {cell_list, true};