replace reaction with fuel, only DT

desc/compute/_equil.py

@@ -858,18 +858,12 @@ def _P_ISS04(params, transforms, profiles, data, **kwargs):
     coordinates="",
     data=["rho", "ni", "<sigma*nu>", "sqrt(g)"],
     resolution_requirement="rtz",
-    reaction="str: Fusion reaction. One of {'T(d,n)4He', 'D(d,p)T', 'D(d,n)3He'}. "
-    + "Default is 'T(d,n)4He'.",
+    fuel="str: Fusion fuel, assuming a 50/50 mix. One of {'DT'}. Default is 'DT'.",
 )
 def _P_fusion(params, transforms, profiles, data, **kwargs):
-    reaction = kwargs.get("fuel", "T(d,n)4He")
-    match reaction:
-        case "T(d,n)4He":
-            energy = 3.52e6 + 14.06e6  # eV
-        case "D(d,p)T":
-            energy = 1.01e6 + 3.02e6  # eV
-        case "D(d,n)3He":
-            energy = 0.82e6 + 2.45e6  # eV
+    energies = {"DT": 3.52e6 + 14.06e6}  # eV
+    fuel = kwargs.get("fuel", "DT")
+    energy = energies.get(fuel)
 
     reaction_rate = jnp.sum(
         data["ni"] ** 2

desc/compute/_profiles.py

@@ -1923,51 +1923,37 @@ def _shear(params, transforms, profiles, data, **kwargs):
     profiles=["ion_temperature"],
     coordinates="r",
     data=["Ti"],
-    reaction="str: Fusion reaction. One of {'T(d,n)4He', 'D(d,p)T', 'D(d,n)3He'}. "
-    + "Default is 'T(d,n)4He'.",
+    fuel="str: Fusion fuel, assuming a 50/50 mix. One of {'DT'}. Default is 'DT'.",
 )
 def _reactivity(params, transforms, profiles, data, **kwargs):
     # Bosch and Hale. “Improved Formulas for Fusion Cross-Sections and Thermal
     # Reactivities.” Nuclear Fusion 32 (April 1992): 611-631.
     # https://doi.org/10.1088/0029-5515/32/4/I07.
-    reaction = kwargs.get("fuel", "T(d,n)4He")
-    match reaction:
-        case "T(d,n)4He":
-            B_G = 34.382
-            mc2 = 1124656
-            C1 = 1.17302e-9
-            C2 = 1.51361e-2
-            C3 = 7.51886e-2
-            C4 = 4.60643e-3
-            C5 = 1.35000e-2
-            C6 = -1.06750e-4
-            C7 = 1.36600e-5
-        case "D(d,p)T":
-            B_G = 31.3970
-            mc2 = 937814
-            C1 = 5.65718e-12
-            C2 = 3.41267e-3
-            C3 = 1.99167e-3
-            C4 = 0
-            C5 = 1.05060e-5
-            C6 = 0
-            C7 = 0
-        case "D(d,n)3He":
-            B_G = 31.3970
-            mc2 = 937814
-            C1 = 5.43360e-12
-            C2 = 5.85778e-3
-            C3 = 7.68222e-3
-            C4 = 0
-            C5 = -2.96400e-6
-            C6 = 0
-            C7 = 0
+    coefficients = {
+        "DT": {
+            "B_G": 34.382,
+            "mc2": 1124656,
+            "C1": 1.17302e-9,
+            "C2": 1.51361e-2,
+            "C3": 7.51886e-2,
+            "C4": 4.60643e-3,
+            "C5": 1.35000e-2,
+            "C6": -1.06750e-4,
+            "C7": 1.36600e-5,
+        }
+    }
+    fuel = kwargs.get("fuel", "DT")
+    coeffs = coefficients.get(fuel)
 
     T = data["Ti"] / 1e3  # keV
     theta = T / (
-        1 - (T * (C2 + T * (C4 + T * C6))) / (1 + T * (C3 + T * (C5 + T * C7)))
+        1
+        - (T * (coeffs["C2"] + T * (coeffs["C4"] + T * coeffs["C6"])))
+        / (1 + T * (coeffs["C3"] + T * (coeffs["C5"] + T * coeffs["C7"])))
     )
-    xi = (B_G**2 / (4 * theta)) ** (1 / 3)
-    sigma_nu = C1 * theta * jnp.sqrt(xi / (mc2 * T**3)) * jnp.exp(-3 * xi)  # cm^3/s
+    xi = (coeffs["B_G"] ** 2 / (4 * theta)) ** (1 / 3)
+    sigma_nu = (
+        coeffs["C1"] * theta * jnp.sqrt(xi / (coeffs["mc2"] * T**3)) * jnp.exp(-3 * xi)
+    )  # cm^3/s
     data["<sigma*nu>"] = sigma_nu / 1e6  # m^3/s
     return data

desc/objectives/_power_balance.py

@@ -49,9 +49,8 @@ class FusionPower(_Objective):
         "auto" selects forward or reverse mode based on the size of the input and output
         of the objective. Has no effect on self.grad or self.hess which always use
         reverse mode and forward over reverse mode respectively.
-    reaction : str, optional
-        Fusion reaction. One of {'T(d,n)4He', 'D(d,p)T', 'D(d,n)3He'}.
-        Default = 'T(d,n)4He'.
+    fuel : str, optional
+        Fusion fuel, assuming a 50/50 mix. One of {'DT'}. Default = 'DT'.
     grid : Grid, optional
         Collocation grid used to compute the intermediate quantities.
         Defaults to ``QuadratureGrid(eq.L_grid, eq.M_grid, eq.N_grid, eq.NFP)``.
