Open shimwell opened 8 months ago
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()
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()
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')
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]
Thanks Remi
instead of a 14MeV or a 2.5MeV Muir distribution we could have a distribution that accounts for DD, TT and DT reactions