diff --git a/pyproject.toml b/pyproject.toml index 9c318f1..2f26878 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,7 +23,7 @@ keywords = ["python", "neutron", "fusion", "source", "openmc", "energy", "tokama dependencies = [ "numpy>=1.9", "matplotlib>=3.2.2", - "NeSST" + "NeSST@git+https://github.com/aidancrilly/NeSST", ] [project.urls] diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 3612e0c..c48e83a 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -112,8 +112,6 @@ def get_neutron_energy_distribution( energy distribution """ - ion_temperature_kev = ion_temperature / 1e3 # convert eV to keV - sum_fuel_isotopes = sum(fuel.values()) if sum_fuel_isotopes > 1.0: raise ValueError( @@ -151,15 +149,15 @@ def get_neutron_energy_distribution( if ["D", "T"] == sorted(set(fuel.keys())): strength_DT = 1.0 strength_DD = nst.yield_from_dt_yield_ratio( - "dd", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"] + "dd", strength_DT, ion_temperature, fuel["D"], fuel["T"] ) strength_TT = nst.yield_from_dt_yield_ratio( - "tt", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"] + "tt", strength_DT, ion_temperature, fuel["D"], fuel["T"] ) total_strength = sum([strength_TT, strength_DD, strength_DT]) - dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) + dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) dd_source = openmc.stats.Normal(mean_value=DDmean, std_dev=DD_std_dev) @@ -184,6 +182,6 @@ def get_neutron_energy_distribution( elif ["T"] == sorted(set(fuel.keys())): strength_TT = 1.0 - dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) + 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] diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 3716f62..cd572c4 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -9,7 +9,7 @@ def fusion_point_source( coordinate: Tuple[float, float, float] = (0.0, 0.0, 0.0), temperature: float = 20000.0, fuel: dict = {"D": 0.5, "T": 0.5}, -): +) -> list[openmc.IndependentSource]: """Creates a list of openmc.IndependentSource objects representing an ICF source. Resulting ICF (Inertial Confinement Fusion) source will have an energy diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index cbf3e9d..097494b 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -118,8 +118,7 @@ def tokamak_source( a = np.random.random(sample_size) * minor_radius alpha = np.random.random(sample_size) * 2 * np.pi - # compute densities, temperatures, neutron source densities and - # convert coordinates + # compute densities, temperatures densities = tokamak_ion_density( mode=mode, ion_density_centre=ion_density_centre, @@ -130,6 +129,8 @@ def tokamak_source( ion_density_separatrix=ion_density_separatrix, r=a, ) + + # compute temperatures temperatures = tokamak_ion_temperature( r=a, mode=mode, @@ -142,10 +143,7 @@ def tokamak_source( major_radius=major_radius, ) - neutron_source_density = tokamak_neutron_source_density(densities, temperatures) - - strengths = neutron_source_density / sum(neutron_source_density) - + # convert coordinates RZ = tokamak_convert_a_alpha_to_R_Z( a=a, alpha=alpha, @@ -156,6 +154,14 @@ def tokamak_source( elongation=elongation, ) + #TODO loop through the fuel reactions + + # compute neutron source densities + neutron_source_density = tokamak_neutron_source_density(densities, temperatures) + + strengths = neutron_source_density / sum(neutron_source_density) + + sources = tokamak_make_openmc_sources( strengths=strengths, angles=angles, @@ -313,8 +319,12 @@ def tokamak_make_openmc_sources( is based on .strengths. Args: + strengths angles ((float, float), optional): rotation of the ring source. - Defaults to (0, 2*np.pi). + Defaults to (0, 2*np.pi). + temperatures + fuel + RZ Returns: list: list of openmc.IndependentSource() @@ -357,7 +367,7 @@ def tokamak_make_openmc_sources( return sources -def tokamak_neutron_source_density(ion_density, ion_temperature): +def tokamak_neutron_source_density(ion_density, ion_temperature, reaction='DT'): """Computes the neutron source density given ion density and ion temperature. @@ -371,23 +381,26 @@ def tokamak_neutron_source_density(ion_density, ion_temperature): ion_density = np.asarray(ion_density) ion_temperature = np.asarray(ion_temperature) - return ion_density**2 * _DT_xs(ion_temperature) + #TODO account for reaction + if reaction == 'DT': + nst.reac_DT(ion_temperature) + if reaction == 'DD' + nst.reac_DD(ion_temperature) + if reaction == 'TT' + nst.reac_TT(ion_temperature) +#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 - Args: - ion_temperature (float, ndarray): ion temperature in keV - + ion_temperature (float, ndarray): ion temperature in eV Returns: float, ndarray: the DT cross section at the given temperature """ - - ion_temperature = np.asarray(ion_temperature) - + ion_temperature_kev = np.asarray(ion_temperature/1e3) c = [ 2.5663271e-18, 19.983026, @@ -397,14 +410,12 @@ def _DT_xs(ion_temperature): 6.6024089e-2, 8.1215505e-3, ] - - U = 1 - ion_temperature * ( - c[2] + ion_temperature * (c[3] - c[4] * ion_temperature) - ) / (1.0 + ion_temperature * (c[5] + c[6] * ion_temperature)) - + U = 1 - ion_temperature_kev * ( + c[2] + ion_temperature_kev * (c[3] - c[4] * ion_temperature_kev) + ) / (1.0 + ion_temperature_kev * (c[5] + c[6] * ion_temperature_kev)) val = ( c[0] - * np.exp(-c[1] * (U / ion_temperature) ** (1 / 3)) - / (U ** (5 / 6) * ion_temperature ** (2 / 3)) + * np.exp(-c[1] * (U / ion_temperature_kev) ** (1 / 3)) + / (U ** (5 / 6) * ion_temperature_kev ** (2 / 3)) ) return val