Skip to content

Commit

Permalink
format
Browse files Browse the repository at this point in the history
  • Loading branch information
shimwell committed Apr 2, 2024
2 parents 4c906fc + d1695ef commit 73eef16
Show file tree
Hide file tree
Showing 4 changed files with 41 additions and 34 deletions.
2 changes: 1 addition & 1 deletion examples/point_source_example.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@

plot = plot_source_energy(
this=settings,
n_samples=1000000, # increase this value for a smoother plot
n_samples=1000000, # increase this value for a smoother plot
energy_bins=np.linspace(0, 16e6, 1000),
yaxis_type="log",
)
Expand Down
5 changes: 1 addition & 4 deletions examples/tokamak_source_example.py
Original file line number Diff line number Diff line change
Expand Up @@ -52,9 +52,6 @@

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()
24 changes: 12 additions & 12 deletions src/openmc_plasma_source/fuel_types.py
Original file line number Diff line number Diff line change
Expand Up @@ -93,18 +93,18 @@ def neutron_energy_std_dev(ion_temperature: float, reaction: str) -> float:
return std_dev



def get_reactions_from_fuel(fuel):
if ["D", "T"] == sorted(set(fuel.keys())):
return ['DT', 'DD', 'TT']
return ["DT", "DD", "TT"]
elif ["D"] == sorted(set(fuel.keys())):
return ['DD']
return ["DD"]
elif ["T"] == sorted(set(fuel.keys())):
return ['TT']
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,
Expand Down Expand Up @@ -160,16 +160,16 @@ def get_neutron_energy_distribution(

reactions = get_reactions_from_fuel(fuel)

if reactions == ['TT']:
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': tt_source}, {'TT': 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': dd_source}, {'DD' :strength_DD}
return {"DD": dd_source}, {"DD": strength_DD}

# DT, DD and TT reaction
else:
Expand All @@ -186,9 +186,9 @@ def get_neutron_energy_distribution(
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')
dNdE_TT = np.trim_zeros(dNdE_TT, "b")
# making array lengths match
E_pspec = E_pspec[:len(dNdE_TT)]
E_pspec = E_pspec[: len(dNdE_TT)]

tt_source = openmc.stats.Tabular(E_pspec, dNdE_TT)

Expand All @@ -202,7 +202,7 @@ def get_neutron_energy_distribution(

# todo look into combining distributions openmc.data.combine_distributions()
return {
'TT': [tt_source, strength_TT / total_strength],
'DD': [dd_source, strength_DD / total_strength],
'DT': [dt_source, strength_DT / total_strength],
"TT": [tt_source, strength_TT / total_strength],
"DD": [dd_source, strength_DD / total_strength],
"DT": [dt_source, strength_DT / total_strength],
}
44 changes: 27 additions & 17 deletions src/openmc_plasma_source/tokamak_source.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
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,
minor_radius: float,
Expand Down Expand Up @@ -132,18 +133,18 @@ def tokamak_source(

fuel_densities = {}
for key, value in fuel.items():
fuel_densities[key] = densities*value
fuel_densities[key] = densities * value

# compute temperatures
temperatures = tokamak_ion_temperature(
r=a,
mode=mode,
pedestal_radius=pedestal_radius,
ion_temperature_pedestal=ion_temperature_pedestal/1e3,
ion_temperature_centre=ion_temperature_centre/1e3,
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/1e3,
ion_temperature_separatrix=ion_temperature_separatrix / 1e3,
major_radius=major_radius,
)

Expand All @@ -159,16 +160,20 @@ def tokamak_source(
)

reactions = get_reactions_from_fuel(fuel)
neutron_source_density={}

neutron_source_density = {}
total_source_density = 0
for reaction in reactions:
neutron_source_density[reaction] = tokamak_neutron_source_density(fuel_densities, temperatures, reaction)
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
neutron_source_density[reaction] = (
neutron_source_density[reaction] / total_source_density
)

sources = tokamak_make_openmc_sources(
strengths=neutron_source_density,
Expand Down Expand Up @@ -353,9 +358,12 @@ def tokamak_make_openmc_sources(
)

# now we have potentially 3 distributions (DT, DD, TT)
for reaction, (energy_distribution, dist_strength) in energy_distributions_and_dist_strengths.items():
for reaction, (
energy_distribution,
dist_strength,
) in energy_distributions_and_dist_strengths.items():

if dist_strength * strengths[reaction][i] > 0.:
if dist_strength * strengths[reaction][i] > 0.0:
my_source = openmc.IndependentSource()

# create a ring source
Expand Down Expand Up @@ -387,19 +395,21 @@ def tokamak_neutron_source_density(ion_density, ion_temperature, reaction):
float, ndarray: Neutron source density (neutron/s/m3)
"""

ion_density = np.asarray(ion_density[reaction[0]]) * np.asarray(ion_density[reaction[1]])
ion_density = np.asarray(ion_density[reaction[0]]) * np.asarray(
ion_density[reaction[1]]
)
ion_temperature = np.asarray(ion_temperature)

if reaction == ['DD']:
if reaction == ["DD"]:
return ion_density * reac_DD(ion_temperature)
elif reaction == ['TT']:
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
return ion_density * reac_DT(ion_temperature) # could use _DT_xs instead


#TODO consider replace with NeSST or getting DD version as well
# 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
Expand All @@ -408,7 +418,7 @@ def _DT_xs(ion_temperature):
Returns:
float, ndarray: the DT cross section at the given temperature
"""
ion_temperature_kev = np.asarray(ion_temperature/1e3)
ion_temperature_kev = np.asarray(ion_temperature / 1e3)
c = [
2.5663271e-18,
19.983026,
Expand Down

0 comments on commit 73eef16

Please sign in to comment.