added reactions func

pyproject.toml

@@ -23,7 +23,7 @@ keywords = ["python", "neutron", "fusion", "source", "openmc", "energy", "tokama
 dependencies = [
     "numpy>=1.9",
     "matplotlib>=3.2.2",
-
+    "NeSST>=1.0.0"
 ]
 
 [project.urls]
@@ -37,7 +37,7 @@ write_to = "src/_version.py"
 tests = [
     "pytest>=5.4.3",
     "hypothesis",
-    "NeSST"
+    "NeSST>=1.0.0"
 ]
 
 [tool.setuptools]

src/openmc_plasma_source/fuel_types.py

@@ -38,8 +38,8 @@ def neutron_energy_mean(ion_temperature: float, reaction: str) -> float:
     ion_temperature_kev = ion_temperature / 1e3  # Ballabio equation accepts KeV units
     mean_delta = (
         a_1
-        * ion_temperature_kev ** (2./3.)
-        / (1. + a_2 * ion_temperature_kev**a_3)
+        * ion_temperature_kev ** (2.0 / 3.0)
+        / (1.0 + a_2 * ion_temperature_kev**a_3)
         + a_4 * ion_temperature_kev
     )
     mean_delta *= 1e3  # converting back to eV
@@ -81,18 +81,30 @@ def neutron_energy_std_dev(ion_temperature: float, reaction: str) -> float:
     ion_temperature_kev = ion_temperature / 1e3  # Ballabio equation accepts KeV units
     delta = (
         a_1
-        * ion_temperature_kev ** (2./3.)
-        / (1. + a_2 * ion_temperature_kev**a_3)
+        * ion_temperature_kev ** (2.0 / 3.0)
+        / (1.0 + a_2 * ion_temperature_kev**a_3)
         + a_4 * ion_temperature_kev
     )
 
     # 2.3548200450309493 on the line below comes from equation 2* math.sqrt(math.log(2)*2)
-    variance = ((w_0 * (1 + delta))**2 * ion_temperature_kev) / 2.3548200450309493**2
+    variance = ((w_0 * (1 + delta)) ** 2 * ion_temperature_kev) / 2.3548200450309493**2
     variance *= 1e6  # converting keV^2 back to eV^2
     std_dev = np.sqrt(variance)
     return std_dev
 
 
+
+def get_reactions_from_fuel(fuel):
+    if ["D", "T"] == sorted(set(fuel.keys())):
+        return ['DT', 'DD', 'TT']
+    elif ["D"] == sorted(set(fuel.keys())):
+        return ['DD']
+    elif ["T"] == sorted(set(fuel.keys())):
+        return ['TT']
+    else:
+        msg = 'reactions of fuel {fuel} could not be found. Supported fuel keys are "T" and "D"'
+        raise ValueError(msg)
+
 def get_neutron_energy_distribution(
     ion_temperature: float,
     fuel: dict,
@@ -112,8 +124,6 @@ def get_neutron_energy_distribution(
         energy distribution
     """
 
-    ion_temperature_kev = ion_temperature / 1e3  # convert eV to keV
-
     sum_fuel_isotopes = sum(fuel.values())
     if sum_fuel_isotopes > 1.0:
         raise ValueError(
@@ -141,26 +151,46 @@ def get_neutron_energy_distribution(
     # 1.0 neutron yield, all reactions scaled by this value
     num_of_vals = 100
     # single grid for TT neutrons
-    E_pspec = np.linspace(0, 12, num_of_vals)  # E_pspec is exspected in MeV units
+    E_pspec = np.linspace(0, 12e6, num_of_vals)  # E_pspec is exspected in MeV units
 
     DDmean = neutron_energy_mean(ion_temperature=ion_temperature, reaction="DD")
     DTmean = neutron_energy_mean(ion_temperature=ion_temperature, reaction="DT")
     DD_std_dev = neutron_energy_std_dev(ion_temperature=ion_temperature, reaction="DD")
     DT_std_dev = neutron_energy_std_dev(ion_temperature=ion_temperature, reaction="DT")
 
-    if ["D", "T"] == sorted(set(fuel.keys())):
+    reactions = get_reactions_from_fuel(fuel)
+
+    if reactions == ['TT']:
+        strength_TT = 1.0
+        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature)
+        tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT)
+        return [tt_source], [strength_TT]
+
+    elif reactions == ["DD"]:
+        strength_DD = 1.0
+        dd_source = openmc.stats.Normal(mean_value=DDmean, std_dev=DD_std_dev)
+        return [dd_source], [strength_DD]
+
+    # DT, DD and TT reaction
+    else:
         strength_DT = 1.0
         strength_DD = nst.yield_from_dt_yield_ratio(
-            "dd", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"]
+            "dd", strength_DT, ion_temperature, fuel["D"], fuel["T"]
         )
         strength_TT = nst.yield_from_dt_yield_ratio(
-            "tt", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"]
+            "tt", strength_DT, ion_temperature, fuel["D"], fuel["T"]
         )
 
         total_strength = sum([strength_TT, strength_DD, strength_DT])
 
-        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev)
-        tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT)
+        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature)
+
+        # removing any zeros from the end of the array
+        dNdE_TT = np.trim_zeros(dNdE_TT, 'b')
+        # making array lengths match
+        E_pspec = E_pspec[:len(dNdE_TT)]
+
+        tt_source = openmc.stats.Tabular(E_pspec, dNdE_TT)
 
         dd_source = openmc.stats.Normal(mean_value=DDmean, std_dev=DD_std_dev)
         # normal could be done with Muir but in this case we have the mean and std dev from NeSST
@@ -176,14 +206,3 @@ def get_neutron_energy_distribution(
             strength_DD / total_strength,
             strength_DT / total_strength,
         ]
-
-    elif ["D"] == sorted(set(fuel.keys())):
-        strength_DD = 1.0
-        dd_source = openmc.stats.muir(e0=DDmean * 1e6, m_rat=4, kt=ion_temperature)
-        return [dd_source], [strength_DD]
-
-    elif ["T"] == sorted(set(fuel.keys())):
-        strength_TT = 1.0
-        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev)
-        tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT)
-        return [tt_source], [strength_TT]

src/openmc_plasma_source/point_source.py

@@ -9,7 +9,7 @@ def fusion_point_source(
     coordinate: Tuple[float, float, float] = (0.0, 0.0, 0.0),
     temperature: float = 20000.0,
     fuel: dict = {"D": 0.5, "T": 0.5},
-):
+) -> list[openmc.IndependentSource]:
     """Creates a list of openmc.IndependentSource objects representing an ICF source.
 
     Resulting ICF (Inertial Confinement Fusion) source will have an energy
@@ -42,14 +42,6 @@ def fusion_point_source(
         ion_temperature=temperature, fuel=fuel
     )
 
-    # if isinstance(energy_distributions, openmc.stats.Normal) or isinstance(energy_distributions, openmc.stats.Discrete) or isinstance(energy_distributions, openmc.stats.Tabular):
-    #     source = openmc.Source()
-    #     source.energy = energy_distributions
-    #     source.space = openmc.stats.Point(coordinate)
-    #     source.angle = openmc.stats.Isotropic()
-    #     return source
-
-    # else:
     for energy_distribution, strength in zip(energy_distributions, strengths):
         source = openmc.IndependentSource()
         source.energy = energy_distribution

src/openmc_plasma_source/tokamak_source.py

@@ -3,7 +3,7 @@
 import numpy as np
 import openmc
 import openmc.checkvalue as cv
-
+import NeSST as nst
 from .fuel_types import get_neutron_energy_distribution
 
 
@@ -118,8 +118,7 @@ def tokamak_source(
     a = np.random.random(sample_size) * minor_radius
     alpha = np.random.random(sample_size) * 2 * np.pi
 
-    # compute densities, temperatures, neutron source densities and
-    # convert coordinates
+    # compute densities, temperatures
     densities = tokamak_ion_density(
         mode=mode,
         ion_density_centre=ion_density_centre,
@@ -130,6 +129,12 @@ def tokamak_source(
         ion_density_separatrix=ion_density_separatrix,
         r=a,
     )
+
+    fuel_densities = {}
+    for key, value in fuel.items():
+        fuel_densities[key] = densities*value
+
+    # compute temperatures
     temperatures = tokamak_ion_temperature(
         r=a,
         mode=mode,
@@ -142,10 +147,7 @@ def tokamak_source(
         major_radius=major_radius,
     )
 
-    neutron_source_density = tokamak_neutron_source_density(densities, temperatures)
-
-    strengths = neutron_source_density / sum(neutron_source_density)
-
+    # convert coordinates
     RZ = tokamak_convert_a_alpha_to_R_Z(
         a=a,
         alpha=alpha,
@@ -156,6 +158,11 @@ def tokamak_source(
         elongation=elongation,
     )
 
+    neutron_source_density = tokamak_neutron_source_density(fuel_densities, temperatures, reaction)
+
+    strengths = neutron_source_density / sum(neutron_source_density)
+
+
     sources = tokamak_make_openmc_sources(
         strengths=strengths,
         angles=angles,
@@ -313,8 +320,12 @@ def tokamak_make_openmc_sources(
     is based on .strengths.
 
     Args:
+        strengths
         angles ((float, float), optional): rotation of the ring source.
-        Defaults to (0, 2*np.pi).
+            Defaults to (0, 2*np.pi).
+        temperatures
+        fuel
+        RZ
 
     Returns:
         list: list of openmc.IndependentSource()
@@ -357,7 +368,7 @@ def tokamak_make_openmc_sources(
     return sources
 
 
-def tokamak_neutron_source_density(ion_density, ion_temperature):
+def tokamak_neutron_source_density(ion_density, ion_temperature, reaction):
     """Computes the neutron source density given ion density and ion
     temperature.
 
@@ -369,25 +380,27 @@ def tokamak_neutron_source_density(ion_density, ion_temperature):
         float, ndarray: Neutron source density (neutron/s/m3)
     """
 
-    ion_density = np.asarray(ion_density)
+    ion_density = np.asarray(ion_density[reaction[0]]) * np.asarray(ion_density[reaction[1]])
     ion_temperature = np.asarray(ion_temperature)
 
-    return ion_density**2 * _DT_xs(ion_temperature)
+    if reaction == 'DT':
+        return ion_density * nst.reac_DT(ion_temperature) # could use _DT_xs instead
+    if reaction == 'DD':
+        return ion_density * nst.reac_DD(ion_temperature)
+    if reaction == 'TT':
+        return ion_density * nst.reac_TT(ion_temperature)
 
 
+#TODO consider replace with NeSST or getting DD version as well
 def _DT_xs(ion_temperature):
     """Sadler–Van Belle formula
     Ref : https://doi.org/10.1016/j.fusengdes.2012.02.025
-
     Args:
-        ion_temperature (float, ndarray): ion temperature in keV
-
+        ion_temperature (float, ndarray): ion temperature in eV
     Returns:
         float, ndarray: the DT cross section at the given temperature
     """
-
-    ion_temperature = np.asarray(ion_temperature)
-
+    ion_temperature_kev = np.asarray(ion_temperature/1e3)
     c = [
         2.5663271e-18,
         19.983026,
@@ -397,14 +410,12 @@ def _DT_xs(ion_temperature):
         6.6024089e-2,
         8.1215505e-3,
     ]
-
-    U = 1 - ion_temperature * (
-        c[2] + ion_temperature * (c[3] - c[4] * ion_temperature)
-    ) / (1.0 + ion_temperature * (c[5] + c[6] * ion_temperature))
-
+    U = 1 - ion_temperature_kev * (
+        c[2] + ion_temperature_kev * (c[3] - c[4] * ion_temperature_kev)
+    ) / (1.0 + ion_temperature_kev * (c[5] + c[6] * ion_temperature_kev))
     val = (
         c[0]
-        * np.exp(-c[1] * (U / ion_temperature) ** (1 / 3))
-        / (U ** (5 / 6) * ion_temperature ** (2 / 3))
+        * np.exp(-c[1] * (U / ion_temperature_kev) ** (1 / 3))
+        / (U ** (5 / 6) * ion_temperature_kev ** (2 / 3))
     )
     return val