@@ -74,13 +73,16 @@ def __init__(
         normalize_target=True,
         loss_function=None,
         deriv_mode="auto",
-        reaction="T(d,n)4He",
+        fuel="DT",
         grid=None,
         name="fusion power",
     ):
+        errorif(
+            fuel not in ["DT"], ValueError, f"fuel must be one of ['DT'], got {fuel}."
+        )
         if target is None and bounds is None:
             target = 0
-        self._reaction = reaction
+        self._fuel = fuel
         self._grid = grid
         super().__init__(
             things=eq,
@@ -107,9 +109,9 @@ def build(self, use_jit=True, verbose=1):
         """
         eq = self.things[0]
         errorif(
-            eq.ion_density is None,
+            eq.electron_density is None,
             ValueError,
-            "Equilibrium must have an ion density profile.",
+            "Equilibrium must have an electron density profile.",
         )
         if self._grid is None:
             self._grid = QuadratureGrid(
@@ -129,7 +131,7 @@ def build(self, use_jit=True, verbose=1):
         self._constants = {
             "profiles": get_profiles(self._data_keys, obj=eq, grid=self._grid),
             "transforms": get_transforms(self._data_keys, obj=eq, grid=self._grid),
-            "reaction": self._reaction,
+            "fuel": self._fuel,
         }
 
         timer.stop("Precomputing transforms")
@@ -168,18 +170,21 @@ def compute(self, params, constants=None):
             params=params,
             transforms=constants["transforms"],
             profiles=constants["profiles"],
-            reaction=constants["reaction"],
+            fuel=constants["fuel"],
         )
         return data["P_fusion"]
 
     @property
-    def reaction(self):
-        """str: Fusion reaction. One of {'T(d,n)4He', 'D(d,p)T', 'D(d,n)3He'}."""
-        return self._reaction
+    def fuel(self):
+        """str: Fusion fuel, assuming a 50/50 mix. One of {'DT'}. Default = 'DT'."""
+        return self._fuel
 
-    @reaction.setter
-    def reaction(self, new):
-        self._reaction = new
+    @fuel.setter
+    def fuel(self, new):
+        errorif(
+            new not in ["DT"], ValueError, f"fuel must be one of ['DT'], got {new}."
+        )
+        self._fuel = new
 
 
 class HeatingPowerISS04(_Objective):

tests/test_objective_funs.py

@@ -48,6 +48,7 @@
     Energy,
     ForceBalance,
     ForceBalanceAnisotropic,
+    FusionPower,
     GenericObjective,
     HeatingPowerISS04,
     Isodynamicity,
@@ -2031,6 +2032,7 @@ class TestComputeScalarResolution:
         CoilLength,
         CoilSetMinDistance,
         CoilTorsion,
+        FusionPower,
         GenericObjective,
         HeatingPowerISS04,
         Omnigenity,
@@ -2092,6 +2094,28 @@ def test_compute_scalar_resolution_bootstrap(self):
             f[i] = obj.compute_scalar(obj.x())
         np.testing.assert_allclose(f, f[-1], rtol=5e-2)
 
+    @pytest.mark.regression
+    def test_compute_scalar_resolution_fusion_power(self):
+        """FusionPower."""
+        eq = self.eq.copy()
+        eq.electron_density = PowerSeriesProfile([1e19, 0, -1e19])
+        eq.electron_temperature = PowerSeriesProfile([1e3, 0, -1e3])
+        eq.ion_temperature = PowerSeriesProfile([1e3, 0, -1e3])
+        eq.atomic_number = 1.0
+
+        f = np.zeros_like(self.res_array, dtype=float)
+        for i, res in enumerate(self.res_array):
+            grid = QuadratureGrid(
+                L=int(self.eq.L * res),
+                M=int(self.eq.M * res),
+                N=int(self.eq.N * res),
+                NFP=self.eq.NFP,
+            )
+            obj = ObjectiveFunction(FusionPower(eq=eq, grid=grid))
+            obj.build(verbose=0)
+            f[i] = obj.compute_scalar(obj.x())
+        np.testing.assert_allclose(f, f[-1], rtol=5e-2)
+
     @pytest.mark.regression
     def test_compute_scalar_resolution_heating_power(self):
         """HeatingPowerISS04."""
@@ -2383,6 +2407,7 @@ class TestObjectiveNaNGrad:
         CoilSetMinDistance,
         CoilTorsion,
         ForceBalanceAnisotropic,
+        FusionPower,
         HeatingPowerISS04,
         Omnigenity,
         PlasmaCoilSetMinDistance,
@@ -2432,6 +2457,22 @@ def test_objective_no_nangrad_bootstrap(self):
         g = obj.grad(obj.x(eq))
         assert not np.any(np.isnan(g)), "redl bootstrap"
 
+    @pytest.mark.unit
+    def test_objective_no_nangrad_fusion_power(self):
+        """FusionPower."""
+        eq = Equilibrium(
+            L=2,
+            M=2,
+            N=2,
+            electron_density=PowerSeriesProfile([1e19, 0, -1e19]),
+            electron_temperature=PowerSeriesProfile([1e3, 0, -1e3]),
+            current=PowerSeriesProfile([1, 0, -1]),
+        )
+        obj = ObjectiveFunction(FusionPower(eq))
+        obj.build()
+        g = obj.grad(obj.x(eq))
+        assert not np.any(np.isnan(g)), "fusion power"
+
     @pytest.mark.unit
     def test_objective_no_nangrad_heating_power(self):
         """HeatingPowerISS04."""