Implement calculation of sqrt(g), B^X, B_X, and |B| from R, Z, lambda, iota, phi'

[lazersos]

Currently VMEC only output the correct values from sqrt(g), B^X, B_X, and |B| when LNYQIST=.true. is set in wrout.f90. When this is done these arrays have a different Fourier spectrum than the R, Z and Lambda arrays. This is going to be addressed in two ways.

1. Code which can use the larger arrays will be modified
2. Where possible such quantities will be calcualted directly from R, Z, lambda, dphi/ds and iota.

For reference
sqrt(g) = R*(RuZs-RsZu) where Ru=dR/du
B^u = phiprime*(iota-Lv)/sqrt(g) where Lv=dlambda/dv
B^v = phiprime(1+Lu)/sqrt(g)
|B| = sqrt(B^u2guu+2B^uB^vguv+b^v2*gvv)

Here's what I suggest be worked on:

• BOOZ_XFORM: Already implemented to handle LNQUIST but we should benchmark this since BOOTSJ, DKES, and NEO rely on these values. (volunteer?)
• LIBSTELL:
  • STEL_TOOLS: Can handle LNQUIST if you re-sort the R, Z, and L arrays. Will upgrade to use R, Z and lambda to calc the above values on internal grid (@lazersos)
  • VMEC_TOOLS: Straightforward to modify (@lazersos)
  • VIRTUAL_CASING_MOD: Similar situation to STEL_TOOLS (@lazersos)
• BNORM: Needs to be implemented and checked. (vounteer?)
• COBRAVMEC: Need to be implemented (@lazersos) - Done
• STELLOPT
  • The code uses STEL_TOOLS for most everything so once that's been done we should be fine.
• BEAMS3D/DIAGNO/FIELDLINES: These rely on the LIBSTELL codes so once that's done they'll be fixed as well. (@lazersos)

Did I miss anything?

[krystophny]

Perhaps related: the output of xn_nyq modes in wout.nc has opposite sign from xn , no matter if LNYQUIST is true or false. One should check the sign convention of zeta once more, and which of xn or xn_nyq produces the intended one. Correction: seems like this affects only matlabVMEC, where xn switches sign but not xn_nyq

[landreman]

? xn and xn_nyq always seems to have the same sign, both in wout.nc files I've generated by running stellopt's vmec myself and in wout files I've received from IPP.

…

On Fri, Feb 7, 2020 at 6:27 AM Christopher Albert ***@***.***> wrote: Perhaps related: the output of xn_nyq modes in wout.nc has opposite sign from xn , no matter if LNYQUIST is true or false. One should check the sign convention of zeta once more, and which of xn or xn_nyq produces the intended one. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#42?email_source=notifications&email_token=ABH4ZA4OYI7EGQNGELC3NBDRBVASDA5CNFSM4KQ5ZS42YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELCTWCQ#issuecomment-583351050>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABH4ZA475HIFCLYIOJEDWR3RBVASDANCNFSM4KQ5ZS4Q> .

-- ======================================= Dr. Matt Landreman Associate Research Scientist Institute for Research in Electronics & Applied Physics University of Maryland 8223 Paint Branch Drive, College Park MD 20742, USA (+1) 651-366-9306 mattland@umd.edu

[krystophny]

Sorry for the confusion - it was really just in matlabVMEC where this occured, and otherwise the sign was the same directly from wout.nc . No idea why - @lazersos what do you think?

[lazersos]

So there was an issue in read_vmec.m where fields of the structure (xnnyq) were being returned when reading a netCDF file, but we were checking for xn_nyq to change the kernel to (mu+nv) from (mu-nv). This has been updated on the git page for matlabVMEC. Unfortunately I no longer have access to my Princeton account so I can't update the files on the Mathworks File exchange. I have posted a comment however pointing to the repo on GitHub.

[krystophny]

Thanks! This whole story became a bit more confusing, and I'll try to summarize it here. We made two runs on the same equilibrium, one time with ntor=mpol=6 and one time with ntor=mpol=12 harmonics. Here are the results for |B| on the flux surface s=0.3 .

Column 1: truncated spectrum (LNYQUIST=.false.) from VMEC .nc (xn_nyq=xn, xm_nyq=xm)
Column 2: full spectrum (LNYQUIST=.true.) from VMEC .nc
Column 3: from geometry via R, Z, phi and lambda and iota (as in SIMPLE alpha tracer)

