fusion power objective

desc/compute/_equil.py

@@ -10,7 +10,7 @@
 """
 
 from interpax import interp1d
-from scipy.constants import mu_0
+from scipy.constants import elementary_charge, mu_0
 
 from desc.backend import jnp
 
@@ -843,3 +843,39 @@ def _P_ISS04(params, transforms, profiles, data, **kwargs):
         )
     ) ** (1 / 0.39)
     return data
+
+
+@register_compute_fun(
+    name="P_fusion",
+    label="P_{fusion}",
+    units="W",
+    units_long="Watts",
+    description="Fusion power",
+    dim=0,
+    params=[],
+    transforms={"grid": []},
+    profiles=[],
+    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'.",
+)
+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
+
+    reaction_rate = jnp.sum(
+        *data["ni"] ** 2
+        * data["<sigma*nu>"]
+        * data["sqrt(g)"]
+        * transforms["grid"].weights
+    )  # reactions/s
+    data["P_fusion"] = reaction_rate * energy * elementary_charge  # J/s
+    return data

desc/compute/_profiles.py

@@ -1909,3 +1909,65 @@ def _shear(params, transforms, profiles, data, **kwargs):
         None,
     )
     return data
+
+
+@register_compute_fun(
+    name="<sigma*nu>",
+    label="\\langle\\sigma\\nu\\rangle",
+    units="m^3 \\cdot s^{-1}",
+    units_long="cubic meters / second",
+    description="Thermal reactivity from Bosch-Hale parameterization",
+    dim=1,
+    params=["Ti_l"],
+    transforms={"grid": []},
+    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'.",
+)
+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
+
+    T = data["Ti"] / 1e3  # keV
+    theta = T / (
+        1 - (T * (C2 + T * (C4 + T * C6))) / (1 + T * (C3 + T * (C5 + T * 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
+    data["<sigma*nu>"] = sigma_nu / 1e6  # m^3/s
+    return data

desc/objectives/_confinement.py

@@ -9,6 +9,179 @@
 from .objective_funs import _Objective
 
 
+class FusionPower(_Objective):
+    """Fusion power.
+
+    P = e E ∫ n^2 ⟨σν⟩ dV (W)
+
+    References
+    ----------
+    https://doi.org/10.1088/0029-5515/32/4/I07.
+    Improved Formulas for Fusion Cross-Sections and Thermal Reactivities.
+    H.-S. Bosch and G.M. Hale. Nucl. Fusion April 1992; 32 (4): 611-631.
+
+    Parameters
+    ----------
+    eq : Equilibrium
+        Equilibrium that will be optimized to satisfy the Objective.
+    target : {float, ndarray}, optional
+        Target value(s) of the objective. Only used if bounds is None.
+        Must be broadcastable to Objective.dim_f. Defaults to ``target=0``.
+    bounds : tuple of {float, ndarray}, optional
+        Lower and upper bounds on the objective. Overrides target.
+        Both bounds must be broadcastable to to Objective.dim_f.
+        Defaults to ``target=0``.
+    weight : {float, ndarray}, optional
+        Weighting to apply to the Objective, relative to other Objectives.
+        Must be broadcastable to to Objective.dim_f
+    normalize : bool, optional
+        Whether to compute the error in physical units or non-dimensionalize.
+    normalize_target : bool, optional
+        Whether target and bounds should be normalized before comparing to computed
+        values. If `normalize` is `True` and the target is in physical units,
+        this should also be set to True.
+    loss_function : {None, 'mean', 'min', 'max'}, optional
+        Loss function to apply to the objective values once computed. This loss function
+        is called on the raw compute value, before any shifting, scaling, or
+        normalization. Note: Has no effect for this objective.
+    deriv_mode : {"auto", "fwd", "rev"}
+        Specify how to compute jacobian matrix, either forward mode or reverse mode AD.
+        "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'.
+    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)``.
+    name : str, optional
+        Name of the objective function.
+
+    """
+
+    _scalar = True
+    _units = "(W)"
+    _print_value_fmt = "Fusion power: {:10.3e} "
+
+    def __init__(
+        self,
+        eq,
+        target=None,
+        bounds=None,
+        weight=1,
+        normalize=True,
+        normalize_target=True,
+        loss_function=None,
+        deriv_mode="auto",
+        reaction="T(d,n)4He",
+        grid=None,
+        name="fusion power",
+    ):
+        if target is None and bounds is None:
+            target = 0
+        self._reaction = reaction
+        self._grid = grid
+        super().__init__(
+            things=eq,
+            target=target,
+            bounds=bounds,
+            weight=weight,
+            normalize=normalize,
+            normalize_target=normalize_target,
+            loss_function=loss_function,
+            deriv_mode=deriv_mode,
+            name=name,
+        )
+
+    def build(self, use_jit=True, verbose=1):
+        """Build constant arrays.
+
+        Parameters
+        ----------
+        use_jit : bool, optional
+            Whether to just-in-time compile the objective and derivatives.
+        verbose : int, optional
+            Level of output.
+
+        """
+        eq = self.things[0]
+        errorif(
+            eq.ion_density is None,
+            ValueError,
+            "Equilibrium must have an ion density profile.",
+        )
+        if self._grid is None:
+            self._grid = QuadratureGrid(
+                L=eq.L_grid,
+                M=eq.M_grid,
+                N=eq.N_grid,
+                NFP=eq.NFP,
+            )
+        self._dim_f = 1
+        self._data_keys = ["P_fusion"]
+
+        timer = Timer()
+        if verbose > 0:
+            print("Precomputing transforms")
+        timer.start("Precomputing transforms")
+
+        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,
+        }
+
+        timer.stop("Precomputing transforms")
+        if verbose > 1:
+            timer.disp("Precomputing transforms")
+
+        if self._normalize:
+            scales = compute_scaling_factors(eq)
+            self._normalization = scales["W_p"]
+
+        super().build(use_jit=use_jit, verbose=verbose)
+
+    def compute(self, params, constants=None):
+        """Compute fusion power.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary of equilibrium or surface degrees of freedom, eg
+            Equilibrium.params_dict
+        constants : dict
+            Dictionary of constant data, eg transforms, profiles etc. Defaults to
+            self.constants
+
+        Returns
+        -------
+        P : float
+            Fusion power (W).
+
+        """
+        if constants is None:
+            constants = self.constants
+        data = compute_fun(
+            "desc.equilibrium.equilibrium.Equilibrium",
+            self._data_keys,
+            params=params,
+            transforms=constants["transforms"],
+            profiles=constants["profiles"],
+            reaction=constants["reaction"],
+        )
+        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
+
+    @reaction.setter
+    def reaction(self, new):
+        self._reaction = new
+
+
 class HeatingPowerISS04(_Objective):
     """Heating power required by the ISS04 energy confinement time scaling.