improved neutron energy

[shimwell]

instead of a 14MeV or a 2.5MeV Muir distribution we could have a distribution that accounts for DD, TT and DT reactions

import NeSST as nst

def create_neutron_source_term(
        ion_temperature: float,
        fuel: dict = {'D':0.5, 'T':0.5},
) -> openmc.stats.Discrete:
    """Finds the energy distribution
"""

    ion_temperature = ion_temperature / 1e3  # convert eV to keV

    sum_fuel_isotopes = sum(fuel.values())
    if sum_fuel_isotopes > 1.:
        raise ValueError(f'isotope fractions within the fuel must sum to be below 1. Not {sum_fuel_isotopes}')

    if sum_fuel_isotopes < 0.:
        raise ValueError(f'isotope must sum to be above 0. Not {sum_fuel_isotopes}')

    for k, v in fuel.dict:
        if k not in ['D', 'T']:
            raise ValueError('Fuel dictionary keys must be either "D" or "T" not "{k}".)
        if v < 0
            raise ValueError('Fuel dictionary values must be above 0 not "{k}".)
        if v > 1
            raise ValueError('Fuel dictionary values must be below 1 not "{k}".)

    #Set atomic fraction of D and T in scattering medium and source
    if 'D' in fuel.keys():
        nst.frac_D_default = fuel['D']
    else:
        nst.frac_D_default = 0

    if 'T' in fuel.keys():
        nst.frac_T_default = fuel['T']
    else:
        nst.frac_T_default = 0

    # 1.0 neutron yield, all reactions scaled by this value
    num_of_vals = 500
    # single grid for DT, DD and TT grid
    E_pspec = np.linspace(0, 20, num_of_vals)  # accepts MeV units

    dNdE_DT_DD_TT = np.zeros(num_of_vals)
    if 'D' in fuel.keys() and 'T' in fuel.keys():
        DTmean, DTvar = nst.DTprimspecmoments(ion_temperature)
        Y_DT = 1.0
        dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar)  # Brysk shape i.e. Gaussian
        dNdE_DT_DD_TT= dNdE_DT_DD_TT + dNdE_DT
    if 'D' in fuel.keys()
        DDmean, DDvar = nst.DDprimspecmoments(ion_temperature)
        Y_DD = nst.yield_from_dt_yield_ratio("dd", 1.0, ion_temperature)
        dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar)  # Brysk shape i.e. Gaussian
        dNdE_DT_DD_TT= dNdE_DT_DD_TT + dNdE_DD
    if 'T' in fuel.keys()
        Y_TT = nst.yield_from_dt_yield_ratio("tt", 1.0, ion_temperature)
        dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature)
        dNdE_DT_DD_TT= dNdE_DT_DD_TT + dNdE_TT

    return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT)

[shimwell]

mixture can be used to combine continuous distributions with probabilities

import openmc
import openmc_source_plotter
import numpy as np

dt_source = openmc.IndependentSource()
dt_source.space = openmc.stats.Point((0, 0, 0))
dt_source.angle = openmc.stats.Isotropic()

# dd_reaction_energy = openmc.stats.muir(e0=2500000.0, m_rat=4.0, kt=20000.0)
# dt_reaction_energy = openmc.stats.muir(e0=14080000.0, m_rat=5.0, kt=20000.0)

dd_reaction_energy = openmc.stats.Normal(mean_value=2500000.0, std_dev=1e6)
dt_reaction_energy = openmc.stats.Normal(mean_value=14080000.0, std_dev=1e6)

both = openmc.stats.Mixture([0.2,0.8], [dd_reaction_energy,dt_reaction_energy])

dt_source.energy = both

settings = openmc.Settings()
settings.particles = 1
settings.batches = 1
settings.source = dt_source

settings.export_to_xml()

materials = openmc.Materials()
sph = openmc.Sphere(r=1000000, boundary_type="vacuum")
cell = openmc.Cell(region=-sph)
geometry = openmc.Geometry([cell])
model = openmc.Model(geometry, materials, settings)

