eschnett/cactusamrex - Change XU5HOJREK4XY4NBCJINLZPKQNKSYOLUDTWR47REFSNQKDOSNDXLQC

CarpetX: Implement reflective boundary conditions (not working yet)

Created by Erik Schnetter on January 8, 2020

XU5HOJREK4XY4NBCJINLZPKQNKSYOLUDTWR47REFSNQKDOSNDXLQC

Dependencies

In channels

main

master

Change contents

Insertion in CarpetX/param.ccl at line 81 [6.17918]
[4.48]
[5.230]

Insertion in CarpetX/param.ccl at line 85 [6.17918]

[4.58]

[7.2]


BOOLEAN reflection_x "Reflection symmetry at the lower x boundary"
{
} "no"
BOOLEAN reflection_y "Reflection symmetry at the lower y boundary"
{
} "no"

Insertion in CarpetX/param.ccl at line 94 [6.17918]

[7.3]

[5.231]

BOOLEAN reflection_z "Reflection symmetry at the lower z boundary"
{
} "no"
BOOLEAN reflection_upper_x "Reflection symmetry at the upper x boundary"
{
} "no"
BOOLEAN reflection_upper_y "Reflection symmetry at the upper y boundary"
{
} "no"
BOOLEAN reflection_upper_z "Reflection symmetry at the upper z boundary"
{
} "no"

Insertion in CarpetX/src/driver.cxx at line 25 [6.15972]
[6.16127]
[6.16169]
```
#include <utility>
```

Insertion in CarpetX/src/driver.cxx at line 190 [6.15972]

[2.3544]

vector<array<int, dim> > get_group_parities(const int gi) {
  DECLARE_CCTK_PARAMETERS;
  assert(gi >= 0);
  const int tags = CCTK_GroupTagsTableI(gi);
  assert(tags >= 0);
  const int nelems = Util_TableGetIntArray(tags, 0, nullptr, "parities");
  if (nelems == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
    // unset (will use index type)
    return {};
  } else if (nelems >= 0) {
    // do nothing
  } else {
    assert(0);
  }
  vector<CCTK_INT> parities1(nelems);
  const int iret =
      Util_TableGetIntArray(tags, nelems, parities1.data(), "parities");
  assert(iret == nelems);
  assert(nelems % dim == 0);
  vector<array<int, dim> > parities(nelems / dim);
  for (size_t n = 0; n < parities.size(); ++n)
    for (int d = 0; d < dim; ++d)
      parities.at(n)[d] = parities1.at(dim * n + d);
  return parities;
}

Insertion in CarpetX/src/driver.cxx at line 389 [6.15972]

[8.1487]

[9.2191]

    groupdata.parities = get_group_parities(gi);
    if (groupdata.parities.empty()) {
      array<int, dim> parity;
      for (int d = 0; d < dim; ++d)
        parity[d] = groupdata.indextype[d] == 0 ? +1 : -1;
      groupdata.parities.resize(groupdata.numvars, parity);
    }
    assert(int(groupdata.parities.size()) == groupdata.numvars);
    for (int vi = 0; vi < groupdata.numvars; ++vi)
      for (int d = 0; d < dim; ++d)
        assert(abs(groupdata.parities.at(vi)[d]) == 1);

Insertion in CarpetX/src/driver.cxx at line 695 [6.15972]

