Internal grid in Miller geometry (leq) is one-point larger than [-pi, pi]

[Peter Hill]

See here: https://bitbucket.org/gyrokinetics/gs2/src/6fbd2a6bc6ac5e5be0e0737f1860a0e0f78026ed/geo/leq.f90#lines-118:124

and here: https://bitbucket.org/gyrokinetics/gs2/src/6fbd2a6bc6ac5e5be0e0737f1860a0e0f78026ed/geo/leq.f90#lines-339

The grids are [-nt:nt] which means they have 2*nt + 1 points, and the grid spacing should be 2*pi / (2*nt + 1 - 1) == pi/nt to be the range [-pi, pi]. As it is, the actual range is [-(pi + dt), (pi + dt)] with dt = pi/(nt - 1)

I’m not sure why there’s this extra point – possibly to avoid using periodic splines?

At any rate, “fixing” it breaks some regression tests, unsurprisingly.

[David Dickinson]

See https://sourceforge.net/p/gyrokinetics/bugs/68/ for some historical discussion of this.

[Joseph Parker]

This leads to segfaults for me when using low values of ntheta (see the attached input file). We think the loop should read:

    do j=-nt,nt
       do i=1,nr
          r = surf%r + dr(i)

          ! this
          t = (j)*pi/real(nt)
          ! not
          ! t = (j)*pi/real(nt-1)

          R_psi(i,j) = Rpos(r, t)
          Z_psi(i,j) = Zpos(r, t)
          ...
       enddo
    enddo   

    ...
    ! calculate derivatives of R_psi, Z_psi
    call derm(R_psi,  drm,  'E')
    call derm(Z_psi,  dzm,  'O')

    ...
    ! grad(psi) in cartesian form
    call eqdcart(dpm, dpcart)

so that the endpoints of t are -pi and pi. With my input file, using bad values of t result in derivatives of R_psi and Z-psi being simultaneously zero, causing division by zero in the subroutine eqdcart. The problem is apparently fixed by dividing by nt instead of nt-1… but as Peter says, this breaks the regression tests.

[David Dickinson]

Just to expand on the issue that Joseph has seen -- in thederm routine we calculate the derivatives in r and theta of the passed data. There’s special handling for the derivatives at the end of the theta grid which assumes there’s some form of periodicity (even, odd etc.) that allows us to use this symmetry to calculate the derivative. If the end of the theta grid isn’t at the +/- pi point this special handling is unlikely to be correct (and indeed in Joseph’s case seems to lead to a point where the derivative of both R and Z wrt theta are essentially 0).

[Joseph Parker]

Another quick note: while the regression test fails with the proposed change, the test that Miller and CHEASE specifications of the geometry give the same result still passes. (The tolerances are pretty large in this test, but the agreement doesn’t get any worse.)

Pinging @Colin Malcolm Roach and @Bill Dorland

[Colin Malcolm Roach]

In leqinit: theta(-nt) =-pi-D, theta(nt)=pi+D I guess this change was introduced by Justin Ball to handle U-D asymmetries, but I am not sure.

These endpoints in the theta grid seem to conflict with assumptions made in derm for the special handling of theta derivatives at theta=+-pi.

dermassumes that the pi boundaries are at -nt, nt, and uses special conditions to evaluate:dfm(:,-nt,2), and dfm(:,-nt,2)

At best the conflict between these assumptions in derm and leqinit may not matter (and just retain redundant code that should be removed) but at worst the conflict may introduce a bug.

Can we document in the code exactly what subroutine derm is doing, and resolve whether this apparent conflict actually matters

[David Dickinson]

@Colin Malcolm Roach I think Joseph's example (attached at the top of the issue) does demonstrate a bug i.e. wrong behaviour for valid input. Peter's already had a go a documenting what derm seems to be trying to do, but this is derived from trying to interpret intent from whatever code is there rather than being based on any other existing documentation.

The centred difference derivative is essentially df_i = f_{i+1}-f_{i-1} when we're at the end of the grid i=nt and we can't access f_{i+1} directly from an array of data, so derm uses symmetry to determine what value f_{i+1} _should_ have.

For char = 'o' (odd?) derm effectively says f_{i+1} = -f_{i-1} and hence df_nt = -2*f_{i-1}.

For char = 'e' (even?) derm effectively says f_{i+1} = f_{i-1} and hence df_nt = f_{i-1} - f_{i-1} = 0

If the end of the grid i=nt isn't at the expected symmetry point theta=pi then this symmetry doesn't hold at this point and hence the derivatives calculated using this are incorrect (consider cos or sin evaluated at theta = pi + pi/2 the symmetry about this point is very different than at theta = pi).

This may have had minor impact on most realistic runs as the shift in theta(nt) from pi will reduce as ntheta increases but it’s really apparent in some of the reduced resolution cases used for testing.

