From 4d93d1bceac255b92939d448099b6cdf69e47eec Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Wed, 19 Mar 2025 11:23:49 +0100 Subject: [PATCH 01/11] fix: repare MPI in 1D - Previously, in multirank MPI, m_cells was empty due to a bad n_cells per subdomains computation - Then I add a specific partitioning for 1D MPI : - We compute the number of cells in the original mesh - Each subdomain get a part of the nb_cells --- include/samurai/mesh.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index cf8f74e24..5970428da 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -783,6 +783,7 @@ namespace samurai this->m_cells[mesh_id_t::cells][start_level] = {start_level, subdomain_box}; */ +<<<<<<< HEAD std::size_t subdomain_start = 0; std::size_t subdomain_end = 0; lcl_type subdomain_cells(start_level, m_domain.origin_point(), m_domain.scaling_factor()); From bf1f099f52191c0bf9c6cf46c242266cdf72d61d Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Tue, 1 Apr 2025 16:54:24 +0200 Subject: [PATCH 02/11] doc: add partition_mesh comments --- include/samurai/mesh.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 5970428da..cf8f74e24 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -783,7 +783,6 @@ namespace samurai this->m_cells[mesh_id_t::cells][start_level] = {start_level, subdomain_box}; */ -<<<<<<< HEAD std::size_t subdomain_start = 0; std::size_t subdomain_end = 0; lcl_type subdomain_cells(start_level, m_domain.origin_point(), m_domain.scaling_factor()); From fe38bf714bd4482debb62871846b6d512930cfc9 Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Thu, 20 Mar 2025 13:47:09 +0100 Subject: [PATCH 03/11] feat: Detect nighbourhood - Not working FOR NOW - We need to reduce the neighbourhood size. Then, we need to construct it based on each subdomain - We construct a first extended nogihbourhood to get all the potential intersection - And then we intersect each subdomain - **BUT**, we need a kind of "expand" feature. - This expand allows us to detect neighbourhood - Nevertheless, maybe a expand of one is not enough. --- include/samurai/mesh.hpp | 61 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index cf8f74e24..3920604a0 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -833,15 +833,76 @@ namespace samurai this->m_cells[mesh_id_t::cells][start_level] = subdomain_cells; + // std::vector m_mpi_neighbourhood_temp; m_mpi_neighbourhood.reserve(static_cast(size) - 1); + // This code find the neighbourhood of each mpi rank + // 1) we have to exchange the local mesh with all the neighbourhood. It is costly and not scalable but is it a easy way to do it + // 2) we find intersection between each subdomain cells_and_ghosts. + // 3) be careful : for periodic conditions, how to treat neighbourhood ?? maybe in another mpi neighbourhood. + for (int ir = 0; ir < size; ++ir) { if (ir != rank) { m_mpi_neighbourhood.push_back(ir); + ; + } + } + std::cout << " nombre voisins du rank " << rank << " before better neighbour : " << m_mpi_neighbourhood.size() << std::endl; + // send + + // Please think about if they are all necessary + construct_subdomain(); + construct_union(); + update_sub_mesh(); + renumbering(); + update_mesh_neighbour(); + + std::vector m_mpi_neighbourhood_temp; + // std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); + + // m_mpi_neighboor is empty lol... + + for (auto& neighbour : m_mpi_neighbourhood) + { + int est_voisin_global = 0; + for (int level = this->min_level(); level < this->max_level(); level++) + { + // I use cells and ghosts becose cell only is not enough + // in contrary, cells and ghost are overkill + // maybe we can do cells for mesh 1 and cells_and_ghosts for mesh 2 + + // We need a kind of ""expand"" operator. + auto overlap = intersection(this->m_cells[mesh_id_t::reference][level], neighbour.mesh.m_cells[mesh_id_t::reference][level] + //,this->m_cells.subdomain() + ) + .on(level); + // tester si voisin + int est_voisin_niveau = 0; + // maybe it exist a method like "overlap.is_empty() so we should use it" + + overlap( + [&](const auto& i, const auto& index) + { + est_voisin_niveau = 1; + }); + + std::cout << " est voisin local : " << est_voisin_niveau << std::endl; + + if (est_voisin_niveau == 1) + { + est_voisin_global = 1; + } + } + if (est_voisin_global) + { + m_mpi_neighbourhood_temp.push_back(neighbour); } } + std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); + + std::cout << " nombre voisins du rank " << rank << " : " << m_mpi_neighbourhood.size() << std::endl; // // Neighbours // m_mpi_neighbourhood.reserve(static_cast(pow(3, dim) - 1)); // auto neighbour = [&](xt::xtensor_fixed> shift) From 777a2caca9c106dd9260dbf87a1eeaacbe95f2df Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Thu, 20 Mar 2025 15:13:24 +0100 Subject: [PATCH 04/11] fix: translate intersection mesh - Previsouly, we use per-level ghost-and-cells - It is not the right method, we need to use subdomain mehs - Then we use translation to operate as an "expand" --- include/samurai/mesh.hpp | 50 ++++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 22 deletions(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 3920604a0..a4568d0b4 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -860,38 +860,44 @@ namespace samurai std::vector m_mpi_neighbourhood_temp; // std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); - // m_mpi_neighboor is empty lol... - for (auto& neighbour : m_mpi_neighbourhood) { int est_voisin_global = 0; - for (int level = this->min_level(); level < this->max_level(); level++) { // I use cells and ghosts becose cell only is not enough // in contrary, cells and ghost are overkill // maybe we can do cells for mesh 1 and cells_and_ghosts for mesh 2 - // We need a kind of ""expand"" operator. - auto overlap = intersection(this->m_cells[mesh_id_t::reference][level], neighbour.mesh.m_cells[mesh_id_t::reference][level] - //,this->m_cells.subdomain() - ) - .on(level); - // tester si voisin - int est_voisin_niveau = 0; - // maybe it exist a method like "overlap.is_empty() so we should use it" - - overlap( - [&](const auto& i, const auto& index) - { - est_voisin_niveau = 1; - }); - - std::cout << " est voisin local : " << est_voisin_niveau << std::endl; - - if (est_voisin_niveau == 1) + // We need a kind of ""expand"" operator. but for now we use bruteforce translation + xt::xtensor_fixed> stencils{ + {-1, 0 }, + {1, 0 }, + {0, -1}, + {0, 1 } + {-1, 1 }, + {1, 1 }, + {-1, -1}, + {1, 1 } + }; + + for (std::size_t is = 0; is < stencils.shape()[0]; ++is) { - est_voisin_global = 1; + auto s = xt::view(stencils, is); + auto translation = translate(this->m_subdomain, s); + auto overlap = intersection(translation, neighbour.mesh.m_subdomain); + // tester si voisin + int est_voisin_niveau = 0; + // maybe it exist a method like "overlap.is_empty() so we should use it" + overlap( + [&](const auto& i, const auto& index) + { + est_voisin_niveau = 1; + }); + if (est_voisin_niveau == 1) + { + est_voisin_global = 1; + } } } if (est_voisin_global) From ffc0e0752502ee5bd030cafb9879cf56bb4db658 Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Thu, 20 Mar 2025 15:20:32 +0100 Subject: [PATCH 05/11] fix: Use the right stencil array --- include/samurai/mesh.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index a4568d0b4..77e1cd120 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -870,7 +870,7 @@ namespace samurai // maybe we can do cells for mesh 1 and cells_and_ghosts for mesh 2 // We need a kind of ""expand"" operator. but for now we use bruteforce translation - xt::xtensor_fixed> stencils{ + xt::xtensor_fixed> stencils{ {-1, 0 }, {1, 0 }, {0, -1}, From 799463cc5caed50e0edc69bf7c9530de0fe38742 Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Thu, 20 Mar 2025 15:28:27 +0100 Subject: [PATCH 06/11] fix: add comma --- include/samurai/mesh.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 77e1cd120..5b52554c1 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -874,7 +874,7 @@ namespace samurai {-1, 0 }, {1, 0 }, {0, -1}, - {0, 1 } + {0, 1 }, {-1, 1 }, {1, 1 }, {-1, -1}, From dcccec8f004a491131aab2492ee960ce3cea717f Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Thu, 20 Mar 2025 16:16:56 +0100 Subject: [PATCH 07/11] fix: remove cout --- include/samurai/mesh.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 5b52554c1..894b1067f 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -848,7 +848,7 @@ namespace samurai ; } } - std::cout << " nombre voisins du rank " << rank << " before better neighbour : " << m_mpi_neighbourhood.size() << std::endl; + // send // Please think about if they are all necessary @@ -908,7 +908,6 @@ namespace samurai std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); - std::cout << " nombre voisins du rank " << rank << " : " << m_mpi_neighbourhood.size() << std::endl; // // Neighbours // m_mpi_neighbourhood.reserve(static_cast(pow(3, dim) - 1)); // auto neighbour = [&](xt::xtensor_fixed> shift) From 984500b672cd698c29f98636b70c5a95135e991c Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Tue, 1 Apr 2025 16:22:56 +0200 Subject: [PATCH 08/11] fix: add find_subdomains in 3D - This is not really beautiful, it needs more work - Todo : - arbitrary nD - create a new special function --- include/samurai/mesh.hpp | 77 ++++++++++++++++++++++++++++++++++------ 1 file changed, 67 insertions(+), 10 deletions(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 894b1067f..4e047a57d 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -714,6 +714,69 @@ namespace samurai } } + template + struct StencilProvider + { + // Optionnel : Déclencher une erreur si une dimension non supportée est utilisée + static_assert(Dim == 2 || Dim == 3, "Dimension non supportée pour les stencils"); + // Ou laisser vide pour obtenir une erreur de symbole non défini si non spécialisé + }; + + // Spécialisation pour dim = 2 + template <> + struct StencilProvider<2> + { + static constexpr std::size_t stencil_size = 8; + inline static const xt::xtensor_fixed> stencils{ + {-1, 0 }, + {1, 0 }, + {0, -1}, + {0, 1 }, + {-1, 1 }, + {1, 1 }, + {-1, -1}, + {1, -1} // Corrigé + }; + // Autres constantes/types dépendant de dim=2 pourraient aller ici + }; + + // Spécialisation pour dim = 3 + template <> + struct StencilProvider<3> + { + static constexpr std::size_t stencil_size = 27; // Ou 26 si on exclut le centre + inline static const xt::xtensor_fixed> stencils{ + {-1, -1, -1}, + {-1, -1, 0 }, + {-1, -1, 1 }, + {-1, 0, -1}, + {-1, 0, 0 }, + {-1, 0, 1 }, + {-1, 1, -1}, + {-1, 1, 0 }, + {-1, 1, 1 }, + {0, -1, -1}, + {0, -1, 0 }, + {0, -1, 1 }, + {0, 0, -1}, + {0, 0, 0 }, + {0, 0, 1 }, // Inclut le centre + {0, 1, -1}, + {0, 1, 0 }, + {0, 1, 1 }, + {1, -1, -1}, + {1, -1, 0 }, + {1, -1, 1 }, + {1, 0, -1}, + {1, 0, 0 }, + {1, 0, 1 }, + {1, 1, -1}, + {1, 1, 0 }, + {1, 1, 1 } + }; + // Autres constantes/types dépendant de dim=3 pourraient aller ici + }; + template void Mesh_base::partition_mesh([[maybe_unused]] std::size_t start_level, [[maybe_unused]] const Box& global_box) { @@ -861,6 +924,9 @@ namespace samurai std::vector m_mpi_neighbourhood_temp; // std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); // m_mpi_neighboor is empty lol... + // + const auto& stencils = StencilProvider::stencils; + for (auto& neighbour : m_mpi_neighbourhood) { int est_voisin_global = 0; @@ -870,16 +936,6 @@ namespace samurai // maybe we can do cells for mesh 1 and cells_and_ghosts for mesh 2 // We need a kind of ""expand"" operator. but for now we use bruteforce translation - xt::xtensor_fixed> stencils{ - {-1, 0 }, - {1, 0 }, - {0, -1}, - {0, 1 }, - {-1, 1 }, - {1, 1 }, - {-1, -1}, - {1, 1 } - }; for (std::size_t is = 0; is < stencils.shape()[0]; ++is) { @@ -897,6 +953,7 @@ namespace samurai if (est_voisin_niveau == 1) { est_voisin_global = 1; + break; } } } From b883ef7cba31451f8a69e99d196ea30747836a3d Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Tue, 1 Apr 2025 16:38:22 +0200 Subject: [PATCH 09/11] fix: Create new function - Now we have a dedicated function to find neighbourhood --- include/samurai/mesh.hpp | 157 ++++++++++++++++++++------------------- 1 file changed, 81 insertions(+), 76 deletions(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 4e047a57d..d0bd237ba 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -170,6 +170,7 @@ namespace samurai void construct_union(); void update_sub_mesh(); void renumbering(); + void find_neighbour(); void partition_mesh(std::size_t start_level, const Box& global_box); void load_balancing(); void load_transfer(const std::vector& load_fluxes); @@ -717,9 +718,7 @@ namespace samurai template struct StencilProvider { - // Optionnel : Déclencher une erreur si une dimension non supportée est utilisée static_assert(Dim == 2 || Dim == 3, "Dimension non supportée pour les stencils"); - // Ou laisser vide pour obtenir une erreur de symbole non défini si non spécialisé }; // Spécialisation pour dim = 2 @@ -735,16 +734,15 @@ namespace samurai {-1, 1 }, {1, 1 }, {-1, -1}, - {1, -1} // Corrigé + {1, -1} }; - // Autres constantes/types dépendant de dim=2 pourraient aller ici }; // Spécialisation pour dim = 3 template <> struct StencilProvider<3> { - static constexpr std::size_t stencil_size = 27; // Ou 26 si on exclut le centre + static constexpr std::size_t stencil_size = 27; inline static const xt::xtensor_fixed> stencils{ {-1, -1, -1}, {-1, -1, 0 }, @@ -759,8 +757,7 @@ namespace samurai {0, -1, 0 }, {0, -1, 1 }, {0, 0, -1}, - {0, 0, 0 }, - {0, 0, 1 }, // Inclut le centre + {0, 0, 1 }, {0, 1, -1}, {0, 1, 0 }, {0, 1, 1 }, @@ -774,9 +771,84 @@ namespace samurai {1, 1, 0 }, {1, 1, 1 } }; - // Autres constantes/types dépendant de dim=3 pourraient aller ici }; + template + void Mesh_base::find_neighbour() + { + mpi::communicator world; + auto rank = world.rank(); + auto size = world.size(); + + m_mpi_neighbourhood.reserve(static_cast(size) - 1); + // This code find the neighbourhood of each mpi rank + // 1) we have to exchange the local mesh with all the neighbourhood. It is costly and not scalable but is it a easy way to do it + // 2) we find intersection between each subdomain cells_and_ghosts. + // 3) be careful : for periodic conditions, how to treat neighbourhood ?? maybe in another mpi neighbourhood. + + for (int ir = 0; ir < size; ++ir) + { + if (ir != rank) + { + m_mpi_neighbourhood.push_back(ir); + ; + } + } + + // send + + // Please think about if they are all necessary + construct_subdomain(); + construct_union(); + update_sub_mesh(); + renumbering(); + update_mesh_neighbour(); + + std::vector m_mpi_neighbourhood_temp; + // std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); + // m_mpi_neighboor is empty lol... + // + const auto& stencils = StencilProvider::stencils; + + for (auto& neighbour : m_mpi_neighbourhood) + { + int est_voisin_global = 0; + { + // I use cells and ghosts becose cell only is not enough + // in contrary, cells and ghost are overkill + // maybe we can do cells for mesh 1 and cells_and_ghosts for mesh 2 + + // We need a kind of ""expand"" operator. but for now we use bruteforce translation + + for (std::size_t is = 0; is < stencils.shape()[0]; ++is) + { + auto s = xt::view(stencils, is); + auto translation = translate(this->m_subdomain, s); + auto overlap = intersection(translation, neighbour.mesh.m_subdomain); + // tester si voisin + int est_voisin_niveau = 0; + // maybe it exist a method like "overlap.is_empty() so we should use it" + overlap( + [&](const auto& i, const auto& index) + { + est_voisin_niveau = 1; + }); + if (est_voisin_niveau == 1) + { + est_voisin_global = 1; + break; + } + } + } + if (est_voisin_global) + { + m_mpi_neighbourhood_temp.push_back(neighbour); + } + } + + std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); + } + template void Mesh_base::partition_mesh([[maybe_unused]] std::size_t start_level, [[maybe_unused]] const Box& global_box) { @@ -896,74 +968,7 @@ namespace samurai this->m_cells[mesh_id_t::cells][start_level] = subdomain_cells; - // std::vector m_mpi_neighbourhood_temp; - m_mpi_neighbourhood.reserve(static_cast(size) - 1); - // This code find the neighbourhood of each mpi rank - // 1) we have to exchange the local mesh with all the neighbourhood. It is costly and not scalable but is it a easy way to do it - // 2) we find intersection between each subdomain cells_and_ghosts. - // 3) be careful : for periodic conditions, how to treat neighbourhood ?? maybe in another mpi neighbourhood. - - for (int ir = 0; ir < size; ++ir) - { - if (ir != rank) - { - m_mpi_neighbourhood.push_back(ir); - ; - } - } - - // send - - // Please think about if they are all necessary - construct_subdomain(); - construct_union(); - update_sub_mesh(); - renumbering(); - update_mesh_neighbour(); - - std::vector m_mpi_neighbourhood_temp; - // std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); - // m_mpi_neighboor is empty lol... - // - const auto& stencils = StencilProvider::stencils; - - for (auto& neighbour : m_mpi_neighbourhood) - { - int est_voisin_global = 0; - { - // I use cells and ghosts becose cell only is not enough - // in contrary, cells and ghost are overkill - // maybe we can do cells for mesh 1 and cells_and_ghosts for mesh 2 - - // We need a kind of ""expand"" operator. but for now we use bruteforce translation - - for (std::size_t is = 0; is < stencils.shape()[0]; ++is) - { - auto s = xt::view(stencils, is); - auto translation = translate(this->m_subdomain, s); - auto overlap = intersection(translation, neighbour.mesh.m_subdomain); - // tester si voisin - int est_voisin_niveau = 0; - // maybe it exist a method like "overlap.is_empty() so we should use it" - overlap( - [&](const auto& i, const auto& index) - { - est_voisin_niveau = 1; - }); - if (est_voisin_niveau == 1) - { - est_voisin_global = 1; - break; - } - } - } - if (est_voisin_global) - { - m_mpi_neighbourhood_temp.push_back(neighbour); - } - } - - std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); + find_neighbour(); // // Neighbours // m_mpi_neighbourhood.reserve(static_cast(pow(3, dim) - 1)); From 3f987513ddb41b125797f6b8e8e80c5bcd8b2bc1 Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Tue, 1 Apr 2025 16:41:37 +0200 Subject: [PATCH 10/11] fix: clean code --- include/samurai/mesh.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index d0bd237ba..868a2c405 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -805,11 +805,8 @@ namespace samurai update_mesh_neighbour(); std::vector m_mpi_neighbourhood_temp; - // std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood); - // m_mpi_neighboor is empty lol... - // - const auto& stencils = StencilProvider::stencils; + const auto& stencils = StencilProvider::stencils; for (auto& neighbour : m_mpi_neighbourhood) { int est_voisin_global = 0; @@ -968,6 +965,7 @@ namespace samurai this->m_cells[mesh_id_t::cells][start_level] = subdomain_cells; + // in theory, at this step, we can construct the neighbourhood by construction. find_neighbour(); // // Neighbours From e0c12375e72dbc1ad50a7e2be61636ccebace065 Mon Sep 17 00:00:00 2001 From: sbstndb/sbstndbs Date: Tue, 1 Apr 2025 16:47:12 +0200 Subject: [PATCH 11/11] fix: use Constexpr - We can use constexrp here. --- include/samurai/mesh.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 868a2c405..add400cb8 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -919,7 +919,7 @@ namespace samurai std::size_t subdomain_end = 0; lcl_type subdomain_cells(start_level, m_domain.origin_point(), m_domain.scaling_factor()); // in 1D MPI, we need a specific partitioning - if (dim == 1) + if constexpr (dim == 1) { std::size_t n_cells = m_domain.nb_cells(); std::size_t n_cells_per_subdomain = n_cells / static_cast(size); @@ -942,7 +942,7 @@ namespace samurai } }); } - else if (dim >= 2) + else if constexpr (dim >= 2) { auto subdomain_nb_intervals = m_domain.nb_intervals() / static_cast(size); subdomain_start = static_cast(rank) * subdomain_nb_intervals;