Merge pull request #591 from PrincetonUniversity/hotfix/mhd_amr_shear…

…_warning2

Add better exclusion for MHD+shearing box+SMR with refinement in x3 dir at shearing boundary

src/bvals/bvals.cpp

@@ -60,6 +60,7 @@
 BoundaryValues::BoundaryValues(MeshBlock *pmb, BoundaryFlag *input_bcs,
                                ParameterInput *pin)
     : BoundaryBase(pmb->pmy_mesh, pmb->loc, pmb->block_size, input_bcs), pmy_block_(pmb),
+      is_shear{}, loc_shear{0, pmy_mesh_->nrbx1*(1L << pmb->loc.level) - 1},
       sb_data_{}, sb_flux_data_{} {
   // Check BC functions for each of the 6 boundaries in turn ---------------------
   for (int i=0; i<6; i++) {
@@ -145,7 +146,6 @@ BoundaryValues::BoundaryValues(MeshBlock *pmb, BoundaryFlag *input_bcs,
           << "<orbital_advection> shboxcoord must be 1 or 2." << std::endl;
       ATHENA_ERROR(msg);
     }
-
     if (pmb->block_size.nx3>1) { // 3D
       if (shearing_box == 2) {
         std::stringstream msg;
@@ -204,15 +204,10 @@ BoundaryValues::BoundaryValues(MeshBlock *pmb, BoundaryFlag *input_bcs,
     }
 
     int level = pmb->loc.level - pmy_mesh_->root_level;
-    // nblx2 is only used for allocating SimpleNeighborBlock arrays; nblx1 for loc_shear
-    // TODO(felker): initialize loc_shear{0, pmy_mesh_->nrbx2*(1L << pmb->loc.level - ..)}
-    // in ctor member initializer list and update as refinement occurs. And nblx2
+    // nblx2 is the only var used in the ctor for alllocating SimpleNeighborBlock arrays
+    nblx1 = pmy_mesh_->nrbx1*(1L << level);
     nblx2 = pmy_mesh_->nrbx2*(1L << level);
-    // is_shear{} in member init_list
-    is_shear[0] = false;
-    is_shear[1] = false;
-    loc_shear[0] = 0;
-    loc_shear[1] = pmy_mesh_->nrbx1*(1L << level) - 1;
+    nblx3 = pmy_mesh_->nrbx3*(1L << level);
 
     if (shearing_box == 1) {
       if (NGHOST+xgh_ > pmb->block_size.nx2) {
@@ -258,18 +253,26 @@ BoundaryValues::~BoundaryValues() {
 //! initializes the shearing block lists
 
 void BoundaryValues::SetupPersistentMPI() {
+  // TODO(KGF): given the fn name, it is confusing that the fn contents are not entirely
+  // wrapped in #ifdef MPI_PARALLEL
   for (auto bvars_it = bvars_main_int.begin(); bvars_it != bvars_main_int.end();
        ++bvars_it) {
     (*bvars_it)->SetupPersistentMPI();
   }
-
   // KGF: begin exclusive shearing-box section in BoundaryValues::SetupPersistentMPI()
+
+  // TODO(KGF): the shearing box metadata is required regardless if MPI is used, so better
+  // to encapsulate it in a separate function
+
   // initialize the shearing block lists
   if (shearing_box != 0) {
     MeshBlock *pmb = pmy_block_;
     int *ranklist = pmy_mesh_->ranklist;
     int *nslist = pmy_mesh_->nslist;
     LogicalLocation *loclist = pmy_mesh_->loclist;
+    // this error needs to be downgraded to a non-fatal warning if the AMR+SMR workaround
+    // to the MHD+shear+refinement restriction is employed (#569)--- or remove entirely,
+    // hoping that RefinementCondition() does not violate below restrictions
     if (pmy_mesh_->adaptive) {
       std::stringstream msg;
       msg << "### FATAL ERROR in BoundaryValues Class" << std::endl
@@ -281,18 +284,34 @@ void BoundaryValues::SetupPersistentMPI() {
           LogicalLocation loc;
           loc.level = pmb->loc.level;
           loc.lx1   = loc_shear[upper];
+          loc.lx2   = pmb->loc.lx2;
+          for (int64_t lx3=0; lx3<nblx3; lx3++) {
+            loc.lx3   = lx3;
+            MeshBlockTree *mbt = pmy_mesh_->tree.FindMeshBlock(loc);
+            // see discussion #569: hydro+shear with SMR (not AMR) works with x1, x3
+            // refinement but MHD with x3 refinement will NOT work without adjusting
+            // corner E-field flux correction to also take into account adjacent x3 block
+            // emf (currently just averages emf values between inner/outermost x1 MBs)
+            if ((mbt == nullptr || mbt->GetGid() == -1) && MAGNETIC_FIELDS_ENABLED) {
+              std::stringstream msg;
+              msg << "### FATAL ERROR in BoundaryValues Class" << std::endl
+                  << "shear_periodic boundaries do NOT work with MHD"
+                  << "if there is refinment in the x3 direction." << std::endl;
+              ATHENA_ERROR(msg);
+            }
+          }
           loc.lx3   = pmb->loc.lx3;
           for (int64_t lx2=0; lx2<nblx2; lx2++) {
             loc.lx2 = lx2;
             MeshBlockTree *mbt = pmy_mesh_->tree.FindMeshBlock(loc);
-            int gid = mbt->GetGid();
-            if (mbt == nullptr || gid == -1) {
+            if (mbt == nullptr || mbt->GetGid() == -1) {
               std::stringstream msg;
               msg << "### FATAL ERROR in BoundaryValues Class" << std::endl
                   << "shear_periodic boundaries do NOT work "
                   << "if there is refinment in the x2 direction." << std::endl;
               ATHENA_ERROR(msg);
             }
+            int gid = mbt->GetGid();
             shbb_[upper][lx2].gid = gid;
             shbb_[upper][lx2].lid = gid - nslist[ranklist[gid]];
             shbb_[upper][lx2].rank = ranklist[gid];
@@ -473,7 +492,6 @@ void BoundaryValues::ApplyPhysicalBoundaries(const Real time, const Real dt,
     ps = pmb->pscalars;
   }
 
-
   NRRadiation *prad = nullptr;
   RadBoundaryVariable *pradbvar = nullptr;
   if (NR_RADIATION_ENABLED || IM_RADIATION_ENABLED) {
@@ -488,8 +506,6 @@ void BoundaryValues::ApplyPhysicalBoundaries(const Real time, const Real dt,
     //pcrbvar = &(pcr->cr_bvar);
   }
 
-
-
   // Apply boundary function on inner-x1 and update W,bcc (if not periodic)
   if (apply_bndry_fn_[BoundaryFace::inner_x1]) {
     DispatchBoundaryFunctions(pmb, pco, time, dt,

src/bvals/bvals.hpp

@@ -183,20 +183,26 @@ class BoundaryValues : public BoundaryBase, //public BoundaryPhysics,
 
   // Shearing box (shared with Field and Hydro)
   // KGF: remove the redundancies in these variables:
-  int shearing_box; // flag for shearing box: 0 = none, 1: xy, 2: xz
-  int joverlap_, joverlap_flux_; // # of cells the shear runs over one block
-  Real ssize_;                  // # of ghost cells in x-z plane
-  Real eps_, eps_flux_;          // fraction part of the shear
+  int shearing_box;               // flag for shearing box: 0 = none, 1: xy, 2: xz
+  int joverlap_, joverlap_flux_;  // # of cells the shear runs over one block
+  Real ssize_;                    // # of ghost cells in x-z plane
+  Real eps_, eps_flux_;           // fraction part of the shear
   Real qomL_;
   int xorder_, xgh_;
-  AthenaArray<Real> pflux_;    // pencil buffer for remapping
+  AthenaArray<Real> pflux_;        // pencil buffer for remapping
+
+
+  std::int64_t nblx1, nblx2, nblx3; // max number of MeshBlocks in each dir
+                                    // on this level of refinement
+
+  //! note: it is possible for a MeshBlock to have is_shear={true, true},
+  //! if it is the only block along entire x1 range of domain
+  bool is_shear[2];                 // inner_x1=0, outer_x1=1
+
+  std::int64_t loc_shear[2];        // x1 LogicalLocation of block(s)
+                                    // on inner/outer shear bndry
 
-  std::int64_t nblx2;
-  //! it is possible for a MeshBlock to have is_shear={true, true},
-  //! if it is the only block along x1
-  bool is_shear[2]; // inner_x1=0, outer_x1=1
   SimpleNeighborBlock *shbb_[2];
-  std::int64_t loc_shear[2];  // x1 LogicalLocation of block(s) on inner/outer shear bndry
 
   // tomo-ono: 3x arrays and 4x arrays are required for int and fc, respectively
   ShearNeighborData<4> sb_data_[2];