Skip to content

fix: mpi and periodic#420

Open
gouarin wants to merge 40 commits intohpc-maths:mainfrom
gouarin:fix-periodic-2
Open

fix: mpi and periodic#420
gouarin wants to merge 40 commits intohpc-maths:mainfrom
gouarin:fix-periodic-2

Conversation

@gouarin
Copy link
Contributor

@gouarin gouarin commented Feb 6, 2026

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_impl method in the mr/mesh.hpp file 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/FiniteVolume called burgers_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

  • I agree to follow this project's Code of Conduct

@gouarin gouarin changed the title Fix: periodic 2 fix: mpi and periodic Feb 6, 2026
Copy link

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

micromamba install -y petsc pkg-config

- name: Install cxx compiler
- name: Installl cxx compiler
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo in step name: "Installl cxx compiler" has an extra 'l'.

Suggested change
- name: Installl cxx compiler
- name: Install cxx compiler

Copilot uses AI. Check for mistakes.
Comment on lines 1083 to 1091
#ifdef SAMURAI_WITH_MPI
mpi::communicator world;
if (m_mpi_neighbourhood.empty() && world.size() > 1)
{
find_neighbourhood();
}
#endif
update_neighbour_subdomain();
}
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Comment on lines 245 to 263
// 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);
});
});
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Comment on lines 323 to 369
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);
});
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Comment on lines +421 to +440
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;
});
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Comment on lines 180 to 187
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;
}
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Comment on lines 57 to 99
@@ -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")
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Copilot uses AI. Check for mistakes.
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
Copy link

Copilot AI Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Suggested change
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

Copilot uses AI. Check for mistakes.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant