Conversation
There was a problem hiding this comment.
Pull request overview
This PR targets correctness of ghost-cell construction for multiresolution meshes in MPI and periodic configurations, and adds an MPI demo/CI coverage to prevent regressions.
Changes:
- Reworks
MRMesh::update_sub_mesh_impl()to build ghost/interface cells using a unified “expand + intersect” approach for single-domain, periodic, and MPI-neighbour cases. - Extends the graduation logic to account for MPI neighbour meshes when determining refinements near boundaries.
- Adds an MPI Burgers OSMP 2D demo and wires it into the build/CI for MPI coverage.
Reviewed changes
Copilot reviewed 14 out of 14 changed files in this pull request and generated 8 comments.
Show a summary per file
| File | Description |
|---|---|
include/samurai/mr/mesh.hpp |
Major rewrite of sub-mesh/ghost construction for MPI + periodic handling, with timers added. |
include/samurai/mesh.hpp |
API/constructor adjustments (non-const cfg(), find_neighbourhood() access), and neighbour exchange behavior changes. |
include/samurai/algorithm/graduation.hpp |
Updates graduation refinement discovery to incorporate MPI neighbour mesh data. |
include/samurai/algorithm/update.hpp |
Reorders periodic/subdomain ghost updates during MR ghost update. |
include/samurai/schemes/fv/FV_scheme.hpp |
Changes behavior for non-box periodic BC application (now silent return). |
demos/FiniteVolume/burgers_os_2d_mpi.cpp |
New MPI demo used as regression/weak-scaling check. |
demos/FiniteVolume/convection_nonlinear_osmp.hpp |
Adds OSMP convection helper used by the new demo. |
demos/FiniteVolume/CMakeLists.txt |
Adds conditional MPI demo target when MPI is enabled. |
.github/workflows/ci.yml |
CI matrix tweaks and adds MPI demo execution; also alters MPI job runner setup. |
python/compare.py |
Removes --tol option and hardcodes tolerance; changes failure reporting. |
include/samurai/static_algorithm.hpp |
Adds compile-time nestedExpand<width>() helper. |
include/samurai/io/hdf5.hpp |
Adjusts HighFive include placement and removes macro block. |
include/samurai/utils.hpp |
Include order tweak. |
include/samurai/field/field_base.hpp |
Fixes dependent type usage (typename). |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
.github/workflows/ci.yml
Outdated
| micromamba install -y petsc pkg-config | ||
|
|
||
| - name: Install cxx compiler | ||
| - name: Installl cxx compiler |
There was a problem hiding this comment.
Typo in step name: "Installl cxx compiler" has an extra 'l'.
| - name: Installl cxx compiler | |
| - name: Install cxx compiler |
| #ifdef SAMURAI_WITH_MPI | ||
| mpi::communicator world; | ||
| if (m_mpi_neighbourhood.empty() && world.size() > 1) | ||
| { | ||
| find_neighbourhood(); | ||
| } | ||
| #endif | ||
| update_neighbour_subdomain(); | ||
| } |
There was a problem hiding this comment.
construct_subdomain() now only exchanges subdomain boxes via update_neighbour_subdomain(), but later code (e.g., MRMesh::update_sub_mesh_impl) relies on neighbour.mesh[mesh_id_t::cells] being populated. Consider also exchanging mesh_id_t::cells here (after find_neighbourhood()) so neighbour meshes have the actual cell layout before ghost/interface computations.
| // Manage neighbourhood for MPI | ||
| // ============================== | ||
| // We do the same with the subdomains of the neighbouring processes | ||
| // to be sure that we have all the ghost cells and the cells at the interface | ||
| // with the other processes. | ||
| for (auto& neighbour : this->mpi_neighbourhood()) | ||
| { | ||
| for (std::size_t level = max_level; level >= ((min_level == 0) ? 1 : min_level); --level) | ||
| for_each_level(neighbour.mesh[mesh_id_t::cells], | ||
| [&](std::size_t level) | ||
| { | ||
| lcl_type& lcl = cell_list[level]; | ||
| auto set = intersection(nestedExpand(neighbour.mesh[mesh_id_t::cells][level], ghost_width()), | ||
| nestedExpand(self(this->subdomain()).on(level), ghost_width())); | ||
| set( | ||
| [&](const auto& interval, const auto& index) | ||
| { | ||
| lcl[index].add_interval(interval); | ||
| }); | ||
| }); |
There was a problem hiding this comment.
In update_sub_mesh_impl(), the MPI ghost construction iterates over neighbour.mesh[cells], but neighbour meshes are not guaranteed to have their cells mesh id exchanged/populated (update_mesh_neighbour() calls were removed and update_meshid_neighbour(cells) is invoked before the neighbourhood is computed in constructors). This makes the intersection sets empty and can leave interface ghost cells missing. Consider exchanging mesh_id_t::cells with neighbours after find_neighbourhood() (or at the start of update_sub_mesh_impl) so neighbour.mesh[cells] is valid before it is used here.
| auto add_prediction_ghosts = [&](auto&& subset_cells, auto&& subset_cells_and_ghosts, auto level) | ||
| { | ||
| lcl_type& lcl_m1 = cell_list[level - 1]; | ||
|
|
||
| expr_2( | ||
| subset_cells_and_ghosts( | ||
| [&](const auto& interval, const auto& index_yz) | ||
| { | ||
| if (level - 1 > 0) | ||
| { | ||
| lcl_type& lcl = cell_list[level - 2]; | ||
|
|
||
| static_nested_loop<dim - 1, -config::prediction_stencil_radius, config::prediction_stencil_radius + 1>( | ||
| [&](auto stencil) | ||
| { | ||
| auto new_interval = interval >> 2; | ||
| lcl[(index_yz >> 2) + stencil].add_interval({new_interval.start - config::prediction_stencil_radius, | ||
| new_interval.end + config::prediction_stencil_radius}); | ||
| }); | ||
| } | ||
| lcl_m1[index_yz].add_interval(interval); | ||
| }); | ||
| } | ||
| this->cells()[mesh_id_t::all_cells] = {cell_list, false}; | ||
|
|
||
| this->update_meshid_neighbour(mesh_id_t::cells_and_ghosts); | ||
| this->update_meshid_neighbour(mesh_id_t::reference); | ||
|
|
||
| for (auto& neighbour : this->mpi_neighbourhood()) | ||
| { | ||
| for (std::size_t level = 0; level <= this->max_level(); ++level) | ||
| if (level - 1 > 0) | ||
| { | ||
| auto expr = intersection(this->subdomain(), neighbour.mesh[mesh_id_t::reference][level]).on(level); | ||
| expr( | ||
| lcl_type& lcl_m2 = cell_list[level - 2]; | ||
| subset_cells( | ||
| [&](const auto& interval, const auto& index_yz) | ||
| { | ||
| lcl_type& lcl = cell_list[level]; | ||
| lcl[index_yz].add_interval(interval); | ||
|
|
||
| if (level > neighbour.mesh[mesh_id_t::reference].min_level()) | ||
| { | ||
| lcl_type& lclm1 = cell_list[level - 1]; | ||
|
|
||
| static_nested_loop<dim - 1, -config::prediction_stencil_radius, config::prediction_stencil_radius + 1>( | ||
| [&](auto stencil) | ||
| { | ||
| auto new_interval = interval >> 1; | ||
| lclm1[(index_yz >> 1) + stencil].add_interval({new_interval.start - config::prediction_stencil_radius, | ||
| new_interval.end + config::prediction_stencil_radius}); | ||
| }); | ||
| } | ||
| lcl_m2[index_yz].add_interval(interval); | ||
| }); | ||
| } | ||
| } | ||
| this->cells()[mesh_id_t::all_cells] = {cell_list, false}; | ||
|
|
||
| using box_t = Box<typename interval_t::value_t, dim>; | ||
|
|
||
| const auto& domain = this->domain(); | ||
| const auto& domain_min_indices = domain.min_indices(); | ||
| const auto& domain_max_indices = domain.max_indices(); | ||
|
|
||
| const auto& subdomain = this->subdomain(); | ||
| const auto& subdomain_min_indices = subdomain.min_indices(); | ||
| const auto& subdomain_max_indices = subdomain.max_indices(); | ||
|
|
||
| // add ghosts for periodicity | ||
| xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> shift; | ||
| xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> min_corner; | ||
| xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> max_corner; | ||
|
|
||
| shift.fill(0); | ||
|
|
||
| #ifdef SAMURAI_WITH_MPI | ||
| std::size_t reference_max_level = mpi::all_reduce(world, | ||
| this->cells()[mesh_id_t::reference].max_level(), | ||
| mpi::maximum<std::size_t>()); | ||
| std::size_t reference_min_level = mpi::all_reduce(world, | ||
| this->cells()[mesh_id_t::reference].min_level(), | ||
| mpi::minimum<std::size_t>()); | ||
|
|
||
| std::vector<ca_type> neighbourhood_extended_subdomain(this->mpi_neighbourhood().size()); | ||
| for (size_t neighbor_id = 0; neighbor_id != neighbourhood_extended_subdomain.size(); ++neighbor_id) | ||
| { | ||
| const auto& neighbor_subdomain = this->mpi_neighbourhood()[neighbor_id].mesh.subdomain(); | ||
| if (not neighbor_subdomain.empty()) | ||
| { | ||
| const auto& neighbor_min_indices = neighbor_subdomain.min_indices(); | ||
| const auto& neighbor_max_indices = neighbor_subdomain.max_indices(); | ||
| for (std::size_t level = reference_min_level; level <= reference_max_level; ++level) | ||
| { | ||
| const std::size_t delta_l = subdomain.level() - level; | ||
| box_t box; | ||
| for (std::size_t d = 0; d < dim; ++d) | ||
| { | ||
| box.min_corner()[d] = (neighbor_min_indices[d] >> delta_l) - ghost_width(); | ||
| box.max_corner()[d] = (neighbor_max_indices[d] >> delta_l) + ghost_width(); | ||
| } | ||
| neighbourhood_extended_subdomain[neighbor_id][level] = {level, box}; | ||
| } | ||
| } | ||
| } | ||
| #endif // SAMURAI_WITH_MPI | ||
| const auto& mesh_ref = this->cells()[mesh_id_t::reference]; | ||
| for (std::size_t level = 0; level <= this->max_level(); ++level) | ||
| { | ||
| const std::size_t delta_l = subdomain.level() - level; | ||
| lcl_type& lcl = cell_list[level]; | ||
| }; | ||
|
|
||
| for (std::size_t d = 0; d < dim; ++d) | ||
| for_each_level( | ||
| this->cells()[mesh_id_t::cells], | ||
| [&](std::size_t level) | ||
| { | ||
| min_corner[d] = (subdomain_min_indices[d] >> delta_l) - ghost_width(); | ||
| max_corner[d] = (subdomain_max_indices[d] >> delta_l) + ghost_width(); | ||
| } | ||
| lca_type lca_extended_subdomain(level, box_t(min_corner, max_corner)); | ||
| for (std::size_t d = 0; d < dim; ++d) | ||
| { | ||
| if (this->is_periodic(d)) | ||
| // own part | ||
| add_prediction_ghosts( | ||
| nestedExpand<config::prediction_stencil_radius>(self(this->cells()[mesh_id_t::cells][level]).on(level - 2)), | ||
| intersection(nestedExpand<config::prediction_stencil_radius>( | ||
| self(this->cells()[mesh_id_t::cells_and_ghosts][level]).on(level - 1)), | ||
| nestedExpand<config::prediction_stencil_radius>(self(this->subdomain()).on(level - 1))), | ||
| level); | ||
|
|
||
| // periodic part | ||
| const int delta_l = int(this->domain().level() - level); | ||
| for (const auto& d : directions) | ||
| { | ||
| shift[d] = (domain_max_indices[d] - domain_min_indices[d]) >> delta_l; | ||
|
|
||
| min_corner[d] = (domain_min_indices[d] >> delta_l) - ghost_width(); | ||
| max_corner[d] = (domain_min_indices[d] >> delta_l) + ghost_width(); | ||
|
|
||
| lca_type lca_min(level, box_t(min_corner, max_corner)); | ||
|
|
||
| min_corner[d] = (domain_max_indices[d] >> delta_l) - ghost_width(); | ||
| max_corner[d] = (domain_max_indices[d] >> delta_l) + ghost_width(); | ||
|
|
||
| lca_type lca_max(level, box_t(min_corner, max_corner)); | ||
|
|
||
| auto set1 = intersection(translate(intersection(mesh_ref[level], lca_min), shift), | ||
| intersection(lca_extended_subdomain, lca_max)); | ||
| set1( | ||
| [&](const auto& i, const auto& index_yz) | ||
| { | ||
| lcl[index_yz].add_interval(i); | ||
| }); | ||
| auto set2 = intersection(translate(intersection(mesh_ref[level], lca_max), -shift), | ||
| intersection(lca_extended_subdomain, lca_min)); | ||
| set2( | ||
| [&](const auto& i, const auto& index_yz) | ||
| { | ||
| lcl[index_yz].add_interval(i); | ||
| }); | ||
| #ifdef SAMURAI_WITH_MPI | ||
| //~ for (const auto& mpi_neighbor : this->mpi_neighbourhood()) | ||
| for (size_t neighbor_id = 0; neighbor_id != this->mpi_neighbourhood().size(); ++neighbor_id) | ||
| { | ||
| const auto& mpi_neighbor = this->mpi_neighbourhood()[neighbor_id]; | ||
| const auto& neighbor_mesh_ref = mpi_neighbor.mesh[mesh_id_t::reference]; | ||
|
|
||
| auto set1_mpi = intersection(translate(intersection(neighbor_mesh_ref[level], lca_min), shift), | ||
| intersection(lca_extended_subdomain, lca_max)); | ||
| set1_mpi( | ||
| [&](const auto& i, const auto& index_yz) | ||
| { | ||
| lcl[index_yz].add_interval(i); | ||
| }); | ||
| auto set2_mpi = intersection(translate(intersection(neighbor_mesh_ref[level], lca_max), -shift), | ||
| intersection(lca_extended_subdomain, lca_min)); | ||
| set2_mpi( | ||
| [&](const auto& i, const auto& index_yz) | ||
| { | ||
| lcl[index_yz].add_interval(i); | ||
| }); | ||
| } | ||
| #endif // SAMURAI_WITH_MPI | ||
| this->cells()[mesh_id_t::all_cells][level] = {lcl}; | ||
|
|
||
| /* reset variables for next iterations. */ | ||
| shift[d] = 0; | ||
| min_corner[d] = (subdomain_min_indices[d] >> delta_l) - ghost_width(); | ||
| max_corner[d] = (subdomain_max_indices[d] >> delta_l) + ghost_width(); | ||
| add_prediction_ghosts( | ||
| intersection(nestedExpand<config::prediction_stencil_radius>( | ||
| translate(this->cells()[mesh_id_t::cells][level], d >> delta_l).on(level - 2)), | ||
| self(this->subdomain()).on(level - 2)), | ||
| intersection(nestedExpand<config::prediction_stencil_radius>( | ||
| translate(this->cells()[mesh_id_t::cells_and_ghosts][level], d >> delta_l).on(level - 1)), | ||
| self(this->subdomain()).on(level - 1)), | ||
| level); | ||
| } | ||
| } | ||
| } | ||
| for (std::size_t level = 0; level < max_level; ++level) | ||
| { | ||
| lcl_type& lcl = cell_list[level + 1]; | ||
| lcl_type lcl_proj{level}; | ||
| auto expr = intersection(this->cells()[mesh_id_t::all_cells][level], this->get_union()[level]); | ||
|
|
||
| expr( | ||
| [&](const auto& interval, const auto& index_yz) | ||
| { | ||
| static_nested_loop<dim - 1, 0, 2>( | ||
| [&](auto s) | ||
| { | ||
| lcl[(index_yz << 1) + s].add_interval(interval << 1); | ||
| }); | ||
| lcl_proj[index_yz].add_interval(interval); | ||
| }); | ||
| this->cells()[mesh_id_t::all_cells][level + 1] = lcl; | ||
| this->cells()[mesh_id_t::proj_cells][level] = lcl_proj; | ||
| } | ||
| this->update_neighbour_subdomain(); | ||
| this->update_meshid_neighbour(mesh_id_t::all_cells); | ||
| }); |
There was a problem hiding this comment.
The prediction-ghost logic uses level - 1 and level - 2 (e.g., cell_list[level - 1], .on(level - 2)) inside a for_each_level() over cells. Since for_each_level() starts at cells.min_level(), this will underflow when the mesh contains level 0/1 (common in MPI/periodic runs) and can index invalid levels. Please guard these operations (e.g., only run this block for level >= 1 / level >= 2) or change the iteration bounds to match the previous code’s safe lower limit.
| for_each_level(this->cells()[mesh_id_t::reference], | ||
| [&](auto level) | ||
| { | ||
| lcl_type& lcl = cell_list[level + 1]; | ||
| lcl_type lcl_proj{level}; | ||
| auto expr = intersection(this->cells()[mesh_id_t::reference][level], this->get_union()[level]); | ||
|
|
||
| expr( | ||
| [&](const auto& interval, const auto& index_yz) | ||
| { | ||
| static_nested_loop<dim - 1, 0, 2>( | ||
| [&](auto s) | ||
| { | ||
| lcl[(index_yz << 1) + s].add_interval(interval << 1); | ||
| }); | ||
| lcl_proj[index_yz].add_interval(interval); | ||
| }); | ||
| this->cells()[mesh_id_t::reference][level + 1] = lcl; | ||
| this->cells()[mesh_id_t::proj_cells][level] = lcl_proj; | ||
| }); |
There was a problem hiding this comment.
The projection/prolongation loop iterates with for_each_level(this->cells()[reference], ...) and then writes to cell_list[level + 1] / reference[level + 1]. If reference.max_level() can reach the maximum allowed refinement level, level + 1 can exceed the valid level range. It would be safer to iterate only up to max_level - 1 (or check level + 1 <= this->max_level() / container max) before accessing/writing level + 1.
| if (!input_field.mesh().domain().is_box()) | ||
| { | ||
| if (input_field.mesh().is_periodic()) | ||
| { | ||
| std::cerr << "Error: apply_directional_bc() not implemented for non-box domains with periodic directions." << std::endl; | ||
| assert(false); | ||
| // std::cerr << "Error: apply_directional_bc() not implemented for non-box domains with periodic directions." << | ||
| // std::endl; assert(false); | ||
| return; | ||
| } |
There was a problem hiding this comment.
apply_directional_bc(): for non-box domains with periodic directions, the previous hard failure was replaced by a silent return. This can hide unsupported configurations and leave boundary/ghost values uninitialized, producing incorrect results without any diagnostic. Please keep a clear runtime error (exception/assert) or implement the missing behavior so callers don’t proceed silently in an invalid state.
| @@ -67,7 +64,7 @@ def compare_meshes(file1, file2, tol): | |||
| index2 = np.argsort(np.asarray([c.tobytes() for c in cells2])) | |||
|
|
|||
| if np.any(cells1.shape != cells2.shape): | |||
| print("shape are not compatibles") | |||
| print("shape are not compatible") | |||
| print(f"{cells1.shape} vs {cells2.shape}") | |||
| sys.exit(f"files {file1} and {file2} are different") | |||
|
|
|||
| @@ -78,6 +75,7 @@ def compare_meshes(file1, file2, tol): | |||
| field1 = construct_fields(mesh1) | |||
| field2 = construct_fields(mesh2) | |||
|
|
|||
| tol = 1e-14 | |||
| for field in field1.keys(): | |||
| if not field in field2.keys(): | |||
| print(f"{field} is not in second file") | |||
| @@ -89,25 +87,23 @@ def compare_meshes(file1, file2, tol): | |||
| np.abs(field1[field][:][index1] - field2[field][:][index2]) > tol | |||
| ) | |||
| # ind = np.where(field1[field][:][index1] != field2[field][:][index2]) | |||
| print(f"-- Values {file1}") | |||
| print(field1[field][:][index1[ind]]) | |||
| print(f"-- Values {file2}") | |||
| print(field2[field][:][index2[ind]]) | |||
| print("-- Difference:") | |||
| print(np.abs(field1[field][:][index1[ind]] - field2[field][:][index2[ind]])) | |||
| print("-- Coordinates (x,y,z):") | |||
| print( | |||
| field1[field][:][index1[ind]], | |||
| field2[field][:][index2[ind]], | |||
| np.abs(field1[field][:][index1[ind]] - field2[field][:][index2[ind]]), | |||
| ) | |||
| print(cells1[index1[ind]], cells2[index2[ind]]) | |||
| # print(np.abs(field1[field][:][index1[ind]]-field2[field][:][index2[ind]])) | |||
| sys.exit(f"{field} is not the same between {file1} and {file2}") | |||
| print(f"{field} is not the same") | |||
| sys.exit() | |||
| print(f"files {file1} and {file2} are the same") | |||
There was a problem hiding this comment.
compare.py hardcodes tol = 1e-14 and removes the --tol CLI option. This makes CI/debug workflows less flexible and also changes failure reporting: sys.exit() is called without a message, losing the file/field context. Consider keeping --tol (with the current default) and preserving an informative non-zero exit message (e.g., include field name and file paths).
.github/workflows/ci.yml
Outdated
| python ../python/compare.py ldc_ink_size_1 ldc_ink_size_2 --tol 1e-11 | ||
|
|
||
| mpiexec -n 2 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=2 --npy=1 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 | ||
| mpiexec -n 2 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=2 --npy=1 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 |
There was a problem hiding this comment.
The MPI weak-scaling test runs the exact same mpiexec -n 2 ... --npx=2 --npy=1 ... command twice in a row. If this duplication isn’t intentional (e.g., warm-up), removing one invocation will reduce CI time without changing coverage.
| mpiexec -n 2 ./demos/FiniteVolume/finite-volume-burgers-os-2d-mpi --npx=2 --npy=1 --min-level=0 --max-level=8 --Tf=6 --nfiles 1 |
Description
This PR fixes a bug that was introduced earlier for the MPI usage and the periodic. The ghost cells were not set properly.
The new implementation of the
update_sub_mesh_implmethod in themr/mesh.hppfile correctly computes all the ghost cells used during the adaptation step and the integration of the numerical scheme.We can see this by observing that we use exactly the same algorithm to add the ghost cells in cases with one domain without periodic conditions, one domain with periodic conditions, and several domains with or without periodic conditions. This makes sense, as for the subdomain or the periodic parts, we simply translate a domain with true cells and add the ghost cells to the subdomain accordingly. The procedure is the same in both cases.
In a future version, AMR mesh will be updated using the same idea.
One of the key elements of this PR is the intense use of the expand operator, which adds one or more cells in all directions of a subset. Unfortunately, the subset implementation is not optimized for this kind of operator, as many unions are involved. The subset implementation is currently being rewritten and should reduce the time required to compute the expanded subset.
Related issue
The MPI version, both with and without periodic, was broken. The graduation method was incorrect because it did not take into account all the levels that could potentially add new cells. And finally, the constraint defined in the PR #320 was not parallel. This PR fixes this issue.
How has this been tested?
To validate the MPI version of samurai and to prevent any regression, a test case is added in
demos/FiniteVolumecalledburgers_os_2d_mpi.cpp. The main idea is to replicate the initial solution on each subdomain and to add periodic conditions. Then, we solve the Burgers equation using an OSMP scheme, checking that the mesh on each subdomain is consistent. This example can also be used to measure the weak scaling in samurai.Code of Conduct
By submitting this PR, you agree to follow our Code of Conduct