Skip to content

Commit

Permalink
feat: Portion for a slice (#136)
Browse files Browse the repository at this point in the history
  • Loading branch information
@gouarin
gouarin authored Oct 18, 2023
1 parent 36b1244 commit 19e19ab
Show file tree
Hide file tree
Showing 5 changed files with 643 additions and 118 deletions.
85 changes: 84 additions & 1 deletion include/samurai/algorithm/update.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -447,6 +447,18 @@ namespace samurai
update_fields_with_old(new_mesh, old_fields, fields.elements(), std::make_index_sequence<sizeof...(T)>{});
}

template <class Mesh, class Fields, std::size_t... Is>
void update_fields(Mesh& new_mesh, Fields& fields, std::index_sequence<Is...>)
{
(update_fields(new_mesh, std::get<Is>(fields)), ...);
}

template <class Mesh, class... T>
void update_fields(Mesh& new_mesh, Field_tuple<T...>& fields)
{
update_fields(new_mesh, fields.elements(), std::make_index_sequence<sizeof...(T)>{});
}

template <class Mesh, class Field, class... Fields>
void update_fields(Mesh& new_mesh, Field& field, Fields&... fields)
{
Expand Down Expand Up @@ -540,7 +552,7 @@ namespace samurai
}

template <class Tag, class Field, class Old_field, class... Fields>
bool update_field_mr(const Tag& tag, Field& field, Old_field& old_field, Fields&... other_fields)
bool update_field_mr(std::integral_constant<bool, true>, const Tag& tag, Field& field, Old_field& old_field, Fields&... other_fields)
{
using mesh_t = typename Field::mesh_t;
static constexpr std::size_t dim = mesh_t::dim;
Expand Down Expand Up @@ -609,4 +621,75 @@ namespace samurai

return false;
}

template <class Tag, class Field, class... Fields>
bool update_field_mr(std::integral_constant<bool, false>, const Tag& tag, Field& field, Fields&... other_fields)
{
using mesh_t = typename Field::mesh_t;
static constexpr std::size_t dim = mesh_t::dim;
using mesh_id_t = typename Field::mesh_t::mesh_id_t;
using interval_t = typename mesh_t::interval_t;
using value_t = typename interval_t::value_t;
using cl_type = typename Field::mesh_t::cl_type;

auto& mesh = field.mesh();
cl_type cl;

for_each_interval(mesh[mesh_id_t::cells],
[&](std::size_t level, const auto& interval, const auto& index)
{
auto itag = interval.start + interval.index;
for (value_t i = interval.start; i < interval.end; ++i)
{
if (tag[itag] & static_cast<int>(CellFlag::refine))
{
if (level < mesh.max_level())
{
static_nested_loop<dim - 1, 0, 2>(
[&](const auto& stencil)
{
auto new_index = 2 * index + stencil;
cl[level + 1][new_index].add_interval({2 * i, 2 * i + 2});
});
}
else
{
cl[level][index].add_point(i);
}
}
else if (tag[itag] & static_cast<int>(CellFlag::keep))
{
cl[level][index].add_point(i);
}
else
{
if (level > mesh.min_level())
{
cl[level - 1][index >> 1].add_point(i >> 1);
}
else
{
cl[level][index].add_point(i);
}
}
itag++;
}
});

mesh_t new_mesh = {cl, mesh.min_level(), mesh.max_level(), mesh.periodicity()};

if (mesh == new_mesh)
{
return true;
}

detail::update_fields(new_mesh, other_fields...);
detail::update_fields(new_mesh, field);

// detail::swap_mesh(new_mesh, old_field, field);
field.mesh().swap(new_mesh);
// std::get<0>(old_field).mesh().swap(new_mesh);

return false;
}
}
51 changes: 37 additions & 14 deletions include/samurai/mr/adapt.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ namespace samurai
}
}

template <class TField, class... TFields>
template <bool keep_previous, class TField, class... TFields>
class Adapt
{
public:
Expand Down Expand Up @@ -152,17 +152,17 @@ namespace samurai
tag_t m_tag;
};

template <class TField, class... TFields>
inline Adapt<TField, TFields...>::Adapt(TField& field, TFields&... fields)
template <bool keep_previous, class TField, class... TFields>
inline Adapt<keep_previous, TField, TFields...>::Adapt(TField& field, TFields&... fields)
: m_fields(field, fields...)
, m_detail("detail", field.mesh())
, m_tag("tag", field.mesh())
{
}

template <class TField, class... TFields>
template <bool keep_previous, class TField, class... TFields>
template <class... Fields>
void Adapt<TField, TFields...>::operator()(double eps, double regularity, Fields&... other_fields)
void Adapt<keep_previous, TField, TFields...>::operator()(double eps, double regularity, Fields&... other_fields)
{
auto& mesh = m_fields.mesh();
std::size_t min_level = mesh.min_level();
Expand All @@ -174,8 +174,12 @@ namespace samurai
}
update_ghost_mr(m_fields);

auto mesh_old = mesh;
old_fields_t old_fields = detail::copy_fields(mesh_old, m_fields);
auto mesh_old = mesh;
old_fields_t old_fields;
if constexpr (keep_previous)
{
old_fields = detail::copy_fields(mesh_old, m_fields);
}

