From 420398ff90b5dfc5e2eb2fe5eaf4d2ae48f0b831 Mon Sep 17 00:00:00 2001 From: Kyle Gerard Felker Date: Tue, 7 May 2024 02:29:40 -0500 Subject: [PATCH 1/6] Add exclusion for SMR+MHD+shearing box --- src/bvals/bvals.cpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/bvals/bvals.cpp b/src/bvals/bvals.cpp index 185f68e931..66dc4a090c 100644 --- a/src/bvals/bvals.cpp +++ b/src/bvals/bvals.cpp @@ -122,6 +122,13 @@ BoundaryValues::BoundaryValues(MeshBlock *pmb, BoundaryFlag *input_bcs, // set parameters for shearing box bc and allocate buffers shearing_box = 0; if (pmy_mesh_->shear_periodic) { + // see discussion #569 + if (pmy_mesh_->multilevel && MAGNETIC_FIELDS_ENABLED) { + std::stringstream msg; + msg << "### FATAL ERROR in BoundaryValues Class" << std::endl + << "SMR is incompatible with shearing box with MHD" << std::endl; + ATHENA_ERROR(msg); + } // It is required to know the reconstruction scheme before pmb->precon is defined. // If a higher-order scheme than PPM is implemented, xgh_ must be larger than 2. std::string input_recon = pin->GetOrAddString("time", "xorder", "2"); @@ -145,7 +152,6 @@ BoundaryValues::BoundaryValues(MeshBlock *pmb, BoundaryFlag *input_bcs, << " 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; @@ -264,12 +270,18 @@ void BoundaryValues::SetupPersistentMPI() { } // KGF: begin exclusive shearing-box section in BoundaryValues::SetupPersistentMPI() + // TODO(KGF): it is confusing that it is not wrapped in #ifdef MPI_PARALLEL + // The shearing box metadata is required regardless if MPI is used, so might be better + // encapsulated 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. See discussion #569 if (pmy_mesh_->adaptive) { std::stringstream msg; msg << "### FATAL ERROR in BoundaryValues Class" << std::endl @@ -473,7 +485,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 +499,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, From aba3865c052c7dc362fb2eae048380e74df20986 Mon Sep 17 00:00:00 2001 From: Kyle Gerard Felker Date: Tue, 4 Jun 2024 10:25:22 -0500 Subject: [PATCH 2/6] Cleanup comments --- src/bvals/bvals.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/bvals/bvals.hpp b/src/bvals/bvals.hpp index 0f1b0ceade..f5a19268ed 100644 --- a/src/bvals/bvals.hpp +++ b/src/bvals/bvals.hpp @@ -183,17 +183,17 @@ 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 pflux_; // pencil buffer for remapping std::int64_t nblx2; - //! it is possible for a MeshBlock to have is_shear={true, true}, - //! if it is the only block along x1 + //! 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 SimpleNeighborBlock *shbb_[2]; std::int64_t loc_shear[2]; // x1 LogicalLocation of block(s) on inner/outer shear bndry From 5555b9a6a951ef7458a4119477404c2b9fc612b0 Mon Sep 17 00:00:00 2001 From: Kyle Gerard Felker Date: Tue, 4 Jun 2024 10:31:49 -0500 Subject: [PATCH 3/6] Avoid potential nullptr deref --- src/bvals/bvals.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/bvals/bvals.cpp b/src/bvals/bvals.cpp index 66dc4a090c..142f2b70cc 100644 --- a/src/bvals/bvals.cpp +++ b/src/bvals/bvals.cpp @@ -297,14 +297,14 @@ void BoundaryValues::SetupPersistentMPI() { for (int64_t lx2=0; lx2tree.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]; From 087df1f6264403b16da60026e2a0060bd2211cf2 Mon Sep 17 00:00:00 2001 From: Kyle Gerard Felker Date: Tue, 4 Jun 2024 13:35:57 -0500 Subject: [PATCH 4/6] Move shearing box variable init to constructor init_list --- src/bvals/bvals.hpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/bvals/bvals.hpp b/src/bvals/bvals.hpp index f5a19268ed..f61343aa54 100644 --- a/src/bvals/bvals.hpp +++ b/src/bvals/bvals.hpp @@ -189,14 +189,20 @@ class BoundaryValues : public BoundaryBase, //public BoundaryPhysics, Real eps_, eps_flux_; // fraction part of the shear Real qomL_; int xorder_, xgh_; - AthenaArray pflux_; // pencil buffer for remapping + AthenaArray pflux_; // pencil buffer for remapping + + + std::int64_t nblx1, nblx2, nblx3; // max number of MeshBlocks in each dir + // on this level of refinement - std::int64_t nblx2; //! 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 + 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 + 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]; From 811ea3c6aad7917506e01f4c57e58b5d55c3de76 Mon Sep 17 00:00:00 2001 From: Kyle Gerard Felker Date: Tue, 4 Jun 2024 13:37:17 -0500 Subject: [PATCH 5/6] Replace simple runtime check for SMR+shear+MHD with nuanced check --- src/bvals/bvals.cpp | 47 ++++++++++++++++++++++++++------------------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/src/bvals/bvals.cpp b/src/bvals/bvals.cpp index 142f2b70cc..85070bf91c 100644 --- a/src/bvals/bvals.cpp +++ b/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++) { @@ -122,13 +123,6 @@ BoundaryValues::BoundaryValues(MeshBlock *pmb, BoundaryFlag *input_bcs, // set parameters for shearing box bc and allocate buffers shearing_box = 0; if (pmy_mesh_->shear_periodic) { - // see discussion #569 - if (pmy_mesh_->multilevel && MAGNETIC_FIELDS_ENABLED) { - std::stringstream msg; - msg << "### FATAL ERROR in BoundaryValues Class" << std::endl - << "SMR is incompatible with shearing box with MHD" << std::endl; - ATHENA_ERROR(msg); - } // It is required to know the reconstruction scheme before pmb->precon is defined. // If a higher-order scheme than PPM is implemented, xgh_ must be larger than 2. std::string input_recon = pin->GetOrAddString("time", "xorder", "2"); @@ -210,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) { @@ -264,15 +253,16 @@ 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): it is confusing that it is not wrapped in #ifdef MPI_PARALLEL - // The shearing box metadata is required regardless if MPI is used, so might be better - // encapsulated in a separate function + + // 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) { @@ -281,7 +271,8 @@ void BoundaryValues::SetupPersistentMPI() { 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. See discussion #569 + // 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 @@ -293,6 +284,22 @@ 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; lx3tree.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 Date: Thu, 6 Jun 2024 10:51:07 -0500 Subject: [PATCH 6/6] Fix operator precedence --- src/bvals/bvals.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bvals/bvals.cpp b/src/bvals/bvals.cpp index 85070bf91c..6433ebd9da 100644 --- a/src/bvals/bvals.cpp +++ b/src/bvals/bvals.cpp @@ -60,7 +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)}, + 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++) {