Skip to content

feat: ennforce contiguous cells at the boundaries #306

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
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

feat: ennforce contiguous cells at the boundaries #306

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
126 changes: 110 additions & 16 deletions include/samurai/algorithm/graduation.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -172,12 +172,13 @@ namespace samurai
}

template <size_t dim, typename TInterval, typename MeshType, size_t max_size, typename TCoord>
void list_intervals_to_remove(const size_t grad_width,
const CellArray<dim, TInterval, max_size>& ca,
[[maybe_unused]] const std::vector<MPI_Subdomain<MeshType>>& mpi_neighbourhood,
const std::array<bool, dim>& is_periodic,
const std::array<int, dim>& nb_cells_finest_level,
std::array<ArrayOfIntervalAndPoint<TInterval, TCoord>, CellArray<dim, TInterval, max_size>::max_size>& out)
void list_interval_to_refine_for_graduation(
const size_t grad_width,
const CellArray<dim, TInterval, max_size>& ca,
[[maybe_unused]] const std::vector<MPI_Subdomain<MeshType>>& mpi_neighbourhood,
const std::array<bool, dim>& is_periodic,
const std::array<int, dim>& nb_cells_finest_level,
std::array<ArrayOfIntervalAndPoint<TInterval, TCoord>, CellArray<dim, TInterval, max_size>::max_size>& out)
{
const size_t max_level = ca.max_level();
const size_t min_level = ca.min_level();
Expand Down Expand Up @@ -218,7 +219,7 @@ namespace samurai
}
}
}
xt::xtensor_fixed<int, xt::xshape<dim>> translation = xt::xscalar(0);
DirectionVector<dim> translation = xt::xscalar(0);
for (size_t d = 0; d != dim; ++d)
{
if (is_periodic[d])
Expand Down Expand Up @@ -260,6 +261,87 @@ namespace samurai
#endif // SAMURAI_WITH_MPI
}

template <size_t dim, typename TInterval, size_t max_size, typename TCoord>
void list_interval_to_refine_for_contiguous_boundary_cells(
const size_t half_stencil_width,
const CellArray<dim, TInterval, max_size>& ca,
const LevelCellArray<dim, TInterval>& domain,
const std::array<bool, dim>& is_periodic,
std::array<ArrayOfIntervalAndPoint<TInterval, TCoord>, CellArray<dim, TInterval, max_size>::max_size>& out)
{
const size_t max_level = ca.max_level();
const size_t min_level = ca.min_level();

// We want to avoid that a flux is computed with ghosts outside of the domain if the cell doesn't touch the boundary,
// because we only want to apply the B.C. on the cells that touch the boundary.
// 1. Case where the boundary is at level L and the jump is going down to L-1:
// We want to have enough contiguous boundary cells to ensure that the stencil at the lower level
// won't go outside the domain.
// To ensure half_stencil_width at L-1, we need 2*half_stencil_width at level L.
// However, since we project the B.C. in the first outside ghost at level L-1, we can reduce the number of contiguous
// cells by 1 at level L-1. This makes, at level L, 2*(half_stencil_width - 1) contiguous cells.
const int n_contiguous_boundary_cells = std::max(int(half_stencil_width), 2 * (int(half_stencil_width) - 2));

// 2. Case where the boundary is at level L and jump is going up:
// Then ensuring half_stencil_width contiguous cells at level L automatically ensures half_stencil_width
// at level L+1.

for_each_cartesian_direction<dim>(
[&](const auto direction_idx, const auto& translation)
{
if (not is_periodic[direction_idx])
{
// Jump level --> level-1
for (size_t level = max_level; level != min_level; --level)
{
auto boundaryCells = difference(ca[level], translate(self(domain).on(level), -translation)).on(level);

for (int i = 2; i <= n_contiguous_boundary_cells; i += 2)
{
// auto refine_subset = intersection(translate(boundaryCells, -i*translation), ca[level-1]).on(level-1);
LevelCellArray<dim, TInterval> translated_boundary(translate(boundaryCells, -i * translation));
auto refine_subset = intersection(translated_boundary, ca[level - 1]).on(level - 1);
refine_subset(
[&](const auto& x_interval, const auto& yz)
{
out[level - 1].push_back(x_interval, yz);
});
}
}
// Jump level --> level+1
for (size_t level = max_level - 1; level != min_level - 1; --level)
{
auto boundaryCells = difference(ca[level], translate(self(domain).on(level), -translation));
for (size_t i = 1; i != half_stencil_width; ++i)
{
auto refine_subset = translate(intersection(translate(boundaryCells, -i * translation), ca[level + 1]).on(level),
i * translation)
.on(level);
refine_subset(
[&](const auto& x_interval, const auto& yz)
{
out[level].push_back(x_interval, yz);
});
}
}
}
});
}

template <size_t dim, typename TInterval, typename MeshType, size_t max_size, typename TCoord>
void list_intervals_to_refine(const size_t grad_width,
const size_t half_stencil_width,
const CellArray<dim, TInterval, max_size>& ca,
const LevelCellArray<dim, TInterval>& domain,
[[maybe_unused]] const std::vector<MPI_Subdomain<MeshType>>& mpi_neighbourhood,
const std::array<bool, dim>& is_periodic,
const std::array<int, dim>& nb_cells_finest_level,
std::array<ArrayOfIntervalAndPoint<TInterval, TCoord>, CellArray<dim, TInterval, max_size>::max_size>& out)
{
list_interval_to_refine_for_graduation(grad_width, ca, mpi_neighbourhood, is_periodic, nb_cells_finest_level, out);
list_interval_to_refine_for_contiguous_boundary_cells(half_stencil_width, ca, domain, is_periodic, out);
}

