Meaning of jtwist for nperiod > 1

[Michael Hardman]

Hi everyone,

I am creating this issue to make the developers aware of a point which I believe may be confusing, although it perhaps is not a bug, but a feature.

For nonlinear simulations we typically take nperiod = 1, corresponding to a theta domain of (-pi, pi). GS2 users learn in these simulations that the input parameter "jtwist" is equivalent to delta j in eqn (20) of Beer Phys. Plasmas 2 (7), July 1995, i.e.,

jtwist = 2*pi *shat *dky / dkx

However, for nperiod > 1, the simulation domain spans (2*nperiod - 1) (-pi,pi) segments in theta, and for consistency the definition of jtwist should be

jtwist = 2*pi (_2*_nperiod - 1)*shat *dky / dkx , where (2*nperiod - 1) = N in the Beer paper.

The issue that I would like to raise is that jtwist may be accidentally defined inconsistently between kt_grids.f90 and dist_fn.fpp.

‌

from kt_grids.f90 lns 1613 box_get_grids

 else
       if (jtwist /= 0) then
          dkx = dky * 2.0*pi*abs(shat)/real(jtwist)
       else
          dkx = dky
       end if
    endif

Notice no factor nperiod

‌

but in dist_fn.fpp lns 2121 init_connected_bc

 ng = ntheta/2 + (nperiod-1)*ntheta

    ! jshift0 corresponds to J (not delta j) from Beer thesis (unsure about +0.01) -- MAB
    ! delta j is the number of akxs (i.e. theta0s) 'skipped' when connecting one
    ! parallel segment to the next to satisfy the twist and shift
    ! boundary conditions. delta j = (ik - 1) * jshift0 . EGH

    if (naky > 1 .and. ntheta0 > 1) then
       jshift0 = nint((theta(ng)-theta(-ng))/(theta0(2,2)-theta0(1,2)))
    else if (naky == 1 .and. ntheta0 > 1 .and. aky(1) /= 0.0) then
       jshift0 = nint((theta(ng)-theta(-ng))/(theta0(2,1)-theta0(1,1)))
    else
       jshift0 = 1
    end if

here jshift0 effectively becomes (2*nperiod - 1) * jtwist, because theta grid can be longer than -pi, pi and dkx (theta0*shat*dky) defined above.

As a result, if one wanted to run a gyrokinetic simulation with nperiod > 1, in box mode, one would have to manually supply

jtwist = delta j_intended / (2*nperiod - 1)

edit: Thanks to Jason for pointing out the above should be a division, not multiplication

‌

where delta j_intended is the intended quantisation index of the grids defined in eq. (20) of Beer 1995.

This is not a problem once one is aware of the issue, but for users this is potentially confusing, and possibly a great way to lose computing time if a supervisor suggests going to nperiod > 1 nonlinearly.

I would be interested to discuss further.

[David Dickinson]

Thanks for flagging this. I’ve been interested in nperiod > 1 runs for a while but hadn’t come across this issue.

Just to expand on this for people not familiar with this bit of code

2*pi*(2*nperiod-1) is equivalent to theta(ntgrid)-theta(-ntgrid) and this is the length of a single “cell”/parallel domain which we stitch together to created our extended domains (or supercells).

dky = ky(1) as ky(0)==0

‌

I guess it is relatively easy to show this if we were to dump out kperp2 and plot this along a supercell for nperiod > 1?

[David Dickinson]

Is it as simple to fix as replacing

dkx = dky * 2.0*pi*abs(shat)/real(jtwist)

with

dkx = dky * 2.0*pi*(2*nperiod-1)*abs(shat)/real(jtwist)

?

‌

I think jtwist is also used in some of the diagnostics in order to calculate extended theta domains but I’m not familiar enough with these diagnostics to know if there use of jtwist is consistent with nperiod>1. From a quick look I think they’re ok, but not 100% confident.

[Michael Hardman]

I would be inclined to say yes, this is the right thing to do.

dkx = dky * 2.0*pi*(2*nperiod-1)*abs(shat)/real(jtwist)

