Issues with gridgen

[David Dickinson]

For eik and s-alpha geometry we pass our initial B(theta) into gridgen in order to “optimise” the placement of the theta grid points.

I believe this optimisation is intended to:

1. Ensure we have theta grid points on all local B extrema.
2. We add sets of theta grid points until we have the requested total number of grid points. These sets comprise all theta grid points with the same value of B.

In practice there are some challenges with this.

In finding the local extrema we have the following issues / constraints

1. We assume an extrema at the first theta grid point (usually theta = -pi).
2. For extrema found close to, but not exactly on, a grid point we move the extrema to sit on the grid point. This can lead to duplicate appearances of these extrema in the initial set of points considered.
3. We only detect extrema between grid points by observing a change in sign of B' across two grid points. This means we do not detect double local extrema potentially associated with inflection points.

In adding additional theta points we have the following issues / constraints

1. We currently add all input theta grid points and additionally insert npadd additional points evenly spaced between input grid points. This means we usually end up with a set of starting theta grid points which give duplicate B(theta) values – (e.g. for symmetric B we might have +theta and -theta). This ends up leading to multiple “sets” which contain the same theta grid points and correspond to the same B value.
2. We actually “overpopulate” our theta grid points and then choose to remove sets of points based on those which have the worst “resolution metric”. When we have duplicate sets, or uniformly spaced points, we can find sets with identical resolution metrics and the choice of which set to remove can be machine dependent. This is something Michael Hardman spent some time looking at and led to the implementation of some more careful logic to provide a deterministic choice.
3. Sometimes sets that we remove will take us below the target ntheta value and we do not avoid this. For example suppose we currently have ntheta_target + 2 points available and the set of theta points with the worst resolution metric contains three theta grid points. We will then end up returning a theta grid with ntheta_target - 1 grid points and the code will likely fail.

‌

Possible fixes:

1. We can remove the assumption of an extrema on the first grid point.
2. We can avoid duplicating extrema with some additional logic
3. Detecting inflection points with two local extrema in the spline can be challenging but it might be possible through sampling the spline derivative on a fine mesh between neighbouring grid points and looking for a change in sign.
4. To avoid the duplicate non-extrema sets we can probably just check if B(theta) has already been seen and if so don’t add theta to the set of starting points.
5. Instead of oversampling the passed grid, we could perhaps pass in a geometrical output with higher ntheta. For example with ntheta_geometry > ntheta we perform the eikcoefs calculation with ntheta_geometry points, down-sample to ntheta and then pass to gridgen, which oversamples the input. Instead we could perform the eikcoefs calculation with ntheta_geometry points and pass this output to gridgen which then just decides which points to exclude.
6. We could ignore the resolution metric from sets with too many points when considering which sets to remove.

‌

Interestingly removing the duplicated points (fixes 2 and 4 in the above) can lead to a different choice of theta grid points even in a circular Miller case:

This suggests that implementing any of these fixes is likely to introduce challenges for reproducing earlier results.

‌

A possible “cleaner” approach might be to:

1. Identify all local extrema in B(theta) – keep the values of B and how many extrema found and how many theta values
2. Sort the B values and calculate the spacing between neighbouring B
3. Choose a new B subdividing the B interval with the largest spacing, find the values of theta with this B and increment the number of grid points. If the number of grid points is less than the target go to 2, otherwise continue
4. If the number of grid points is larger than the target identify a B value with an appropriate number to remove and remove it
5. return

In other words here we identify B as the quantity we wish to resolve as uniformly as possible and use this directly, rather than the current approach which is a little more theta grid focused.

[David Dickinson]

See PR #944 for some more discussion and some code changes to avoid duplicates etc.

‌