The two observations are

1. In the 6x6 truncation causes a big inconsistency, while in 12x12 fields are optically the same
2. Artificial smoothing filter by truncating the spectrum in 6x6 leads to a field that matches the more accurate 12x12 calculation (by chance?).

This means that when using low resolution VMEC runs such as 6x6, one might get results matching reality better by artificially truncating / filtering the spectrum of the B field, even though it's not consistent with low-resolution R,Z geometry. It would be interesting to investigate this in detail.

6x6

truncated	full	geometry

12x12

truncated	full	geometry

[lazersos]

@krystophny Interesting. I've been reading through bit more of VMEC to better understand how the GMNC array is constructed. One thing that VMEC does internally is to separate the even and odd modes of R,Z, and Lambda. They then divide the odd modes by sqrt(s), then transform. So to compose a value in real space it looks like

r=r_e+sqrt(s)*r_o

where the r_e and r_o are the real space Fourier transforms of the even and odd modes separately. This has an effect that to calculate df/ds (on the half grid) you get something like

dfds(i) = (ns-1)*(f_e(i)-f_e(i-1) + sqrt(s)*(f_o(i)-f_o(i-1))) + 0.25*(f_o(i)+f_o(i-1))/sqrt(s)

where I've taken sqrt(s) to be on the half grid for simplicity. This third term is a combination of the chain rule applied to sqrt(s)*f_o and an averaging of f_o onto the half grid (hence the 0.25 coefficient instead of 0.5).

This has implication for the computation of the cross product. Essentially it has extra terms as can be seen in jacobian.f. Using some matlab pseudo-code I'm able to calculate sqrt(g), B^u, and B^v in real space form R, Z, and Lambda. There is a small error when comparing to the LNYQUIST=T version of the GMNC, BSUPUMNC, and BSUPVMNC arrays but I suspect this is roundoff and not significant.

I've raised the issue with the ORNL/VMEC GitHub site and am tracking it there with their developers. The end goal will be to only use the basic quantities from VMEC in the future so as to be LNYQUIST ambivalent.

[lazersos]

I just wanted to give a quick explanation of why this is an issue. The VMEC jacobian looks like
R*(dR/du*dZ/ds-dR/ds*dZ/du). Because R and Z have opposite parity (cosine vs sine) the convolutions of dR/du*dZ/ds give terms which go like cos^2 and sin^2. However, to properly Fourier transform these terms back into Fourier space you need twice the harmonic content as before. Hence our need for the Nyquist sized array for these values.

[jonathanschilling]

This third term is a combination of the chain rule applied to sqrt(s)*f_o and an averaging of f_o onto the half grid (hence the 0.25 coefficient instead of 0.5).

@lazersos Thanks for pointing this out! I always wondered about the 0.25 in Compute_Currents in read_wout_mod.f90.

[zhucaoxiang]

@lazersos You also changed TERPSICHORE to use mnmax_nyq. There is a bug of not broadcasting the mnmax_nyq and it is now fixed in #90 . But I am not quite sure if it is necessary to use the Nyquist number in TERPSCHORE.

[lazersos]

@zhucaoxiang it is because TERPSICHORE needs the Jacobian which is NYQUIST sized.

[lazersos]

Just an update to this issue. The even and odd modes of R, Z, and Lambda do need to be interpolated differently. The even modes can be linearly interpolated in flux but the Odd modes must be interpolated in rho, so the harmonics need to be multiplied by 1/sqrt(s). Doing this one can properly interpolate R, Z, and lambda. However, the radial derivatives still have issues as we approach the magnetic axis. Essentially the Jacobian becomes singular which messes everything up. So I wasn't able to switch GetBcyl over to using R,Z and lambda solely. However, I did fix the interpolation of the Jacobian which improves calculation of B near the axis.

Looking at VMEC, internally the code recomputes lambda on the half grid before outputting (for back-compatibility with old version). However, this may not be done correctly near the axis (or may be hiding special treatment).

[jcschmitt]

Similar to using TERPSICHORE, when using stellopt or vmec with cobravmec, LNYQUIST=F must be specified in the vmec-portion of the input namelist.