Jason P was interested in going to nperiod > 1, so hopefully he will be able to provide a test case for us to confirm that the output makes sense with this change. I am not aware of any diagnostics which use jtwist explicitly (except perhaps one of the unused parallel correlation length diagnostics). I can take a look if required.

[Michael Hardman]

@Bill Dorland @Michael Barnes @Jason Parisi Do you have opinions on this?

[Jason Parisi]

I agree with David’s suggestion. In this implementation, where Delta k_x is multiplied by (2*nperiod -1), we still have iky*jtwist independent ballooning mode chains. Some of our post processing scripts (@Michael Hardman) do use jtwist as an input, so we’ll need to modify that. Checking the diagnostics now.

[Jason Parisi]

The only diagnostic I see issues with is write_correlation_extend and correlation_extend in diagnostics_turbulence.f90, which I believe is largely unused.

[Michael Hardman]

I don't think that we will need to edit post processing scripts. I (we?) want the meaning of jtwist to be the integer spacing in the kx grid between connected modes. In the current definition this is not the case for nperiod > 1. In fact we need to make this fix in order to make a post processing scripts (and init_connected_bc) work as intended

[Jason Parisi]

Hi all,

After discussion with Michael Hardman, I have made a branch of GS2 with the correct dkx for all nperiod values, which I have confirmed works on linear simulations for nperiod > 1. When testing on linear runs, there is still a disagreement in the growth rates for a ‘range’-like run and a ‘box’-like run, but I guess we’ve known this for a long time. I am currently testing with nonlinear simulations, but I don’t foresee any strange behavior. The branch name is origin/feature/fast_response_flux_dist_basic_hyper_correct_nperiod. It also has some additional features I have tested fairly extensively, so I will create a pull request.

Jason

[David Dickinson]

Thanks for implementing the fix. I’m curious how different you are finding box and range linear runs. This is something I looked at a long time ago during my PhD work and I thought the agreement was generally rather good once we’d dealt with things like ensuring the duplicate theta points had a unique value.

[Jason Parisi]

I found that the growth rate varied by a lot (a factor of 4), but I was using ntheta = 8. I have heard that the agreement gets better once ntheta increases in value. Is this consistent with what you have seen?

[David Dickinson]

I think I’ve only tested agreement in cases where the ballooning results are converged in grid size. I’d be very concerned about even a 10% difference when comparing range and box results when the setup should be equivalent. If you have a test case which shows poor agreement I’d be interested to explore this further.

[Jason Parisi]

@David Dickinson Just emailed you the input files.

[David Dickinson]

Just to record the results of this investigation here are comparisons of the mode structure and complex frequency for different ntheta

‌

Ntheta = 8

Ntheta = 32

Ntheta = 128

‌

[Michael Hardman]

Thanks David for making this comparison and sharing the plots here. This case does not seem so pathological for "ordinary" resolution. Were you able to get any insight into why code does a slightly different calculation for the two different grid layouts?

[David Dickinson]

I’ve not really dug into it for this case but the disagreement appears largest at the flux tube cell boundaries. The most obvious place I can think of there being a difference is in the invert_rhs_linked routine where we adjust the local solutions to ensure continuity across the linked boundaries, which we don’t need to do in ballooning space. I suppose this could lead to small ~floating point round off level differences. I don’t really expect this to be the cause of any significant discrepancy.

[Michael Hardman]

From these plots of the potential it does not look like that roundoff error is cumulative. I suppose one could plot the distribution function h_{+} of only forward going particles, and normalise to the outgoing value to see if there was a diverging difference between h_{+} for the two different methods.

[David Dickinson]

Yes that sounds sensible. I’ll dig into this a bit more to see if I can see any obvious cause. Originally I would have expected even the ntheta=8 case to have reasonable agreement.

[Michael Hardman]

Might be worth setting wfbbc_option = “passing” to see if that fixes it. Wfb is the only bit of vspace that seems problematic.

[David Dickinson]

Trying the different wfb treatments doesn’t appear to change things significantly. I’ve also compared the response matrix for this supercell and for ntheta=8 the relative difference between FT and BS is around 2.5%. This doesn’t mean the issue is in the response matrix (the two field modules agree near exactly) as the response matrix calculation mostly calls timeadv. It does help make it clear that there is a difference.

