Skip to content

fix: boundary conditions for 3d case #441

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
gouarin wants to merge 4 commits into
base: main
Choose a base branch
from
Open

fix: boundary conditions for 3d case #441

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8bcc930
fix
gouarin May 4, 2026
09039ed
...
gouarin May 4, 2026
44ed113
update lid driven cavity results
gouarin May 4, 2026
e91024b
fix pre-commit
gouarin May 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 37 additions & 8 deletions include/samurai/algorithm/update.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -361,11 +361,20 @@ namespace samurai

auto& mesh = field.mesh();

std::size_t nnz = 0;
for (std::size_t d = 0; d < dim; ++d)
{
if (direction[d] != 0)
{
nnz++;
}
}

for (std::size_t delta_l = 1; delta_l <= 2; ++delta_l) // lower level (1 or 2)
{
auto proj_level = level - delta_l;

auto fine_inner_corner = self(mesh.corner(direction)).on(level);
auto fine_inner_corner = intersection(self(mesh.corner(direction)).on(level), mesh[mesh_id_t::cells][level]);
auto fine_outer_corner = intersection(translate(fine_inner_corner, direction), mesh[mesh_id_t::reference][level]);
auto projection_ghost = intersection(fine_outer_corner.on(proj_level), mesh[mesh_id_t::reference][proj_level]);

Expand All @@ -374,15 +383,25 @@ namespace samurai
{
using index_t = std::decay_t<decltype(index)>;

auto i_child = (1 << delta_l) * i;
i_child.start += direction[0] == -1 ? ((1 << delta_l) - 1) : 0;
i_child.end = i_child.start + 1;
i_child.step = 1;
index_t index_child = (1 << delta_l) * index;
auto i_child = (i << delta_l) + (direction[0] == -1 ? ((1 << delta_l) - 1) : 0); // this is the interval of the child
// cell in the fine level
if (nnz == dim) // || direction[0] != 0)
{
i_child.end = i_child.start + 1; // if we are projecting a corner ghost, we want only 1 child, so end = start + 1
i_child.step = 1;
}
Comment thread
gouarin marked this conversation as resolved.
else
{
i_child.step = 1 << delta_l;
}

index_t index_child = index << delta_l;
Comment thread
gouarin marked this conversation as resolved.

for (std::size_t d = 0; d < dim - 1; ++d)
{
index_child[d] += direction[d + 1] == -1 ? ((1 << delta_l) - 1) : 0;
}

field(proj_level, i, index) = field(level, i_child, index_child);
});
if (proj_level == 0)
Expand All @@ -409,8 +428,18 @@ namespace samurai
for_each_diagonal_direction<dim>(
[&](auto& direction)
{
auto d = find_direction_index(direction);
if (!mesh.is_periodic(d))
// Skip if any non-zero component of the direction is periodic:
// a periodic direction has no real boundary, so no corner ghost to fill.
bool any_periodic = false;
for (std::size_t d = 0; d < dim; ++d)
{
if (direction[d] != 0 && mesh.is_periodic(d))
{
any_periodic = true;
break;
}
}
if (!any_periodic)
{
update_outer_corners_by_polynomial_extrapolation(level, direction, field);
project_corner_below(level, direction, field); // project to level-1 and level-2
Expand Down
117 changes: 112 additions & 5 deletions include/samurai/bc/apply_field_bc.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -322,14 +322,121 @@ namespace samurai
return; // No outer corners in 1D
}

static constexpr std::size_t extrap_stencil_size = 2;
static constexpr std::size_t max_stencil_size_PE = PolynomialExtrapolation<Field, 2>::max_stencil_size_implemented_PE;

auto& domain = detail::get_mesh(field.mesh());
PolynomialExtrapolation<Field, extrap_stencil_size> bc(domain, ConstantBc<Field>(), true);
int ghost_width = field.mesh().ghost_width();
const auto& domain = detail::get_mesh(field.mesh());
const auto& corner_lca = field.mesh().corner(direction);

auto corner = self(field.mesh().corner(direction)).on(level);
// Step 1: Fill the diagonal ghost cells layer by layer using stencil sizes 2, 4, ..., 2*ghost_width
for (int ghost_layer = 1; ghost_layer <= ghost_width; ++ghost_layer)
{
int stencil_s = 2 * ghost_layer;
static_for<2, max_stencil_size_PE + 1>::apply(
[&](auto stencil_size_)
Comment thread
gouarin marked this conversation as resolved.
{
static constexpr int stencil_size = static_cast<int>(stencil_size_());
if constexpr (stencil_size % 2 == 0)
{
if (stencil_s == stencil_size)
{
PolynomialExtrapolation<Field, stencil_size> bc(domain, ConstantBc<Field>(), true);
auto corner = self(corner_lca).on(level);
__apply_extrapolation_bc__cells<stencil_size>(bc, level, field, direction, corner);
}
}
});
}

// Step 2: Fill off-diagonal ghost cells using Cartesian constant extrapolation
// For each non-zero component of the direction, apply extrapolation from the diagonal ghost cells
using mesh_id_t = typename Field::mesh_t::mesh_id_t;

// Count non-zero direction components
std::size_t num_nonzero = 0;
std::array<std::size_t, Field::dim> nonzero_dirs;
for (std::size_t d = 0; d < Field::dim; ++d)
{
if (direction[d] != 0)
{
nonzero_dirs[num_nonzero] = d;
++num_nonzero;
}
}

// Check if the corner cells exist at this level
// If they don't, skip the off-diagonal fill entirely
auto corner_at_level = self(corner_lca).on(level);

PolynomialExtrapolation<Field, 2> bc_e(domain, ConstantBc<Field>(), true);
auto stencil_0_e = bc_e.get_stencil(std::integral_constant<std::size_t, 2>());

for (std::size_t d = 0; d < Field::dim; ++d)
{
if (direction[d] == 0)
{
continue;
}

DirectionVector<Field::dim> e_dir;
e_dir.fill(0);
e_dir[d] = direction[d];

auto stencil_e = convert_for_direction(stencil_0_e, e_dir);
auto analyzer_e = make_stencil_analyzer(stencil_e);

auto& mesh = field.mesh();
auto has_outer_neighbour = translate(self(mesh[mesh_id_t::reference][level]).on(level), -e_dir);

// For each layer, extrapolate from diagonal/off-diagonal ghost cells in the e_dir direction
for (int g = 1; g <= ghost_width; ++g)
{
// Source: translate(corner, direction + (g-1)*e_dir)
auto source_set = translate(corner_at_level, direction + (g - 1) * e_dir);
auto valid_source = intersection(source_set, has_outer_neighbour).on(level);
__apply_bc_on_subset(bc_e, field, valid_source, analyzer_e, e_dir);
}

// For true 3D corners (3 non-zero components), apply secondary sweeps
// to fill cells at edge-pair positions like (N+1,N+1,N)
if (num_nonzero == 3)
{
// For each pair of the OTHER non-zero directions, apply sweep in e_dir
// from the cells filled by those directions
for (std::size_t i = 0; i < num_nonzero; ++i)
{
if (nonzero_dirs[i] == d)
{
continue; // skip the current direction
}
for (std::size_t j = i + 1; j < num_nonzero; ++j)
{
if (nonzero_dirs[j] == d)
{
continue; // skip the current direction
}

// Create unit direction vectors for the other two directions
DirectionVector<Field::dim> e_dir_i;
e_dir_i.fill(0);
e_dir_i[nonzero_dirs[i]] = direction[nonzero_dirs[i]];

__apply_extrapolation_bc__cells<extrap_stencil_size>(bc, level, field, direction, corner);
DirectionVector<Field::dim> e_dir_j;
e_dir_j.fill(0);
e_dir_j[nonzero_dirs[j]] = direction[nonzero_dirs[j]];

// Source: cells filled by sweeps in directions i and j
// = translate(corner, direction + e_dir_i) and translate(corner, direction + e_dir_j)
auto source_i = translate(self(corner_lca).on(level), direction + e_dir_i);
auto source_j = translate(self(corner_lca).on(level), direction + e_dir_j);
auto combined_source = union_(source_i, source_j);

auto valid_source = intersection(combined_source, has_outer_neighbour).on(level);
__apply_bc_on_subset(bc_e, field, valid_source, analyzer_e, e_dir);
}
}
}
}
}

template <class Field>
Expand Down
109 changes: 77 additions & 32 deletions include/samurai/mesh.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -946,41 +946,86 @@ namespace samurai
template <class D, class Config>
SAMURAI_INLINE void Mesh_base<D, Config>::construct_corners()
{
using direction_t = DirectionVector<dim>;

// Generic implementation for any dim >= 2.
// For a diagonal direction d (each component ±1), the "corner" subset is:
// domain \ union_{i=0}^{dim-1}( translate(domain, e_i * -d[i]) )
// where e_i is the canonical unit vector along axis i.
// We build the union_ call generically using an index sequence.
m_corners.clear();
for_each_diagonal_direction<dim>(
[&](const auto& direction)
{
// Build a tuple of dim translates, one per axis
auto translates = [&]<std::size_t... Is>(std::index_sequence<Is...>)
static_assert(
dim <= 6,
"Corner construction is currently only implemented up to 6D due to the combinatorial number of cases. Please implement the general ND version if you need higher dimensions.");
if constexpr (dim > 1)
Comment thread
gouarin marked this conversation as resolved.
{
using direction_t = DirectionVector<dim>;

// For a diagonal direction d, the "corner" (inner cells at the boundary corner/edge) is:
// domain \ union_{i: d[i] != 0}( translate(domain, e_i * -d[i]) )
// Only non-zero components are included: a zero component would give
// translate(domain, 0) = domain, making the union contain the domain itself
// and the difference empty. This is the correct formula for both true corner
// directions (all components ±1) and edge directions in 3D (one zero component).
m_corners.clear();
for_each_diagonal_direction<dim>(
[&](const auto& direction)
{
// For axis I: shift is -direction[I] along axis I, 0 elsewhere
auto make_translate = [&]<std::size_t I>(std::integral_constant<std::size_t, I>)
// Collect shifts only for non-zero components
std::array<direction_t, dim> nonzero_shifts;
std::size_t n = 0;
for (std::size_t I = 0; I < dim; ++I)
{
direction_t shift{};
shift.fill(0);
shift[I] = -direction[I];
return translate(m_domain, shift);
};
return std::make_tuple(make_translate(std::integral_constant<std::size_t, Is>{})...);
}(std::make_index_sequence<dim>{});

// Apply union_ to the tuple of translates
auto u = std::apply(
[](const auto&... ts)
if (direction[I] != 0)
{
nonzero_shifts[n].fill(0);
nonzero_shifts[n][I] = -direction[I];
++n;
}
}

// Apply union_ to non-zero-component translates only.
// for_each_diagonal_direction guarantees n >= 2.
// In 3D: n == 2 for edge directions, n == 3 for true corner directions.
// TODO: make a ND version
if (n == 2)
{
return union_(ts...);
},
translates);

m_corners.push_back(difference(m_domain, u).to_lca());
});
m_corners.push_back(
difference(m_domain, union_(translate(m_domain, nonzero_shifts[0]), translate(m_domain, nonzero_shifts[1])))
.to_lca());
}
else if (n == 3)
{
m_corners.push_back(difference(m_domain,
union_(translate(m_domain, nonzero_shifts[0]),
translate(m_domain, nonzero_shifts[1]),
translate(m_domain, nonzero_shifts[2])))
.to_lca());
}
else if (n == 4)
{
m_corners.push_back(difference(m_domain,
union_(translate(m_domain, nonzero_shifts[0]),
translate(m_domain, nonzero_shifts[1]),
translate(m_domain, nonzero_shifts[2]),
translate(m_domain, nonzero_shifts[3])))
.to_lca());
}
else if (n == 5)
{
m_corners.push_back(difference(m_domain,
union_(translate(m_domain, nonzero_shifts[0]),
translate(m_domain, nonzero_shifts[1]),
translate(m_domain, nonzero_shifts[2]),
translate(m_domain, nonzero_shifts[3]),
translate(m_domain, nonzero_shifts[4])))
.to_lca());
}
else if (n == 6)
{
m_corners.push_back(difference(m_domain,
union_(translate(m_domain, nonzero_shifts[0]),
translate(m_domain, nonzero_shifts[1]),
translate(m_domain, nonzero_shifts[2]),
translate(m_domain, nonzero_shifts[3]),
translate(m_domain, nonzero_shifts[4]),
translate(m_domain, nonzero_shifts[5])))
.to_lca());
}
});
}
}

template <class D, class Config>
Expand Down
Binary file modified tests/reference/finite_volume/test_finite_volume_demo_lid_driven_cavity.h5
View file
Open in desktop
Binary file not shown.
Loading