Internal wfb-like bounce point treatment issues/inconsistencies

In PR #698 we introduced the use of is_lower_bounce_point(ig, il) and similar methods to detect bounce points instead of the previous approaches built off of forbid(ig,il) and forbid(ig+1, il). For example this allows

if (forbid(ig, il) .and. .not. forbid(ig+1, il) to be replaced with if is_lower_bounce_point(ig + 1, il)

This is intended to both improve code clarity and prepares the ground for treating the wfb as a trapped particle.

The is_lower_bounce_point method is built of two checks:

Is the point considered a bounce point – here we check that 1-al(il) * bmag(ig) is sufficiently close to zero.
We check if we are allowed to visit the next theta grid point to the right, i.e. .not. forbid(ig+1, il).

‌

This works well for simple geometries. In geometries with local bmag minima/maxima we end up with internal wfb-like points. By this I mean a place where 1 - al(il) * bmag(ig) approaches zero, but where the adjacent theta grid points are accessible. As it stands this means that is_lower_bounce_point(ig + 1, il) can return true whilst (forbid(ig, il) .and. .not. forbid(ig+1, il) would evaluate to .false. As a result, since PR #698 geometries with local maxima can return different results (and may be subject to numerical instability).

A simple fix for this is to replace step 2 of is_lower_point_method with the inverse logic – test if the point to the left is forbidden. This will reproduce the behaviour prior to PR #698, but means that we essentially do not decouple either side of the internal wfb bounce point. I think our code is not able to handle this decoupling currently and some of this is discussed in issue #148.

‌

As such I’d propose that we adopt the simple fix for now whilst we work on improving our ability to handle these internal wfb points.

‌

Comments (2)