Repository: plasmacontrol/desc
Issue: BootstrapRedlConsistency default grid

The default grid for the `BootstrapRedlConsistency` objective includes the $\\rho=1$ surface. Typical density and pressure profiles vanish at the boundary surface; this causes division by 0 in the Redl formula, and results in the objective returning a `NaN` value. 

The easy fix is to change the default grid to exclude the $\\rho=1$ surface. We may also want to do some checking when the objective is built to ensure the profiles and grid are consistent with each other.  

Changed files:
tests/test_bootstrap.py, desc/profiles.py, tests/test_profiles.py, desc/objectives/_bootstrap.py

diff --git a/desc/objectives/_bootstrap.py b/desc/objectives/_bootstrap.py
index ce84f45b..e02ab1ab 100644
--- a/desc/objectives/_bootstrap.py
+++ b/desc/objectives/_bootstrap.py
@@ -77,7 +77,6 @@ class BootstrapRedlConsistency(_Objective):
         helicity=(1, 0),
         name="Bootstrap current self-consistency (Redl)",
     ):
-
         if target is None and bounds is None:
             target = 0
         assert (len(helicity) == 2) and (int(helicity[1]) == helicity[1])
@@ -114,7 +113,7 @@ class BootstrapRedlConsistency(_Objective):
                 N=eq.N_grid,
                 NFP=eq.NFP,
                 sym=eq.sym,
-                rho=np.linspace(1 / 5, 1, 5),
+                rho=np.linspace(1 / eq.L, 1, eq.L) - 1 / (2 * eq.L),
             )
         else:
             grid = self._grid
diff --git a/desc/profiles.py b/desc/profiles.py
index ca7efe67..c6db8b40 100644
--- a/desc/profiles.py
+++ b/desc/profiles.py
@@ -32,7 +32,6 @@ class Profile(IOAble, ABC):
     _io_attrs_ = ["_name", "_grid", "_params"]
 
     def __init__(self, grid=None, name=""):
-
         self.name = name
         self.grid = grid if grid is not None else LinearGrid(L=20)
 
@@ -475,7 +474,6 @@ class ProductProfile(Profile):
     _io_attrs_ = Profile._io_attrs_ + ["_profiles"]
 
     def __init__(self, *profiles, **kwargs):
-
         self._profiles = []
         for profile in profiles:
             assert isinstance(profile, Profile), (
@@ -768,6 +766,9 @@ class PowerSeriesProfile(Profile):
             profile in power series basis fit to given data.
 
         """
+        if sym and sym != "auto":
+            x = x**2
+            order = order // 2
         params = np.polyfit(x, y, order, rcond=rcond, w=w, full=False)[::-1]
         return cls(params, grid=grid, sym=sym, name=name)
 
@@ -801,7 +802,6 @@ class SplineProfile(Profile):
     def __init__(
         self, values=[0, 0, 0], knots=None, grid=None, method="cubic2", name=""
     ):
-
         super().__init__(grid, name)
 
         values = np.atleast_1d(values)
@@ -915,7 +915,6 @@ class MTanhProfile(Profile):
     """
 
     def __init__(self, params=[0, 0, 1, 1, 0], grid=None, name=""):
-
         super().__init__(grid, name)
 
         self._params = params
diff --git a/tests/test_bootstrap.py b/tests/test_bootstrap.py
index e9327c96..455c9060 100644
--- a/tests/test_bootstrap.py
+++ b/tests/test_bootstrap.py
@@ -1620,9 +1620,9 @@ def test_bootstrap_objective_build():
         BootstrapRedlConsistency(eq=eq).build()
 
     eq = Equilibrium(
-        L=3,
-        M=3,
-        N=3,
+        L=4,
+        M=4,
+        N=4,
         NFP=2,
         electron_temperature=1e3,
         electron_density=1e21,
@@ -1631,10 +1631,10 @@ def test_bootstrap_objective_build():
     obj = BootstrapRedlConsistency(eq=eq)
     obj.build()
     # make sure default grid has the right nodes
-    assert obj._transforms["grid"].num_theta == 13
-    assert obj._transforms["grid"].num_zeta == 13
-    assert obj._transforms["grid"].num_rho == 5
+    assert obj._transforms["grid"].num_theta == 17
+    assert obj._transforms["grid"].num_zeta == 17
+    assert obj._transforms["grid"].num_rho == 4
     np.testing.assert_allclose(
         obj._transforms["grid"].nodes[obj._transforms["grid"].unique_rho_idx, 0],
-        np.array([0.2, 0.4, 0.6, 0.8, 1.0]),
+        np.array([0.125, 0.375, 0.625, 0.875]),
     )
diff --git a/tests/test_profiles.py b/tests/test_profiles.py
index 16d4c1d8..3433237e 100644
--- a/tests/test_profiles.py
+++ b/tests/test_profiles.py
@@ -49,7 +49,6 @@ class TestProfiles:
         np.testing.assert_allclose(theta_err, 0, atol=2e-11)
 
     @pytest.mark.unit
-    @pytest.mark.slow
     def test_close_values(self):
         """Test that different forms of the same profile give similar values."""
         pp = PowerSeriesProfile(
@@ -81,6 +80,21 @@ class TestProfiles:
         sp4 = mp.to_spline()
         np.testing.assert_allclose(sp3(x), sp4(x), rtol=1e-5, atol=1e-2)
 
+    @pytest.mark.unit
+    def test_PowerSeriesProfile_even_sym(self):
+        """Test that even symmetry is enforced properly in PowerSeriesProfile."""
+        pp = PowerSeriesProfile(params=np.array([4, 0, -2, 0, 3, 0, -1]), sym="auto")
+        assert pp.sym == "even"  # auto symmetry detects it is even
+        sp = pp.to_spline()
+
+        pp_o = sp.to_powerseries(sym="auto")
+        assert not pp_o.sym  # default conversion from spline is not symmetric
+
+        pp_e = sp.to_powerseries(sym="even")
+        assert pp_e.sym == "even"  # check that even symmetry is enforced
+        # check that this matches the original parameters
+        np.testing.assert_allclose(pp.params, pp_e.params, rtol=1e-3)
+
     @pytest.mark.unit
     def test_SplineProfile_methods(self):
         """Test that all methods of SplineProfile work as intended without errors."""
@@ -346,16 +360,6 @@ class TestProfiles:
             a = sp * pp
             a.params = sp.params
 
-    @pytest.mark.unit
-    def test_profile_conversion(self):
-        """Test converting to FourierZernikeProfile."""
-        pp = PowerSeriesProfile(
-            modes=np.array([0, 1, 2, 4]), params=np.array([1, 0, -2, 1]), sym="auto"
-        )
-        zp = pp.to_fourierzernike(L=6, M=6, N=0)
-        x = np.linspace(0, 1, 100)
-        np.testing.assert_allclose(pp(x), zp(x), atol=1e-10)
-
     @pytest.mark.unit
     def test_default_profiles(self):
         """Test that default profiles are just zeros."""