MHD shearing box with SMR

[jonosquire]

I have been trying to run an MHD shearing box with static mesh refinement and been running into issues. I add SMR in the vertical direction only for the usual shearing domain with shear_periodic boundary conditions in x. After a short time (~0.5 orbits), the integrator goes unstable with the timestep dropping to zero because of a very large radial magnetic field that develops in the 2-3 cells around the refinement location at the radial boundaries. Just 1 level of refinement causes the issue, but the blowup/instability occurs faster if higher levels are added.

In order to diagnose the issue I have tried (in various combinations): the vl2, rk2, and rk3 time integrators with various conservative CFLs; order 2, 2c and 3 reconstruction; different initial conditions (from the strat.cpp pgen) and different field strengths; different orbital advection schemes (including no orbital advection); setting qshear=0 (still with the shear_periodic BCs); different mesh refinement locations, meshblock configurations and resolutions; different sound speeds and density floors. The problem persists with all of these. Turning off the magnetic field (pure hydro shearing box), or using periodic (rather than shear_periodic) boundary conditions in x1 removes the issue, and the grid setup appears to be as expected with all of the meshblocks at the same z facing the x1 shear-periodic boundaries having the same level.

Has anybody come across similar problems? I'm at a loss for what to try at this point – is there an issue with with shear_periodic BCs and SMR? I've attached an input file that should reproduce the issue with the strat.cpp pgen in the current athena version (this is effectively the default input file with different units, a finer grid, and the SMR added). I've also attached a screenshot showing the radial field blowup at the boundary (the colormap is saturated; the points right on the boundary at z=+-2.5 continue growing to infinity).

Thanks in advance for the help.
athinput.strat.txt

[tomidakn]

@tomo-ono is the best person to ask.

[tomo-ono]

I think that this is not surprising. At the shear-periodic boundaries, we have two treatments. One of these is the flux correction where we adjust values belonging to cells facing on the x-boundaries. For cell-centered values like hydrodynamic quantities, it is OK if we calculate average between values in the innermost meshblock and outermost meshblock. However, it is required to take into account the next meshblock in the z direction when we adjust values of the electric field at meshblock corners. Without refinement, this treatment would be skipped. On the other hand, it is required to update this process for your runs.

[jonosquire]

Thanks for the thoughts. So, you're saying that in order to have this type of setup work stably, it requires changing how the electric field is calculated between the two refinement levels at the shear_periodic boundaries? Do you know of any potential simple workarounds to try, as this sounds rather involved. Or, perhaps this could be listed as a known issue in the documentation?

[tomo-ono]

We have not supported the refinement in shearing box officially yet. I guess that it is required to give changes to "athena/src/bvals/fc/flux_correction_fc.cpp". However, it could be not very easy.

[jonosquire]

Ok, thanks, I'll take a look. Should I flag this as an Issue?

[c-white]

I used the automatic tool to elevate this to an issue (#571), given that I think there are many users out there who are using or who want to use refinement with their MHD shearing boxes.

[trwaters]

It's currently possible to run a nested SMR grid in the shearing box, as long as the meshblock size is chosen so that the shearing boundary is not refined. To do this, we had to define an AMR criterion that refines by the desired number of levels within a given radius. We run this for 0 cycles with refinement=adaptive to let the code automatically determine the meshblock structure. Then we restart the simulation with refinement=static.

[felker]

clever workaround! seems like this could be integrated in Athena++ as an automatic feature as a temporary patch until flux_correction_fc.cpp is fully generalized

[trwaters]

Credit goes to Adam Dempsey, who used the shearing box with SMR to study accretion onto a binary black hole embedded in an AGN disk
https://ui.adsabs.harvard.edu/abs/2022ApJ...940..155D/abstract

To use the pure SMR setup he employed (see Sec 3.1), there shouldn't be anything to patch. For the AMR+SMR approach I describe above, you'd need to change this check in src/bvals/bvals.cpp from an ATHENA_ERROR to a warning:

 273     if (pmy_mesh_->adaptive) { 
 274       std::stringstream msg; 
 275       msg << "### FATAL ERROR in BoundaryValues Class" << std::endl 
 276           << "shear_periodic boundaries do NOT work with AMR at present." << std::endl; 
 277       ATHENA_ERROR(msg); 
 278     } else {

[felker]

Thanks; it is confusing that this entire shearing box-related neighbor metadata setup is in a function named BoundaryValues::SetupPersistentMPI() and isnt wrapped in #ifdef MPI_PARALLEL like all the related variable functions e.g.

athena/src/bvals/cc/bvals_cc.cpp

Lines 614 to 615 in 574a21d

	void CellCenteredBoundaryVariable::SetupPersistentMPI() {
	#ifdef MPI_PARALLEL