Skip to content
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

Jump to bottom

feat: Portion for a slice #136

Merged
merged 6 commits into from
Oct 18, 2023
Merged

feat: Portion for a slice #136

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
78522b5
add hash for tuple and interval
gouarin Oct 18, 2023
c7d15bd
add a flag to keep previous solution or not in MRA
gouarin Oct 18, 2023
ba2bacb
fix criteria for 3d
gouarin Oct 18, 2023
dae7c4a
add portion for a slice defined by intervals
gouarin Oct 18, 2023
9101844
fix adapt
gouarin Oct 18, 2023
94fed48
fix
gouarin Oct 18, 2023
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
85 changes: 84 additions & 1 deletion include/samurai/algorithm/update.hpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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