Fix bug in desc_output_to_input utility function for asymmetric equ…

…ilibria profiles (#1587)

CHANGELOG.md

@@ -47,6 +47,7 @@ Bug Fixes
 - Fixes the coil currents in ``desc.coils.initialize_modular_coils`` to now give the correct expected linking current.
 - ``desc.objectives.PlasmaVesselDistance`` now correctly accounts for multiple field periods on both the equilibrium and the vessel surface. Previously it only considered distances within a single field period.
 - Sets ``os.environ["JAX_PLATFORMS"] = "cpu"`` instead of ``os.environ["JAX_PLATFORM_NAME"] = "cpu"`` when doing ``set_device("cpu")``.
+- Fixes bug in ``desc.input_reader.desc_output_to_input`` utility function for asymmetric equilibria profiles, where the full profile resolution was not being saved.
 - Fixes bug when passing only `sym` into `.change_resolution` for ``FourierRZToroidalSurface``, ``ZernikeRZToroidalSection`` and ``FourierRZCurve``.
 
 Performance Improvements

desc/input_reader.py

@@ -906,12 +906,12 @@ def desc_output_to_input(  # noqa: C901 - fxn too complex
         pres_profile.change_resolution(L=L_profile)
         profile.change_resolution(L=L_profile)
 
-        prof_modes = np.zeros((L_profile, 3))
-        prof_modes[:, 0] = np.arange(L_profile)
+        prof_modes = np.zeros((L_profile + 1, 3))
+        prof_modes[:, 0] = np.arange(L_profile + 1)
         p1 = copy_coeffs(pres_profile.params, pres_profile.basis.modes, prof_modes)
         p2 = copy_coeffs(profile.params, profile.basis.modes, prof_modes)
         f.write("\n# pressure and rotational transform/current profiles\n")
-        for l in range(L_profile):
+        for l in range(L_profile + 1):
             f.write(
                 "l: {:3d}  p = {:15.8E}  {} = {:15.8E}\n".format(
                     int(l), p1[l], char, p2[l]

tests/inputs/iotest_HELIOTRON

@@ -1,45 +1,47 @@
-# This is the DESC input file for a fixed iota heliotron
+# DESC input file generated from the output file:
+# tests/inputs/iotest_HELIOTRON.h5
 
 # global parameters
 sym = 1
 NFP = 19
-Psi = 1.0
+Psi = 1.00000000E+00
 
 # spectral resolution
-L_rad  = 9
-M_pol  =  14
-N_tor  =   3
-
-# continuation parameters
-bdry_ratio = 1
-pres_ratio = 1
-pert_order = 2
+L_rad = 9
+M_pol = 14
+N_tor = 3
+L_grid = 28
+M_grid = 28
+N_grid = 6
 
 # solver tolerances
-ftol = 1e-3
-xtol = 1e-6
-gtol = 1e-6
-maxiter = 50
+ftol = 1.000E-02
+xtol = 1.000E-06
+gtol = 1.000E-06
+maxiter = 100
 
 # solver methods
-optimizer         = lsq-exact
-objective         = force
+optimizer = lsq-exact
+objective = force
 spectral_indexing = ansi
 
 # pressure and rotational transform/current profiles
-l:   0	p =   1E2    i = 0.5
-l:   1	p =   0      i = 0.1
-l:   2	p =   1E2    i = 0
-l:   4	p =   0      i = 0
-
-# magnetic axis initial guess
-n:   0  R0 =  1.0E+1  Z0 =  0.0E+0
+l:   0  p =  1.00000000E+02  i =  5.00000000E-01
+l:   1  p =  1.00000000E-05  i =  1.00000000E-01
+l:   2  p =  1.00000000E+02  i =  0.00000000E+00
+l:   3  p =  0.00000000E+00  i =  0.00000000E+00
+l:   4  p =  0.00000000E+00  i =  0.00000000E+00
+l:   5  p =  0.00000000E+00  i =  0.00000000E+00
+l:   6  p =  0.00000000E+00  i =  0.00000000E+00
+l:   7  p =  0.00000000E+00  i =  0.00000000E+00
+l:   8  p =  0.00000000E+00  i =  0.00000000E+00
+l:   9  p =  0.00000000E+00  i =  0.00000000E+00
 
 # fixed-boundary surface shape
-m:   0  n:   0  R1 =  1.0E+1
-m:   1  n:   0  R1 = -1.0E+0
-m:   1  n:   1  R1 = -3.0E-1
-m:  -1  n:  -1  R1 =  3.0E-1
-m:  -1  n:   0  Z1 =  1.0E+0
-m:  -1  n:   1  Z1 = -3.0E-1
-m:   1  n:  -1  Z1 = -3.0E-1
+l:   0  m:  -1  n:  -1  R1 =  3.00000000E-01  Z1 =  0.00000000E+00
+l:   0  m:   0  n:   0  R1 =  1.00000000E+01  Z1 =  0.00000000E+00
+l:   0  m:   1  n:   0  R1 = -1.00000000E+00  Z1 =  0.00000000E+00
+l:   0  m:   1  n:   1  R1 = -3.00000000E-01  Z1 =  0.00000000E+00
+l:   0  m:   1  n:  -1  R1 =  0.00000000E+00  Z1 = -3.00000000E-01
+l:   0  m:  -1  n:   0  R1 =  0.00000000E+00  Z1 =  1.00000000E+00
+l:   0  m:  -1  n:   1  R1 =  0.00000000E+00  Z1 = -3.00000000E-01

tests/test_input_output.py

@@ -178,12 +178,14 @@ def test_desc_output_to_input(tmpdir_factory):
     tmp_path = tmpdir.join("input_iotest_HELIOTRON")
     tmpout_path = tmpdir.join("iotest_HELIOTRON.h5")
     shutil.copyfile(outfile_path, tmpout_path)
+    eq = load(outfile_path)[-1]
 
     ir1 = InputReader()
     ir1.desc_output_to_input(str(tmp_path), str(tmpout_path))
     ir1 = InputReader(cl_args=[str(tmp_path)])
     arr1 = ir1.parse_inputs()[-1]["surface"]
     pres1 = ir1.parse_inputs()[-1]["pressure"]
+    iota1 = ir1.parse_inputs()[-1]["iota"]
 
     arr1 = arr1[arr1[:, 1].argsort()]
     arr1mneg = arr1[arr1[:, 1] < 0]
@@ -193,6 +195,8 @@ def test_desc_output_to_input(tmpdir_factory):
     ir2 = InputReader(cl_args=[str(desc_input_truth)])
     arr2 = ir2.parse_inputs()[-1]["surface"]
     pres2 = ir2.parse_inputs()[-1]["pressure"]
+    iota2 = ir2.parse_inputs()[-1]["iota"]
+
     arr2 = arr2[arr2[:, 1].argsort()]
     arr2mneg = arr2[arr2[:, 1] < 0]
     arr2mpos = arr2[arr2[:, 1] >= 0]
@@ -215,7 +219,10 @@ def test_desc_output_to_input(tmpdir_factory):
         atol=1e-8,
     )
 
-    np.testing.assert_allclose(pres1[pres1[:, 1] > 0], pres2[pres2[:, 1] > 0])
+    np.testing.assert_allclose(pres1, pres2)
+    np.testing.assert_allclose(pres1[:, 1], eq.pressure.params)
+    np.testing.assert_allclose(iota1, iota2)
+    np.testing.assert_allclose(iota1[:, 1], eq.iota.params)
 
 
 @pytest.mark.unit