Open lazersos opened 4 years ago
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
? 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 notifications@github.com 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 https://github.com/PrincetonUniversity/STELLOPT/issues/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
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?
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.
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 |
---|---|---|
@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.
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.
This third term is a combination of the chain rule applied to
sqrt(s)*f_o
and an averaging off_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
.
@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.
@zhucaoxiang it is because TERPSICHORE needs the Jacobian which is NYQUIST sized.
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).
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 98b166764088cdc179e0768b2ece2b2df293191c (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:
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)
@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.
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.
@jcschmitthttps://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 GitHubhttps://github.com/PrincetonUniversity/STELLOPT/issues/42#issuecomment-1065944365, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACUG23NTSQLTL3MCG2TSZILU7TTEVANCNFSM4KQ5ZS4Q. Triage notifications on the go with GitHub Mobile for iOShttps://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Androidhttps://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: @.***>
@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.
@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.
@lazersos Hmm... Can we still trust all the previous results?
Currently VMEC only output the correct values from sqrt(g), B^X, B_X, and |B| when
LNYQIST=.true.
is set inwrout.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.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:
Did I miss anything?