Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update VMECIO to specify Nyquist spectrum and fix issues with asymmetric wouts #1211

Merged
merged 24 commits into from
Sep 1, 2024

Conversation

dpanici
Copy link
Collaborator

@dpanici dpanici commented Aug 21, 2024

  • plot_vmec_comparison sneakily only worked for sym=True wout files, this changes it to work correctly
  • allow M_nyq and N_nyq to be passed as keywords to VMECIO.save
  • fix typo in axis loading from an asymmetric VMEC equilibrium
  • fix sneaky bug in asymmetric wout saving where some of the derived quantities were being fit with the transform directly instead of using the fullfit function which did the fit AND negated the correct modes to account for poloidal angle reversal... this will help in Add Ability to Write Boozmn.nc style output files #680 and also resolve Saving asymmetric eq as a wout has some errors #1062 I believe, as the magnetic field coeffs being saved before were wrt the wrong poloidal angle (and the opposite one as the jacobian fourier coefficients saved).
  • make test with an asymmetric stellarator

hopefully resolves #1062

Notes

For an NFP=4 stellarator, VMEC with NZETA=200 yields a max nyquist N of 400 while NZETA=100 yields a max nyquist N of 200, which leads me to believe the max nyquist toroidal modenumber is set as NFP * NZETA/2 or something along those lines (I support with a line from educational VMEC : https://github.com/jonathanschilling/educational_VMEC/blob/ad4db0a56b23c3324df3d6101781b69b2e3b7e13/src/fixaray.f90#L28)

default values of ntheta, nzeta are here

defaults:
NTHETA = 2M+6
NZETA = 2N+4

then nyquist freq determined by
mnyq = MAX(0, NTHETA, 2*(M-1))
nnyq = MAX(0, NZETA, 2*ntor)

translates to default of M+3 (in practice it seems to be M+4... don't understand completely but this is what the wout has)) and N+2 for the nyquist max spectrum defaults

When VMEC is given a right-handed bdry (such that sqrt(g)>0 ), it seems to change the modes such that the new bdry is left-handed (but the way of changing the bdry handedness changes where the theta=0 point is on the bdry, and is NOT how we enforce the change in handedness in DESC) (this is bc they change [theta -> pi - theta] (https://github.com/jonathanschilling/educational_VMEC/blob/ad4db0a56b23c3324df3d6101781b69b2e3b7e13/src/flip_theta.f90#L29).)
For example, the asymmetric rotating ellipse bdry given by

image

will be, in the resulting wout file, replaced by

image

This is a pitfall if trying to compare the wout from a VMEC run to a DESC equilibrium saved as a wout, as the output bdry modes in VMEC are not what you used as the input modes, AND the theta definition is changed not only in its direction but also with a pi shift.

@dpanici dpanici changed the title Fix some VMECIO issues with asymmetric equilibria Update VMECIO to specify Nyquist spectrum and fix issues with asymmetric equilibria Aug 21, 2024
Copy link
Contributor

github-actions bot commented Aug 21, 2024

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     +1.18 +/- 6.76     | +6.19e-03 +/- 3.55e-02 |  5.31e-01 +/- 3.4e-02  |  5.25e-01 +/- 1.1e-02  |
 test_equilibrium_init_medres            |     -6.58 +/- 4.76     | -2.86e-01 +/- 2.07e-01 |  4.06e+00 +/- 2.8e-02  |  4.34e+00 +/- 2.0e-01  |
 test_equilibrium_init_highres           |     +1.53 +/- 4.30     | +8.30e-02 +/- 2.34e-01 |  5.52e+00 +/- 2.1e-01  |  5.44e+00 +/- 1.1e-01  |
 test_objective_compile_dshape_current   |     +0.47 +/- 3.07     | +1.79e-02 +/- 1.17e-01 |  3.82e+00 +/- 6.6e-02  |  3.80e+00 +/- 9.6e-02  |
 test_objective_compute_dshape_current   |     -0.64 +/- 2.65     | -2.23e-05 +/- 9.15e-05 |  3.43e-03 +/- 4.6e-05  |  3.46e-03 +/- 7.9e-05  |
 test_objective_jac_dshape_current       |     -2.09 +/- 7.59     | -8.72e-04 +/- 3.17e-03 |  4.08e-02 +/- 2.2e-03  |  4.17e-02 +/- 2.3e-03  |
 test_perturb_2                          |     -0.91 +/- 4.06     | -1.63e-01 +/- 7.25e-01 |  1.77e+01 +/- 4.5e-01  |  1.79e+01 +/- 5.7e-01  |
 test_proximal_freeb_jac                 |     +0.87 +/- 1.86     | +6.54e-02 +/- 1.40e-01 |  7.60e+00 +/- 1.2e-01  |  7.54e+00 +/- 7.3e-02  |
 test_solve_fixed_iter                   |     +0.44 +/- 60.99    | +2.23e-02 +/- 3.08e+00 |  5.07e+00 +/- 2.2e+00  |  5.05e+00 +/- 2.1e+00  |
 test_build_transform_fft_midres         |    +10.72 +/- 3.71     | +6.40e-02 +/- 2.21e-02 |  6.61e-01 +/- 2.0e-02  |  5.97e-01 +/- 8.8e-03  |
 test_build_transform_fft_highres        |     +5.70 +/- 3.86     | +5.70e-02 +/- 3.86e-02 |  1.06e+00 +/- 3.1e-02  |  1.00e+00 +/- 2.4e-02  |
 test_equilibrium_init_lowres            |     +1.12 +/- 3.66     | +4.23e-02 +/- 1.38e-01 |  3.81e+00 +/- 1.1e-01  |  3.77e+00 +/- 8.3e-02  |
 test_objective_compile_atf              |     -1.09 +/- 4.04     | -8.57e-02 +/- 3.17e-01 |  7.75e+00 +/- 2.2e-01  |  7.83e+00 +/- 2.3e-01  |
 test_objective_compute_atf              |     -2.48 +/- 2.62     | -2.58e-04 +/- 2.72e-04 |  1.01e-02 +/- 1.1e-04  |  1.04e-02 +/- 2.5e-04  |
 test_objective_jac_atf                  |     -2.46 +/- 2.43     | -4.76e-02 +/- 4.71e-02 |  1.89e+00 +/- 4.3e-02  |  1.94e+00 +/- 2.0e-02  |
 test_perturb_1                          |     +1.36 +/- 5.05     | +1.69e-01 +/- 6.28e-01 |  1.26e+01 +/- 2.0e-01  |  1.24e+01 +/- 5.9e-01  |
 test_proximal_jac_atf                   |     +1.15 +/- 0.90     | +9.17e-02 +/- 7.20e-02 |  8.08e+00 +/- 5.0e-02  |  7.98e+00 +/- 5.2e-02  |
 test_proximal_freeb_compute             |     -1.24 +/- 0.90     | -2.27e-03 +/- 1.65e-03 |  1.81e-01 +/- 8.6e-04  |  1.84e-01 +/- 1.4e-03  |

@dpanici
Copy link
Collaborator Author

dpanici commented Aug 21, 2024

check NZETA does what I expect

Copy link

codecov bot commented Aug 23, 2024

Codecov Report

Attention: Patch coverage is 75.00000% with 6 lines in your changes missing coverage. Please review.

Project coverage is 95.31%. Comparing base (2ebdf6a) to head (95a3292).
Report is 25 commits behind head on master.

Files with missing lines Patch % Lines
desc/vmec.py 75.00% 6 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1211      +/-   ##
==========================================
- Coverage   95.33%   95.31%   -0.03%     
==========================================
  Files          90       90              
  Lines       22702    22713      +11     
==========================================
+ Hits        21644    21649       +5     
- Misses       1058     1064       +6     
Files with missing lines Coverage Δ
desc/vmec.py 91.91% <75.00%> (-0.41%) ⬇️

... and 1 file with indirect coverage changes

@dpanici dpanici marked this pull request as ready for review August 28, 2024 04:15
@dpanici dpanici requested review from a team, rahulgaur104, f0uriest and ddudt and removed request for a team August 29, 2024 15:00
@dpanici dpanici changed the title Update VMECIO to specify Nyquist spectrum and fix issues with asymmetric equilibria Update VMECIO to specify Nyquist spectrum and fix issues with asymmetric wouts Aug 29, 2024
if lmnc is None:
lmbda_mnc = lambda s: 0
else:
lmbda_mnc = interpolate.CubicSpline(si, lmnc)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is cubic spline the best way to interpolate fourier coefficients in the radial direction?
A fourier coefficient may not always be smooth.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might not be, esp. at edge, but I don't feel the need to change as this is not really used anywhere important in the code

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think vmec itself implicitly assumes the fourier coeffs at at least c1 smooth

atol=1e-4, # nearly zero bc vacuum
)
np.testing.assert_allclose(
vmec.variables["DWell"][20:100], desc.variables["DWell"][20:100], rtol=1e-2
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why exclude points close to the axis? Because of strong variations?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

b/c VMEC is notorious for being inaccurate near-axis. But in this case the two solutions (the DESC and VMEC) one are not exactly similar (even visually) near-axis, so we also want to exclude it for that reason.

I dont need this test to be a validation of DESC v VMEC, I just want it to tell me if I am way off in the asymmetric quantities basically

@@ -242,8 +248,9 @@ def save(cls, eq, path, surfs=128, verbose=1): # noqa: C901 - FIXME - simplify
NFP = eq.NFP
M = eq.M
N = eq.N
M_nyq = M + 4
N_nyq = N + 2 if N > 0 else 0
M_nyq = M + 4 if M_nyq is None else M_nyq
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are these still the correct defaults? It sounded like N_nyq = NZETA*2?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also you can use setdefault for a cleaner version of this

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are still correct, I think. At least VMEC by default gives M nyq of Mpol-1+4, and Mpol-1 = our Mpol

f0uriest
f0uriest previously approved these changes Aug 30, 2024
@@ -208,7 +208,9 @@ def load(
return eq

@classmethod
def save(cls, eq, path, surfs=128, verbose=1): # noqa: C901 - FIXME - simplify
def save( # noqa: C901 - FIXME - simplify
cls, eq, path, surfs=128, verbose=1, M_nyq=None, N_nyq=None
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I won't force you to change this, but I would prefer to consistently have verbose always be the last kwarg

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I won't just because changing it means changing positional args that others might be using and that would break code. I know we have not committed to an API but for this I will keep as is

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is cool! In a future PR I would be in favor of requiring stuff like verbose to be keyword arguments. That would help with backward-compatibility of the code

desc/vmec.py Show resolved Hide resolved
desc/vmec.py Show resolved Hide resolved
@dpanici dpanici merged commit b683e85 into master Sep 1, 2024
23 of 24 checks passed
@dpanici dpanici deleted the dp/vmecio-asym branch September 1, 2024 15:23
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Saving asymmetric eq as a wout has some errors
5 participants