plot = model.plot_source_energy(n_samples=2000000, energy_bins=np.linspace(0,20e6, 1000))
plot.show()

[shimwell]

combine distributions can be used to combine discrete and tabular

import openmc
import openmc_source_plotter
import numpy as np

dt_source = openmc.IndependentSource()
dt_source.space = openmc.stats.Point((0, 0, 0))
dt_source.angle = openmc.stats.Isotropic()

dd_reaction_energy = openmc.stats.Discrete([2e6,3e6,4e6], [1,0.7,0.5])
dt_reaction_energy = openmc.stats.Discrete([12e6,13e6,14e6], [2,1.4,0.5])
both = openmc.data.combine_distributions([dd_reaction_energy,dt_reaction_energy], [0.5,0.5])

dt_source.energy = both

settings = openmc.Settings()
settings.particles = 1
settings.batches = 1
settings.source = dt_source

settings.export_to_xml()

materials = openmc.Materials()
sph = openmc.Sphere(r=1000000, boundary_type="vacuum")
cell = openmc.Cell(region=-sph)
geometry = openmc.Geometry([cell])
model = openmc.Model(geometry, materials, settings)

plot = model.plot_source_energy(n_samples=2000000, energy_bins=np.linspace(0,20e6, 1000))
plot.show()

[shimwell]

Ballabio method

a_DD_s= [ 4.69515, -0.040729, 0.47, 0.81844, 18.225, 2.1525 ]
a_DD_w = [1.7013e-3, 0.16888,  0.49,  7.9460e-4, 8.4619e-3, 8.3241e-4]

a_DT_s = [5.30509, 2.4736e-3, 1.84, 1.3818, 37.771, 0.92181]
a_DT_w = [5.1068e-4, 7.6223e-3, 1.78,  8.7691e-5, 2.0199e-3, 5.9501e-5 ]

Ti = np.linspace(0.5, 100, 99)

a= a_DD_s

E_shift_DD = []
E_shift_DT = []
for t in Ti:
    a = a_DD_s
    if t<40:
        e = (a[0]/ (1+a[1]*(t**a[2])))*t**(2/3) + a[3]*t
    else:
        e = a[4] + a[5]*t
    E_shift_DD.append(e)

    a = a_DT_s
    if t < 40:
        e = (a[0] / (1 + a[1] * (t ** a[2]))) * t ** (2 / 3) + a[3] * t
    else:
        e = a[4] + a[5] * t
    E_shift_DT.append(e)

E_w_DD = []
E_w_DT = []
for t in Ti:
    a = a_DD_w
    if t<40:
        s = (a[0]/ (1+a[1]*(t**a[2])))*t**(2/3) + a[3]*t
    else:
        s = a[4] + a[5]*t
    E_w_DD.append(82.542*(1+s)*math.sqrt(t))

    a = a_DT_w
    if t < 40:
        s = (a[0] / (1 + a[1] * (t ** a[2]))) * t ** (2 / 3) + a[3] * t
    else:
        s = a[4] + a[5] * t
    E_w_DT.append(177.259*(1+s)*math.sqrt(t))

fig, ax = plt.subplots(2)
ax[0].plot(Ti, E_shift_DD, label = 'DD')
ax[0].plot(Ti, E_shift_DT, label = 'DT')

[RemDelaporteMathurin]

Ballabio method

This won't scale well as the length of Ti increases. I would avoid iterating through numpy arrays like this.
Consider for instance:

Ti = np.linspace(0.5, 100, 99)

a = a_DD_s

E_shift_DT = np.zeros_like(Ti)

idx = np.where(Ti < 40)
E_shift_DT[idx] = (a[0] / (1 + a[1] * (Ti[idx] ** a[2]))) * Ti[idx] ** (2 / 3) + a[3] * Ti[idx]

idx = np.where(Ti >= 40)
E_shift_DT[idx] = a[4] + a[5] * Ti[idx]

[shimwell]

Thanks Remi