[10.466]

}
void apply_physbcs(const Box &, const FArrayBox &, int, int, const Geometry &,
                   CCTK_REAL, const Vector<BCRec> &, int, int) { // do nothing

Insertion in CarpetX/src/driver.cxx at line 700 [6.15972]

[10.468]


tuple<CarpetXPhysBCFunct, Vector<BCRec> >
get_boundaries(const GHExt::LevelData &leveldata,
               const GHExt::LevelData::GroupData &groupdata) {
  DECLARE_CCTK_PARAMETERS;

Insertion in CarpetX/src/driver.cxx at line 706 [6.15972]

[10.469]

[9.2450]

  const array<array<bool, 3>, 2> is_periodic{{
      {{
          periodic || periodic_x,
          periodic || periodic_y,
          periodic || periodic_z,
      }},
      {{
          periodic || periodic_x,
          periodic || periodic_y,
          periodic || periodic_z,
      }},
  }};
  const array<array<bool, 3>, 2> is_reflect{{
      {{
          !!reflection_x,
          !!reflection_y,
          !!reflection_z,
      }},
      {{
          !!reflection_upper_x,
          !!reflection_upper_y,
          !!reflection_upper_z,
      }},
  }};
  const auto makebc{[&](const int vi, const int dir, const int face) {
    assert(dir >= 0 && dir < dim);
    assert(face >= 0 && face < 2);
    if (is_periodic[face][dir])
      return BCType::int_dir;
    if (is_reflect[face][dir])
      return groupdata.parities.at(vi)[dir] > 0 ? BCType::reflect_even
                                                : BCType::reflect_odd;
    return BCType::ext_dir;
  }};
  Vector<BCRec> bcs(groupdata.numvars);
  for (int vi = 0; vi < groupdata.numvars; ++vi)
    bcs.at(vi) = BCRec(makebc(vi, 0, 0), makebc(vi, 1, 0), makebc(vi, 2, 0),
                       makebc(vi, 0, 1), makebc(vi, 1, 1), makebc(vi, 2, 1));
  const auto apply_physbc{[](const Box &, const FArrayBox &, int, int,
                             const Geometry &, CCTK_REAL, const Vector<BCRec> &,
                             int, int) {}};
  CarpetXPhysBCFunct physbc(ghext->amrcore->Geom(leveldata.level), bcs,
                            apply_physbcs);
  return {move(physbc), move(bcs)};
}

Deletion in CarpetX/src/driver.cxx at line 799 [6.15972]
∅:D[11.697] → [9.3180:3238]
∅:D[2.4224] → [9.3180:3238]
B:BD[9.3180] → [9.3180:3238]
```
    PhysBCFunctNoOp cphysbc;
    PhysBCFunctNoOp fphysbc;
```

Replacement in CarpetX/src/driver.cxx at line 800 [6.15972]

B:BD[9.3274] → [12.1165:1192]

B:BD[12.1192] → [7.109:541]

∅:D[4.576] → [9.3448:3503]

∅:D[12.1637] → [9.3448:3503]

∅:D[7.541] → [9.3448:3503]

B:BD[9.3448] → [9.3448:3503]

    // boundary conditions
    const BCRec bcrec(
        periodic || periodic_x ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_y ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_z ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_x ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_y ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_z ? BCType::int_dir : BCType::ext_dir);
    const Vector<BCRec> bcs(groupdata.numvars, bcrec);

[9.3274]

[9.3503]


    auto physbc_bcs = get_boundaries(leveldata, groupdata);
    CarpetXPhysBCFunct &physbc = get<0>(physbc_bcs);
    const Vector<BCRec> &bcs = get<1>(physbc_bcs);

Replacement in CarpetX/src/driver.cxx at line 834 [6.15972]

∅:D[13.64] → [14.1354:1428]

B:BD[9.3931] → [14.1354:1428]

            ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0, reffact,

[13.64]

[10.528]

            ghext->amrcore->Geom(level), physbc, 0, physbc, 0, reffact,

Deletion in CarpetX/src/driver.cxx at line 890 [6.15972]
B:BD[9.5502] → [15.2321:2379]
```
    PhysBCFunctNoOp cphysbc;
    PhysBCFunctNoOp fphysbc;
```

Replacement in CarpetX/src/driver.cxx at line 891 [6.15972]

B:BD[15.2415] → [15.2415:2442]

B:BD[15.2442] → [7.542:974]

∅:D[4.1094] → [15.2887:2942]

∅:D[7.974] → [15.2887:2942]

B:BD[15.2887] → [15.2887:2942]

    // boundary conditions
    const BCRec bcrec(
        periodic || periodic_x ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_y ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_z ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_x ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_y ? BCType::int_dir : BCType::ext_dir,
        periodic || periodic_z ? BCType::int_dir : BCType::ext_dir);
    const Vector<BCRec> bcs(groupdata.numvars, bcrec);

[15.2415]

[2.4901]


    auto physbc_bcs = get_boundaries(leveldata, groupdata);
    CarpetXPhysBCFunct &physbc = get<0>(physbc_bcs);
    const Vector<BCRec> &bcs = get<1>(physbc_bcs);

Replacement in CarpetX/src/driver.cxx at line 952 [6.15972]

B:BD[2.6836] → [2.6836:6916]

                           ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0,

[2.6836]

[2.6916]

                           ghext->amrcore->Geom(level), physbc, 0, physbc, 0,

Insertion in CarpetX/src/driver.hxx at line 15 [6.15275]
[16.3199]
[17.2709]
```
#include <AMReX_FillPatchUtil.H>
```
Insertion in CarpetX/src/driver.hxx at line 278 [6.15275]
[8.2636]
[17.2742]
```
      vector<array<int, dim> > parities;
```

Insertion in CarpetX/src/driver.hxx at line 297 [6.15275]

[11.894]

[18.4387]


typedef void apply_physbcs_t(const Box &, const FArrayBox &, int, int,
                             const Geometry &, CCTK_REAL, const Vector<BCRec> &,
                             int, int);
typedef PhysBCFunct<apply_physbcs_t *> CarpetXPhysBCFunct;
tuple<CarpetXPhysBCFunct, Vector<BCRec> >
get_boundaries(const GHExt::LevelData &leveldata,
               const GHExt::LevelData::GroupData &groupdata);

Replacement in CarpetX/src/schedule.cxx at line 401 [20.531]

∅:D[3.2328] → [19.254:332]

B:BD[19.254] → [19.254:332]

               msg().c_str(), CCTK_FullVarName(groupdata.firstvarindex + vi),

[3.2328]

[3.2329]

        msg().c_str(), CCTK_FullVarName(groupdata.firstvarindex + vi),

Replacement in CarpetX/src/schedule.cxx at line 406 [20.531]

∅:D[3.2605] → [19.368:431]

B:BD[19.368] → [19.368:431]

               string(groupdata.valid.at(tl).at(vi)).c_str());

