Handling of duplicate v|| = 0 points in Lorentz collisions

[David Dickinson]

Due to our use of pitch angle and sign of v|| (sigma) GS2 ends up with duplicate mesh points corresponding to the v||=0 point. For numerical convenience / symmetry we currently duplicate the distribution function at this point, holding the value of g for both v_{||} = 0_+ and v_{||} = 0_-. The algorithm in invert_rhs should ensure that these duplicate points always have a unique value (ignoring the WFB for now, where this isn’t true).

Non-wfb considerations

In velocity space integrations the pitch angle weight associated with a bouncing pitch angle only represents half of the v|| volume associated with the bounce point (e.g. the positive or negative side of the v|| domain). This works out nicely for our duplicate scheme where we end up summing over sigma at this point, g(ig_bounce, sigma = 1, il) * wl(ig_bounce, il) + g(ig_bounce, sigma = 2, il) * wl(ig_bounce, il), Given g(ig_bounce, sigma = 1, il)==g(ig_bounce, sigma = 2, il) for a bouncing particle this is equivalent to 2 * g(ig_bounce, sigma = 1, il) * wl(ig_bounce, il) – in other words this accounts for the full v|| volume associated with the bounce point.

As the conservative scheme in the Lorentz collision operator uses the pitch angle weights but we don’t include both g(ig_bounce, sigma = {1,2}) values it is important that we double the effective weight at this point in order to truly reflect the volume associated with the bounce point. @Michael Hardman introduced this fix in commit 5e0a21d2b

As we don’t include both g(ig_bounce, sigma = {1,2}) we must choose how we select the value of g to use at the bounce point in the collision operator. We could either pick a particular sigma or we could average the values. This choice is influenced by the input vpar_zero_mean, which defaults to true meaning we take the mean. As the invert_rhs algorithm is expected to guarantee that the values are identical for either sigma (within ~floating precision) one might expect that either option for vpar_zero_mean should lead to the same results and hence we could consider this to be redundant for non-wfb pitch angles.

‌

WFB considerations

The WFB (the pitch angle bouncing at B = Bmax) is slightly more complicated for (at least) two reasons:

1. There is not yet consistent treatment of this pitch angle as either passing or trapped (Note: I think of the wfb as trapped and this might influence how I view the issues outlined here). Whilst some steps in this direction have been made there is outstanding work left in order to ensure that is treated consistently as a trapped particle.
2. As a trapped particle the pitch angle weight here can be exactly zero.

Considering the WFB as a passing particle it is clear that we can have g(ig_vpar0, sigma = 1, il) /= g(ig_vpar0, sigma = 2, il) and hence it is not clear how we should select g to use for the vpar = 0 point. Currently this is dealt with in the same way as trapped particles – we either just take one sigma or take the mean. Neither of these appear entirely appropriate for a passing particle.

Considering the WFB as a trapped particle we should be able to take the same approach as regular trapped particles to select the value of g to use for the v|| = 0 point [Note: Currently we do not correctly enforce the bounce condition for trapped WFB at internal boundaries, i.e. for nperiod > 1 the WFB is not correctly bouncing]. However, the fact that the pitch angle weight is zero remains a challenge for a trapped WFB as the conservative differencing operator (eq. 42 of Barnes et al) includes a 1/weight factor. This had been treated by removing the WFB from the Lorentz system at its bounce point and then setting the output value of g_WFB by a linear interpolation between the neighbouring passing pitch angles (or more recently replacing with the pre-collision value). This is the behaviour for special_wfb_lorentz = .true. (not the default). More recently the default of special_wfb_lorentz was changed and instead special handing was added to the base scheme where the pitch angle weight was replaced for the WFB with twice the v|| spacing between 0 and the nearest passing pitch. By manipulating eq. 42 of Barnes et al at j = wfb and making use of the symmetries of the coefficients at this point (G_{j+1/2} == G_{j-1/2} , x_j == 0 and x_{j+1} = - x_{j-1}) one can show that this reduces to [h_{j+1} + h_{j-1} - 2 h_j] / w_j. If w_j tends to zero the only way for this to remain consistent is if the numerator also goes to zero. One might recognise that the term in square brackets shows the central difference approximation of the second derivative of h evaluated at x_j. In other words this scheme may require the second derivative of h to be zero across the v_|| = 0 point for the WFB. We don’t currently enforce this. One way to achieve this would be to set g_wfb = 0.5 * (g_least_passing(sigma = 1) + g_least_passing(sigma=2)) prior to the collisional solve.

‌