Incorrect handling of two-sided bounce points (wfb and internal wfb-like)

[David Dickinson]

The numerical instability @Michael Hardman reports in issue #147 has been tracked down to inconsistent handling of the duplicate v||=0 points for the wfb in the parallel solve (invert_rhs_1) and the Lorentz collision operator.

Our current treatment of the wfb does not ensure that the distribution function is single valued for both signs of parallel velocity at the points where v|| → 0 (i.e. when bmag = bmax). The Lorentz operator enforces that the distribution function is single valued for the v|| = 0 point after the application of collisions. This inconsistency can result in numerical instabilities such as reported in issue #147 where we have demonstrated that applying the Lorentz collision operator with zero collision frequency is sufficient to trigger the instability. @Michael Hardman has demonstrated that excluding the wfb from the calculation removes the instability.

This same issue can also arise for internal wfb-like particles. These arise when the magnetic field has internal maxima with bmag < bmax. Again in issue #147 we have demonstrated that this scenario leads to a similar numerical instability.

Our current handling of trapped particles only identifies bounce points where there is a grid point corresponding to the forbidden region so we do not see the “internal” bounce points of these wfb modes. Going beyond this, our algorithm for ensuring a single valued distribution function at the bounce point currently requires a forbidden grid point to the left of lower bounce points (and possibly requires this at upper bounce points as well).

These are issues which could be addressed by changing how we detect bounce points and some modification to our algorithm (less elegant, more logic etc.). The primary challenge with addressing this is that these internal bounce points are both lower and upper bounce points and the distribution function for these two bounce orbits need not be single valued. We therefore would need to change the data-structures to store this extra data. To help make this clearer consider an internal bounce point – at a lower bounce point we would usually store g(bounce_right, sigma=1), g(bounce_right, sigma=2). At an upper bounce point we store g(bounce_left, sigma=1), g(bounce_left, sigma=2). Hence for an internal bounce point we would expect to have to store four values but we only have two storage elements typically as bounce_right==bounce_left.

‌

One way to avoid a change in data structure is to note that g(bounce_right, sigma=1) == g(bounce_right, sigma=2) and g(bounce_left, sigma=1)== g(bounce_left, sigma=2), so whilst we have four values there are only two unique values. This could allow us to choose to avoid storing the duplicated v||=0 data and instead use sigma at the bounce point to determine which side of the bounce point we are storing. For internal wfb-like bounce points (and the wfb) this would allow us to store the unique values from either side of the bounce point. For regular bounce points we would find that one sigma is unchanged whilst the other sigma is zero. There would be no change to passing particles.

This change could impact on many areas of the code so is not to be taken lightly but should allow us to consistently treat regular and wfb-like trapped particles and hopefully avoid bugs associated with the duplicate v||=0 point.

An alternative fix is to ensure that we also avoid including wfb-like pitch angles within our grid. This would break the condition that there is a pitch angle which bounces at each theta grid point, but this should be consistently handled throughout most of the code I think. We would need to address places where there are wfb specific logic (checks like il == ng2+1 and similar).

[David Dickinson]

The parallel solve algorithm would undoubtedly become more complicated than it currently is if we opt for the approach of storing just in one sigma to handle the case of the internal bounce points. Before we add in the homogeneous solution to ensure the distribution function is single valued for the two signs of v||=0 on one side of the bounce point we do have an additional value which we need to keep track of (so three values total). One approach may be to treat each bounce orbit in turn for both v|| rather than each sign of v|| in turn for all theta.

[Colin Malcolm Roach]

Thanks @David Dickinson for our interesting discussion of this last week, and for this post. Looks like we have a lot of progress here and that we may be getting close to cracking this old chestnut! Here are some of my comments/thoughts on the problem.

Classes of Trapped Particle

GS2 needs to treat more carefully at least three different classes of trapped particle (which it handles now), and generalise this to possible equilibria with non-monotonic B(theta).

