Skip to content

Commit

Permalink
CarpetX: Correct handling symmetry boundaries
Browse files Browse the repository at this point in the history
  • Loading branch information
eschnett committed Jul 27, 2023
1 parent 3f93f42 commit 69030c7
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 61 deletions.
3 changes: 0 additions & 3 deletions CarpetX/src/boundaries_impl.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,6 @@ void BoundaryCondition::apply_on_face() const {
const boundary_t boundary_x = groupdata.boundaries[f][0];

if ((symmetry_x == symmetry_t::none && boundary_x == boundary_t::none) ||
(boundary_x == boundary_t::none) ||
(symmetry_x == symmetry_t::interpatch &&
boundary_x == boundary_t::none) ||
symmetry_x == symmetry_t::periodic)
Expand Down Expand Up @@ -113,7 +112,6 @@ void BoundaryCondition::apply_on_face_symbcx(
const boundary_t boundary_y = groupdata.boundaries[f][1];

if ((symmetry_y == symmetry_t::none && boundary_y == boundary_t::none) ||
(boundary_y == boundary_t::none) ||
(symmetry_y == symmetry_t::interpatch &&
boundary_y == boundary_t::none) ||
symmetry_y == symmetry_t::periodic)
Expand Down Expand Up @@ -158,7 +156,6 @@ void BoundaryCondition::apply_on_face_symbcxy(
const boundary_t boundary_z = groupdata.boundaries[f][2];

if ((symmetry_z == symmetry_t::none && boundary_z == boundary_t::none) ||
(boundary_z == boundary_t::none) ||
(symmetry_z == symmetry_t::interpatch &&
boundary_z == boundary_t::none) ||
symmetry_z == symmetry_t::periodic)
Expand Down
95 changes: 44 additions & 51 deletions CarpetX/src/driver.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -114,21 +114,23 @@ std::ostream &operator<<(std::ostream &os, const boundary_t boundary) {
return os;
}

array<array<symmetry_t, dim>, 2> get_symmetries() {
array<array<symmetry_t, dim>, 2> get_symmetries(const int patch) {
// patch < 0 return symmetries without taking interpatch boundaries into
// account
DECLARE_CCTK_PARAMETERS;

const array<array<bool, 3>, 2> is_outer_boundary{{
{{
bool(!CCTK_EQUALS(boundary_x, "none")),
bool(!CCTK_EQUALS(boundary_y, "none")),
bool(!CCTK_EQUALS(boundary_z, "none")),
}},
{{
bool(!CCTK_EQUALS(boundary_upper_x, "none")),
bool(!CCTK_EQUALS(boundary_upper_y, "none")),
bool(!CCTK_EQUALS(boundary_upper_z, "none")),
}},
}};
array<array<bool, 3>, 2> is_interpatch{
{{{false, false, false}}, {{false, false, false}}}};
if (patch >= 0 &&
CCTK_IsFunctionAliased("MultiPatch_GetBoundarySpecification2")) {
CCTK_INT is_interpatch_boundary[2 * dim];
const int ierr = MultiPatch_GetBoundarySpecification2(
patch, 2 * dim, is_interpatch_boundary);
assert(!ierr);
for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
is_interpatch[f][d] = is_interpatch_boundary[2 * d + f];
}
const array<array<bool, 3>, 2> is_periodic{{
{{bool(periodic_x), bool(periodic_y), bool(periodic_z)}},
{{bool(periodic_x), bool(periodic_y), bool(periodic_z)}},
Expand All @@ -138,33 +140,31 @@ array<array<symmetry_t, dim>, 2> get_symmetries() {
{{bool(reflection_upper_x), bool(reflection_upper_y),
bool(reflection_upper_z)}},
}};

for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
assert(is_outer_boundary[f][d] + is_periodic[f][d] +
is_reflection[f][d] <=
assert(is_interpatch[f][d] + is_periodic[f][d] + is_reflection[f][d] <=
1);

array<array<symmetry_t, dim>, 2> symmetries;
for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
symmetries[f][d] = is_outer_boundary[f][d] ? symmetry_t::none
: is_periodic[f][d] ? symmetry_t::periodic
: is_reflection[f][d] ? symmetry_t::reflection
: symmetry_t::none;
symmetries[f][d] = is_interpatch[f][d] ? symmetry_t::interpatch
: is_periodic[f][d] ? symmetry_t::periodic
: is_reflection[f][d] ? symmetry_t::reflection
: symmetry_t::none;

return symmetries;
}

array<array<boundary_t, dim>, 2> get_default_boundaries() {
DECLARE_CCTK_PARAMETERS;

const array<array<symmetry_t, dim>, 2> symmetries = get_symmetries();
const array<array<symmetry_t, dim>, 2> symmetries = get_symmetries(-1);
array<array<bool, 3>, 2> is_symmetry;
for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
is_symmetry[f][d] = symmetries[f][d] != symmetry_t::none
// && symmetries[f][d] != symmetry_t::outer_boundary
;
is_symmetry[f][d] = symmetries[f][d] != symmetry_t::none;

const array<array<bool, 3>, 2> is_dirichlet{{
{{
Expand Down Expand Up @@ -238,7 +238,7 @@ array<array<boundary_t, dim>, 2> get_default_boundaries() {
array<array<boundary_t, dim>, 2> get_group_boundaries(const int gi) {
DECLARE_CCTK_PARAMETERS;

const array<array<symmetry_t, dim>, 2> symmetries = get_symmetries();
const array<array<symmetry_t, dim>, 2> symmetries = get_symmetries(-1);
array<array<bool, 3>, 2> is_symmetry;
for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
Expand Down Expand Up @@ -1187,20 +1187,7 @@ GHExt::PatchData::PatchData(const int patch) : patch(patch) {
const amrex::Vector<amrex::IntVect> reffacts{}; // empty

// Symmetries
if (CCTK_IsFunctionAliased("MultiPatch_GetBoundarySpecification2")) {
CCTK_INT is_interpatch_boundary[2 * dim];
const int ierr = MultiPatch_GetBoundarySpecification2(
patch, 2 * dim, is_interpatch_boundary);
assert(!ierr);
// TODO: Set this in get_symmetries()
for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
symmetries[f][d] = is_interpatch_boundary[2 * d + f]
? symmetry_t::interpatch
: symmetry_t::none;
} else {
symmetries = get_symmetries();
}
symmetries = get_symmetries(patch);

{
const std::array faces{"lower", "upper"};
Expand Down Expand Up @@ -1262,14 +1249,6 @@ GHExt::PatchData::PatchData(const int patch) : patch(patch) {
}
}

bool GHExt::PatchData::all_faces_have_symmetries() const {
bool res = true;
for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
res &= symmetries.at(f).at(d) != symmetry_t::none;
return res;
}

GHExt::PatchData::LevelData::LevelData(const int patch, const int level,
const amrex::BoxArray &ba,
const amrex::DistributionMapping &dm,
Expand Down Expand Up @@ -1487,6 +1466,18 @@ void GHExt::PatchData::LevelData::GroupData::free_tmp_mfabs() const {

////////////////////////////////////////////////////////////////////////////////

bool GHExt::PatchData::LevelData::GroupData::
all_faces_have_symmetries_or_boundaries() const {
const auto &symmetries = ghext->patchdata.at(patch).symmetries;
bool res = true;
for (int f = 0; f < 2; ++f)
for (int d = 0; d < dim; ++d)
res &= symmetries.at(f).at(d) != symmetry_t::none ||
(symmetries.at(f).at(d) == symmetry_t::none &&
boundaries.at(f).at(d) != boundary_t::none);
return res;
}

void GHExt::PatchData::LevelData::GroupData::apply_boundary_conditions(
amrex::MultiFab &mfab) const {
DECLARE_CCTK_PARAMETERS;
Expand Down Expand Up @@ -1834,9 +1825,10 @@ void CactusAmrCore::MakeNewLevelFromCoarse(
groupdata, *groupdata.mfab.at(tl), *coarsegroupdata.mfab.at(tl),
patchdata.amrcore->Geom(level - 1), patchdata.amrcore->Geom(level),
interpolator, groupdata.bcrecs);
const auto outer_valid = patchdata.all_faces_have_symmetries()
? make_valid_outer()
: valid_t();
const auto outer_valid =
groupdata.all_faces_have_symmetries_or_boundaries()
? make_valid_outer()
: valid_t();
assert(outer_valid == make_valid_outer());
for (int vi = 0; vi < groupdata.numvars; ++vi) {
groupdata.valid.at(tl).at(vi).set_all(
Expand Down Expand Up @@ -1956,8 +1948,9 @@ void CactusAmrCore::RemakeLevel(const int level, const amrex::Real time,
amrex::Interpolater *const interpolator =
get_interpolator(groupdata.indextype);

const auto outer_valid =
patchdata.all_faces_have_symmetries() ? make_valid_outer() : valid_t();
const auto outer_valid = groupdata.all_faces_have_symmetries_or_boundaries()
? make_valid_outer()
: valid_t();
assert(outer_valid == make_valid_outer());

const nan_handling_t nan_handling = groupdata.do_checkpoint
Expand Down
2 changes: 1 addition & 1 deletion CarpetX/src/driver.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,6 @@ struct GHExt {
int patch;

array<array<symmetry_t, dim>, 2> symmetries;
bool all_faces_have_symmetries() const;

// AMReX grid structure
// TODO: convert this from unique_ptr to optional
Expand Down Expand Up @@ -264,6 +263,7 @@ struct GHExt {
array<int, dim> nghostzones;

array<array<boundary_t, dim>, 2> boundaries;
bool all_faces_have_symmetries_or_boundaries() const;
vector<array<int, dim> > parities;
vector<CCTK_REAL> dirichlet_values;
vector<CCTK_REAL> robin_values;
Expand Down
10 changes: 4 additions & 6 deletions CarpetX/src/schedule.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -2227,7 +2227,7 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups,
std::vector<std::function<void()> > tasks1;
std::vector<std::function<void()> > tasks2;

const bool have_multipatch_boundaries =
static const bool have_multipatch_boundaries =
CCTK_IsFunctionAliased("MultiPatch_Interpolate");

active_levels->loop([&](auto &restrict leveldata) {
Expand Down Expand Up @@ -2285,8 +2285,7 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups,
return "SyncGroupsByDirI after syncing: "
"Mark ghost zones as valid";
});
if (ghext->patchdata.at(leveldata.patch)
.all_faces_have_symmetries())
if (groupdata.all_faces_have_symmetries_or_boundaries())
groupdata.valid.at(tl).at(vi).set_outer(true, []() {
return "SyncGroupsByDirI after syncing: "
"Mark outer boundaries as valid";
Expand Down Expand Up @@ -2359,8 +2358,7 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups,
return "SyncGroupsByDirI after prolongation: "
"Mark ghost zones as valid";
});
if (ghext->patchdata.at(leveldata.patch)
.all_faces_have_symmetries())
if (groupdata.all_faces_have_symmetries_or_boundaries())
groupdata.valid.at(tl).at(vi).set_outer(true, []() {
return "SyncGroupsByDirI after prolongation: "
"Mark outer boundaries as valid";
Expand Down Expand Up @@ -2391,7 +2389,7 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups,
tasks2.clear();
synchronize();

if (CCTK_IsFunctionAliased("MultiPatch_Interpolate")) {
if (have_multipatch_boundaries) {
std::vector<CCTK_INT> cactusvarinds;
for (int group : groups) {
const auto &groupdata =
Expand Down

0 comments on commit 69030c7

Please sign in to comment.