[3.2605]

[19.431]

        string(groupdata.valid.at(tl).at(vi)).c_str());

Replacement in CarpetX/src/schedule.cxx at line 519 [20.531]

∅:D[3.2957] → [21.8286:8348]

B:BD[22.924] → [21.8286:8348]

        string(scalargroupdata.valid.at(tl).at(vi)).c_str());

[3.2957]

[23.2]

                string(scalargroupdata.valid.at(tl).at(vi)).c_str());

Insertion in CarpetX/src/schedule.cxx at line 1956 [20.531]

[24.475]

      auto physbc_bcs = get_boundaries(leveldata, groupdata);
      CarpetXPhysBCFunct &physbc = get<0>(physbc_bcs);
      const Vector<BCRec> &bcs = get<1>(physbc_bcs);

Replacement in CarpetX/src/schedule.cxx at line 1976 [20.531]

∅:D[14.8103] → [24.621:668]

∅:D[25.2267] → [24.621:668]

∅:D[26.11086] → [24.621:668]

B:BD[24.621] → [24.621:668]

B:BD[24.668] → [16.5544:5612]

          groupdata.mfab.at(tl)->FillBoundary(
              ghext->amrcore->Geom(leveldata.level).periodicity());

[25.2267]

[26.11087]

          FillPatchSingleLevel(
              *groupdata.mfab.at(tl), 0.0, {&*groupdata.mfab.at(tl)}, {0.0}, 0,
              0, groupdata.numvars, ghext->amrcore->Geom(level), physbc, 0);

Replacement in CarpetX/src/schedule.cxx at line 1997 [20.531]
∅:D[11.8726] → [24.1081:1147]
B:BD[24.1081] → [24.1081:1147]
```
        PhysBCFunctNoOp cphysbc;
        PhysBCFunctNoOp fphysbc;
```
[11.8726]
[24.1147]

Deletion in CarpetX/src/schedule.cxx at line 1999 [20.531]

B:BD[24.1187] → [12.2424:2455]

B:BD[12.2455] → [27.10877:10904]

B:BD[27.10904] → [7.1311:1744]

∅:D[27.11361] → [24.1402:1461]

∅:D[4.1746] → [24.1402:1461]

∅:D[12.2924] → [24.1402:1461]

∅:D[7.1744] → [24.1402:1461]

B:BD[24.1402] → [24.1402:1461]

        // boundary conditions
        const BCRec bcrec(
            periodic || periodic_x ? BCType::int_dir : BCType::ext_dir,
            periodic || periodic_y ? BCType::int_dir : BCType::ext_dir,
            periodic || periodic_z ? BCType::int_dir : BCType::ext_dir,
            periodic || periodic_x ? BCType::int_dir : BCType::ext_dir,
            periodic || periodic_y ? BCType::int_dir : BCType::ext_dir,
            periodic || periodic_z ? BCType::int_dir : BCType::ext_dir);
        const Vector<BCRec> bcs(groupdata.numvars, bcrec);

Deletion in CarpetX/src/schedule.cxx at line 2001 [20.531]

B:BD[26.11313] → [28.0:1]

B:BD[28.1] → [2.10486:10844]

∅:D[21.14145] → [28.366:367]

∅:D[2.10844] → [28.366:367]

B:BD[28.366] → [28.366:367]

B:BD[28.367] → [2.10845:11199]

∅:D[2.11199] → [28.396:397]

B:BD[28.396] → [28.396:397]

B:BD[28.397] → [2.11200:11222]

∅:D[21.14482] → [28.598:599]

∅:D[2.11222] → [28.598:599]

B:BD[28.598] → [28.598:599]


          // // Only prolongate valid grid functions
          // bool all_invalid = true;
          // for (int vi = 0; vi < groupdata.numvars; ++vi)
          //   all_invalid &=
          //       !coarsegroupdata.valid.at(tl).at(vi).get().valid_any() &&
          //       !groupdata.valid.at(tl).at(vi).get().valid_int;
          // if (all_invalid) {
          //   for (int vi = 0; vi < groupdata.numvars; ++vi)
          //     groupdata.valid.at(tl).at(vi).set(valid_t(), [] {
          //       return "SyncGroupsByDirI skipping prolongation: Mark as "
          //              "invalid because neither coarse grid nor fine grid "
          //              "interior are valid";
          //     });
          // } else {

Insertion in CarpetX/src/schedule.cxx at line 2002 [20.531]

[2.11282]

#warning "TODO: interpolation does not require boundaries (nor ghosts?)"
#warning "TODO: also for regridding?"

Replacement in CarpetX/src/schedule.cxx at line 2029 [20.531]

B:BD[2.12607] → [2.12607:12677]

              cphysbc, 0, fphysbc, 0, reffact, interpolator, bcs, 0);

[2.12607]

[2.12677]

              physbc, 0, physbc, 0, reffact, interpolator, bcs, 0);