NaNs for small radius with RKL

[aborderies]

Hi all,

I'm having trouble with the RKL integration method. The setup I'm using consists of a spherical 2D axisymmetric grid, logarithmically spaced in radius. I've just noticed that when I activate non-ideal MHD terms (ambipolar or ohmic scattering) using RKL, with rotation, and start my grid at a sufficiently small radius, NaNs appear in the magnetic field after only a few cycles.

Here is my setup that reproduce that: Collapse_Setup.zip, which consists of a sphere of homogeneous density with solid rotation and a constant magnetic field.

From what I've been able to test in this particular case, the crash appears for the X1-grid going from 1e-4 to 1e4, but not for the grid going from 1e-3 or more to 1e4. Nor does the crash appear if I remove the rotation (parameter beta=0 in the ini file), nor if I use explicit integration for non-ideal effects.
The NaNs appears near the axis and the inner boundary, but changing the boudary conditions doesn't seem to influence the apparition of the NaNs from my tests.

Is this a known problem with the RKL method? And if yes, is there a solution to prevent that ?

Thanks in advance for the help !

[glesur]

Hi @Anto6453 , yes this happens sometimes with the RKL scheme, hard to know why in this precise case, but my experience is that Ambipolar diffusion can sometime render the RKL scheme unstable (as it is the case here). It could also be an issue with your boundary conditions. I'm having a look at your setup, I'll keep you posted.

[aborderies]

Hi @glesur, thanks for your answer and for your time !
According to my tests, omic resistivity alone also causes NaN to appear. Concerning the boundary conditions, I've also tested the one implemented in idefix (reflective and outflow) on the X1-beg boundary and it doesn't seem to have much influence on the appearance of NaNs.
However, one thing that seems to prevent or at least reduce the instability is the cell size in the X2-dir (theta direction) near the axis. I tested a stretched spacing near the axis as follows:
X2-grid 3 0.0 16 s+ 0.5 350 u 2.641592653589793 16 s- 3.141592653589793
which gives the same resolution between $0.5$ rad and $\pi-0.5$ rad as a uniformly spaced 512 grid between $0$ and $\pi$, but with much larger cells around the axis, and the NaNs don't seem to appear (at least not directly) this time.

[glesur]

Ok, I think I found the culprit: an old regularisation procedure of the pole that was fixed in the non-rkl version of the RHS, but not fixed in the RKL. This uses a "small_number" to regularise the area around the pole, that turns out to be large for your choice of grid (with delta r=1e-7!)

The easy fix is just to comment out the #ifdef block there:

idefix/src/rkl/rkl.hpp

Line 766 in 54482d5

#if GEOMETRY != CARTESIAN

(not needed anymore I believe)

A PR will be on its way soon.

[aborderies]

That seems to solve my problem, thanks !