Set up reading ky from n0 in gs2 input file

[bpatel2107]

This PR lets you read in GS2 input files where the toroidal mode number is specified with rhostar where ky = n0 * rhostar * drho_dpsi

A slight issue here is that GS2 calculates drho_dpsi internally and to calculate something equivalent in pyro we do $\frac{r}{q}\frac{B_{unit}}{B0}$. The $\frac{B{unit}}{B_0}$ is calculated from a loop integral and I imagine GS2 does something similar. The main issue is that GS2 and Pyro will only really agree if that integral is done with ~high enough~ similar resolution. Pyro does this with 256 theta points whilst GS2 would use ntheta_geometry if specified, otherwise ntheta and when ntheta is relatively low they won't really agree.

I could force ntheta_geometry= 256 but that seem arbitrary. Alternatively I could do the integral in pyro at lower resolution but that doesn't seem nice either though leads to better agreement with the GS2 calc.

This doesn't impact reading the outputs which takes ky from the netCDF, but it would impact a conversion from a GS2 input file to another code (which is how I got here...)

Maybe the real solution is to allow for n to be specified...

[ZedThree]

@d7919 and I talked about potentially pulling out the GS2 geometry calculations into a standalone project which we could then add a Python wrapper to. Is that something that would be useful here?

[d7919]

This PR lets you read in GS2 input files where the toroidal mode number is specified with rhostar where ky = n0 * rhostar * drho_dpsi

A slight issue here is that GS2 calculates drho_dpsi internally and to calculate something equivalent in pyro we do rqBunitB0. The BunitB0 is calculated from a loop integral and I imagine GS2 does something similar. The main issue is that GS2 and Pyro will only really agree if that integral is done with ~high enough~ similar resolution. Pyro does this with 256 theta points whilst GS2 would use ntheta_geometry if specified, otherwise ntheta and when ntheta is relatively low they won't really agree.

I could force ntheta_geometry= 256 but that seem arbitrary. Alternatively I could do the integral in pyro at lower resolution but that doesn't seem nice either though leads to better agreement with the GS2 calc.

This doesn't impact reading the outputs which takes ky from the netCDF, but it would impact a conversion from a GS2 input file to another code (which is how I got here...)

Maybe the real solution is to allow for n to be specified...

Presumably one only needs to worry about the conversion from n to ky when trying to convert the input to another format? In this case it shouldn't matter* if GS2 and pyro have different drho_dpsi calculations as we are trying to ask pyro to produce an input file to study a particular toroidal mode number -- as long as pyro calculates drho_dpsi accurately this is achieved. Of course, this may mean that running GS2 on the original input and the generated input with the other code will result in the codes using different ky and hence not agreeing if compared directly. This should be less of an issue when producing plots of spectra (i.e. growth rate vs ky).

*Of course the user may actually be expecting pyro to produce an input file which should reproduce the GS2 results exactly.

If you use the input file ky for anything other than conversion then that seems slightly surprising, at least when you have an output file you should get ky from there.

[bpatel2107]

That is the main issue I'm facing converting as I have a GS2 case where n0 was specified and I am trying to run a CGYRO case with this as a benchmark. I could of course look at the ky in the output but I am not keen on having the GKInput object having dependencies on the GKOutput object. It also used different reference values meaning which resulted in #311.

I may be being slightly pedantic here as the "ky". For the case I'm looking at, if I force pyro to calculate drho_dpsi at the same resolution as GS2 they are within 1.4% (ntheta=32). If I use a higher ntheta_geometry=256 it drops down to 0.8% so we're not miles off.

At the moment we don't support specifying a toroidal mode number so pyro will always set "aky" in the input file. Of course GS2 doesn't care about that but it may be confusing as a user could think they are changing ky when they actually aren't.

I guess there's a little concern in that drho_dpsi can depend on the ntheta/ntheta_geometry used (both for pyro and gs2). You could offer the option to specify drho_dpsi but at that point you may as well just specify ky

[d7919]

I wasn't trying to suggest that we use the output as a part of the conversion, but rather we shouldn't use the input when comparing outputs. One options here is GS2_n --> CGYRO --> GS2_ky which would then at least give two consistent inputs (CGYRO and GS2_ky), although it's less helpful if you've already run GS2_n. Alternatively might just need to emit a warning when converting from GS2_n?

If you just want consistent inputs then you could use ntheta/ntheta_geometry to do the pyro calculation, sure it's not converged but presumably neither is GS2 there.

[bpatel2107]

