RZ values are nested?

examples/point_source_example.py

@@ -34,12 +34,12 @@
 # optionally if you would like to plot the energy of particles then another package can be used
 # https://github.com/fusion-energy/openmc_source_plotter
 
-# from openmc_source_plotter import plot_source_energy
-
-# plot = plot_source_energy(
-#     this=settings,
-#     n_samples=1000000, # increase this value for a smoother plot
-#     energy_bins=np.linspace(0, 16e6, 1000),
-#     yaxis_type="log",
-# )
-# plot.show()
+from openmc_source_plotter import plot_source_energy
+
+plot = plot_source_energy(
+    this=settings,
+    n_samples=1000000, # increase this value for a smoother plot
+    energy_bins=np.linspace(0, 16e6, 1000),
+    yaxis_type="log",
+)
+plot.show()

examples/tokamak_source_example.py

@@ -17,20 +17,20 @@
 my_sources = tokamak_source(
     elongation=1.557,
     ion_density_centre=1.09e20,
-    ion_density_peaking_factor=1,
     ion_density_pedestal=1.09e20,
+    ion_density_peaking_factor=1,
     ion_density_separatrix=3e19,
-    ion_temperature_centre=45.9,
+    ion_temperature_centre=45.9e3,
+    ion_temperature_pedestal=6.09e3,
+    ion_temperature_separatrix=0.1e3,
     ion_temperature_peaking_factor=8.06,
-    ion_temperature_pedestal=6.09,
-    ion_temperature_separatrix=0.1,
+    ion_temperature_beta=6,
     major_radius=906,
     minor_radius=292.258,
     pedestal_radius=0.8 * 292.258,
     mode="H",
     shafranov_factor=0.44789,
     triangularity=0.270,
-    ion_temperature_beta=6,
     fuel={"D": 0.5, "T": 0.5},
 )
 
@@ -50,11 +50,11 @@
 # optionally if you would like to plot the direction of particles then another package can be used
 # https://github.com/fusion-energy/openmc_source_plotter
 
-# from openmc_source_plotter import plot_source_direction
+from openmc_source_plotter import plot_source_direction
 
-# plot = plot_source_direction(
-#     this=settings,
-#     n_samples=200
-# )
+plot = plot_source_direction(
+    this=settings,
+    n_samples=200
+)
 
-# plot.show()
+plot.show()

src/openmc_plasma_source/fuel_types.py

@@ -164,12 +164,12 @@ def get_neutron_energy_distribution(
         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]
+        return {'TT': tt_source}, {'TT': 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]
+        return {'DD': dd_source}, {'DD' :strength_DD}
 
     # DT, DD and TT reaction
     else:
@@ -201,8 +201,8 @@ def get_neutron_energy_distribution(
         # dt_source = openmc.stats.muir(e0=DTmean * 1e6, m_rat=5, kt=ion_temperature)
 
         # todo look into combining distributions openmc.data.combine_distributions()
-        return [tt_source, dd_source, dt_source], [
-            strength_TT / total_strength,
-            strength_DD / total_strength,
-            strength_DT / total_strength,
-        ]
+        return {
+            'TT': [tt_source, strength_TT / total_strength],
+            'DD': [dd_source, strength_DD / total_strength],
+            'DT': [dt_source, strength_DT / total_strength],
+        }

src/openmc_plasma_source/tokamak_source.py

@@ -4,8 +4,8 @@
 import openmc
 import openmc.checkvalue as cv
 import NeSST as nst
-from .fuel_types import get_neutron_energy_distribution
-
+from .fuel_types import get_neutron_energy_distribution, get_reactions_from_fuel
+from NeSST.spectral_model import reac_DD, reac_TT, reac_DT
 
 def tokamak_source(
     major_radius: float,
@@ -54,14 +54,14 @@ def tokamak_source(
         ion_density_pedestal (float): Ion density at pedestal (m-3)
         ion_density_separatrix (float): Ion density at separatrix (m-3)
         ion_temperature_centre (float): Ion temperature at the plasma
-            centre (keV)
+            centre (eV)
         ion_temperature_peaking_factor (float): Ion temperature peaking
             factor (referred in [1] as ion temperature exponent alpha_T)
         ion_temperature_beta (float): Ion temperature beta exponent
             (referred in [1] as ion temperature exponent beta_T)
-        ion_temperature_pedestal (float): Ion temperature at pedestal (keV)
+        ion_temperature_pedestal (float): Ion temperature at pedestal (eV)
         ion_temperature_separatrix (float): Ion temperature at separatrix
-            (keV)
+            (eV)
         pedestal_radius (float): Minor radius at pedestal (cm)
         shafranov_factor (float): Shafranov factor (referred in [1] as esh)
             also known as outward radial displacement of magnetic surfaces
@@ -139,11 +139,11 @@ def tokamak_source(
         r=a,
         mode=mode,
         pedestal_radius=pedestal_radius,
-        ion_temperature_pedestal=ion_temperature_pedestal,
-        ion_temperature_centre=ion_temperature_centre,
+        ion_temperature_pedestal=ion_temperature_pedestal/1e3,
+        ion_temperature_centre=ion_temperature_centre/1e3,
         ion_temperature_beta=ion_temperature_beta,
         ion_temperature_peaking_factor=ion_temperature_peaking_factor,
-        ion_temperature_separatrix=ion_temperature_separatrix,
+        ion_temperature_separatrix=ion_temperature_separatrix/1e3,
         major_radius=major_radius,
     )
 
@@ -158,19 +158,27 @@ 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,
-        temperatures=temperatures,
-        fuel=fuel,
-        RZ=RZ,
-    )
-    return sources
+    reactions = get_reactions_from_fuel(fuel)
+
+    neutron_source_density={}
+    total_source_density = 0
+    for reaction in reactions:
+        neutron_source_density[reaction] = tokamak_neutron_source_density(fuel_densities, temperatures, reaction)
+        total_source_density += sum(neutron_source_density[reaction])
+
+    all_sources = []
+    for reaction in reactions:
+        neutron_source_density[reaction] = neutron_source_density[reaction]/total_source_density
+
+        sources = tokamak_make_openmc_sources(
+            strengths=neutron_source_density,
+            angles=angles,
+            temperatures=temperatures,
+            fuel=fuel,
+            RZ=RZ,
+        )
+        all_sources = all_sources + sources
+    return all_sources
 
 
 def tokamak_ion_density(
@@ -333,38 +341,37 @@ def tokamak_make_openmc_sources(
 
     sources = []
     # create a ring source for each sample in the plasma source
-    for i in range(len(strengths)):
+    for i, (RZ_val, temperature) in enumerate(zip(RZ, temperatures)):
         # extract the RZ values accordingly
-        radius = openmc.stats.Discrete([RZ[0][i]], [1])
-        z_values = openmc.stats.Discrete([RZ[1][i]], [1])
+        radius = openmc.stats.Discrete([RZ_val[0]], [1])
+        z_values = openmc.stats.Discrete([RZ_val[1]], [1])
         angle = openmc.stats.Uniform(a=angles[0], b=angles[1])
-        strength = strengths[i]
 
-        energy_distributions, dist_strengths = get_neutron_energy_distribution(
-            ion_temperature=temperatures[i],
+        energy_distributions_and_dist_strengths = get_neutron_energy_distribution(
+            ion_temperature=temperature,
             fuel=fuel,
         )
 
         # now we have potentially 3 distributions (DT, DD, TT)
-        for energy_distribution, dist_strength in zip(
-            energy_distributions, dist_strengths
-        ):
-            my_source = openmc.IndependentSource()
+        for reaction, (energy_distribution, dist_strength) in energy_distributions_and_dist_strengths.items():
 
-            # create a ring source
-            my_source.space = openmc.stats.CylindricalIndependent(
-                r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0)
-            )
-            my_source.angle = openmc.stats.Isotropic()
+            if dist_strength * strengths[reaction][i] > 0.:
+                my_source = openmc.IndependentSource()
 
-            my_source.energy = energy_distribution
+                # create a ring source
+                my_source.space = openmc.stats.CylindricalIndependent(
+                    r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0)
+                )
+                my_source.angle = openmc.stats.Isotropic()
 
-            # the strength of the source (its probability) is given by the
-            # strength of the energy distribution and the location distribution
-            my_source.strength = dist_strength * strength
+                my_source.energy = energy_distribution
 
-            # append to the list of sources
-            sources.append(my_source)
+                # the strength of the source (its probability) is given by the
+                # strength of the energy distribution and the location distribution
+                my_source.strength = dist_strength * strengths[reaction][i]
+
+                # append to the list of sources
+                sources.append(my_source)
     return sources
 
 
@@ -383,12 +390,13 @@ def tokamak_neutron_source_density(ion_density, ion_temperature, reaction):
     ion_density = np.asarray(ion_density[reaction[0]]) * np.asarray(ion_density[reaction[1]])
     ion_temperature = np.asarray(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)
+    if reaction == ['DD']:
+        return ion_density * reac_DD(ion_temperature)
+    elif reaction == ['TT']:
+        return ion_density * reac_TT(ion_temperature)
+    # ['DT', 'DD', 'TT']
+    else:
+        return ion_density * reac_DT(ion_temperature) # could use _DT_xs instead
 
 
 #TODO consider replace with NeSST or getting DD version as well