On branch Lazerson with commit 98b1667 (HEAD -> Lazerson, origin/Lazerson)

I did a benchmarking test of vmec, stellopt and cobravmec through the various ways that the code(s) handle the LNYQUIST flag, and through the various ways that that codes can be run: vmec stand-alone, stellopt in 'single_iter' mode, and stellopt in standard 'lmdiff_bounded' mode with niter=1, (similar to lmdiff, but with bounds), and 'lmdiff_bounded' mode with niter>1.

So, I have 12 types of code paths:

1. A standalone vmec run, with LNYQUIST=F, followed by a standalone cobravmec run.
2. A standalone vmec run, with LNYQUIST=F, followed by a standalone cobravmec run.
3. A standalone vmec run, with LNYQUIST UNSPECIFIED (commented out), followed by a standalone cobravmec run.
4. A stellopt run with opt_type = 'one_iter', followed by a standalone cobravmec run using the wout file from the stellopt run with LNYQUIST=F in the vmec portion of the input namelist
5. Same as 4, with LNYQUIST=T
6. Same as 4 with LNYQUIST UNSPECIFIED (commented out)
7. A stellopt run with opt_type = 'lmdiff_bounded' and NITER=1, and I inspect the cost function output file (stellopt.whatever) to determine the actually ballooning growth rates. LNYQUIST=F in the vmec portion of input namelist.
8. Same as 7, with LNYQUIST=T
9. Same as 7, LNYQUIST UNSPECIFIED (commented out)
10. Same as 7, with NITER = 10 (just to make sure that stellopt runs with multiple iterations)
11. Same as 8, with NITER = 10
12. Same as 9, with NITER = 10

Only the runs with LNYQUIST=False produce the correct results.
This is the same situation as using TERPSICHORE.
So: LNYQUIST must be specified and it must be set to FALSE (or f, or 0, or .false., or whatever logical flag works in that spot)

cobravmec nyquist checks.pdf

[lazersos]

@jcschmitt I can see in COBRAVMEC at order_input the change which needs to happen. However, I disagree that any result obtained with LYQUIST=F is 'correct'. When LNYQUIST=F is set the modes for anything other than the R, Z and Lambda harmonics are artificially truncated. I will make the change to order_input with the correct behavior. This change will propagate to STELLOPT as well.

[jcschmitt]

Sounds good. On Monday, I’ll go over this with Aaron and Ben (they wrote a julia version that agrees with what I called the ‘correct’ results. )

…

________________________________ From: Samuel Lazerson ***@***.***> Sent: Saturday, March 12, 2022 1:04 PM To: PrincetonUniversity/STELLOPT ***@***.***> Cc: John Schmitt ***@***.***>; Mention ***@***.***> Subject: [EXT] Re: [PrincetonUniversity/STELLOPT] Implement calculation of sqrt(g), B^X, B_X, and |B| from R, Z, lambda, iota, phi' (#42) CAUTION: Email Originated Outside of Auburn. @jcschmitt<https://github.com/jcschmitt> I can see in COBRAVMEC at order_input the change which needs to happen. However, I disagree that any result obtained with LYQUIST=F is 'correct'. When LNYQUIST=F is set the modes for anything other than the R, Z and Lambda harmonics are artificially truncated. I will make the change to order_input with the correct behavior. This change will propagate to STELLOPT as well. — Reply to this email directly, view it on GitHub<#42 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACUG23NTSQLTL3MCG2TSZILU7TTEVANCNFSM4KQ5ZS4Q>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[lazersos]

@jcschmitt If they're not using the full nyquist sized results then it's isn't 'correct'. I've pushed a branch which should fix the issue for COBRAVMEC, but initial tests indicate it's slower since the searches over Lambda now have to deal with the larger mode structure.

[lazersos]

@jcschmitt In the branch I've pushed #147 the code now properly handles LNYQUIST sized arrays. However, I've noticed that inside of COBRAVMEC all arrays are linearly interpolated from the half to full mesh. This isn't entirely correct as the odd modes should go as sqrt(s). I'm going to open this as a new issue since I think the proper handling of interpolation half/full and calculation of df/ds affects many codes.

[zhucaoxiang]

@lazersos Hmm... Can we still trust all the previous results?