Numerical instability at high theta resolution due to pitch angle collision operator

[David Dickinson]

A number of users have come across an issue arising at high theta resolution for collisional linear simulations. This manifests in either a numerical instability with very high growth rate or values becoming NaN after the first timesteps.

‌

The attached gives an example of this case. With ntheta = 64 the simulation is well behaved, whilst at ntheta = 128 the reported phi2 is NaN at the first output step with t>0.

‌

This issue arises only in cases which include the pitch angle scattering collision operator. It appears to to be due to the trapped pitch angle grid spacing. This influences both the weights wl and the xi grid spacing, the product of which appears in the denominator of the off-diagonal coefficients in the Lorentz matrix (I think this product is probably the relevant quantity determining if this issue arises or not). I expect this is making the matrix ill-conditioned/diagonally non-dominant.

‌

One workaround available in next/master is to set theta_grid_gridgen_knobs:alknob large (e.g. 1e5) rather than the default of 0. This ensures that gridgen pays attention to the pitch angle grid spacing. The downside of this is that it can lead to potentially worse theta grid spacing etc.

An alternative workaround is to limit the number of pitch angle grid points (I think a version of this has been implemented in a branch).

<s>A more robust fix which allows high pitch angle resolution may require a change to the numerical treatment in the existing collision operator.</s> Edit, see comment below regarding jend – whilst the matrix was becoming diagonally non-dominant this was due to incorrect bounds calculations, fixing this seems to avoid the issue in the matrix.

[David Dickinson]

It seems that this issue arises due to incorrect values of jend near at points near the inboard side. This arises due to the way we determine if a particular pitch angle is the one which bounces at this theta grid point:

if( (1 - al(il) * bmag(ig) > -bouncefuzz) .and. (1 - al(il + 1) * bmag(ig) > -bouncefuzz)) then ”Consider index 'il' does not bounce at 'ig'

Here bouncefuzz is an input which defaults to 1e-5. In high resolution cases the change in al between neighbouring il points is quite small. This means that for a pitch angle which does bounce at this ig (i.e. such that 1-al(il)*bmag(ig) = 0 ) the next pitch angle might return a -ve but small magnitude value, e.g. 1-al(il+1)*bmag(ig) = -1e-6 . This would lead to il + 1 being considered as the bouncing pitch angle index rather than il. This then leads to issues in the collision operator as it tries to use pitch angle weights that are zero (due to being outside the allowed range of xi).

‌

I think an appropriate fix is to reduce the default of bouncefuzz to something much smaller, say epsilon(0.0) , which should avoid these issues whilst still allowing for some small floating round off in the calculation 1 - al(il) * bmag(ig) returning slightly less than zero (which I think is the purpose of bouncefuzz /= 0).

[David Dickinson]

See PR #313 for a proposed fix.

[David Dickinson]

Resolved in next through PR #313