I can also note that flux tube with no links agrees exactly with the equivalent ballooning space – not surprising, but at least somewhat reassuring.

[David Dickinson]

Interestingly the initial field seems to already be in disagreement. I think this may help narrow it down to something in gamtot/antot or our Maxwell field calculation. I’m using the default ginit option for a width of -1 the two approaches are in good agreement and the initial field magnitude has dropped to almost 1e-308 by the end of the domain. If I change width to the default of -3.5 or -10.0 the initial field more obviously disagrees. Even when the magnitude agrees it seems that the phase may vary.

The choice of width doesn’t change the converged omega so this may be an unrelated point.

[David Dickinson]

I’ve started to try working with a case with most of the source term disabled by forcing the drifts to zero and setting tprim=fprim=0. In this case I find both setups are strongly damped from collisions but the damping rate is different. I therefore experimented with different collision models leading to the following comparisons for the ntheta=8 case

collision_model = ‘default’:

collision_model = ‘ediffuse’:

collision_model = ‘lorentz’:

It seems that when the energy diffusion is included there is a more significant different between FT and BS and the mode structure becomes much more jagged. The collisionless case is also quite jagged and not in good agreement, so it’s not just a consequence of energy diffusion.

Note the omega given in the legends for these cases is wrong.

[David Dickinson]

By looking at the source term calculating at the initial time it seems clear that there are differences between FT and BS even when forcing phi=phinew=1 and tprim=fprim=cvdriftknob=gbdriftknob=0. In this limit the only field dependent source term is v||.grad phi. The field independent part of the source term is very close between the two cases. This suggests that there is a difference in the array vpar between FT and BS. This involves v|| and some geometrical coefficients. It’s possible to show that gradpar , which enters vpar, is slightly different between the two setups. Switching from Miller to s-alpha I find that the field dependent source term is identical between FT and BS, further indicating that at least part of the difference is arising from the geometry handling. There are several places this could come in and I’ll dig further.

[David Dickinson]

Here’s a comparison of the source at ntheta=8

and ntheta= 128

‌

This is when forcing phi=phinew=1 and tprim=fprim=cvdriftknob=gbdriftknob=0

[David Dickinson]

It seems there is certainly something happening leading to somewhat different geometrical quantities being calculated. Here’s a plot showing gds2, gds21 and gds22 for the -pi to pi range of FT and BS runs:

‌

[David Dickinson]

I’ve managed to finally track this down to a suspected bug in our use of Miller. When we initialise Miller we call leqin and pass ntgrid to indicate the number of theta grid points in use. Unfortunately we passed ntgrid which is the number of points in the extended domain, i.e. (2*nperiod-1)*(2*nth+1) where 2*nth+1 is the number of theta grid points for the range -pi to pi (i.e. ntheta+1). This meant that the theta grid calculated internally to the local equilibrium moduleleq was higher resolution than intended and therefore when we calculated index space derivatives (i.e. differences of theta on the grid) and used these in calls to routines like bgradient we weren’t using a consistent theta grid so only actually used the 2*nth+1 values in the centre of the high resolution grid.

I think this bug has been present since at least release 2.0 and I think will have caused mistakes for any Miller run with nperiod > 1. The impact of this will likely be more severe with higher nperiod. I’ll hopefully push a bug fix tomorrow once I’ve had chance to check this over again and understand why increasing ntheta seemed to help this problem go away. With this bug fix I find that the geometrical quantities now match very well between flux tube and ballooning space runs and the source term and mode structure match closely. For the ntheta=8 case we now find omega = 0.32977119+0.12224306i for flux tube and omega = 0.33026171+0.12309022i for ballooning space. Not a perfect match, but much better than without this fix.

[David Dickinson]

See PR #372 for proposed fix. Note that even with this fix it appears geometrical quantities can disagree slightly between runs with nperiod = 1 and nperiod = 2.

[David Dickinson]

Fixed in release 8.0.6