One-sided anomaly near polar axis in GRMHD with SMR

[pomelo-akiyama]

Hi all,

I've been running a 3D GRMHD MAD simulation with Athena++ (Kerr-Schild, a = 0.98, tilt_angle = 0, full θ coverage, polar boundary condition, rk2/xorder = 2) and I'm seeing a numerical anomaly. I'd appreciate any thoughts or pointers from the community.

Starting around t ≈ 23,000–24,000, the north hemisphere (z > 0) begins to develop anomalous behavior near the polar axis — density starts climbing in a column along +z, and the Lorentz factor begins rising toward gamma_max. Over the next ~2,000–3,000 time units, this spreads from the axis outward. By t ≈ 26,000, essentially the entire north hemisphere is contaminated: $\rho$ is anomalously high and $\gamma$ is saturating against the gamma_max = 50 ceiling across large areas.

Meanwhile, the south hemisphere stays fine.

The comparison panel below shows $\rho$ (top) and $\gamma$ (bottom) in the x–z plane at three times:

And here's the evolution sampled every Δt = 500:

$\rho$:

$\gamma$:

I'm using refinement = static with multiple Level 2 boxes placed near both poles (r ∈ [30, 1200], θ/π ∈ [~0.005, ~0.17] for the north, symmetric counterpart for the south). The anomaly starts in the region covered by these refinement boxes and spreads outward from there. The refinement layout is symmetric between north and south, but only the north develops the problem.

The following relevant input parameters are used:

Parameter	Value
Problem	gr_torus
Coordinates	Kerr-Schild
Black hole spin	a = 0.98
Tilt angle	0.0°
Domain	r ∈ [1.1, 1200], θ ∈ [0, π], φ ∈ [0, 2π]
Root grid	88 × 32 × 16
Meshblock	22 × 16 × 8
Boundary Conditions	radial outflow, polar polar, azimuthal periodic
Refinement	static, Level 2 near both poles
gamma_max	50
sigma_max	100
beta_min	0.1
Density floor	max(10⁻² r⁻¹·⁵, 10⁻⁶)
Integrator	rk2, xorder = 2

Configure line:

> python configure.py --prob=gr_torus -g -b --coord=kerr-schild --flux=hlle --nghost=2 --cxx icpx -mpi -hdf5

However, this simulation costs computational resources. I am therefore unable to provide a quick reproducible example.

Any suggestions for diagnostics that might help pin down the root cause?

Thanks!

[tomidakn]

In your configuration, the refined grid should be contacting the pole, and I have to admit that Mesh refinement touching the pole is not very well tested. Does it run smoothly without refinement? If so, I appreciate if you could look in to the ghost cell values on the pole. You can either print them in the code, or use VTK output with ghost_zone enabled. If some values (particularly B-fields) there are not appropriate, we have to fix it.

I am not very familiar with GR. And unfortunately, the most knowledgeable person left the field, so we are on our own.

[pomelo-akiyama]

Thanks for the response and the pointers.

I have run the no-refinement version before, but only for a short duration — it seemed fine at early times. I'll extend it to check whether the issue appears at late times without SMR. This will take a while.

In the meantime, for the simulation with SMR, I have .rst files from around the onset of the anomaly, so I'll restart from there with VTK ghost zone output enabled. Will report back with findings.

[pomelo-akiyama]

Hi @tomidakn, thanks for the suggestions. As you requested, I examined the ghost cell values on the pole (with ghost-zone VTK output enabled). I also ran a no-refinement version and, additionally, a global uniform Level-1 version (no SMR level boundaries) to further isolate the SMR factor.

The analysis below was prepared with AI assistance, which I've tried my best to understand and verify. Apologies in advance if anything is off. The raw data files (VTK + athdf) are available here (Google Drive) if you'd like to verify directly.

Ghost-Cell Diagnostics

Short answer: the polar BC implementation itself is fine. Ghost-cell values were compared against the corresponding physical cells in the $\varphi + \pi$ partner block (with the appropriate sign flips for $v^\theta$, $v^\varphi$, $B^\theta$, $B^\varphi$). The $B$-field components match to machine precision (diff = 0.0), and hydro variables agree to $\mathcal{O}(10^{-5})$.

Here's a concrete example at $r \approx 110 r_g$ (block 520 ghost vs. block 664 physical, $\varphi$ offset = $\pi$):

