diff --git a/include/samurai/algorithm/update.hpp b/include/samurai/algorithm/update.hpp index 43b4c8268..4c0b6888f 100644 --- a/include/samurai/algorithm/update.hpp +++ b/include/samurai/algorithm/update.hpp @@ -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]); @@ -374,15 +383,25 @@ namespace samurai { using index_t = std::decay_t; - 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; + } + else + { + i_child.step = 1 << delta_l; + } + + index_t index_child = index << delta_l; + 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) @@ -409,8 +428,18 @@ namespace samurai for_each_diagonal_direction( [&](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 diff --git a/include/samurai/bc/apply_field_bc.hpp b/include/samurai/bc/apply_field_bc.hpp index 9b28137c6..439ec62b3 100644 --- a/include/samurai/bc/apply_field_bc.hpp +++ b/include/samurai/bc/apply_field_bc.hpp @@ -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::max_stencil_size_implemented_PE; - auto& domain = detail::get_mesh(field.mesh()); - PolynomialExtrapolation bc(domain, ConstantBc(), 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_) + { + static constexpr int stencil_size = static_cast(stencil_size_()); + if constexpr (stencil_size % 2 == 0) + { + if (stencil_s == stencil_size) + { + PolynomialExtrapolation bc(domain, ConstantBc(), true); + auto corner = self(corner_lca).on(level); + __apply_extrapolation_bc__cells(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 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 bc_e(domain, ConstantBc(), true); + auto stencil_0_e = bc_e.get_stencil(std::integral_constant()); + + for (std::size_t d = 0; d < Field::dim; ++d) + { + if (direction[d] == 0) + { + continue; + } + + DirectionVector 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 e_dir_i; + e_dir_i.fill(0); + e_dir_i[nonzero_dirs[i]] = direction[nonzero_dirs[i]]; - __apply_extrapolation_bc__cells(bc, level, field, direction, corner); + DirectionVector 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 diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp index 77f1f91ca..213c68d97 100644 --- a/include/samurai/mesh.hpp +++ b/include/samurai/mesh.hpp @@ -946,41 +946,86 @@ namespace samurai template SAMURAI_INLINE void Mesh_base::construct_corners() { - using direction_t = DirectionVector; - - // 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( - [&](const auto& direction) - { - // Build a tuple of dim translates, one per axis - auto translates = [&](std::index_sequence) + 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) + { + using direction_t = DirectionVector; + + // 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( + [&](const auto& direction) { - // For axis I: shift is -direction[I] along axis I, 0 elsewhere - auto make_translate = [&](std::integral_constant) + // Collect shifts only for non-zero components + std::array 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::make_index_sequence{}); - - // 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 diff --git a/tests/reference/finite_volume/test_finite_volume_demo_lid_driven_cavity.h5 b/tests/reference/finite_volume/test_finite_volume_demo_lid_driven_cavity.h5 index 73e467613..9164103cc 100644 Binary files a/tests/reference/finite_volume/test_finite_volume_demo_lid_driven_cavity.h5 and b/tests/reference/finite_volume/test_finite_volume_demo_lid_driven_cavity.h5 differ