diff --git a/examples/point_source_example.py b/examples/point_source_example.py index 46d15c2..740406b 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -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", ) diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index 4711476..ebc286c 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -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() diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index ef3ed03..e0c5ce8 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -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, @@ -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: @@ -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) @@ -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], } diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 3b30b94..b1983e1 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -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, @@ -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, ) @@ -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, @@ -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 @@ -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 @@ -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,