Add one to `nxi`

[David Dickinson]

I've noticed that nxi is almost always used as nxi+1 everywhere (particularly after @Michael Hardman introduced some fixes in collisions). I'd propose we redefine nxi -> nxi + 1to avoid a lot of the confusion/noise in the code regarding this.

‌

The only time we don't use nxi+1 seems to be when 2*ng2>nxi, however note that we define nxi = max(2*nlambda -1, 2*ng2) and that 2*nlambda >= 2*ng2 by definition so in these cases nxi == 2*ng2 == 2*nlambda. If we were to simply say nxi = 2*nlambda I think we can avoid the need for adding on 1 in most uses whilst still being correct for those places where we don’t currently add on one. I know @Colin Malcolm Roach has looked at this before and noted the importance of including the nxi+1 point in velocity space integrations .

On a similar theme, we also allocate quite a few arrays in collisions in le layout with energy size negrid + 1 but we don’t ever seem to use this extra energy grid point in any calculation (we just set it to 0 and then don’t look at it again). I think we can remove this from the allocations.

I can push a branch with these changes at some point but though I would see if anyone has any arguments to counter this (or in support of this).

[Michael Hardman]

Hi David,

In my opinion it is better to keep the meaning of nxi as the number of unique points on the xi grid. From this perspective nxi + 1 appears entirely naturally as the maximum upper bound on the velocity space integrals. Besides the naming convention, I would note that many users and developers now will be used to the definition of nxi in its current form and peculiarity. Moving to a different, and equally peculiar, definition of nxi will likely lead to further confusion (and bugs!). I do not think it is worth pursuing this change unless there is a clear physics benefit derived from doing so. Finally, unless you're confident that the collision operator is tested to a very high level, I think there are good reasons for not changing low-level code for cosmetic purposes.

The historical reason for why le layout has the size you point out may be due to the need to make gs2 e grid compatible with sfincs or neo in the past. In any case, as I argue above for nxi, I believe this to should be left alone unless there is a positive physics reason to change it.

Best, Michael

[David Dickinson]

Thanks @Michael Hardman that's a fair point that a change to the meaning on nxi may initially be confusing -- I don't think this is always a strong reason for not trying to improve things but I agree this means that there has to be a good reason for the change. The motivation for this proposal isn't to improve physics capability but it isn't really just cosmetic either -- my motivation is to try to improve the maintainability, clarity and consistency of the code in order to help reduce the likelihood of bugs being introduced. It's very easy to forget to manually add the +1 that's required in these loops for example and there's more overhead for new users/develops to understand why we have +1 everywhere for this loop type and no other loop type. I think the approach you adopted in collisions of introducing a different named variable is probably the best balance between maintaining the existing meaning of nxi whilst avoiding the manual +1 everywhere, but we should probably roll it out to the reset of the relevant code (probably just le_grids) so that it is used consistently.

The change for the energy dimension extent is again not motivated by physics (as this point isn’t on our energy grid!) but to reduce time spent by people trying to understand why we have this unusual limit and to reduce the scope for bugs (e.g. accidentally using this extra point in an average or similar). The reason this extra point exists is, I think, to make some of the matrix solves (and similar sweeping type loops) “self-starting”; we have loops like:

ged(negrid+1,ielo) = 0.0
do ie = negrid, 1, -1
  ged(ie,ielo) = (delta(ie) - ec1(ie,ielo)*ged(ie+1,ielo))*ebetaa(ie,ielo)
end do

where for the first iteration ie = negrid and we access ged at ie+1, having just set ged(negrid+1, ielo) = 0 we could equally just do

ged(negrid,ielo) = delta(negrid) * ebetaa(negrid, ielo)
do ie = negrid-1, 1, -1
  ged(ie,ielo) = (delta(ie) - ec1(ie,ielo)*ged(ie+1,ielo))*ebetaa(ie,ielo)
end do

which does away with the need for this extra point that is not on our grid.

[Michael Hardman]

Hi David,

I think that defining a new variable, perhaps called something like nxi_ulim, to be used in the loops as the upper limit is reasonable. I think there is a good reason for why nxi +1 should be understood as the upper limit in integrations. Physically there is only one vpar =0 on the xi grid, whereas we are forced to keep 2 such points on the lambda grid.

Regarding the point about the size of the energy dimension in the le layout, your observation that the 2 methods of starting the tri diagonal solve are equivalent is correct. However I observe that the same "self starting" behaviour is used in the tri-diagonal solve used to invert the pitch angle scattering matrix. Unless you are also proposing to change this loop, then it would be best to keep the "self starting" behaviour in the energy scattering for consistency. Personally I think there is some elegance in a tri-diagonal solve which has non-trivial content on only one line. Changing the size of the le layout array will surely require changes to the redistributes. Do these registries have unit tests already? If not then I think there is a strong conservative reason for keeping things the way they are.

Finally, I will note that as a present developer I greatly value the familiar code. Said another way, whilst changing how the code is written may make it more sustainable according to modern accepted standards of good practice, minor changes like this may have the effect of alienating the researcher developer base which gs2 relies on to survive.

I'm aware that I may sound like some fuddy-duddy who doesn't appreciate the software support which you are providing more generally. That is not all the case. I'm only trying to convey my concern that changes to the base level code in gs2 should only be made when there is a physics reason to do so.

Perhaps we could ask @Michael Barnes @Bill Dorland to comment?

[David Dickinson]

A large part of the motivation for moving the code from sourceforge and adopting new workflows etc. was exactly to ensure this sort of input and discussion about proposals was possible, so I do really welcome your views here – we don’t want to get into a situation where there’s no discussion and changes are made against the will of the community. The more input and discussion the better.

Regarding the energy dimension – yes I would ensure that both dimensions are treated consistently. The redistributes shouldn’t need changing – these points aren’t on the grid so there’s no distribution function information to move about corresponding to these points. There are tests for these but they haven’t been pulled out into proper unit tests yet, they still live in the main code.