Numerical instability at large kx in nonlinear simulations

[David Dickinson]

We occasionally observe sudden rapid growth in phi2_by_ky during nonlinear simulations. Attached is a very simple input file which shows this sort of behaviour (collisionless, electrostatic, adiabatic electrons, 2 ky modes) and can run on a couple of cores. Below is an example of phi2_by_ky for this case with two different nx. We see that the lower nx displays the instability a bit later.

‌

‌

This instability appears to be avoided here with either very large bakdif or somewhat higher (but still reasonable) ntheta, from 8 to 32. Whilst it looks like the zonal phi2 takes off first, it is unclear if it is this mode “going bad” first or not.

This happens with both the split (on grid) and unsplit (bakdif decentered) nonlinear approaches.

‌

Whilst it seems there is an obvious cause and reasonable solution (low ntheta), run failures with very similar characteristics have been observed where a more reasonable ntheta (>= 32) has been used so simply increasing ntheta may not always be a feasible solution.

‌

‌

[David Dickinson]

As discussed in issue #144, we know the homogeneous solution that we add in to ensure consistency at bounce points and periodicity of zonals can be highly oscillatory and essentially grid scale. In that issue it looks like this may cancel out grid scale in the inhomogeneous solution, but clearly this could be a source for grid scale (like we’re seeing here) if the “cancellation” is not exact.

Interestingly forcing the homogeneous part to zero (i.e. g_h = 0) seems to give a simulation which avoids the “blow up” seen here and also gives a zonal phi2 which closely follows the full case with higher ntheta. It is not sufficient to simply avoid calling self_periodic and it is not sufficient to simply zero g_h for the trapped and wfb. This leaves only the linked BC which uses g_h. Switching to zero and zeroing g_h for zonals, trapped and wfb avoids the issue (this is the same as zeroing g_h entirely so no surprise). Zero bc and g_h = 0 for trapped and wfb (i.e. only using for periodicity) leads to the issue, as does Zero bc and g_h = 0 for zonal modes only (i.e. g_h for trapped + wfb still used). This seems to therefore indicate that any use of g_h leads us to the issue.