Merge pull request #711 from PlasmaControl/dp/io-hotfix

Continuation _solve_axisym and Surface from_input_file hotfix

desc/continuation.py

@@ -39,7 +39,7 @@ def _solve_axisym(
     L, M, N, L_grid, M_grid, N_grid = eq.L, eq.M, eq.N, eq.L_grid, eq.M_grid, eq.N_grid
     spectral_indexing = eq.spectral_indexing
 
-    Mi = min(M // 2, mres_step) if mres_step > 0 else M
+    Mi = min(M, mres_step) if mres_step > 0 else M
     Li = int(np.ceil(L / M) * Mi)
     Ni = 0
     L_gridi = np.ceil(L_grid / L * Li).astype(int)

desc/geometry/surface.py

@@ -272,7 +272,12 @@ def from_input_file(cls, path):
 
         """
         f = open(path)
-        if "&INDATA" in f.readlines()[0].upper():  # vmec input, convert to desc
+        isVMEC = False
+        for line in f.readlines():
+            if "&INDATA" in line.upper():
+                isVMEC = True
+                break
+        if isVMEC:  # vmec input, convert to desc
             inputs = InputReader.parse_vmec_inputs(f)[-1]
         else:
             inputs = InputReader().parse_inputs(f)[-1]

tests/inputs/input.DSHAPE

@@ -1,3 +1,4 @@
+! Comment line before INDATA start
 &indata
   ! DSHAPE Tokamak from Hirshman and Whitson 1983
   MGRID_FILE = 'none'

tests/test_examples.py

@@ -10,7 +10,7 @@
 from qsc import Qsc
 
 import desc.examples
-from desc.continuation import solve_continuation_automatic
+from desc.continuation import _solve_axisym, solve_continuation_automatic
 from desc.equilibrium import EquilibriaFamily, Equilibrium
 from desc.geometry import FourierRZToroidalSurface
 from desc.grid import LinearGrid
@@ -52,6 +52,15 @@ def test_SOLOVEV_vacuum(SOLOVEV_vac):
     np.testing.assert_allclose(data["iota"], 0, atol=1e-16)
     np.testing.assert_allclose(data["|J|"], 0, atol=3e-3)
 
+    # test that solving with the continuation method works correctly
+    # when eq resolution is lower than the mres_step
+    eq.change_resolution(L=3, M=3)
+    eqf = _solve_axisym(eq, mres_step=6)
+    assert len(eqf) == 1
+    assert eqf[-1].L == eq.L
+    assert eqf[-1].M == eq.M
+    assert eqf[-1].N == eq.N
+
 
 @pytest.mark.regression
 @pytest.mark.solve