[Colin Malcolm Roach]

@David Dickinson @Joseph Parker @Peter Hill

• I agree that the endpoint theta derivatives in derm are garbage, as leq_initmoved the endpoint so it is NO LONGER +-pi. Nevertheless this still may not matter… are endpoints derivatives used later in the code?
• theta derivatives evaluated using centered differences for all grid points in range [-nt+1:nt-1]. These will evaluate the derivatives at +- pi, which corresponds to theta(-nt+1), theta(nt-1)
• Justin’s new grid runs from [-nt+1:nt-1]: i.e. theta has moved two of its points to beyond the pi boundaries. It might have been cleaner to ADD two points outside the pi boundaries, if the external points are needed to cope with his tilted equilibria. Do the lost gridpoints get added back later to [-ntgrid:ntgrid], e.g. by interpolation?

[David Dickinson]

@Colin Malcolm Roach I think they do get used later (in some form) as with the current code Joseph’s test gives Nans but with the “fix” the test gives sensible behaviour. I don’t think those points outside pi are needed for the tilted case – I think the change happened in this routine to allow up-down asymmetry, originally the theta loop was from 1 tont , i.e. it only did half the range and assumed symmetry for the rest, Justin had to relax this assumption by changing the lower limit of the loop to -nt. This changed the number of points in the loop from nt to 2*nt+1, I think this means there are now nt+1 points to cover the range (0,pi) rather than nt points for the same range in the original code. In order to keep the same resolution of theta as the original code the range was changed to (0,pi+delta) so there were still nt points in the range (0,pi) and then one extra after (at least this is what I think the motivation was). This meant that one could compare grids with and without his changes to ensure things were in agreement as everything was on the same resolution grid. Here the fix will change (increase) the resolution of the theta grid used internally to leq rather than changing the range of theta covered. This means regression tests would need updating but I think that’s ok as we understand why.

[Colin Malcolm Roach]

I assumed extra points were exactly to get theta derivatives at +-pi in the tilted case, where none of the symmetry options work. The code achieved this by creating the extra points outside -pi:pi. We can calculate this derivative on the usual grid (with endpoints at+-pi) WITHOUT using any symmetry flags simply as follows:

df/dtheta(-nt) = df/dtheta(nt) = 0.5*(f(-nt+1)-f(nt-1))

This should work without needing symmetry flags for any equilibrium quantity (that has to be periodic in theta). I have not thought this through, but if we modify derm to do this more general theta derivative, could we then revert to the more usual (and perhaps less dangerous) theta grid where theta(nt)=pi andtheta(-nt)=-pi?

[Colin Malcolm Roach]

New branch created to resolve this issue https://bitbucket.org/gyrokinetics/gs2/branch/bugfix/leqtheta.

• Changes leq local theta grid to allocate all gridpoints in range [-pi:pi].
• Removes redundant symmetry assumptions used to evaluate theta derivatives of equm objects at theta=pi in subroutine derm.
• More general approach in derm uses values at pi-dtheta, and -pi+dtheta, and applies to any periodic single valued function in theta.
• New code should handle Justin Ball’s tilted equilbria without need for extra theta grid points outside [-pi:pi].
• New code is significantly more accurate at low ntheta, because all locally allocated theta grid points are used INSIDE [-pi:pi].

• Attached figure shows equilibrium metrics as functions of theta from original and “fixed” codes, for various choices of ntheta. Metrics were calculated for an leq equilbrium taken from the theta_grid unit tests. The “fixed” code starts to be a more obvious improvement for ntheta<=16.

  * Note that the fixed code appears to “fix” the mysterious segfault reported above by Joseph Parker in an earlier comment. * Note that the new code “FAILS” the theta_grid unit test, probably because it changes the result slightly. This will obviously have to be resolved before the code is merged into next.

See pull request https://bitbucket.org/gyrokinetics/gs2/pull-requests/196

[Peter Hill]

Thanks @Colin Malcolm Roach , those graphs clearly demonstrate that the two approaches converge to the same result, with the "fixed" version converging faster, which is what we expected. Because the failing tests in test_asym_geo*.pf are all regression tests, that just means we need to update the expected answers with the new results.

[David Dickinson]

The difference seems quite significant for low ntheta, with the fixed version generally doing a better job of matching the high resolution cases (as we’d expect as the fix effectively increases the resolution by using two extra points).

[David Dickinson]

I guess this also demonstrates that this issue shouldn’t have impacted any runs that were converged in ntheta

[Colin Malcolm Roach]

Yes… nothing was absolutely broken, the code was just slightly less accurate than it could have been.

[David Dickinson]

Update golden numbers in asymGeo pfunit test to account for recent bugfix to leq (see bug #77)

→ <<cset 7af6b1f6197d>>

[David Dickinson]

Fixed in next, will be in 8.1