Meaning of jtwist for nperiod > 1
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.
Comments (29)
-
-
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.
-
reporter 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.
-
reporter @William Dorland @Michael Barnes @Jason Parisi Do you have opinions on this?
-
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.
-
reporter - edited description
-
The only diagnostic I see issues with is write_correlation_extend and correlation_extend in diagnostics_turbulence.f90, which I believe is largely unused.
-
reporter 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
-
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
-
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.
-
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?
-
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.
-
@David Dickinson Just emailed you the input files.
-
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
-
reporter 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?
-
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. -
reporter 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.
-
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.
-
reporter Might be worth setting wfbbc_option = “passing” to see if that fixes it. Wfb is the only bit of vspace that seems problematic.
-
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.
-
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.
-
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 casecollision_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.
-
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
andtprim=fprim=cvdriftknob=gbdriftknob=0
. In this limit the only field dependent source term isv||.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 arrayvpar
between FT and BS. This involvesv||
and some geometrical coefficients. It’s possible to show thatgradpar
, which entersvpar
, 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. -
Here’s a comparison of the source at ntheta=8
and ntheta= 128
This is when forcing
phi=phinew=1
andtprim=fprim=cvdriftknob=gbdriftknob=0
-
It seems there is certainly something happening leading to somewhat different geometrical quantities being calculated. Here’s a plot showing
gds2
,gds21
andgds22
for the -pi to pi range of FT and BS runs:
-
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 passntgrid
to indicate the number of theta grid points in use. Unfortunately we passedntgrid
which is the number of points in the extended domain, i.e.(2*nperiod-1)*(2*nth+1)
where2*nth+1
is the number of theta grid points for the range-pi
topi
(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 likebgradient
we weren’t using a consistent theta grid so only actually used the2*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 andomega = 0.33026171+0.12309022i
for ballooning space. Not a perfect match, but much better than without this fix. -
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.
-
-
- changed status to resolved
Fixed in release 8.0.6
- Log in to comment
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 totheta(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)
asky(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?