Variable	Ghost cell	Partner physical	Absolute diff
$B^r$	$-1.93194 \times 10^{-3}$	$-1.93194 \times 10^{-3}$	0.0
$B^\theta$	$-3.94532 \times 10^{-5}$	$-3.94532 \times 10^{-5}$	0.0
$B^\varphi$	$+5.82414 \times 10^{-3}$	$+5.82414 \times 10^{-3}$	0.0
$\rho$	$2.29935 \times 10^{-4}$	$2.29944 \times 10^{-4}$	$8.3 \times 10^{-9}$
$p_\text{gas}$	$3.22667 \times 10^{-4}$	$3.22712 \times 10^{-4}$	$4.4 \times 10^{-8}$

However, a more revealing diagnostic is the gradient between the ghost cell and its adjacent physical cell right at the polar axis.

The polar BC fills ghost cells from the $\varphi + \pi$ partner. So ghost[1] at $(\theta \approx -0.7°, \varphi)$ gets its value from physical[2] at $(\theta \approx +0.7°, \varphi + \pi)$. Meanwhile, the actual physical[2] at the local $\varphi$ sits right next to this ghost cell. If the flow is axisymmetric near the pole, these two should be nearly identical. In the MAD turbulent flow, they're not — and the discrepancy is extreme:

$r$ range ($r_g$)	$\rho$ phys/ghost (max)	$\lvert B \rvert$ phys/ghost (max)
$[1, 30)$	$3\times$	$10\times$
$[30, 100)$	$\mathbf{15\times}$	$\mathbf{83\times}$
$[100, 300)$	$\mathbf{215\times}$	$\mathbf{1862\times}$
$[300, 600)$	$\mathbf{577\times}$	$\mathbf{16698\times}$

The PLM reconstruction at the $\theta = 0$ interface sees jumps of $10^2 \sim 10^4$ $\times$ in adjacent cells. This likely drives a positive feedback loop:

1. Turbulent azimuthal asymmetry near the pole → ghost cell (from $\varphi + \pi$) and physical cell (at local $\varphi$) differ
2. PLM reconstruction across this jump overshoots → distorted interface values → bad fluxes
3. Distorted physical values get written into the partner's ghost cells next timestep
4. Gradient grows → repeat

The polar BC is doing exactly what it should — the problem is that in a non-axisymmetric flow, "copy from $\varphi + \pi$" inherently creates $\theta$-direction gradients at the pole that don't correspond to any real physical gradient (since $\theta = 0$ is a coordinate singularity where all $\varphi$ directions converge to a single point). And the PLM reconstructor has no way to know this gradient is artificial.

Additional Runs with Different SMR Configurations

As you suggested, I also ran the simulation without refinement, and additionally with global uniform Level-1 refinement (no level boundaries). Here's what happened:

No refinement (base grid $88 \times 32 \times 16$): No numerical contamination appeared, but the jet physics never really developed either — $\gamma_\text{max} \approx 1.2$ throughout the run, likely because the resolution is too coarse to resolve the relevant processes.

Global uniform Level-1 (effective $176 \times 64 \times 32$, no level boundaries, no prolongation/restriction): Numerical contamination still appeared, but with a different pattern — both hemispheres were affected, and instead of the one-sided pole-outward propagation seen in the original case, it manifested more as block-shaped artifacts.

Videos ($\rho$ left panel, $\gamma$ right panel):

• No refinement:
  MAD98_no-refinement.webm

• Global uniform L1:
  MAD98_global-L1.webm

N/S Asymmetry Profile (Original Case, $t = 25001$)

For completeness, here's the quantitative N/S comparison at the polar axis (first/last physical $\theta$-cell, $\varphi$-averaged):

$r$ range ($r_g$)	$\langle\rho\rangle$ North	$\langle\rho\rangle$ South	$\rho$ N/S	$\langle B^2\rangle$ North	$\langle B^2\rangle$ South	$B^2$ N/S	$\gamma_\text{max}$ N	$\gamma_\text{max}$ S
[1, 10)	8.73e-03	7.87e-03	1.1×	2.29e+01	4.07e+00	5.6×	1.3	1.2
[10, 30)	5.05e-04	4.23e-04	1.2×	1.44e-01	7.49e-02	1.9×	1.4	1.5
[30, 100)	3.06e-04	1.85e-04	1.7×	7.94e-03	1.05e-02	0.8×	11	2.4
[100, 300)	6.38e-05	1.80e-05	3.5×	4.01e-04	9.62e-05	4.2×	16	7.0
[300, 600)	3.51e-04	5.76e-06	61×	3.69e-03	8.92e-06	413×	9.0	4.6
[600, 1200)	1.37e-04	1.92e-06	72×	1.24e-03	1.01e-06	1235×	8.2	3.0

