Refactor Region Size

CHANGELOG.md

@@ -26,6 +26,7 @@
 - [[PR 890]](https://github.com/parthenon-hpc-lab/parthenon/pull/890) Fix bugs in sparse communication and prolongation
 
 ### Infrastructure (changes irrelevant to downstream codes)
+- [[PR 918]](https://github.com/parthenon-hpc-lab/parthenon/pull/918) Refactor RegionSize 
 - [[PR 901]](https://github.com/parthenon-hpc-lab/parthenon/pull/901) Implement shared element ownership model
 
 ### Removed (removing behavior/API/varaibles/...)

example/advection/parthenon_app_inputs.cpp

@@ -204,8 +204,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, SimTime &tm) {
   if (Globals::my_rank == 0) {
     // normalize errors by number of cells
     auto mesh_size = mesh->mesh_size;
-    Real vol = (mesh_size.x1max - mesh_size.x1min) * (mesh_size.x2max - mesh_size.x2min) *
-               (mesh_size.x3max - mesh_size.x3min);
+    Real vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) *
+               (mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) *
+               (mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR));
     l1_err /= vol;
     // compute rms error
     max_max_over_l1 = std::max(max_max_over_l1, (max_err / l1_err));
@@ -237,8 +238,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, SimTime &tm) {
     }
 
     // write errors
-    std::fprintf(pfile, "%d  %d", mesh_size.nx1, mesh_size.nx2);
-    std::fprintf(pfile, "  %d  %d", mesh_size.nx3, tm.ncycle);
+    std::fprintf(pfile, "%d  %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR));
+    std::fprintf(pfile, "  %d  %d", mesh_size.nx(X3DIR), tm.ncycle);
     std::fprintf(pfile, "  %e ", l1_err);
     std::fprintf(pfile, "  %e  %e  ", max_max_over_l1, max_err);
     std::fprintf(pfile, "\n");

example/calculate_pi/calculate_pi.cpp

@@ -49,10 +49,10 @@ void SetInOrOut(MeshBlockData<Real> *rc) {
   if (use_sparse) {
     auto &bs = pmb->block_size;
     // check if block falls on radius.
-    Real coords[4][2] = {{bs.x1min, bs.x2min},
-                         {bs.x1min, bs.x2max},
-                         {bs.x1max, bs.x2min},
-                         {bs.x1max, bs.x2max}};
+    Real coords[4][2] = {{bs.xmin(X1DIR), bs.xmin(X2DIR)},
+                         {bs.xmin(X1DIR), bs.xmax(X2DIR)},
+                         {bs.xmax(X1DIR), bs.xmin(X2DIR)},
+                         {bs.xmax(X1DIR), bs.xmax(X2DIR)}};
 
     bool fully_outside = true;
 

example/particle_tracers/particle_tracers.cpp

@@ -345,12 +345,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   const Real &z_max = pmb->coords.Xf<3>(kb.e + 1);
 
   const auto mesh_size = pmb->pmy_mesh->mesh_size;
-  const Real x_min_mesh = mesh_size.x1min;
-  const Real y_min_mesh = mesh_size.x2min;
-  const Real z_min_mesh = mesh_size.x3min;
-  const Real x_max_mesh = mesh_size.x1max;
-  const Real y_max_mesh = mesh_size.x2max;
-  const Real z_max_mesh = mesh_size.x3max;
+  const Real x_min_mesh = mesh_size.xmin(X1DIR);
+  const Real y_min_mesh = mesh_size.xmin(X2DIR);
+  const Real z_min_mesh = mesh_size.xmin(X3DIR);
+  const Real x_max_mesh = mesh_size.xmax(X1DIR);
+  const Real y_max_mesh = mesh_size.xmax(X2DIR);
+  const Real z_max_mesh = mesh_size.xmax(X3DIR);
 
   const Real kwave = 2. * M_PI / (x_max_mesh - x_min_mesh);
 

example/stochastic_subgrid/parthenon_app_inputs.cpp

@@ -171,9 +171,9 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, SimTime &tm) {
     if (Globals::my_rank == 0) {
       // normalize errors by number of cells
       auto mesh_size = mesh->mesh_size;
-      Real vol = (mesh_size.x1max - mesh_size.x1min) *
-                 (mesh_size.x2max - mesh_size.x2min) *
-                 (mesh_size.x3max - mesh_size.x3min);
+      Real vol = (mesh_size.xmax(X1DIR) - mesh_size.xmin(X1DIR)) *
+                 (mesh_size.xmax(X2DIR) - mesh_size.xmin(X2DIR)) *
+                 (mesh_size.xmax(X3DIR) - mesh_size.xmin(X3DIR));
       l1_err /= vol;
       // compute rms error
       max_max_over_l1 = std::max(max_max_over_l1, (max_err / l1_err));
@@ -205,8 +205,8 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, SimTime &tm) {
       }
 
       // write errors
-      std::fprintf(pfile, "%d  %d", mesh_size.nx1, mesh_size.nx2);
-      std::fprintf(pfile, "  %d  %d", mesh_size.nx3, tm.ncycle);
+      std::fprintf(pfile, "%d  %d", mesh_size.nx(X1DIR), mesh_size.nx(X2DIR));
+      std::fprintf(pfile, "  %d  %d", mesh_size.nx(X3DIR), tm.ncycle);
       std::fprintf(pfile, "  %e ", l1_err);
       std::fprintf(pfile, "  %e  %e  ", max_max_over_l1, max_err);
       std::fprintf(pfile, "\n");

src/bvals/bvals.cpp

@@ -71,14 +71,14 @@ BoundaryValues::BoundaryValues(std::weak_ptr<MeshBlock> wpmb, BoundaryFlag *inpu
   CheckBoundaryFlag(block_bcs[BoundaryFace::outer_x1], CoordinateDirection::X1DIR);
 
   std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
-  if (pmb->block_size.nx2 > 1) {
+  if (!pmb->block_size.symmetry(X2DIR)) {
     nface_ = 4;
     nedge_ = 4;
     CheckBoundaryFlag(block_bcs[BoundaryFace::inner_x2], CoordinateDirection::X2DIR);
     CheckBoundaryFlag(block_bcs[BoundaryFace::outer_x2], CoordinateDirection::X2DIR);
   }
 
-  if (pmb->block_size.nx3 > 1) {
+  if (!pmb->block_size.symmetry(X3DIR)) {
     nface_ = 6;
     nedge_ = 12;
     CheckBoundaryFlag(block_bcs[BoundaryFace::inner_x3], CoordinateDirection::X3DIR);
@@ -151,14 +151,14 @@ BoundarySwarms::BoundarySwarms(std::weak_ptr<MeshBlock> wpmb, BoundaryFlag *inpu
   CheckBoundaryFlag(block_bcs[BoundaryFace::outer_x1], CoordinateDirection::X1DIR);
 
   std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
-  if (pmb->block_size.nx2 > 1) {
+  if (!pmb->block_size.symmetry(X2DIR)) {
     nface_ = 4;
     nedge_ = 4;
     CheckBoundaryFlag(block_bcs[BoundaryFace::inner_x2], CoordinateDirection::X2DIR);
     CheckBoundaryFlag(block_bcs[BoundaryFace::outer_x2], CoordinateDirection::X2DIR);
   }
 
-  if (pmb->block_size.nx3 > 1) {
+  if (!pmb->block_size.symmetry(X3DIR)) {
     nface_ = 6;
     nedge_ = 12;
     CheckBoundaryFlag(block_bcs[BoundaryFace::inner_x3], CoordinateDirection::X3DIR);

src/bvals/bvals_base.cpp

@@ -114,7 +114,7 @@ BoundaryBase::BoundaryBase(Mesh *pm, LogicalLocation iloc, RegionSize isize,
 
   if (pmy_mesh_->multilevel) { // SMR or AMR
     // allocate surface area array
-    int nc1 = block_size_.nx1 + 2 * Globals::nghost;
+    int nc1 = block_size_.nx(X1DIR) + 2 * Globals::nghost;
     sarea_[0] = ParArrayND<Real>(PARARRAY_TEMP, nc1);
     sarea_[1] = ParArrayND<Real>(PARARRAY_TEMP, nc1);
   }
@@ -309,13 +309,13 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
   myfx2 = ((loc.lx2() & 1LL) == 1LL);
   myfx3 = ((loc.lx3() & 1LL) == 1LL);
   myox1 = ((loc.lx1() & 1LL) == 1LL) * 2 - 1;
-  if (block_size_.nx2 > 1) myox2 = ((loc.lx2() & 1LL) == 1LL) * 2 - 1;
-  if (block_size_.nx3 > 1) myox3 = ((loc.lx3() & 1LL) == 1LL) * 2 - 1;
+  if (!block_size_.symmetry(X2DIR)) myox2 = ((loc.lx2() & 1LL) == 1LL) * 2 - 1;
+  if (!block_size_.symmetry(X3DIR)) myox3 = ((loc.lx3() & 1LL) == 1LL) * 2 - 1;
 
   int nf1 = 1, nf2 = 1;
   if (pmy_mesh_->multilevel) {
-    if (block_size_.nx2 > 1) nf1 = 2;
-    if (block_size_.nx3 > 1) nf2 = 2;
+    if (!block_size_.symmetry(X2DIR)) nf1 = 2;
+    if (!block_size_.symmetry(X3DIR)) nf2 = 2;
   }
   int bufid = 0;
   nneighbor = 0;
@@ -367,7 +367,7 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
       nneighbor++;
     }
   }
-  if (block_size_.nx2 == 1) {
+  if (block_size_.nx(X2DIR) == 1) {
     SetNeighborOwnership();
     Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
     return;
@@ -415,7 +415,7 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
   }
 
   // x3 face
-  if (block_size_.nx3 > 1) {
+  if (!block_size_.symmetry(X3DIR)) {
     for (int n = -1; n <= 1; n += 2) {
       neibt = tree.FindNeighbor(loc, 0, 0, n);
       if (neibt == nullptr) {
@@ -502,7 +502,7 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
     }
   }
 
-  if (block_size_.nx3 == 1) {
+  if (block_size_.nx(X3DIR) == 1) {
     SetNeighborOwnership();
     Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
     return;

src/bvals/comms/bnd_info.cpp

@@ -413,41 +413,41 @@ BndInfo BndInfo::GetSetCCFluxCor(std::shared_ptr<MeshBlock> pmb, const NeighborB
     else
       si = ++ei;
     if (nb.ni.fi1 == 0)
-      ej -= pmb->block_size.nx2 / 2;
+      ej -= pmb->block_size.nx(X2DIR) / 2;
     else
-      sj += pmb->block_size.nx2 / 2;
+      sj += pmb->block_size.nx(X2DIR) / 2;
     if (nb.ni.fi2 == 0)
-      ek -= pmb->block_size.nx3 / 2;
+      ek -= pmb->block_size.nx(X3DIR) / 2;
     else
-      sk += pmb->block_size.nx3 / 2;
+      sk += pmb->block_size.nx(X3DIR) / 2;
   } else if (nb.fid == BoundaryFace::inner_x2 || nb.fid == BoundaryFace::outer_x2) {
     out.dir = X2DIR;
     if (nb.fid == BoundaryFace::inner_x2)
       ej = sj;
     else
       sj = ++ej;
     if (nb.ni.fi1 == 0)
-      ei -= pmb->block_size.nx1 / 2;
+      ei -= pmb->block_size.nx(X1DIR) / 2;
     else
-      si += pmb->block_size.nx1 / 2;
+      si += pmb->block_size.nx(X1DIR) / 2;
     if (nb.ni.fi2 == 0)
-      ek -= pmb->block_size.nx3 / 2;
+      ek -= pmb->block_size.nx(X3DIR) / 2;
     else
-      sk += pmb->block_size.nx3 / 2;
+      sk += pmb->block_size.nx(X3DIR) / 2;
   } else if (nb.fid == BoundaryFace::inner_x3 || nb.fid == BoundaryFace::outer_x3) {
     out.dir = X3DIR;
     if (nb.fid == BoundaryFace::inner_x3)
       ek = sk;
     else
       sk = ++ek;
     if (nb.ni.fi1 == 0)
-      ei -= pmb->block_size.nx1 / 2;
+      ei -= pmb->block_size.nx(X1DIR) / 2;
     else
-      si += pmb->block_size.nx1 / 2;
+      si += pmb->block_size.nx(X1DIR) / 2;
     if (nb.ni.fi2 == 0)
-      ej -= pmb->block_size.nx2 / 2;
+      ej -= pmb->block_size.nx(X2DIR) / 2;
     else
-      sj += pmb->block_size.nx2 / 2;
+      sj += pmb->block_size.nx(X2DIR) / 2;
   } else {
     PARTHENON_FAIL("Flux corrections only occur on faces for CC variables.");
   }

src/coordinates/uniform_cartesian.hpp

@@ -28,19 +28,15 @@ class UniformCartesian {
  public:
   UniformCartesian() = default;
   UniformCartesian(const RegionSize &rs, ParameterInput *pin) {
-    dx_[0] = (rs.x1max - rs.x1min) / rs.nx1;
-    dx_[1] = (rs.x2max - rs.x2min) / rs.nx2;
-    dx_[2] = (rs.x3max - rs.x3min) / rs.nx3;
+    for (auto &dir : {X1DIR, X2DIR, X3DIR}) {
+      dx_[dir - 1] = (rs.xmax(dir) - rs.xmin(dir)) / rs.nx(dir);
+      istart_[dir - 1] = (!rs.symmetry(dir) ? Globals::nghost : 0);
+      xmin_[dir - 1] = rs.xmin(dir) - istart_[dir - 1] * dx_[dir - 1];
+    }
     area_[0] = dx_[1] * dx_[2];
     area_[1] = dx_[0] * dx_[2];
     area_[2] = dx_[0] * dx_[1];
     cell_volume_ = dx_[0] * dx_[1] * dx_[2];
-    istart_[0] = Globals::nghost;
-    istart_[1] = (rs.nx2 > 1 ? Globals::nghost : 0);
-    istart_[2] = (rs.nx3 > 1 ? Globals::nghost : 0);
-    xmin_[0] = rs.x1min - istart_[0] * dx_[0];
-    xmin_[1] = rs.x2min - istart_[1] * dx_[1];
-    xmin_[2] = rs.x3min - istart_[2] * dx_[2];
   }
   UniformCartesian(const UniformCartesian &src, int coarsen)
       : istart_(src.GetStartIndex()) {

src/defs.hpp

@@ -57,23 +57,6 @@ static_assert(NDIM >= 3,
 //! \struct RegionSize
 //  \brief physical size and number of cells in a Mesh or a MeshBlock
 
-struct RegionSize { // aggregate and POD type; do NOT reorder member declarations:
-  Real x1min, x2min, x3min;
-  Real x1max, x2max, x3max;
-  Real x1rat, x2rat, x3rat; // ratio of dxf(i)/dxf(i-1)
-  // the size of the root grid or a MeshBlock should not exceed std::int32_t limits
-  int nx1, nx2, nx3; // number of active cells (not including ghost zones)
-};
-
-//----------------------------------------------------------------------------------------
-// enums used everywhere
-// (not specifying underlying integral type (C++11) for portability & performance)
-
-// TODO(felker): C++ Core Guidelines Enum.5: Don’t use ALL_CAPS for enumerators
-// (avoid clashes with preprocessor macros). Enumerated type definitions in this file and:
-// io_wrapper.hpp, bvals.hpp, field_diffusion.hpp,
-// task_list.hpp, ???
-
 //------------------
 // named, weakly typed / unscoped enums:
 //------------------
@@ -86,6 +69,45 @@ struct RegionSize { // aggregate and POD type; do NOT reorder member declaration
 // X3DIR z, phi, etc...
 enum CoordinateDirection { NODIR = -1, X0DIR = 0, X1DIR = 1, X2DIR = 2, X3DIR = 3 };
 
+struct RegionSize {
+  RegionSize() = default;
+  RegionSize(std::array<Real, 3> xmin, std::array<Real, 3> xmax, std::array<Real, 3> xrat,
+             std::array<int, 3> nx)
+      : xmin_(xmin), xmax_(xmax), xrat_(xrat),
+        nx_(nx), symmetry_{nx[0] == 1, nx[1] == 1, nx[2] == 1} {}
+  RegionSize(std::array<Real, 3> xmin, std::array<Real, 3> xmax, std::array<Real, 3> xrat,
+             std::array<int, 3> nx, std::array<bool, 3> symmetry)
+      : xmin_(xmin), xmax_(xmax), xrat_(xrat), nx_(nx), symmetry_(symmetry) {}
+
+  std::array<Real, 3> xmin_, xmax_, xrat_; // xrat is ratio of dxf(i)/dxf(i-1)
+  std::array<int, 3> nx_;
+  std::array<bool, 3> symmetry_;
+
+  Real &xmin(CoordinateDirection dir) { return xmin_[dir - 1]; }
+  const Real &xmin(CoordinateDirection dir) const { return xmin_[dir - 1]; }
+
+  Real &xmax(CoordinateDirection dir) { return xmax_[dir - 1]; }
+  const Real &xmax(CoordinateDirection dir) const { return xmax_[dir - 1]; }
+
+  Real &xrat(CoordinateDirection dir) { return xrat_[dir - 1]; }
+  const Real &xrat(CoordinateDirection dir) const { return xrat_[dir - 1]; }
+
+  int &nx(CoordinateDirection dir) { return nx_[dir - 1]; }
+  const int &nx(CoordinateDirection dir) const { return nx_[dir - 1]; }
+
+  bool &symmetry(CoordinateDirection dir) { return symmetry_[dir - 1]; }
+  const bool &symmetry(CoordinateDirection dir) const { return symmetry_[dir - 1]; }
+};
+
+//----------------------------------------------------------------------------------------
+// enums used everywhere
+// (not specifying underlying integral type (C++11) for portability & performance)
+
+// TODO(felker): C++ Core Guidelines Enum.5: Don’t use ALL_CAPS for enumerators
+// (avoid clashes with preprocessor macros). Enumerated type definitions in this file and:
+// io_wrapper.hpp, bvals.hpp, field_diffusion.hpp,
+// task_list.hpp, ???
+
 // identifiers for all 6 faces of a MeshBlock
 constexpr int BOUNDARY_NFACES = 6;
 enum BoundaryFace {
@@ -98,6 +120,28 @@ enum BoundaryFace {
   outer_x3 = 5
 };
 
+inline BoundaryFace GetInnerBoundaryFace(CoordinateDirection dir) {
+  if (dir == X1DIR) {
+    return BoundaryFace::inner_x1;
+  } else if (dir == X2DIR) {
+    return BoundaryFace::inner_x2;
+  } else if (dir == X3DIR) {
+    return BoundaryFace::inner_x3;
+  }
+  return BoundaryFace::undef;
+}
+
+inline BoundaryFace GetOuterBoundaryFace(CoordinateDirection dir) {
+  if (dir == X1DIR) {
+    return BoundaryFace::outer_x1;
+  } else if (dir == X2DIR) {
+    return BoundaryFace::outer_x2;
+  } else if (dir == X3DIR) {
+    return BoundaryFace::outer_x3;
+  }
+  return BoundaryFace::undef;
+}
+
 //------------------
 // strongly typed / scoped enums (C++11):
 //------------------

src/interface/swarm.cpp

@@ -50,12 +50,12 @@ SwarmDeviceContext Swarm::GetDeviceContext() const {
   context.x_max_ = pmb->coords.Xf<1>(ib.e + 1);
   context.y_max_ = pmb->coords.Xf<2>(jb.e + 1);
   context.z_max_ = pmb->coords.Xf<3>(kb.e + 1);
-  context.x_min_global_ = mesh_size.x1min;
-  context.x_max_global_ = mesh_size.x1max;
-  context.y_min_global_ = mesh_size.x2min;
-  context.y_max_global_ = mesh_size.x2max;
-  context.z_min_global_ = mesh_size.x3min;
-  context.z_max_global_ = mesh_size.x3max;
+  context.x_min_global_ = mesh_size.xmin(X1DIR);
+  context.x_max_global_ = mesh_size.xmax(X1DIR);
+  context.y_min_global_ = mesh_size.xmin(X2DIR);
+  context.y_max_global_ = mesh_size.xmax(X2DIR);
+  context.z_min_global_ = mesh_size.xmin(X3DIR);
+  context.z_max_global_ = mesh_size.xmax(X3DIR);
   context.ndim_ = pmb->pmy_mesh->ndim;
   context.my_rank_ = Globals::my_rank;
   context.coords_ = pmb->coords;