for (std::size_t i = 0; i < max_level - min_level; ++i)
{
Expand Down Expand Up @@ -230,9 +234,13 @@ namespace samurai
}
}

template <class TField, class... TFields>
template <bool keep_previous, class TField, class... TFields>
template <class... Fields>
bool Adapt<TField, TFields...>::harten(std::size_t ite, double eps, double regularity, old_fields_t& old_fields, Fields&... other_fields)
bool Adapt<keep_previous, TField, TFields...>::harten(std::size_t ite,
double eps,
double regularity,
old_fields_t& old_fields,
Fields&... other_fields)
{
auto& mesh = m_fields.mesh();

Expand Down Expand Up @@ -317,8 +325,9 @@ namespace samurai

for (std::size_t is = 0; is < stencil.shape(0); ++is)
{
auto s = xt::view(stencil, is);
auto subset = intersection(translate(mesh[mesh_id_t::cells][level], s), mesh[mesh_id_t::all_cells][level - 1]).on(level);
auto s = xt::view(stencil, is);
auto subset = intersection(translate(mesh[mesh_id_t::cells][level], s), mesh[mesh_id_t::all_cells][level - 1], mesh.domain())
.on(level);

subset.apply_op(make_graduation(m_tag));
}
Expand Down Expand Up @@ -356,14 +365,28 @@ namespace samurai
update_tag_periodic(level, m_tag);
}

update_ghost_mr(old_fields);
update_ghost_mr(other_fields...);
return update_field_mr(m_tag, m_fields, old_fields, other_fields...);

if constexpr (keep_previous)
{
update_ghost_mr(old_fields);
return update_field_mr(std::integral_constant<bool, true>{}, m_tag, m_fields, old_fields, other_fields...);
}
else
{
return update_field_mr(std::integral_constant<bool, false>{}, m_tag, m_fields, other_fields...);
}
}

template <class... TFields>
auto make_MRAdapt(TFields&... fields)
{
return Adapt<TFields...>(fields...);
return Adapt<true, TFields...>(fields...);
}

template <bool keep_previous, class... TFields>
auto make_MRAdapt(TFields&... fields)
{
return Adapt<keep_previous, TFields...>(fields...);
}
} // namespace samurai
60 changes: 30 additions & 30 deletions include/samurai/mr/criteria.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,13 +126,13 @@ namespace samurai
// (xt::abs(detail(level, 2*i+1, 2*j+1))/maxd <
// eps);
auto mask = (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) < eps)
and (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) < eps);
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) < eps)
&& (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) < eps)
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) < eps)
&& (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) < eps)
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) < eps)
&& (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) < eps)
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) < eps);

xt::masked_view(tag(fine_level, 2 * i, 2 * j, 2 * k), mask) = static_cast<int>(CellFlag::coarsen);
xt::masked_view(tag(fine_level, 2 * i + 1, 2 * j, 2 * k), mask) = static_cast<int>(CellFlag::coarsen);
Expand All @@ -154,14 +154,14 @@ namespace samurai
// (xt::abs(detail(level, 2*i+1,
// 2*j+1))/maxd < eps), {1}) > (size-1);
auto mask = xt::sum((xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) < eps)
and (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) < eps)
and (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) < eps)
and (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) < eps),
{1})
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) < eps)
&& (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) < eps)
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) < eps)
&& (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) < eps)
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) < eps)
&& (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) < eps)
&& (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) < eps),
{detail.is_soa ? 0 : 1})
> (size - 1);

xt::masked_view(tag(fine_level, 2 * i, 2 * j, 2 * k), mask) = static_cast<int>(CellFlag::coarsen);
Expand Down Expand Up @@ -381,13 +381,13 @@ namespace samurai
// ((xt::abs(detail(fine_level, 2*i+1,
// 2*j+1))/maxd) > eps);
auto mask = (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) > eps)
or (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) > eps);
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) > eps)
|| (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) > eps)
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) > eps)
|| (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) > eps)
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) > eps)
|| (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) > eps)
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) > eps);

xt::masked_view(tag(fine_level, 2 * i, 2 * j, 2 * k), mask) = static_cast<int>(CellFlag::refine);
xt::masked_view(tag(fine_level, 2 * i + 1, 2 * j, 2 * k), mask) = static_cast<int>(CellFlag::refine);
Expand All @@ -409,14 +409,14 @@ namespace samurai
// ((xt::abs(detail(fine_level, 2*i+1,
// 2*j+1))/maxd) > eps), {1}) > 0;
auto mask = xt::sum((xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) > eps)
or (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) > eps)
or (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) > eps)
or (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) > eps),
{1})
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k)) > eps)
|| (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k)) > eps)
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k)) > eps)
|| (xt::abs(detail(fine_level, 2 * i, 2 * j, 2 * k + 1)) > eps)
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j, 2 * k + 1)) > eps)
|| (xt::abs(detail(fine_level, 2 * i, 2 * j + 1, 2 * k + 1)) > eps)
|| (xt::abs(detail(fine_level, 2 * i + 1, 2 * j + 1, 2 * k + 1)) > eps),
{detail.is_soa ? 0 : 1})
> 0;

xt::masked_view(tag(fine_level, 2 * i, 2 * j, 2 * k), mask) = static_cast<int>(CellFlag::refine);
Expand Down
Loading

0 comments on commit 19e19ab

Please sign in to comment.