Merge branch 'master' into dependabot/pip/interpax-fft-gte-0.0.6-and-lte-0.1.0

Author: dpanici

CHANGELOG.md

@@ -8,6 +8,8 @@ Performance Improvements
 Bug Fixes
 
 - Fixes bug that was always setting NFP=1 in ``to_FourierRZ`` methods.
+- Fixes ``VMECIO.save`` metadata for current-density variables and corrects the
+  asymmetric ``currvmns`` magnetic-axis extrapolation.
 
 
 v0.17.2

desc/vmec.py

@@ -53,10 +53,11 @@ def load(
         boundary is loaded and the DESC Equilibrium R,Z are constrained to
         match the given VMEC boundary.
 
-        NOTE: This is only a fit, so the DESC Equilibrium returned is not
-        expected to be in force balance. It is recommended to solve the
-        Equilibrium once loaded before using the Equilibrium for any
-        analysis.
+        NOTE: This is only a fit of the VMEC R, Z and lambda, so the DESC
+        Equilibrium returned is not expected to be in force balance. It is
+        recommended to solve the Equilibrium once loaded before using the
+        Equilibrium for any analysis or computing any derived quantities in
+        DESC (e.g. anything that is not just R, Z or lambda).
 
         Parameters
         ----------
@@ -734,7 +735,7 @@ def save(  # noqa: C901
         jcuru[0] = 0
 
         jcurv = file.createVariable("jcurv", np.float64, ("radius",))
-        jcuru.long_name = "flux surface average of sqrt(g)*J^zeta, on full mesh"
+        jcurv.long_name = "flux surface average of sqrt(g)*J^zeta, on full mesh"
         jcurv.units = "A/m^3"
         jcurv[:] = surface_averages(
             grid_full,
@@ -1331,7 +1332,7 @@ def fullfit(x):
         # TODO (#1379): evaluate current at rho=0 nodes instead of extrapolation
         if not eq.sym:
             currvmns[:, :] = -s
-            currumns[0, :] = -(
+            currvmns[0, :] = -(
                 s[1, :] - (s[2, :] - s[1, :]) / (s_full[2] - s_full[1]) * s_full[1]
             )
         timer.stop("J^zeta*sqrt(g)")

tests/test_vmec.py

@@ -1027,6 +1027,25 @@ def test_vmec_save_asym(VMEC_save_asym):
     np.testing.assert_allclose(
         vmec.variables["jcurv"][20:100], desc.variables["jcurv"][20:100], rtol=2
     )
+    assert desc.variables["jcuru"].long_name == (
+        "flux surface average of sqrt(g)*J^theta, on full mesh"
+    )
+    assert desc.variables["jcurv"].long_name == (
+        "flux surface average of sqrt(g)*J^zeta, on full mesh"
+    )
+    s_full = np.linspace(0, 1, desc.dimensions["radius"].size)
+    expected_currumns_axis = desc.variables["currumns"][1, :] - (
+        (desc.variables["currumns"][2, :] - desc.variables["currumns"][1, :])
+        / (s_full[2] - s_full[1])
+        * s_full[1]
+    )
+    expected_currvmns_axis = desc.variables["currvmns"][1, :] - (
+        (desc.variables["currvmns"][2, :] - desc.variables["currvmns"][1, :])
+        / (s_full[2] - s_full[1])
+        * s_full[1]
+    )
+    np.testing.assert_allclose(desc.variables["currumns"][0, :], expected_currumns_axis)
+    np.testing.assert_allclose(desc.variables["currvmns"][0, :], expected_currvmns_axis)
     np.testing.assert_allclose(
         vmec.variables["DShear"][20:100], desc.variables["DShear"][20:100], rtol=6e-2
     )