I wasn't trying to suggest that we use the output as a part of the conversion, but rather we shouldn't use the input when comparing outputs. One options here is GS2_n --> CGYRO --> GS2_ky which would then at least give two consistent inputs (CGYRO and GS2_ky), although it's less helpful if you've already run GS2_n. Alternatively might just need to emit a warning when converting from GS2_n?

To be honest even when converting from GS2 to CGYRO we use bunit_over_b0 which inversely proportional to dpsi_drho so that loop integral will be there regardless. The question I guess is which resolution do we use when calculating bunit_over_b0. What you suggest above is effectively the same as using the default resolution (ntheta=256) when calculating bunit_over_b0

If you just want consistent inputs then you could use ntheta/ntheta_geometry to do the pyro calculation, sure it's not converged but presumably neither is GS2 there.

So that is what I've done here https://github.com/pyro-kinetics/pyrokinetics/blob/2150b3531463b6254bcc177a14c8699b51850f3b/src/pyrokinetics/gk_code/gs2.py#L421-L432

which bring that values calculated by both codes to within a 1% or so.

If we leave it to the default then we are doing what you suggest above.

https://github.com/pyro-kinetics/pyrokinetics/blob/2150b3531463b6254bcc177a14c8699b51850f3b/src/pyrokinetics/local_geometry/local_geometry.py#L525-L540

Really not sure there's a clear cut way here. Maybe a warning is best.

[d7919]

Yes I think a warning is probably best. How expensive is the bunit_over_b0 calculation? Could we do it with the input ntheta and the default 256 so that the warning can include some rough estimate of the error?

[bpatel2107]

Not expensive at all. Can include it but it doesn't account for the error between the GS2 and the pyro calculation due to different numerical implementations. To be honest for a simple integral like that I'm surprised they are not closer

[d7919]

Just looking at the integral -- https://github.com/pyro-kinetics/pyrokinetics/blob/2150b3531463b6254bcc177a14c8699b51850f3b/src/pyrokinetics/local_geometry/local_geometry.py#L561C10-L561C15

If dL is defined from 0 up to and including 2 pi then doesn't this sum double count the contribution from the duplicate theta = 0 = 2 * pi?

[bpatel2107]

I think you might be right which is very annoying as that will slightly tweak anytime where a GS2/GENE <-> CGYRO/TGLF conversion occurred. The impact may be very small but still there.

The odd thing is that GS2 and pyro are closer when we include that point. For the case I'm looking at now (ntheta_geometry=256)

GS2 drho_dpsi: 1.8309115626684531 Pyro drho_dpsi (double count): 1.8309140191745965 Pyro drho_dpsi (remove last point): 1.8322749988846894

The percentage errors go from 0.0001% to 0.07% (I was miscalculating my error before they actually were closer)

[bpatel2107]

@d7919 I thought you might be right but I tried a simple test with a unit circle and it seems that we should include it.

import numpy as np

theta = np.linspace(0, 2 * np.pi, 1000000)

R = np.cos(theta)
Z = np.sin(theta)

dR = (np.roll(R, 1) - np.roll(R, -1)) / 2.0
dZ = (np.roll(Z, 1) - np.roll(Z, -1)) / 2.0
dL = np.sqrt(dR**2 + dZ**2)

integral = np.sum(dL) / (2 * np.pi)

print(integral)

Gives me

> 0.9999999999934204

Our integral seems to be under but I have no idea why, if anything it should be over with the double counting.

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 84.84848% with 5 lines in your changes are missing coverage. Please review.

Project coverage is 85.39%. Comparing base (930d32c) to head (87991ea). Report is 163 commits behind head on unstable.

:exclamation: Current head 87991ea differs from pull request most recent head 4145c1c

Please upload reports for the commit 4145c1c to get more accurate results.

Files	Patch %	Lines
src/pyrokinetics/gk_code/gs2.py	82.75%	5 Missing :warning:

Additional details and impacted files

```diff @@ Coverage Diff @@ ## unstable #324 +/- ## ============================================ + Coverage 81.06% 85.39% +4.32% ============================================ Files 47 48 +1 Lines 8124 8771 +647 ============================================ + Hits 6586 7490 +904 + Misses 1538 1281 -257 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[bpatel2107]

I've moved the Bunit/B0 issue to #346 so we are good to go on this end. Need to handle #344 properly but this as is should be fine

[bpatel2107]

@d7919 I think I'm good with this now so will merge.

[bpatel2107]

@mrhardman this also fixes the issue you were facing in STELLA

[bpatel2107]

Now it fixes the issue in STELLA...