• wfb: pitch angles that bounce at any gridpoint that is either a global or local maximum in B on the theta grid
• normal: pitch angles that bounce at any gridpoint that is neither a local/global maximum or minimum
• ttp: pitch angles that bounce at any gridpoint that is either a global or local minimum in B on the theta grid

Velocity Space Weights/Cells in the Trapped Region

More care is required to define the vspace region represented by each class of trapped lambda gridpoint, to ensure that all particles in the cell are consistent with the imposed parallel boundary condition. In particular the associated vspace region should ONLY INCLUDE trapped particles that bounce close to (within dtheta) of the associated bouncepoint! (Otherwise the parallel BC is inappropriate for some of the particles!)

For simplicity, lets start by considering a uniformly spaced lambda grid, 0 < lambda < 1/Bmin, where 1/Bmax < lambda < 1/Bmin is the trapped region, and dlambda is the grid spacing in lambda. Imposing a self-consistent parallel BC in each region requires the following constraints on the vspace weights:

• wfb: lambda_wfb < lambda <lambda_wfb+dlambda/2 : this ensures NO PARTICLES are included with FINITE v_|| at the gridpoint that will NOT BOUNCE within dtheta of theta_wfb.

  • This is why lambda_wfb CANNOT represent both trapped and passing particles.
  • vspace weight is ZERO for a wfb at its bounce point

• normal: lambda_wfb-dlambda/2 < lambda < lambda_wfb+dlambda/2 : this includes trapped particles that bounce withing dtheta/2 on EITHER SIDE of theta_bounce.

• ttp: lambda_ttp - dlambda/2 < lambda <lambda_ttp: ensures we ONLY include particles that can reach theta_ttp.

In nonmonotonic B, the presence of local maxima and minima will require further careful adjustments to vspace weights for neighbouring pitchangles: lambda=lambda_wfb+dlambda, and lamba=lambda_ttp-dlambda

Duplicated Points with Non-Unique Values in g at wfb Bounce Point and Memory Considerations

@David Dickinson commented on the possible need to store 4 values for g at a wfb bouncepoint: g(bounce_right, sigma=1) == g(bounce_right, sigma=2) and g(bounce_left, sigma=1)== g(bounce_left, sigma=2) in order to evolve the distribution function. sigma is obviously meaningless at a bouncepoint, and so only two numbers numbers must be stored. @David Dickinson has made an excellent suggestion to use both sigmas to store g from the right and left bounce orbits!

These two values should indeed be different because they represent different trapped regions of r-vspace : particles that bounce on different sides of the gridpoint, NONE of which REACH theta_wfb! Non-uniqueness here IS NO ISSUE, because the vspace weight for the wfb at its gridpoint should be ZERO. This means that wfb does not contribute to vspace integrals or fields at theta_wfb!

Excluding wfb from the Grid

@David Dickinson suggested removing wfb gridpoint. While it is certainly possible to exclude the wfb pitch angle gridpoint from the grid, it is not physical to exclude particles in vspace region lambda_wfb < lambda < lambda_wfb+dlambda/2. One could accommodate this vspace in the weight corresponding to lambda=lambda_wfb+dlambda, and sacrifice accuracy slightly by imposoing that these particles bounce at theta_max-dtheta. I would prefer to keep the lambda_wfb gridpoint, and just implement the algorithm and vspace weights correctly, now we know how it can be done properly!

[David Dickinson]

@Colin Malcolm Roach thanks for your comments. I think that the velocity weights are handled consistently currently and we have some velocity integration tests that at least suggest the weights have some of the correct basic properties (i.e. we can integrate known functions to a very good degree of accuracy provided we have high enough resolution). I think we also need to consider how the weights are used within the collision operator as this provides a further complication. For example we have new_trap_int which controls how we calculate the pitch angle weights in the trapped region. These can result in very different weights for individual pitch angles and whilst they both lead to correct integration it is not clear to me that they both lead to sensible weights for use in the collision operator (although I haven’t thought about this much until recently).