helicity_n now optional in j_dot_B_Redl. In Quasisymmetry, renamed m and n to helicity_m and helicity_n for consistency with QuasisymmetryRatioResidual and j_dot_B_Redl

Author: landreman

src/simsopt/mhd/bootstrap.py

@@ -165,7 +165,7 @@ def integrand(lambd):
     return Bmin, Bmax, epsilon, fsa_B2, fsa_1overB, f_t
 
 
-def j_dot_B_Redl(ne, Te, Ti, Zeff, helicity_n, s=None, G=None, R=None, iota=None,
+def j_dot_B_Redl(ne, Te, Ti, Zeff, helicity_n=None, s=None, G=None, R=None, iota=None,
                  epsilon=None, f_t=None, psi_edge=None, nfp=None,
                  geom=None, plot=False):
     r"""
@@ -199,6 +199,10 @@ def j_dot_B_Redl(ne, Te, Ti, Zeff, helicity_n, s=None, G=None, R=None, iota=None
     :math:`\left<\vec{J}\cdot\vec{B}\right>` will be computed on this
     same set of flux surfaces.
 
+    If you provide a :obj:`RedlGeomBoozer` object for ``geom``, then
+    it is not necessary to specify the argument ``helicity_n`` here,
+    in which case ``helicity_n`` will be taken from ``geom``.
+
     Args:
         ne: A :obj:`~simsopt.mhd.profiles.Profile` object with the electron density profile.
         Te: A :obj:`~simsopt.mhd.profiles.Profile` object with the electron temperature profile.
@@ -247,6 +251,9 @@ def j_dot_B_Redl(ne, Te, Ti, Zeff, helicity_n, s=None, G=None, R=None, iota=None
         f_t = geom_data.f_t
         nfp = geom_data.nfp
 
+    if helicity_n is None:
+        helicity_n = geom.helicity_n
+
     helicity_N = nfp * helicity_n
     if Zeff is None:
         Zeff = ProfilePolynomial(1.0)

src/simsopt/mhd/boozer.py

@@ -243,10 +243,10 @@ class Quasisymmetry(Optimizable):
     Args:
         boozer: A Boozer object on which the calculation will be based.
         s: The normalized toroidal magnetic flux for the flux surface to analyze. Should be in the range [0, 1].
-        m: The poloidal mode number of the symmetry you want to achive.
-           The departure from symmetry B(m * theta - nfp * n * zeta) will be reported.
-        n: The toroidal mode number of the symmetry you want to achieve.
-           The departure from symmetry B(m * theta - nfp * n * zeta) will be reported.
+        helicity_m: The poloidal mode number of the symmetry you want to achive.
+           The departure from symmetry ``B(helicity_m * theta - nfp * helicity_n * zeta)`` will be reported.
+        helicity_n: The toroidal mode number of the symmetry you want to achieve.
+           The departure from symmetry ``B(helicity_m * theta - nfp * helicity_n * zeta)`` will be reported.
         normalization: A uniform normalization applied to all bmnc harmonics.
            If ``"B00"``, the symmetry-breaking modes will be divided by the m=n=0 mode amplitude
            on the same surface. If ``"symmetric"``, the symmetry-breaking modes will be
@@ -258,17 +258,17 @@ class Quasisymmetry(Optimizable):
     def __init__(self,
                  boozer: Boozer,
                  s: Union[float, Iterable[float]],
-                 m: int,
-                 n: int,
+                 helicity_m: int,
+                 helicity_n: int,
                  normalization: str = "B00",
                  weight: str = "even") -> None:
         """
         Constructor
 
         """
         self.boozer = boozer
-        self.m = m
-        self.n = n
+        self.helicity_m = helicity_m
+        self.helicity_n = helicity_n
         self.normalization = normalization
         self.weight = weight
 
@@ -308,23 +308,23 @@ def J(self) -> np.ndarray:
             xm = self.boozer.bx.xm_b
             xn = self.boozer.bx.xn_b / self.boozer.bx.nfp
 
-            if self.m != 0 and self.m != 1:
+            if self.helicity_m != 0 and self.helicity_m != 1:
                 raise ValueError("m for quasisymmetry should be 0 or 1.")
 
             # Find the indices of the symmetric modes:
-            if self.n == 0:
+            if self.helicity_n == 0:
                 # Quasi-axisymmetry
                 symmetric = (xn == 0)
 
-            elif self.m == 0:
+            elif self.helicity_m == 0:
                 # Quasi-poloidal symmetry
                 symmetric = (xm == 0)
 
             else:
                 # Quasi-helical symmetry
-                symmetric = (xm * self.n + xn * self.m == 0)
-                # Stellopt takes the "and" of this with mod(xm, self.m),
-                # which does not seem necessary since self.m must be 1 to
+                symmetric = (xm * self.helicity_n + xn * self.helicity_m == 0)
+                # Stellopt takes the "and" of this with mod(xm, self.helicity_m),
+                # which does not seem necessary since self.helicity_m must be 1 to
                 # get here.
             nonsymmetric = np.logical_not(symmetric)
 

tests/mhd/test_bootstrap.py

@@ -636,10 +636,11 @@ def test_Redl_sfincs_tokamak_benchmark(self):
         for geom_method in range(2):
             if geom_method == 0:
                 geom = RedlGeomVmec(vmec, surfaces)
+                jdotB, details = j_dot_B_Redl(ne, Te, Ti, Zeff, helicity_n, geom=geom, plot=False)
             else:
                 geom = RedlGeomBoozer(boozer, surfaces, helicity_n)
+                jdotB, details = j_dot_B_Redl(ne, Te, Ti, Zeff, geom=geom, plot=False)
 
-            jdotB, details = j_dot_B_Redl(ne, Te, Ti, Zeff, helicity_n, geom=geom, plot=False)
             jdotB_history.append(jdotB)
 
             # The relative error is a bit larger at s \approx 1, where the