The transition at $r \approx 30 r_g$ lines up with the outer edge of <refinement2> — in the original mesh, that's where the Level-2 polar region transitions to a different AMR structure.

Raw Data

The data archive (Google Drive) contains the snapshot at $t = 25001 r_g/c$ with ghost zones enabled. A note on file numbering: this snapshot was produced by restarting from the original simulation's restart file with a modified input file that enables ghost-zone output. As a result, Athena++ continued each output stream's counter from the restart point, leading to mismatched numbering across formats:

File	Number	Physical time
mad98.block{N}.prim.02501.vtk (× 800)	02501	$t = 25001$
mad98.block{N}.user.02501.vtk (× 800)	02501	$t = 25001$
mad98.prim.00251.athdf	00251	$t = 25001$
mad98.user.00001.athdf	00001	$t = 25001$

The VTK and HDF5 prim counters inherited the original run's numbering, while the HDF5 user output was newly added for this diagnostic run and started from 0. All four correspond to the same physical time.

If you'd like to inspect the polar ghost cells directly, the 128 blocks touching the north polar axis ($\theta = 0$) are:

  0,   1,   4,   5,   8,   9,  12,  13,  32,  33,  36,  37,  40,  41,  44,  45,
 64,  65,  68,  69,  72,  73,  76,  77,  96,  97, 100, 101, 104, 105, 108, 109,
256, 257, 260, 261, 264, 265, 268, 269, 288, 289, 292, 293, 296, 297, 300, 301,
320, 321, 324, 325, 328, 329, 332, 333, 352, 353, 356, 357, 360, 361, 364, 365,
512, 513, 516, 517, 520, 521, 524, 525, 530, 531, 534, 535, 538, 539, 542, 543,
548, 549, 552, 553, 556, 557, 560, 561, 566, 567, 570, 571, 574, 575, 578, 579,
656, 657, 660, 661, 664, 665, 668, 669, 674, 675, 678, 679, 682, 683, 686, 687,
692, 693, 696, 697, 700, 701, 704, 705, 710, 711, 714, 715, 718, 719, 722, 723

These blocks have $\theta \in [-0.049, 0.442]$ rad (2 ghost cells at $\theta < 0$, 18 physical cells). The numbering follows the octree Morton-curve ordering — blocks come in pairs (e.g. 0/1, 4/5, …) that are adjacent in $\varphi$, and the polar BC partner for a given block is the one offset by $\Delta\varphi = \pi$ (e.g. block 260 ↔ block 4).

Would appreciate any thoughts on the cross-pole gradient issue — happy to discuss further.

[tomidakn]

Thank you for the update. I quickly looked at the data, but please understand I am very busy.

I guess the fact it is one-sided is just coincidence (it can happen randomly on either side) and probably not fundamental.

The data shows anomalously different magnetic fields (all the directions) on the pole, which likely means that the polar boundary condition for the magnetic field is not working properly. Yes, it properly copies the data from the other side of the pole, but it does not mean it is correct nor stable.

First of all, the pole is a singular point, and all the faces and edges of the first active cell coincide. The hydrodynamic fluxes through these faces do not matter so much because the surface area is zero and they do not contribute to evolution. However, the EMF_r along the radial edges on the pole (j=js for example) directly appears in the evolution of the magnetic fields in the first active cells along the pole. I think the treatment of this is the root of the issue.

When the polar boundary is in use, we replace the EMF_r with its average over all the edges in the phi-direction (see how flux_north/south_* are processed in flux_correction_fc.cpp), because they are actually the same edge. This is intended to be good for everyone, but it may be bad for everyone, possibly because of the turbulent nature of your problem. If the EMF_r on the pole before averaging are very different, the averaging actually can destroy the smooth distribution.

So, how can we fix it? I do not have a silver bullet, and it can be problem specific. I'm afraid that it takes trials-and-errors. I can suggest, but without any warranty, to try for example...

1. Use the reflecting or outflow boundaries (without averaging). This can break the consistency on the pole, but anyway it is already broken.
2. Use the polar boundary (copy from the other side) but skip averaging. This has the same issue as 1, though.
3. Write your own boundary function based on a sensible assumption.
4. Add (magnetic) diffusion in the phi direction around the pole to keep the polar region smooth.
   Note that you may need to hack the code and/or write your own boundary function to implement these (possible) fixes.

[pomelo-akiyama]

Thank you very much for looking into this. I really appreciate it. Your explanation of how the polar BC works is very helpful. I'll experiment with the approaches and follow up here if I find anything useful.