// if add the intervals in add_m_interval
// if dim = 2 then add_m_interval stores the y coord
// if dim > 2 then add_intercal contains the 'inner_stencil' i.e. the coordinates y+s_x, z+s_z, etc.
Expand Down Expand Up @@ -317,17 +399,29 @@ namespace samurai

template <std::size_t dim, class TInterval, class MeshType, size_t max_size>
size_t make_graduation(CellArray<dim, TInterval, max_size>& ca,
const LevelCellArray<dim, TInterval>& domain,
[[maybe_unused]] const std::vector<MPI_Subdomain<MeshType>>& mpi_neighbourhood,
const std::array<bool, dim>& is_periodic,
const std::array<int, dim>& nb_cells_finest_level,
const size_t grad_width = 1)
const size_t grad_width = 1,
const size_t half_stencil_width = 1 // helf of width of the numerical scheme's stencil.
)
{
using ca_type = CellArray<dim, TInterval, max_size>;
using coord_type = typename ca_type::lca_type::coord_type;

const size_t max_level = ca.max_level();
const size_t min_level = ca.min_level();

const auto& min_indices = domain.min_indices();
const auto& max_indices = domain.max_indices();

std::array<int, dim> nb_cells_finest_level;

for (size_t d = 0; d != max_indices.size(); ++d)
{
nb_cells_finest_level[d] = max_indices[d] - min_indices[d];
}

std::vector<TInterval> add_p_interval;
std::vector<coord_type> add_p_inner_stencil;
std::vector<size_t> add_p_idx;
Expand All @@ -351,18 +445,18 @@ namespace samurai
// Then, if the non-graduated is not taged as keep, we coarsen it
ca_add_p.clear();
ca_remove_p.clear();
list_intervals_to_remove(grad_width, ca, mpi_neighbourhood, is_periodic, nb_cells_finest_level, remove_m_all);
list_intervals_to_refine(grad_width, half_stencil_width, ca, domain, mpi_neighbourhood, is_periodic, nb_cells_finest_level, remove_m_all);

add_p_interval.clear();
add_p_inner_stencil.clear();
add_p_idx.clear();
for (size_t level = min_level; level != max_level + 1; ++level)
{
#ifdef SAMURAI_WITH_MPI
// #ifdef SAMURAI_WITH_MPI
remove_m_all[level].remove_overlapping_intervals();
#else
remove_m_all[level].sort_intervals();
#endif // SAMURAI_WITH_MPI
// #else
// remove_m_all[level].sort_intervals();
// #endif // SAMURAI_WITH_MPI
const size_t imax = remove_m_all[level].size();
for (size_t i = 0; i != imax; ++i)
{
Expand Down Expand Up @@ -424,10 +518,10 @@ namespace samurai

std::vector<MPI_Subdomain<DummyMesh>> mpi_neighbourhood;
std::array<bool, dim> is_periodic;
std::array<int, dim> nb_cells_finest_level;
LevelCellArray<dim, TInterval> domain;

is_periodic.fill(false);
return make_graduation(ca, mpi_neighbourhood, is_periodic, nb_cells_finest_level, grad_width);
return make_graduation(ca, domain, mpi_neighbourhood, is_periodic, grad_width);
}

template <std::size_t dim, class TInterval, size_t max_size, class Tag>
Expand Down
17 changes: 7 additions & 10 deletions include/samurai/mr/adapt.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -311,18 +311,15 @@ namespace samurai
// on test adapt_test/2.mutliple_fields with:
// linux-mamba (clang-18, ubuntu-24.04, clang, clang-18, clang-18, clang++-18)
// while the code bellow do not.
const auto& min_indices = mesh.domain().min_indices();
const auto& max_indices = mesh.domain().max_indices();

std::array<int, mesh_t::dim> nb_cells_finest_level;

for (size_t d = 0; d != max_indices.size(); ++d)
{
nb_cells_finest_level[d] = max_indices[d] - min_indices[d];
}

ca_type new_ca = update_cell_array_from_tag(mesh[mesh_id_t::cells], m_tag);
make_graduation(new_ca, mesh.mpi_neighbourhood(), mesh.periodicity(), nb_cells_finest_level, mesh_t::config::graduation_width);
make_graduation(new_ca,
mesh.domain(),
mesh.mpi_neighbourhood(),
mesh.periodicity(),
mesh_t::config::graduation_width,
mesh_t::config::max_stencil_width);

mesh_t new_mesh{new_ca, mesh};
#ifdef SAMURAI_WITH_MPI
mpi::communicator world;
Expand Down
15 changes: 15 additions & 0 deletions include/samurai/stencil.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,21 @@ namespace samurai
return 0;
}

template <std::size_t dim, class Func>
void for_each_cartesian_direction(Func&& f)
{
DirectionVector<dim> direction;
direction.fill(0);
for (std::size_t d = 0; d < dim; ++d)
{
direction[d] = 1;
f(d, direction);
direction[d] = -1;
f(d, direction);
direction[d] = 0;
}
}

//-----------------------//
// Useful stencils //
//-----------------------//
Expand Down
Loading