diff --git a/desc/compute/_equil.py b/desc/compute/_equil.py index 3f048c34a..c681c530a 100644 --- a/desc/compute/_equil.py +++ b/desc/compute/_equil.py @@ -858,18 +858,12 @@ def _P_ISS04(params, transforms, profiles, data, **kwargs): coordinates="", data=["rho", "ni", "", "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 diff --git a/desc/compute/_profiles.py b/desc/compute/_profiles.py index 0c386e677..290512d85 100644 --- a/desc/compute/_profiles.py +++ b/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 / 1e6 # m^3/s return data diff --git a/desc/objectives/_power_balance.py b/desc/objectives/_power_balance.py index 56a98de2b..14c3e5849 100644 --- a/desc/objectives/_power_balance.py +++ b/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): diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py index ff566c641..2a1915187 100644 --- a/tests/test_objective_funs.py +++ b/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."""