diff --git a/Algo/test/test_roots/.empty-directory b/Algo/test/test_roots/.empty-directory new file mode 100644 index 000000000..e69de29bb diff --git a/CarpetX/src/boundaries.cxx b/CarpetX/src/boundaries.cxx index 1448a25f6..b5a94b9e2 100644 --- a/CarpetX/src/boundaries.cxx +++ b/CarpetX/src/boundaries.cxx @@ -30,9 +30,9 @@ BoundaryCondition::BoundaryCondition( const GHExt::PatchData::LevelData::GroupData &groupdata, const amrex::Box &box, amrex::FArrayBox &dest) : groupdata(groupdata), patchdata(ghext->patchdata.at(groupdata.patch)), - geom(patchdata.amrcore->Geom(groupdata.level)), - dest(dest), imin{geom.Domain().smallEnd(0), geom.Domain().smallEnd(1), - geom.Domain().smallEnd(2)}, + geom(patchdata.amrcore->Geom(groupdata.level)), dest(dest), + imin{geom.Domain().smallEnd(0), geom.Domain().smallEnd(1), + geom.Domain().smallEnd(2)}, imax{geom.Domain().bigEnd(0) + 1 + !groupdata.indextype[0], geom.Domain().bigEnd(1) + 1 + !groupdata.indextype[1], geom.Domain().bigEnd(2) + 1 + !groupdata.indextype[2]}, @@ -252,8 +252,6 @@ void BoundaryCondition::apply() const { const boundary_t boundary = groupdata.boundaries[face][dir]; if ((symmetry == symmetry_t::none && boundary == boundary_t::none) || - (symmetry == symmetry_t::outer_boundary && - boundary == boundary_t::none) || symmetry == symmetry_t::interpatch || symmetry == symmetry_t::periodic) { symmetries[dir] = symmetry_t::none; diff --git a/CarpetX/src/boundaries_impl.hxx b/CarpetX/src/boundaries_impl.hxx index 7d0d1fb12..6c0b47cf0 100644 --- a/CarpetX/src/boundaries_impl.hxx +++ b/CarpetX/src/boundaries_impl.hxx @@ -69,8 +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) || - (symmetry_x == symmetry_t::outer_boundary && - boundary_x == boundary_t::none) || (symmetry_x == symmetry_t::interpatch && boundary_x == boundary_t::none) || symmetry_x == symmetry_t::periodic) @@ -114,8 +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) || - (symmetry_y == symmetry_t::outer_boundary && - boundary_y == boundary_t::none) || (symmetry_y == symmetry_t::interpatch && boundary_y == boundary_t::none) || symmetry_y == symmetry_t::periodic) @@ -160,8 +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) || - (symmetry_z == symmetry_t::outer_boundary && - boundary_z == boundary_t::none) || (symmetry_z == symmetry_t::interpatch && boundary_z == boundary_t::none) || symmetry_z == symmetry_t::periodic) diff --git a/CarpetX/src/driver.cxx b/CarpetX/src/driver.cxx index cfea78b58..ec535a563 100644 --- a/CarpetX/src/driver.cxx +++ b/CarpetX/src/driver.cxx @@ -82,8 +82,6 @@ std::ostream &operator<<(std::ostream &os, const symmetry_t symmetry) { switch (symmetry) { case symmetry_t::none: return os << "none"; - case symmetry_t::outer_boundary: - return os << "outer_boundary"; case symmetry_t::interpatch: return os << "interpatch"; case symmetry_t::periodic: @@ -116,21 +114,23 @@ std::ostream &operator<<(std::ostream &os, const boundary_t boundary) { return os; } -array, 2> get_symmetries() { +array, 2> get_symmetries(const int patch) { + // patch < 0 return symmetries without taking interpatch boundaries into + // account DECLARE_CCTK_PARAMETERS; - const array, 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, 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, 2> is_periodic{{ {{bool(periodic_x), bool(periodic_y), bool(periodic_z)}}, {{bool(periodic_x), bool(periodic_y), bool(periodic_z)}}, @@ -140,19 +140,19 @@ array, 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, 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::outer_boundary - : is_periodic[f][d] ? symmetry_t::periodic - : is_reflection[f][d] ? symmetry_t::reflection - : symmetry_t::outer_boundary; + 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; } @@ -160,12 +160,11 @@ array, 2> get_symmetries() { array, 2> get_default_boundaries() { DECLARE_CCTK_PARAMETERS; - const array, 2> symmetries = get_symmetries(); + const array, 2> symmetries = get_symmetries(-1); array, 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, 2> is_dirichlet{{ {{ @@ -239,12 +238,11 @@ array, 2> get_default_boundaries() { array, 2> get_group_boundaries(const int gi) { DECLARE_CCTK_PARAMETERS; - const array, 2> symmetries = get_symmetries(); + const array, 2> symmetries = get_symmetries(-1); array, 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; array, 2> boundaries = get_default_boundaries(); @@ -1189,20 +1187,7 @@ GHExt::PatchData::PatchData(const int patch) : patch(patch) { const amrex::Vector 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"}; @@ -1264,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, @@ -1489,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; @@ -1836,15 +1825,17 @@ 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(); - assert(outer_valid == make_valid_outer()); + const auto outer_valid = + groupdata.all_faces_have_symmetries_or_boundaries() + ? make_valid_outer() + : valid_t(); for (int vi = 0; vi < groupdata.numvars; ++vi) { groupdata.valid.at(tl).at(vi).set_all( make_valid_int() | make_valid_ghosts() | outer_valid, []() { return "MakeNewLevelFromCoarse after prolongation"; }); - // check_valid(groupdata, vi, tl, nan_handling, []() { + // This cannot be called because it would access the data + // with old metadata + // check_valid(leveldata, groupdata, vi, tl, nan_handling, []() { // return "MakeNewLevelFromCoarse after prolongation"; // }); } @@ -1958,8 +1949,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 diff --git a/CarpetX/src/driver.hxx b/CarpetX/src/driver.hxx index 488e1e9bb..46e10fbdc 100644 --- a/CarpetX/src/driver.hxx +++ b/CarpetX/src/driver.hxx @@ -40,7 +40,6 @@ using rat64 = rational; // Symmetries are domain properties enum class symmetry_t { none, - outer_boundary, interpatch, periodic, reflection, @@ -48,7 +47,7 @@ enum class symmetry_t { std::ostream &operator<<(std::ostream &os, const symmetry_t symmetry); // Boundary conditions are group properties. They are valid only for faces where -// the domain symmetry is `outer_boundary`. +// the domain symmetry is `none`. enum class boundary_t { none, symmetry_boundary, @@ -209,7 +208,6 @@ struct GHExt { int patch; array, 2> symmetries; - bool all_faces_have_symmetries() const; // AMReX grid structure // TODO: convert this from unique_ptr to optional @@ -265,6 +263,7 @@ struct GHExt { array nghostzones; array, 2> boundaries; + bool all_faces_have_symmetries_or_boundaries() const; vector > parities; vector dirichlet_values; vector robin_values; diff --git a/CarpetX/src/schedule.cxx b/CarpetX/src/schedule.cxx index 079209176..ba177668d 100644 --- a/CarpetX/src/schedule.cxx +++ b/CarpetX/src/schedule.cxx @@ -117,11 +117,9 @@ GridDesc::GridDesc(const GHExt::PatchData::LevelData &leveldata, } // Boundaries - const auto &symmetries = ghext->patchdata.at(leveldata.patch).symmetries; for (int d = 0; d < dim; ++d) for (int f = 0; f < 2; ++f) - bbox[f][d] = vbx[orient(d, f)] == domain[orient(d, f)] && - symmetries[f][d] != symmetry_t::none; + bbox[f][d] = vbx[orient(d, f)] == domain[orient(d, f)]; // Thread tile box for (int d = 0; d < dim; ++d) { @@ -2229,7 +2227,7 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups, std::vector > tasks1; std::vector > tasks2; - const bool have_multipatch_boundaries = + static const bool have_multipatch_boundaries = CCTK_IsFunctionAliased("MultiPatch_Interpolate"); active_levels->loop([&](auto &restrict leveldata) { @@ -2270,37 +2268,35 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups, ghext->patchdata.at(leveldata.patch) .amrcore->Geom(leveldata.level)); - tasks1.emplace_back([&tasks2, &leveldata, &groupdata, - have_multipatch_boundaries, nan_handling, tl, - fillpatch_continue = - std::move(fillpatch_continue)]() { - auto fillpatch_finish = fillpatch_continue(); - - tasks2.emplace_back([&leveldata, &groupdata, - have_multipatch_boundaries, nan_handling, tl, - fillpatch_finish = - std::move(fillpatch_finish)]() { - fillpatch_finish(); - - for (int vi = 0; vi < groupdata.numvars; ++vi) { - groupdata.valid.at(tl).at(vi).set_ghosts(true, []() { - return "SyncGroupsByDirI after syncing: " - "Mark ghost zones as valid"; - }); - if (ghext->patchdata.at(leveldata.patch) - .all_faces_have_symmetries()) - groupdata.valid.at(tl).at(vi).set_outer(true, []() { - return "SyncGroupsByDirI after syncing: " - "Mark outer boundaries as valid"; - }); - poison_invalid(leveldata, groupdata, vi, tl); - if (!have_multipatch_boundaries) - check_valid(leveldata, groupdata, vi, tl, nan_handling, []() { - return "SyncGroupsByDirI after syncing"; - }); - } - }); - }); + tasks1.emplace_back( + [&tasks2, &leveldata, &groupdata, nan_handling, tl, + fillpatch_continue = std::move(fillpatch_continue)]() { + auto fillpatch_finish = fillpatch_continue(); + + tasks2.emplace_back( + [&leveldata, &groupdata, nan_handling, tl, + fillpatch_finish = std::move(fillpatch_finish)]() { + fillpatch_finish(); + + for (int vi = 0; vi < groupdata.numvars; ++vi) { + groupdata.valid.at(tl).at(vi).set_ghosts(true, []() { + return "SyncGroupsByDirI after syncing: " + "Mark ghost zones as valid"; + }); + 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"; + }); + poison_invalid(leveldata, groupdata, vi, tl); + if (!have_multipatch_boundaries) + check_valid(leveldata, groupdata, vi, tl, + nan_handling, []() { + return "SyncGroupsByDirI after syncing"; + }); + } + }); + }); } // for tl } else { // if leveldata.level > 0 @@ -2361,8 +2357,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"; @@ -2393,7 +2388,7 @@ int SyncGroupsByDirI(const cGH *restrict cctkGH, int numgroups, tasks2.clear(); synchronize(); - if (CCTK_IsFunctionAliased("MultiPatch_Interpolate")) { + if (have_multipatch_boundaries) { std::vector cactusvarinds; for (int group : groups) { const auto &groupdata = diff --git a/Loop/src/loop.hxx b/Loop/src/loop.hxx index 8abae7d43..453016b67 100644 --- a/Loop/src/loop.hxx +++ b/Loop/src/loop.hxx @@ -338,9 +338,9 @@ public: for (int ni = -1; ni <= +1; ++ni) { if ((ni == 0) + (nj == 0) + (nk == 0) == rank) { - if ((ni != 0 && bbox[ni == -1 ? 0 : 1][0]) || - (nj != 0 && bbox[nj == -1 ? 0 : 1][1]) || - (nk != 0 && bbox[nk == -1 ? 0 : 1][2])) { + if ((ni != 0 && bbox[ni < 0 ? 0 : 1][0]) || + (nj != 0 && bbox[nj < 0 ? 0 : 1][1]) || + (nk != 0 && bbox[nk < 0 ? 0 : 1][2])) { const vect inormal{ni, nj, nk}; @@ -394,9 +394,9 @@ public: for (int ni = -1; ni <= +1; ++ni) { if ((ni == 0) + (nj == 0) + (nk == 0) == rank) { - if ((ni != 0 && !bbox[ni == -1 ? 0 : 1][0]) || - (nj != 0 && !bbox[nj == -1 ? 0 : 1][1]) || - (nk != 0 && !bbox[nk == -1 ? 0 : 1][2])) { + if ((ni != 0 && !bbox[ni < 0 ? 0 : 1][0]) || + (nj != 0 && !bbox[nj < 0 ? 0 : 1][1]) || + (nk != 0 && !bbox[nk < 0 ? 0 : 1][2])) { const vect inormal{ni, nj, nk}; @@ -478,9 +478,9 @@ public: for (int ni = -1; ni <= +1; ++ni) { if ((ni == 0) + (nj == 0) + (nk == 0) == rank) { - if ((ni == 0 || !bbox[ni == -1 ? 0 : 1][0]) && - (nj == 0 || !bbox[nj == -1 ? 0 : 1][1]) && - (nk == 0 || !bbox[nk == -1 ? 0 : 1][2])) { + if ((ni == 0 || !bbox[ni < 0 ? 0 : 1][0]) && + (nj == 0 || !bbox[nj < 0 ? 0 : 1][1]) && + (nk == 0 || !bbox[nk < 0 ? 0 : 1][2])) { const vect inormal{ni, nj, nk}; @@ -1351,6 +1351,8 @@ template struct is_GF3D5 : std::false_type {}; template struct is_GF3D5 > : std::true_type {}; template inline constexpr bool is_GF3D5_v = is_GF3D5::value; +//////////////////////////////////////////////////////////////////////////////// + template struct GF3D5vector { static_assert((std::is_same_v)); typedef T value_type; diff --git a/Loop/src/loop_device.hxx b/Loop/src/loop_device.hxx index da4854a87..1536b435c 100644 --- a/Loop/src/loop_device.hxx +++ b/Loop/src/loop_device.hxx @@ -220,9 +220,9 @@ public: for (int ni = -1; ni <= +1; ++ni) { if ((ni == 0) + (nj == 0) + (nk == 0) == rank) { - if ((ni != 0 && bbox[ni == -1 ? 0 : 1][0]) || - (nj != 0 && bbox[nj == -1 ? 0 : 1][1]) || - (nk != 0 && bbox[nk == -1 ? 0 : 1][2])) { + if ((ni != 0 && bbox[ni < 0 ? 0 : 1][0]) || + (nj != 0 && bbox[nj < 0 ? 0 : 1][1]) || + (nk != 0 && bbox[nk < 0 ? 0 : 1][2])) { const vect inormal{ni, nj, nk}; @@ -278,9 +278,9 @@ public: for (int ni = -1; ni <= +1; ++ni) { if ((ni == 0) + (nj == 0) + (nk == 0) == rank) { - if ((ni != 0 && !bbox[ni == -1 ? 0 : 1][0]) || - (nj != 0 && !bbox[nj == -1 ? 0 : 1][1]) || - (nk != 0 && !bbox[nk == -1 ? 0 : 1][2])) { + if ((ni != 0 && !bbox[ni<0 ? 0 : 1][0]) || + (nj != 0 && !bbox[nj<0 ? 0 : 1][1]) || + (nk != 0 && !bbox[nk<0 ? 0 : 1][2])) { const vect inormal{ni, nj, nk}; @@ -344,9 +344,9 @@ public: for (int ni = -1; ni <= +1; ++ni) { if ((ni == 0) + (nj == 0) + (nk == 0) == rank) { - if ((ni == 0 || !bbox[ni == -1 ? 0 : 1][0]) && - (nj == 0 || !bbox[nj == -1 ? 0 : 1][1]) && - (nk == 0 || !bbox[nk == -1 ? 0 : 1][2])) { + if ((ni == 0 || !bbox[ni < 0 ? 0 : 1][0]) && + (nj == 0 || !bbox[nj < 0 ? 0 : 1][1]) && + (nk == 0 || !bbox[nk < 0 ? 0 : 1][2])) { const vect inormal{ni, nj, nk}; diff --git a/TestOutput/test/output-every0/.empty-directory b/TestOutput/test/output-every0/.empty-directory new file mode 100644 index 000000000..e69de29bb diff --git a/TestOutput/test/output-tsv-every0/.empty-directory b/TestOutput/test/output-tsv-every0/.empty-directory new file mode 100644 index 000000000..e69de29bb