From 77af3f9782f30b38a9d1911d84e40070e61b092e Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 20 Jan 2024 22:35:32 +0000 Subject: [PATCH 01/44] improved energy distribution --- src/openmc_plasma_source/fuel_types.py | 111 +++++++++++++++------ src/openmc_plasma_source/point_source.py | 10 +- src/openmc_plasma_source/ring_source.py | 13 +-- src/openmc_plasma_source/tokamak_source.py | 8 +- 4 files changed, 95 insertions(+), 47 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 27299b9..509d2bd 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -1,37 +1,92 @@ -"""fuel_types.py +import NeSST as nst +import numpy as np +import openmc -Defines dictionary for determining mean energy and mass of reactants -for a given fusion fuel type. -""" +def get_neutron_energy_distribution( + ion_temperature: float, + fuel: dict = {'D':0.5, 'T':0.5}, +) -> openmc.stats.Discrete: + """Finds the energy distribution + Parameters + ---------- + ion_temperature : float + temperature of plasma ions in eV + fuel : dict + isotopes as keys and atom fractions as values -class Fuel: - def __init__(self, mean_energy, mass_of_reactants): - self.mean_energy = mean_energy - self.mass_of_reactants = mass_of_reactants + Returns + ------- + openmc.stats.Discrete + energy distribution + """ - @property - def mean_energy(self): - return self._mean_energy + ion_temperature = ion_temperature / 1e3 # convert eV to keV - @mean_energy.setter - def mean_energy(self, value): - if value <= 0: - raise (ValueError("mean_energy needs to be strictly positive")) - self._mean_energy = value + 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}') - @property - def mass_of_reactants(self): - return self._mass_of_reactants + if sum_fuel_isotopes < 0.: + raise ValueError(f'isotope must sum to be above 0. Not {sum_fuel_isotopes}') - @mass_of_reactants.setter - def mass_of_reactants(self, value): - if value <= 0: - raise (ValueError("mass_of_reactants needs to be strictly positive")) - self._mass_of_reactants = value + for k, v in fuel.dict: + if k not in ['D', 'T']: + raise ValueError(f'Fuel dictionary keys must be either "D" or "T" not "{k}".') + if v < 0: + raise ValueError(f'Fuel dictionary values must be above 0 not "{k}".') + if v > 1: + raise ValueError(f'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'] + max_energy_mev=5 + else: + nst.frac_D_default = 0 -fuel_types = { - "DD": Fuel(mean_energy=2450000.0, mass_of_reactants=4), - "DT": Fuel(mean_energy=14080000.0, mass_of_reactants=5), -} + if 'T' in fuel.keys(): + nst.frac_T_default = fuel['T'] + max_energy_mev=12 + else: + nst.frac_T_default = 0 + + if 'T' in fuel.keys() and 'D' in fuel.keys(): + max_energy_mev=20 + + # 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, max_energy_mev, 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) + DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) + + Y_DT = 1.0 + Y_DD = nst.yield_from_dt_yield_ratio("dd", Y_DT, ion_temperature) + Y_TT = nst.yield_from_dt_yield_ratio("tt", Y_DT, ion_temperature) + + dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar) # Brysk shape i.e. Gaussian + dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian + dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) + dNdE_DT_DD_TT= dNdE_DT + dNdE_DD + dNdE_TT + + if 'D' in fuel.keys() and 'T' not in fuel.keys(): + DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) + DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) + + Y_DD = 1.0 + + dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian + dNdE_DT_DD_TT= dNdE_DD + + if 'D' not in fuel.keys() and 'T' in fuel.keys(): + + Y_TT = 1.0 + + dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) + dNdE_DT_DD_TT= dNdE_TT + + return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index e123d1c..664ccfd 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -1,7 +1,7 @@ import openmc from typing import Tuple -from .fuel_types import fuel_types +from .fuel_types import get_neutron_energy_distribution class FusionPointSource(openmc.IndependentSource): @@ -21,7 +21,7 @@ def __init__( self, coordinate: Tuple[float, float, float] = (0.0, 0.0, 0.0), temperature: float = 20000.0, - fuel: str = "DT", + fuel: dict = {"D": 0.5, "T": 0.5}, ): # Set local attributes self.coordinate = coordinate @@ -35,11 +35,7 @@ def __init__( # performed after the super init as these are Source attributes self.space = openmc.stats.Point(self.coordinate) self.angle = openmc.stats.Isotropic() - self.energy = openmc.stats.muir( - e0=self.fuel.mean_energy, - m_rat=self.fuel.mass_of_reactants, - kt=self.temperature, - ) + self.energy = get_neutron_energy_distribution(ion_temperature=temperature, fuel=fuel) @property def coordinate(self): diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index 9a4c501..bc09169 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -2,7 +2,7 @@ import numpy as np from typing import Tuple -from .fuel_types import fuel_types +from .fuel_types import get_neutron_energy_distribution class FusionRingSource(openmc.IndependentSource): @@ -25,15 +25,14 @@ def __init__( angles: Tuple[float, float] = (0, 2 * np.pi), z_placement: float = 0, temperature: float = 20000.0, - fuel: str = "DT", + fuel: dict = {"D": 0.5, "T": 0.5}, ): # Set local attributes self.radius = radius self.angles = angles self.z_placement = z_placement self.temperature = temperature - self.fuel_type = fuel - self.fuel = fuel_types[self.fuel_type] + self.fuel = fuel # Call init for openmc.Source super().__init__() @@ -46,11 +45,7 @@ def __init__( origin=(0.0, 0.0, 0.0), ) self.angle = openmc.stats.Isotropic() - self.energy = openmc.stats.muir( - e0=self.fuel.mean_energy, - m_rat=self.fuel.mass_of_reactants, - kt=self.temperature, - ) + self.energy = get_neutron_energy_distribution(ion_temperature=temperature, fuel=fuel) @property def radius(self): diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index a594e95..6f9a0fe 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -2,6 +2,8 @@ import numpy as np from typing import Tuple +from .fuel_types import get_neutron_energy_distribution + class TokamakSource: """Plasma neutron source sampling. @@ -68,6 +70,7 @@ def __init__( shafranov_factor: float, angles: Tuple[float, float] = (0, 2 * np.pi), sample_size: int = 1000, + fuel: dict = {"D": 0.5, "T": 0.5}, ) -> None: # Assign attributes self.major_radius = major_radius @@ -88,6 +91,7 @@ def __init__( self.shafranov_factor = shafranov_factor self.angles = angles self.sample_size = sample_size + self.fuel = fuel # Perform sanity checks for inputs not caught by properties if self.minor_radius >= self.major_radius: @@ -385,9 +389,7 @@ def make_openmc_sources(self): ) my_source.angle = openmc.stats.Isotropic() - my_source.energy = openmc.stats.muir( - e0=14080000.0, m_rat=5.0, kt=self.temperatures[i] - ) + my_source.energy = get_neutron_energy_distribution(ion_temperature=self.temperature[i], fuel=self.fuel) # the strength of the source (its probability) is given by # self.strengths From 547154dfc1a7af36b556afa2fb3ec7e8ce1b2cf5 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 20 Jan 2024 22:35:55 +0000 Subject: [PATCH 02/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/fuel_types.py | 56 ++++++++++++---------- src/openmc_plasma_source/point_source.py | 4 +- src/openmc_plasma_source/ring_source.py | 4 +- src/openmc_plasma_source/tokamak_source.py | 4 +- 4 files changed, 39 insertions(+), 29 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 509d2bd..5e8856e 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -2,9 +2,10 @@ import numpy as np import openmc + def get_neutron_energy_distribution( - ion_temperature: float, - fuel: dict = {'D':0.5, 'T':0.5}, + ion_temperature: float, + fuel: dict = {"D": 0.5, "T": 0.5}, ) -> openmc.stats.Discrete: """Finds the energy distribution @@ -24,35 +25,39 @@ def get_neutron_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 > 1.0: + 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}') + if sum_fuel_isotopes < 0.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(f'Fuel dictionary keys must be either "D" or "T" not "{k}".') + if k not in ["D", "T"]: + raise ValueError( + f'Fuel dictionary keys must be either "D" or "T" not "{k}".' + ) if v < 0: raise ValueError(f'Fuel dictionary values must be above 0 not "{k}".') if v > 1: raise ValueError(f'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'] - max_energy_mev=5 + # Set atomic fraction of D and T in scattering medium and source + if "D" in fuel.keys(): + nst.frac_D_default = fuel["D"] + max_energy_mev = 5 else: nst.frac_D_default = 0 - if 'T' in fuel.keys(): - nst.frac_T_default = fuel['T'] - max_energy_mev=12 + if "T" in fuel.keys(): + nst.frac_T_default = fuel["T"] + max_energy_mev = 12 else: nst.frac_T_default = 0 - - if 'T' in fuel.keys() and 'D' in fuel.keys(): - max_energy_mev=20 + + if "T" in fuel.keys() and "D" in fuel.keys(): + max_energy_mev = 20 # 1.0 neutron yield, all reactions scaled by this value num_of_vals = 500 @@ -60,7 +65,7 @@ def get_neutron_energy_distribution( E_pspec = np.linspace(0, max_energy_mev, 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(): + if "D" in fuel.keys() and "T" in fuel.keys(): DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) @@ -71,22 +76,21 @@ def get_neutron_energy_distribution( dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar) # Brysk shape i.e. Gaussian dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) - dNdE_DT_DD_TT= dNdE_DT + dNdE_DD + dNdE_TT + dNdE_DT_DD_TT = dNdE_DT + dNdE_DD + dNdE_TT - if 'D' in fuel.keys() and 'T' not in fuel.keys(): + if "D" in fuel.keys() and "T" not in fuel.keys(): DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) Y_DD = 1.0 dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian - dNdE_DT_DD_TT= dNdE_DD - - if 'D' not in fuel.keys() and 'T' in fuel.keys(): + dNdE_DT_DD_TT = dNdE_DD + if "D" not in fuel.keys() and "T" in fuel.keys(): Y_TT = 1.0 - + dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) - dNdE_DT_DD_TT= dNdE_TT + dNdE_DT_DD_TT = dNdE_TT return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 664ccfd..e0d4ff5 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -35,7 +35,9 @@ def __init__( # performed after the super init as these are Source attributes self.space = openmc.stats.Point(self.coordinate) self.angle = openmc.stats.Isotropic() - self.energy = get_neutron_energy_distribution(ion_temperature=temperature, fuel=fuel) + self.energy = get_neutron_energy_distribution( + ion_temperature=temperature, fuel=fuel + ) @property def coordinate(self): diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index bc09169..31ffb12 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -45,7 +45,9 @@ def __init__( origin=(0.0, 0.0, 0.0), ) self.angle = openmc.stats.Isotropic() - self.energy = get_neutron_energy_distribution(ion_temperature=temperature, fuel=fuel) + self.energy = get_neutron_energy_distribution( + ion_temperature=temperature, fuel=fuel + ) @property def radius(self): diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 6f9a0fe..ace9811 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -389,7 +389,9 @@ def make_openmc_sources(self): ) my_source.angle = openmc.stats.Isotropic() - my_source.energy = get_neutron_energy_distribution(ion_temperature=self.temperature[i], fuel=self.fuel) + my_source.energy = get_neutron_energy_distribution( + ion_temperature=self.temperature[i], fuel=self.fuel + ) # the strength of the source (its probability) is given by # self.strengths From efa65dc214290291048b4450141160a151d5c165 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 20 Jan 2024 22:39:49 +0000 Subject: [PATCH 03/44] doc strings --- src/openmc_plasma_source/point_source.py | 2 +- src/openmc_plasma_source/ring_source.py | 6 +++--- src/openmc_plasma_source/tokamak_source.py | 1 + 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 664ccfd..fd326b8 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -14,7 +14,7 @@ class FusionPointSource(openmc.IndependentSource): coordinate (tuple[float,float,float]): Location of the point source. Each component is measured in metres. temperature (float): Temperature of the source (eV). - fuel_type (str): The fusion fuel mix. Either 'DT' or 'DD'. + fuel (dict): Isotopes as keys and atom fractions as values """ def __init__( diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index bc09169..8229b96 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -8,15 +8,15 @@ class FusionRingSource(openmc.IndependentSource): """An openmc.Source object with some presets to make it more convenient for fusion simulations using a ring source. All attributes can be changed - after initialization if required. Default isotropic ring source with a Muir - energy distribution. + after initialization if required. Default isotropic ring source with a + realistic energy distribution. Args: radius (float): the inner radius of the ring source, in metres angles (iterable of floats): the start and stop angles of the ring in radians z_placement (float): Location of the ring source (m). Defaults to 0. temperature (float): the temperature to use in the Muir distribution in eV, - fuel_type (str): The fusion fuel mix. Either 'DT' or 'DD'. + fuel (dict): Isotopes as keys and atom fractions as values """ def __init__( diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 6f9a0fe..ff57bf0 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -48,6 +48,7 @@ class TokamakSource: radians sample_size (int, optional): number of neutron sources. Defaults to 1000. + fuel (dict): Isotopes as keys and atom fractions as values """ def __init__( From 1c33d0103bd5d8732ff73419fda3de0716781176 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sat, 20 Jan 2024 22:56:41 +0000 Subject: [PATCH 04/44] added nesst --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index c421db0..7fc3d6f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ keywords = ["python", "neutron", "fusion", "source", "openmc", "energy", "tokama dependencies = [ "numpy>=1.9", "matplotlib>=3.2.2", + "NeSST" ] [project.urls] From ff728e64f0641e76042dc32ac898e3a9c3cd3481 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 30 Jan 2024 15:39:46 +0000 Subject: [PATCH 05/44] simpler switch statements --- examples/point_source_example.py | 2 +- examples/ring_source_example.py | 1 + examples/tokamak_source_example.py | 1 + src/openmc_plasma_source/__init__.py | 1 + src/openmc_plasma_source/fuel_types.py | 43 +++++++++++------------- src/openmc_plasma_source/point_source.py | 14 +------- 6 files changed, 24 insertions(+), 38 deletions(-) diff --git a/examples/point_source_example.py b/examples/point_source_example.py index fb040a1..e8f8f14 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -11,7 +11,7 @@ geometry = openmc.Geometry([cell]) # define the source -my_source = FusionPointSource(coordinate=(0, 0, 0), temperature=20000.0, fuel="DT") +my_source = FusionPointSource(coordinate=(0, 0, 0), temperature=20000.0, fuel={"D":0.5, "T":0.5}) # Tell OpenMC we're going to use our custom source settings = openmc.Settings() diff --git a/examples/ring_source_example.py b/examples/ring_source_example.py index f09d570..7c2b4ab 100644 --- a/examples/ring_source_example.py +++ b/examples/ring_source_example.py @@ -16,6 +16,7 @@ radius=700, angles=(0.0, 2 * math.pi), # 360deg source temperature=20000.0, + fuel={"D":0.5,"T":0.5} ) settings = openmc.Settings() diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index 9c98390..d8eb53b 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -29,6 +29,7 @@ shafranov_factor=0.44789, triangularity=0.270, ion_temperature_beta=6, + fuel={"D":0.5,"T":0.5} ) # Tell OpenMC we're going to use our custom source diff --git a/src/openmc_plasma_source/__init__.py b/src/openmc_plasma_source/__init__.py index e7553a0..cd52fa2 100644 --- a/src/openmc_plasma_source/__init__.py +++ b/src/openmc_plasma_source/__init__.py @@ -14,3 +14,4 @@ from .tokamak_source import TokamakSource from .ring_source import FusionRingSource from .point_source import FusionPointSource +from .fuel_types import get_neutron_energy_distribution \ No newline at end of file diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 5e8856e..eb8ff17 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -5,7 +5,7 @@ def get_neutron_energy_distribution( ion_temperature: float, - fuel: dict = {"D": 0.5, "T": 0.5}, + fuel: dict, ) -> openmc.stats.Discrete: """Finds the energy distribution @@ -33,7 +33,7 @@ def get_neutron_energy_distribution( if sum_fuel_isotopes < 0.0: raise ValueError(f"isotope must sum to be above 0. Not {sum_fuel_isotopes}") - for k, v in fuel.dict: + for k, v in fuel.items(): if k not in ["D", "T"]: raise ValueError( f'Fuel dictionary keys must be either "D" or "T" not "{k}".' @@ -43,54 +43,49 @@ def get_neutron_energy_distribution( if v > 1: raise ValueError(f'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"] + if ['D'] == sorted(set(fuel.keys())): max_energy_mev = 5 - else: - nst.frac_D_default = 0 - - if "T" in fuel.keys(): - nst.frac_T_default = fuel["T"] + elif ['T'] == sorted(set(fuel.keys())): max_energy_mev = 12 - else: - nst.frac_T_default = 0 - - if "T" in fuel.keys() and "D" in fuel.keys(): + elif ['D', 'T'] == sorted(set(fuel.keys())): max_energy_mev = 20 + print(max_energy_mev,'MeV') + # 1.0 neutron yield, all reactions scaled by this value - num_of_vals = 500 + num_of_vals = 50 # single grid for DT, DD and TT grid E_pspec = np.linspace(0, max_energy_mev, 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(): + if ['D', 'T'] == sorted(set(fuel.keys())): DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) Y_DT = 1.0 - Y_DD = nst.yield_from_dt_yield_ratio("dd", Y_DT, ion_temperature) - Y_TT = nst.yield_from_dt_yield_ratio("tt", Y_DT, ion_temperature) + Y_DD = nst.yield_from_dt_yield_ratio("dd", Y_DT, ion_temperature, fuel["D"], fuel["T"]) + Y_TT = nst.yield_from_dt_yield_ratio("tt", Y_DT, ion_temperature, fuel["D"], fuel["T"]) dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar) # Brysk shape i.e. Gaussian dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) dNdE_DT_DD_TT = dNdE_DT + dNdE_DD + dNdE_TT - - if "D" in fuel.keys() and "T" not in fuel.keys(): + print('dt') + return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) + elif ['D'] == sorted(set(fuel.keys())): DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) + print('d') Y_DD = 1.0 dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian - dNdE_DT_DD_TT = dNdE_DD + return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DD) - if "D" not in fuel.keys() and "T" in fuel.keys(): + elif ['T'] == sorted(set(fuel.keys())): Y_TT = 1.0 + print('t') dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) - dNdE_DT_DD_TT = dNdE_TT + return openmc.stats.Discrete(E_pspec * 1e6, dNdE_TT) - return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 96f9395..ddfc27f 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -26,8 +26,7 @@ def __init__( # Set local attributes self.coordinate = coordinate self.temperature = temperature - self.fuel_type = fuel - self.fuel = fuel_types[self.fuel_type] + self.fuel = fuel # Call init for openmc.Source super().__init__() @@ -64,14 +63,3 @@ def temperature(self, value): self._temperature = value else: raise ValueError("Temperature must be strictly positive float.") - - @property - def fuel_type(self): - return self._fuel_type - - @fuel_type.setter - def fuel_type(self, value): - if value in fuel_types: - self._fuel_type = value - else: - raise KeyError("Invalid fuel type") From e5ea26f9816d6b1ab8bdad21695ff03837ff8e99 Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 30 Jan 2024 15:41:57 +0000 Subject: [PATCH 06/44] [skip ci] Apply formatting changes --- examples/point_source_example.py | 4 +++- examples/ring_source_example.py | 2 +- examples/tokamak_source_example.py | 2 +- src/openmc_plasma_source/__init__.py | 2 +- src/openmc_plasma_source/fuel_types.py | 29 ++++++++++++++------------ 5 files changed, 22 insertions(+), 17 deletions(-) diff --git a/examples/point_source_example.py b/examples/point_source_example.py index e8f8f14..032d6af 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -11,7 +11,9 @@ geometry = openmc.Geometry([cell]) # define the source -my_source = FusionPointSource(coordinate=(0, 0, 0), temperature=20000.0, fuel={"D":0.5, "T":0.5}) +my_source = FusionPointSource( + coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.5, "T": 0.5} +) # Tell OpenMC we're going to use our custom source settings = openmc.Settings() diff --git a/examples/ring_source_example.py b/examples/ring_source_example.py index 7c2b4ab..a6f7c28 100644 --- a/examples/ring_source_example.py +++ b/examples/ring_source_example.py @@ -16,7 +16,7 @@ radius=700, angles=(0.0, 2 * math.pi), # 360deg source temperature=20000.0, - fuel={"D":0.5,"T":0.5} + fuel={"D": 0.5, "T": 0.5}, ) settings = openmc.Settings() diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index d8eb53b..4bd202b 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -29,7 +29,7 @@ shafranov_factor=0.44789, triangularity=0.270, ion_temperature_beta=6, - fuel={"D":0.5,"T":0.5} + fuel={"D": 0.5, "T": 0.5}, ) # Tell OpenMC we're going to use our custom source diff --git a/src/openmc_plasma_source/__init__.py b/src/openmc_plasma_source/__init__.py index cd52fa2..b3ce3a6 100644 --- a/src/openmc_plasma_source/__init__.py +++ b/src/openmc_plasma_source/__init__.py @@ -14,4 +14,4 @@ from .tokamak_source import TokamakSource from .ring_source import FusionRingSource from .point_source import FusionPointSource -from .fuel_types import get_neutron_energy_distribution \ No newline at end of file +from .fuel_types import get_neutron_energy_distribution diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index eb8ff17..426d9f8 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -43,14 +43,14 @@ def get_neutron_energy_distribution( if v > 1: raise ValueError(f'Fuel dictionary values must be below 1 not "{k}".') - if ['D'] == sorted(set(fuel.keys())): + if ["D"] == sorted(set(fuel.keys())): max_energy_mev = 5 - elif ['T'] == sorted(set(fuel.keys())): + elif ["T"] == sorted(set(fuel.keys())): max_energy_mev = 12 - elif ['D', 'T'] == sorted(set(fuel.keys())): + elif ["D", "T"] == sorted(set(fuel.keys())): max_energy_mev = 20 - print(max_energy_mev,'MeV') + print(max_energy_mev, "MeV") # 1.0 neutron yield, all reactions scaled by this value num_of_vals = 50 @@ -58,34 +58,37 @@ def get_neutron_energy_distribution( E_pspec = np.linspace(0, max_energy_mev, num_of_vals) # accepts MeV units dNdE_DT_DD_TT = np.zeros(num_of_vals) - if ['D', 'T'] == sorted(set(fuel.keys())): + if ["D", "T"] == sorted(set(fuel.keys())): DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) Y_DT = 1.0 - Y_DD = nst.yield_from_dt_yield_ratio("dd", Y_DT, ion_temperature, fuel["D"], fuel["T"]) - Y_TT = nst.yield_from_dt_yield_ratio("tt", Y_DT, ion_temperature, fuel["D"], fuel["T"]) + Y_DD = nst.yield_from_dt_yield_ratio( + "dd", Y_DT, ion_temperature, fuel["D"], fuel["T"] + ) + Y_TT = nst.yield_from_dt_yield_ratio( + "tt", Y_DT, ion_temperature, fuel["D"], fuel["T"] + ) dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar) # Brysk shape i.e. Gaussian dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) dNdE_DT_DD_TT = dNdE_DT + dNdE_DD + dNdE_TT - print('dt') + print("dt") return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) - elif ['D'] == sorted(set(fuel.keys())): + elif ["D"] == sorted(set(fuel.keys())): DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) - print('d') + print("d") Y_DD = 1.0 dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DD) - elif ['T'] == sorted(set(fuel.keys())): + elif ["T"] == sorted(set(fuel.keys())): Y_TT = 1.0 - print('t') + print("t") dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) return openmc.stats.Discrete(E_pspec * 1e6, dNdE_TT) - From 293c544d53e0e1afc92e63caa23e541c7cde5f46 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 1 Mar 2024 18:07:24 +0000 Subject: [PATCH 07/44] [skip ci] refactor to functions to support lists --- examples/point_source_example.py | 24 +++++---- src/openmc_plasma_source/fuel_types.py | 32 ++++++----- src/openmc_plasma_source/point_source.py | 69 ++++++++++-------------- 3 files changed, 58 insertions(+), 67 deletions(-) diff --git a/examples/point_source_example.py b/examples/point_source_example.py index 032d6af..6019814 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -1,6 +1,7 @@ import openmc from openmc_plasma_source import FusionPointSource from pathlib import Path +import openmc_source_plotter # just making use of a local cross section xml file, replace with your own cross sections or comment out openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml" @@ -11,18 +12,19 @@ geometry = openmc.Geometry([cell]) # define the source -my_source = FusionPointSource( - coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.5, "T": 0.5} -) - +plot= my_source = FusionPointSource( + coordinate=(0, 0, 0), temperature=20000.0, fuel={"T": 1.} + # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.5, "T": 1.} +).plot_source_energy(n_samples=10000) +plot.show() # Tell OpenMC we're going to use our custom source -settings = openmc.Settings() -settings.run_mode = "fixed source" -settings.batches = 10 -settings.particles = 1000 -settings.source = my_source +# settings = openmc.Settings() +# settings.run_mode = "fixed source" +# settings.batches = 10 +# settings.particles = 1000 +# settings.source = my_source -model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) +# model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) -model.run() +# model.run() diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 426d9f8..3faf806 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -22,7 +22,7 @@ def get_neutron_energy_distribution( energy distribution """ - ion_temperature = ion_temperature / 1e3 # convert eV to keV + ion_temperature_kev = ion_temperature / 1e3 # convert eV to keV sum_fuel_isotopes = sum(fuel.values()) if sum_fuel_isotopes > 1.0: @@ -55,30 +55,34 @@ def get_neutron_energy_distribution( # 1.0 neutron yield, all reactions scaled by this value num_of_vals = 50 # single grid for DT, DD and TT grid - E_pspec = np.linspace(0, max_energy_mev, num_of_vals) # accepts MeV units + E_pspec = np.linspace(0, 12e6, num_of_vals) # accepts MeV units dNdE_DT_DD_TT = np.zeros(num_of_vals) if ["D", "T"] == sorted(set(fuel.keys())): - DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) - DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) + # DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) + # DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) Y_DT = 1.0 Y_DD = nst.yield_from_dt_yield_ratio( - "dd", Y_DT, ion_temperature, fuel["D"], fuel["T"] + "dd", Y_DT, ion_temperature_kev, fuel["D"], fuel["T"] ) Y_TT = nst.yield_from_dt_yield_ratio( - "tt", Y_DT, ion_temperature, fuel["D"], fuel["T"] + "tt", Y_DT, ion_temperature_kev, fuel["D"], fuel["T"] ) - dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar) # Brysk shape i.e. Gaussian - dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian - dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) - dNdE_DT_DD_TT = dNdE_DT + dNdE_DD + dNdE_TT + # dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar) # Brysk shape i.e. Gaussian + # dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian + dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) + # dNdE_DT_DD_TT = dNdE_DT + dNdE_DD + dNdE_TT + dNdE_DT_DD_TT = dNdE_TT print("dt") - return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) + tt_source = openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) + dd_source = openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) + dt_source = openmc.stats.muir(e0=14.06e6, m_rat=4, kt=ion_temperature) + return [tt_source, dd_source, dt_source] elif ["D"] == sorted(set(fuel.keys())): - DTmean, DTvar = nst.DTprimspecmoments(ion_temperature) - DDmean, DDvar = nst.DDprimspecmoments(ion_temperature) + DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) + DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) print("d") Y_DD = 1.0 @@ -90,5 +94,5 @@ def get_neutron_energy_distribution( Y_TT = 1.0 print("t") - dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature) + dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) return openmc.stats.Discrete(E_pspec * 1e6, dNdE_TT) diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index ddfc27f..02324d3 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -3,8 +3,11 @@ from .fuel_types import get_neutron_energy_distribution - -class FusionPointSource(openmc.IndependentSource): +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}, +): """An openmc.Source object with some presets to make it more convenient for fusion simulations using a point source. All attributes can be changed after initialization if required. Default isotropic point source at the @@ -17,49 +20,31 @@ class FusionPointSource(openmc.IndependentSource): fuel (dict): Isotopes as keys and atom fractions as values """ - def __init__( - self, - coordinate: Tuple[float, float, float] = (0.0, 0.0, 0.0), - temperature: float = 20000.0, - fuel: dict = {"D": 0.5, "T": 0.5}, + if ( + isinstance(coordinate, tuple) + and len(coordinate) == 3 + and all(isinstance(x, (int, float)) for x in coordinate) ): - # Set local attributes - self.coordinate = coordinate - self.temperature = temperature - self.fuel = fuel - - # Call init for openmc.Source - super().__init__() + pass + else: + raise ValueError("coordinate must be a tuple of three floats.") - # performed after the super init as these are Source attributes - self.space = openmc.stats.Point(self.coordinate) - self.angle = openmc.stats.Isotropic() - self.energy = get_neutron_energy_distribution( - ion_temperature=temperature, fuel=fuel - ) + if isinstance(temperature, (int, float)) and temperature > 0: + raise ValueError("Temperature must be a float.") + if temperature > 0: + raise ValueError("Temperature must be positive float.") - @property - def coordinate(self): - return self._coordinate - @coordinate.setter - def coordinate(self, value): - if ( - isinstance(value, tuple) - and len(value) == 3 - and all(isinstance(x, (int, float)) for x in value) - ): - self._coordinate = value - else: - raise ValueError("coordinate must be a tuple of three floats.") + sources = [] - @property - def temperature(self): - return self._temperature + energy_distributions = get_neutron_energy_distribution( + ion_temperature=temperature, fuel=fuel + ) + for energy_distribution in energy_distributions: + source = openmc.Source() + source.energy = energy_distribution + source.space = openmc.stats.Point(coordinate) + source.angle = openmc.stats.Isotropic() + sources.append(source) - @temperature.setter - def temperature(self, value): - if isinstance(value, (int, float)) and value > 0: - self._temperature = value - else: - raise ValueError("Temperature must be strictly positive float.") + return sources From 7633d919483192c03fcebe36ede65e16400f4a1c Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 1 Mar 2024 19:27:29 +0000 Subject: [PATCH 08/44] improved distributions --- examples/point_source_example.py | 36 +++++++++++++++--------- src/openmc_plasma_source/__init__.py | 2 +- src/openmc_plasma_source/fuel_types.py | 32 ++++++++++++--------- src/openmc_plasma_source/point_source.py | 18 +++++++++--- 4 files changed, 56 insertions(+), 32 deletions(-) diff --git a/examples/point_source_example.py b/examples/point_source_example.py index 6019814..adc9586 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -1,5 +1,5 @@ import openmc -from openmc_plasma_source import FusionPointSource +from openmc_plasma_source import fusion_point_source from pathlib import Path import openmc_source_plotter @@ -12,19 +12,29 @@ geometry = openmc.Geometry([cell]) # define the source -plot= my_source = FusionPointSource( - coordinate=(0, 0, 0), temperature=20000.0, fuel={"T": 1.} - # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.5, "T": 1.} -).plot_source_energy(n_samples=10000) -plot.show() +my_source = fusion_point_source( + # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.01, "T": 0.99} + coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.5, "T": 0.5} + # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.99, "T": 0.01} + # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 1.} + # coordinate=(0, 0, 0), temperature=20000.0, fuel={"T": 1.} +) + # Tell OpenMC we're going to use our custom source -# settings = openmc.Settings() -# settings.run_mode = "fixed source" -# settings.batches = 10 -# settings.particles = 1000 -# settings.source = my_source +settings = openmc.Settings() +settings.run_mode = "fixed source" +settings.batches = 10 +settings.particles = 1000 +settings.source = my_source +import numpy as np +plot = settings.plot_source_energy( + n_samples=100000, + energy_bins=np.linspace(0,16e6, 1000), + yaxis_type = 'log', +) +plot.show() -# model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) +model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) -# model.run() +model.run() diff --git a/src/openmc_plasma_source/__init__.py b/src/openmc_plasma_source/__init__.py index b3ce3a6..7da4fd0 100644 --- a/src/openmc_plasma_source/__init__.py +++ b/src/openmc_plasma_source/__init__.py @@ -13,5 +13,5 @@ from .tokamak_source import TokamakSource from .ring_source import FusionRingSource -from .point_source import FusionPointSource +from .point_source import fusion_point_source from .fuel_types import get_neutron_energy_distribution diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 3faf806..87e4162 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -43,19 +43,19 @@ def get_neutron_energy_distribution( if v > 1: raise ValueError(f'Fuel dictionary values must be below 1 not "{k}".') - if ["D"] == sorted(set(fuel.keys())): - max_energy_mev = 5 - elif ["T"] == sorted(set(fuel.keys())): - max_energy_mev = 12 - elif ["D", "T"] == sorted(set(fuel.keys())): - max_energy_mev = 20 + # if ["D"] == sorted(set(fuel.keys())): + # max_energy_mev = 5 + # elif ["T"] == sorted(set(fuel.keys())): + # max_energy_mev = 12 + # elif ["D", "T"] == sorted(set(fuel.keys())): + # max_energy_mev = 20 - print(max_energy_mev, "MeV") + # print(max_energy_mev, "MeV") # 1.0 neutron yield, all reactions scaled by this value - num_of_vals = 50 + num_of_vals = 100 # single grid for DT, DD and TT grid - E_pspec = np.linspace(0, 12e6, num_of_vals) # accepts MeV units + E_pspec = np.linspace(0, 12, num_of_vals) # E_pspec is exspected in MeV units dNdE_DT_DD_TT = np.zeros(num_of_vals) if ["D", "T"] == sorted(set(fuel.keys())): @@ -76,10 +76,12 @@ def get_neutron_energy_distribution( # dNdE_DT_DD_TT = dNdE_DT + dNdE_DD + dNdE_TT dNdE_DT_DD_TT = dNdE_TT print("dt") - tt_source = openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) + # tt_source = openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) + tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_DT_DD_TT) dd_source = openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) - dt_source = openmc.stats.muir(e0=14.06e6, m_rat=4, kt=ion_temperature) - return [tt_source, dd_source, dt_source] + dt_source = openmc.stats.muir(e0=14.06e6, m_rat=5, kt=ion_temperature) + # todo look into combining distributions openmc.data.combine_distributions() + return [tt_source, dd_source, dt_source], [Y_TT, Y_DD, Y_DT] elif ["D"] == sorted(set(fuel.keys())): DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) @@ -88,11 +90,13 @@ def get_neutron_energy_distribution( Y_DD = 1.0 dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian - return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DD) + # return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DD) + return openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) elif ["T"] == sorted(set(fuel.keys())): Y_TT = 1.0 print("t") dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) - return openmc.stats.Discrete(E_pspec * 1e6, dNdE_TT) + return openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) + # return openmc.stats.Discrete(E_pspec * 1e6, dNdE_TT) diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 02324d3..c076e4a 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -29,22 +29,32 @@ def fusion_point_source( else: raise ValueError("coordinate must be a tuple of three floats.") - if isinstance(temperature, (int, float)) and temperature > 0: + if not isinstance(temperature, (int, float)): raise ValueError("Temperature must be a float.") - if temperature > 0: + if temperature <= 0: raise ValueError("Temperature must be positive float.") sources = [] - energy_distributions = get_neutron_energy_distribution( + energy_distributions, strengths = get_neutron_energy_distribution( ion_temperature=temperature, fuel=fuel ) - for energy_distribution in energy_distributions: + + # if isinstance(energy_distributions, openmc.stats.Normal) or isinstance(energy_distributions, openmc.stats.Discrete) or isinstance(energy_distributions, openmc.stats.Tabular): + # source = openmc.Source() + # source.energy = energy_distributions + # source.space = openmc.stats.Point(coordinate) + # source.angle = openmc.stats.Isotropic() + # return source + + # else: + for energy_distribution, strength in zip(energy_distributions, strengths): source = openmc.Source() source.energy = energy_distribution source.space = openmc.stats.Point(coordinate) source.angle = openmc.stats.Isotropic() + source.strength = strength sources.append(source) return sources From 1eca213643235d089d84e32d8f8e50ffe81ef749 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 1 Mar 2024 19:29:40 +0000 Subject: [PATCH 09/44] [skip ci] Apply formatting changes --- examples/point_source_example.py | 9 ++++++--- src/openmc_plasma_source/point_source.py | 4 ++-- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/examples/point_source_example.py b/examples/point_source_example.py index adc9586..7afd154 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -14,7 +14,9 @@ # define the source my_source = fusion_point_source( # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.01, "T": 0.99} - coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.5, "T": 0.5} + coordinate=(0, 0, 0), + temperature=20000.0, + fuel={"D": 0.5, "T": 0.5}, # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.99, "T": 0.01} # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 1.} # coordinate=(0, 0, 0), temperature=20000.0, fuel={"T": 1.} @@ -28,10 +30,11 @@ settings.source = my_source import numpy as np + plot = settings.plot_source_energy( n_samples=100000, - energy_bins=np.linspace(0,16e6, 1000), - yaxis_type = 'log', + energy_bins=np.linspace(0, 16e6, 1000), + yaxis_type="log", ) plot.show() diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index c076e4a..80b3ec0 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -3,6 +3,7 @@ from .fuel_types import get_neutron_energy_distribution + def fusion_point_source( coordinate: Tuple[float, float, float] = (0.0, 0.0, 0.0), temperature: float = 20000.0, @@ -34,13 +35,12 @@ def fusion_point_source( if temperature <= 0: raise ValueError("Temperature must be positive float.") - sources = [] energy_distributions, strengths = get_neutron_energy_distribution( ion_temperature=temperature, fuel=fuel ) - + # if isinstance(energy_distributions, openmc.stats.Normal) or isinstance(energy_distributions, openmc.stats.Discrete) or isinstance(energy_distributions, openmc.stats.Tabular): # source = openmc.Source() # source.energy = energy_distributions From 53b7f985698a8930eb5903998e38ce5cfea3303f Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 5 Mar 2024 01:26:20 +0000 Subject: [PATCH 10/44] making use of new plotter api --- examples/point_source_example.py | 19 +++--- examples/ring_source_example.py | 6 +- examples/tokamak_source_example.py | 2 + src/openmc_plasma_source/__init__.py | 10 +-- src/openmc_plasma_source/fuel_types.py | 62 +++++++------------ .../plotting/plot_tokamak_source.py | 2 +- src/openmc_plasma_source/point_source.py | 3 +- src/openmc_plasma_source/ring_source.py | 5 +- src/openmc_plasma_source/tokamak_source.py | 5 +- tests/test_fuel_types.py | 3 +- tests/test_plotting.py | 7 +-- tests/test_point_source.py | 6 +- tests/test_ring_source.py | 6 +- tests/test_tokamak_source.py | 9 +-- 14 files changed, 65 insertions(+), 80 deletions(-) diff --git a/examples/point_source_example.py b/examples/point_source_example.py index 7afd154..d37d404 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -1,7 +1,10 @@ +from pathlib import Path +import numpy as np + import openmc +from openmc_source_plotter import plot_source_energy + from openmc_plasma_source import fusion_point_source -from pathlib import Path -import openmc_source_plotter # just making use of a local cross section xml file, replace with your own cross sections or comment out openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml" @@ -13,13 +16,9 @@ # define the source my_source = fusion_point_source( - # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.01, "T": 0.99} coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.5, "T": 0.5}, - # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 0.99, "T": 0.01} - # coordinate=(0, 0, 0), temperature=20000.0, fuel={"D": 1.} - # coordinate=(0, 0, 0), temperature=20000.0, fuel={"T": 1.} ) # Tell OpenMC we're going to use our custom source @@ -29,15 +28,11 @@ settings.particles = 1000 settings.source = my_source -import numpy as np -plot = settings.plot_source_energy( +plot = plot_source_energy( + this=settings, n_samples=100000, energy_bins=np.linspace(0, 16e6, 1000), yaxis_type="log", ) plot.show() - -model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) - -model.run() diff --git a/examples/ring_source_example.py b/examples/ring_source_example.py index a6f7c28..ee17cbf 100644 --- a/examples/ring_source_example.py +++ b/examples/ring_source_example.py @@ -1,8 +1,10 @@ -import openmc -from openmc_plasma_source import FusionRingSource import math from pathlib import Path +import openmc + +from openmc_plasma_source import FusionRingSource + # just making use of a local cross section xml file, replace with your own cross sections or comment out openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml" diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index 4bd202b..af973a5 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -1,5 +1,7 @@ from pathlib import Path + import openmc + from openmc_plasma_source import TokamakSource # just making use of a local cross section xml file, replace with your own cross sections or comment out diff --git a/src/openmc_plasma_source/__init__.py b/src/openmc_plasma_source/__init__.py index 7da4fd0..ac1a8b6 100644 --- a/src/openmc_plasma_source/__init__.py +++ b/src/openmc_plasma_source/__init__.py @@ -1,7 +1,7 @@ try: - from importlib.metadata import version, PackageNotFoundError + from importlib.metadata import PackageNotFoundError, version except (ModuleNotFoundError, ImportError): - from importlib_metadata import version, PackageNotFoundError + from importlib_metadata import PackageNotFoundError, version try: __version__ = version("openmc_plasma_source") except PackageNotFoundError: @@ -11,7 +11,7 @@ __all__ = ["__version__"] -from .tokamak_source import TokamakSource -from .ring_source import FusionRingSource -from .point_source import fusion_point_source from .fuel_types import get_neutron_energy_distribution +from .point_source import fusion_point_source +from .ring_source import FusionRingSource +from .tokamak_source import TokamakSource diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 87e4162..82bbaa1 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -43,60 +43,42 @@ def get_neutron_energy_distribution( if v > 1: raise ValueError(f'Fuel dictionary values must be below 1 not "{k}".') - # if ["D"] == sorted(set(fuel.keys())): - # max_energy_mev = 5 - # elif ["T"] == sorted(set(fuel.keys())): - # max_energy_mev = 12 - # elif ["D", "T"] == sorted(set(fuel.keys())): - # max_energy_mev = 20 - - # print(max_energy_mev, "MeV") - # 1.0 neutron yield, all reactions scaled by this value num_of_vals = 100 - # single grid for DT, DD and TT grid + # single grid for TT neutrons E_pspec = np.linspace(0, 12, num_of_vals) # E_pspec is exspected in MeV units - dNdE_DT_DD_TT = np.zeros(num_of_vals) + DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) + print('DTmean', DTmean) + DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) + if ["D", "T"] == sorted(set(fuel.keys())): - # DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) - # DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) - Y_DT = 1.0 - Y_DD = nst.yield_from_dt_yield_ratio( - "dd", Y_DT, ion_temperature_kev, fuel["D"], fuel["T"] + strength_DT = 1.0 + strength_DD = nst.yield_from_dt_yield_ratio( + "dd", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"] ) - Y_TT = nst.yield_from_dt_yield_ratio( - "tt", Y_DT, ion_temperature_kev, fuel["D"], fuel["T"] + strength_TT = nst.yield_from_dt_yield_ratio( + "tt", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"] ) - # dNdE_DT = Y_DT * nst.Qb(E_pspec, DTmean, DTvar) # Brysk shape i.e. Gaussian - # dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian - dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) - # dNdE_DT_DD_TT = dNdE_DT + dNdE_DD + dNdE_TT - dNdE_DT_DD_TT = dNdE_TT - print("dt") - # tt_source = openmc.stats.Discrete(E_pspec * 1e6, dNdE_DT_DD_TT) - tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_DT_DD_TT) + dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) + tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) dd_source = openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) dt_source = openmc.stats.muir(e0=14.06e6, m_rat=5, kt=ion_temperature) # todo look into combining distributions openmc.data.combine_distributions() - return [tt_source, dd_source, dt_source], [Y_TT, Y_DD, Y_DT] - elif ["D"] == sorted(set(fuel.keys())): - DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) - DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) + return [tt_source, dd_source, dt_source], [strength_TT, strength_DD, strength_DT] - print("d") - Y_DD = 1.0 + elif ["D"] == sorted(set(fuel.keys())): - dNdE_DD = Y_DD * nst.Qb(E_pspec, DDmean, DDvar) # Brysk shape i.e. Gaussian - # return openmc.stats.Discrete(E_pspec * 1e6, dNdE_DD) - return openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) + strength_DD = 1.0 + dd_source = openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) + return [dd_source], [strength_DD] elif ["T"] == sorted(set(fuel.keys())): - Y_TT = 1.0 - print("t") - dNdE_TT = Y_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) - return openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) - # return openmc.stats.Discrete(E_pspec * 1e6, dNdE_TT) + strength_TT = 1.0 + dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) + tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) + return [tt_source], [strength_TT] + diff --git a/src/openmc_plasma_source/plotting/plot_tokamak_source.py b/src/openmc_plasma_source/plotting/plot_tokamak_source.py index 74d25f5..107a2e1 100644 --- a/src/openmc_plasma_source/plotting/plot_tokamak_source.py +++ b/src/openmc_plasma_source/plotting/plot_tokamak_source.py @@ -1,6 +1,6 @@ import matplotlib.pyplot as plt -from matplotlib import cm import numpy as np +from matplotlib import cm def scatter_tokamak_source(source, ax=None, quantity=None, aspect="equal", **kwargs): diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 80b3ec0..646f27f 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -1,6 +1,7 @@ -import openmc from typing import Tuple +import openmc + from .fuel_types import get_neutron_energy_distribution diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index 051988f..cb66de4 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -1,7 +1,8 @@ -import openmc -import numpy as np from typing import Tuple +import numpy as np +import openmc + from .fuel_types import get_neutron_energy_distribution diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index bbbcf42..a0d58a1 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -1,7 +1,8 @@ -import openmc -import numpy as np from typing import Tuple +import numpy as np +import openmc + from .fuel_types import get_neutron_energy_distribution diff --git a/tests/test_fuel_types.py b/tests/test_fuel_types.py index 2d68717..110ff02 100644 --- a/tests/test_fuel_types.py +++ b/tests/test_fuel_types.py @@ -1,6 +1,7 @@ -from openmc_plasma_source.fuel_types import Fuel, fuel_types import pytest +from openmc_plasma_source.fuel_types import Fuel, fuel_types + @pytest.mark.parametrize("energy,mass", [(2.5e7, 5), (15, 30)]) def test_fuel_with_correct_inputs(energy, mass): diff --git a/tests/test_plotting.py b/tests/test_plotting.py index 7f1a5f7..4f11a15 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -1,11 +1,10 @@ import matplotlib.pyplot as plt import numpy as np -from openmc_plasma_source import ( - TokamakSource, - plotting as ops_plt, -) import pytest +from openmc_plasma_source import TokamakSource +from openmc_plasma_source import plotting as ops_plt + @pytest.fixture def tokamak_source(): diff --git a/tests/test_point_source.py b/tests/test_point_source.py index 2e4ab13..0c383d3 100644 --- a/tests/test_point_source.py +++ b/tests/test_point_source.py @@ -1,8 +1,8 @@ -from openmc_plasma_source import FusionPointSource - +import numpy as np import openmc import pytest -import numpy as np + +from openmc_plasma_source import FusionPointSource def test_creation(): diff --git a/tests/test_ring_source.py b/tests/test_ring_source.py index 35de338..a8c8158 100644 --- a/tests/test_ring_source.py +++ b/tests/test_ring_source.py @@ -1,8 +1,8 @@ -from openmc_plasma_source import FusionRingSource - +import numpy as np import openmc import pytest -import numpy as np + +from openmc_plasma_source import FusionRingSource def test_creation(): diff --git a/tests/test_tokamak_source.py b/tests/test_tokamak_source.py index 27379aa..8e557db 100644 --- a/tests/test_tokamak_source.py +++ b/tests/test_tokamak_source.py @@ -1,9 +1,10 @@ -from openmc_plasma_source import TokamakSource -from openmc import IndependentSource import numpy as np - import pytest -from hypothesis import given, settings, assume, strategies as st +from hypothesis import assume, given, settings +from hypothesis import strategies as st +from openmc import IndependentSource + +from openmc_plasma_source import TokamakSource @pytest.fixture From 42ae6a55db8e1216137bd6dc9137478978135da6 Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 5 Mar 2024 01:26:39 +0000 Subject: [PATCH 11/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/fuel_types.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 82bbaa1..8a1e26d 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -49,7 +49,7 @@ def get_neutron_energy_distribution( E_pspec = np.linspace(0, 12, num_of_vals) # E_pspec is exspected in MeV units DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) - print('DTmean', DTmean) + print("DTmean", DTmean) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) if ["D", "T"] == sorted(set(fuel.keys())): @@ -67,7 +67,11 @@ def get_neutron_energy_distribution( dd_source = openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) dt_source = openmc.stats.muir(e0=14.06e6, m_rat=5, kt=ion_temperature) # todo look into combining distributions openmc.data.combine_distributions() - return [tt_source, dd_source, dt_source], [strength_TT, strength_DD, strength_DT] + return [tt_source, dd_source, dt_source], [ + strength_TT, + strength_DD, + strength_DT, + ] elif ["D"] == sorted(set(fuel.keys())): @@ -81,4 +85,3 @@ def get_neutron_energy_distribution( dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) return [tt_source], [strength_TT] - From 4e68053eaf1aa95645fb38f7bccfab3f1c27e884 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 5 Mar 2024 01:51:34 +0000 Subject: [PATCH 12/44] converted to method instead of class --- src/openmc_plasma_source/ring_source.py | 139 ++++++++---------------- 1 file changed, 48 insertions(+), 91 deletions(-) diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index cb66de4..28236e9 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -6,7 +6,13 @@ from .fuel_types import get_neutron_energy_distribution -class FusionRingSource(openmc.IndependentSource): +def fusion_ring_source( + radius: float, + angles: Tuple[float, float] = (0, 2 * np.pi), + z_placement: float = 0, + temperature: float = 20000.0, + fuel: dict = {"D": 0.5, "T": 0.5}, +): """An openmc.Source object with some presets to make it more convenient for fusion simulations using a ring source. All attributes can be changed after initialization if required. Default isotropic ring source with a @@ -20,96 +26,47 @@ class FusionRingSource(openmc.IndependentSource): fuel (dict): Isotopes as keys and atom fractions as values """ - def __init__( - self, - radius: float, - angles: Tuple[float, float] = (0, 2 * np.pi), - z_placement: float = 0, - temperature: float = 20000.0, - fuel: dict = {"D": 0.5, "T": 0.5}, - ): - # Set local attributes - self.radius = radius - self.angles = angles - self.z_placement = z_placement - self.temperature = temperature - self.fuel = fuel - - # Call init for openmc.Source - super().__init__() - - # performed after the super init as these are Source attributes - self.space = openmc.stats.CylindricalIndependent( - r=openmc.stats.Discrete([self.radius], [1]), - phi=openmc.stats.Uniform(a=self.angles[0], b=self.angles[1]), - z=openmc.stats.Discrete([self.z_placement], [1]), - origin=(0.0, 0.0, 0.0), + if isinstance(radius, (int, float)) and radius > 0: + pass + else: + raise ValueError("Radius must be a float strictly greater than 0.") + + if ( + isinstance(angles, tuple) + and len(angles) == 2 + and all( + isinstance(angle, (int, float)) and -2 * np.pi <= angle <= 2 * np.pi + for angle in angles ) - self.angle = openmc.stats.Isotropic() - self.energy = get_neutron_energy_distribution( - ion_temperature=temperature, fuel=fuel + ): + pass + else: + raise ValueError( + "Angles must be a tuple of floats between zero and 2 * np.pi" ) - @property - def radius(self): - return self._radius - - @radius.setter - def radius(self, value): - if isinstance(value, (int, float)) and value > 0: - self._radius = value - else: - raise ValueError("Radius must be a float strictly greater than 0.") - - @property - def angles(self): - return self._angles - - @angles.setter - def angles(self, value): - if ( - isinstance(value, tuple) - and len(value) == 2 - and all( - isinstance(angle, (int, float)) and -2 * np.pi <= angle <= 2 * np.pi - for angle in value - ) - ): - self._angles = value - else: - raise ValueError( - "Angles must be a tuple of floats between zero and 2 * np.pi" - ) - - @property - def z_placement(self): - return self._z_placement - - @z_placement.setter - def z_placement(self, value): - if isinstance(value, (int, float)): - self._z_placement = value - else: - raise TypeError("Z placement must be a float.") - - @property - def temperature(self): - return self._temperature - - @temperature.setter - def temperature(self, value): - if isinstance(value, (int, float)) and value > 0: - self._temperature = value - else: - raise ValueError("Temperature must be a float strictly greater than 0.") - - @property - def fuel_type(self): - return self._fuel_type - - @fuel_type.setter - def fuel_type(self, value): - if value in fuel_types.keys(): - self._fuel_type = value - else: - raise KeyError("Invalid fuel type.") + if isinstance(z_placement, (int, float)): + pass + else: + raise TypeError("Z placement must be a float.") + + if isinstance(temperature, (int, float)) and temperature > 0: + pass + else: + raise ValueError("Temperature must be a float strictly greater than 0.") + + source = openmc.IndependentSource() + + # performed after the super init as these are Source attributes + source.space = openmc.stats.CylindricalIndependent( + r=openmc.stats.Discrete([radius], [1]), + phi=openmc.stats.Uniform(a=angles[0], b=angles[1]), + z=openmc.stats.Discrete([z_placement], [1]), + origin=(0.0, 0.0, 0.0), + ) + source.angle = openmc.stats.Isotropic() + source.energy = get_neutron_energy_distribution( + ion_temperature=temperature, fuel=fuel + ) + + return source From 61112ad92dc41da6966a523221b6377aed831cf2 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 5 Mar 2024 02:01:02 +0000 Subject: [PATCH 13/44] started ts --- src/openmc_plasma_source/tokamak_source.py | 92 ++++++++-------------- 1 file changed, 35 insertions(+), 57 deletions(-) diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index a0d58a1..0506bdf 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -6,7 +6,27 @@ from .fuel_types import get_neutron_energy_distribution -class TokamakSource: +def tokamak_source( + major_radius: float, + minor_radius: float, + elongation: float, + triangularity: float, + mode: str, + ion_density_centre: float, + ion_density_peaking_factor: float, + ion_density_pedestal: float, + ion_density_separatrix: float, + ion_temperature_centre: float, + ion_temperature_peaking_factor: float, + ion_temperature_beta: float, + ion_temperature_pedestal: float, + ion_temperature_separatrix: float, + pedestal_radius: float, + shafranov_factor: float, + angles: Tuple[float, float] = (0, 2 * np.pi), + sample_size: int = 1000, + fuel: dict = {"D": 0.5, "T": 0.5}, +): """Plasma neutron source sampling. This class greatly relies on models described in [1] @@ -52,62 +72,20 @@ class TokamakSource: fuel (dict): Isotopes as keys and atom fractions as values """ - def __init__( - self, - major_radius: float, - minor_radius: float, - elongation: float, - triangularity: float, - mode: str, - ion_density_centre: float, - ion_density_peaking_factor: float, - ion_density_pedestal: float, - ion_density_separatrix: float, - ion_temperature_centre: float, - ion_temperature_peaking_factor: float, - ion_temperature_beta: float, - ion_temperature_pedestal: float, - ion_temperature_separatrix: float, - pedestal_radius: float, - shafranov_factor: float, - angles: Tuple[float, float] = (0, 2 * np.pi), - sample_size: int = 1000, - fuel: dict = {"D": 0.5, "T": 0.5}, - ) -> None: - # Assign attributes - self.major_radius = major_radius - self.minor_radius = minor_radius - self.elongation = elongation - self.triangularity = triangularity - self.ion_density_centre = ion_density_centre - self.ion_density_peaking_factor = ion_density_peaking_factor - self.mode = mode - self.pedestal_radius = pedestal_radius - self.ion_density_pedestal = ion_density_pedestal - self.ion_density_separatrix = ion_density_separatrix - self.ion_temperature_centre = ion_temperature_centre - self.ion_temperature_peaking_factor = ion_temperature_peaking_factor - self.ion_temperature_pedestal = ion_temperature_pedestal - self.ion_temperature_separatrix = ion_temperature_separatrix - self.ion_temperature_beta = ion_temperature_beta - self.shafranov_factor = shafranov_factor - self.angles = angles - self.sample_size = sample_size - self.fuel = fuel - - # Perform sanity checks for inputs not caught by properties - if self.minor_radius >= self.major_radius: - raise ValueError("Minor radius must be smaller than major radius") - - if self.pedestal_radius >= self.minor_radius: - raise ValueError("Pedestal radius must be smaller than minor radius") - - if abs(self.shafranov_factor) >= 0.5 * minor_radius: - raise ValueError("Shafranov factor must be smaller than 0.5*minor radius") - - # Create a list of souces - self.sample_sources() - self.sources = self.make_openmc_sources() + + # Perform sanity checks for inputs not caught by properties + if minor_radius >= major_radius: + raise ValueError("Minor radius must be smaller than major radius") + + if pedestal_radius >= minor_radius: + raise ValueError("Pedestal radius must be smaller than minor radius") + + if abs(shafranov_factor) >= 0.5 * minor_radius: + raise ValueError("Shafranov factor must be smaller than 0.5*minor radius") + + # Create a list of souces + self.sample_sources() + self.sources = self.make_openmc_sources() @property def major_radius(self): From 6dfd64c01772f22e810fd4fd4a75450a3d19f8ea Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 5 Mar 2024 02:01:53 +0000 Subject: [PATCH 14/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/ring_source.py | 4 +--- src/openmc_plasma_source/tokamak_source.py | 1 - 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index 28236e9..996a633 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -41,9 +41,7 @@ def fusion_ring_source( ): pass else: - raise ValueError( - "Angles must be a tuple of floats between zero and 2 * np.pi" - ) + raise ValueError("Angles must be a tuple of floats between zero and 2 * np.pi") if isinstance(z_placement, (int, float)): pass diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 0506bdf..05a26b5 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -72,7 +72,6 @@ def tokamak_source( fuel (dict): Isotopes as keys and atom fractions as values """ - # Perform sanity checks for inputs not caught by properties if minor_radius >= major_radius: raise ValueError("Minor radius must be smaller than major radius") From a70b0188b4651a67a10b36b7e22f338b089b8a0b Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 6 Mar 2024 13:44:46 +0000 Subject: [PATCH 15/44] adapted tests for new api --- src/openmc_plasma_source/__init__.py | 4 +- src/openmc_plasma_source/fuel_types.py | 10 ++-- tests/test_fuel_types.py | 49 ++++++++++-------- tests/test_plotting.py | 12 ++--- tests/test_point_source.py | 16 +++--- tests/test_ring_source.py | 22 ++++---- tests/test_tokamak_source.py | 70 ++++++++++++-------------- 7 files changed, 92 insertions(+), 91 deletions(-) diff --git a/src/openmc_plasma_source/__init__.py b/src/openmc_plasma_source/__init__.py index ac1a8b6..8a1ddda 100644 --- a/src/openmc_plasma_source/__init__.py +++ b/src/openmc_plasma_source/__init__.py @@ -13,5 +13,5 @@ from .fuel_types import get_neutron_energy_distribution from .point_source import fusion_point_source -from .ring_source import FusionRingSource -from .tokamak_source import TokamakSource +from .ring_source import fusion_ring_source +from .tokamak_source import tokamak_source diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 8a1e26d..45e6ac6 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -49,8 +49,10 @@ def get_neutron_energy_distribution( E_pspec = np.linspace(0, 12, num_of_vals) # E_pspec is exspected in MeV units DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) - print("DTmean", DTmean) + # print("DTmean", DTmean) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) + # print("DDmean", DDmean) + #todo make use of these DTvar and DDvar values in muir or gaussian distribution if ["D", "T"] == sorted(set(fuel.keys())): @@ -64,8 +66,8 @@ def get_neutron_energy_distribution( dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) - dd_source = openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) - dt_source = openmc.stats.muir(e0=14.06e6, m_rat=5, kt=ion_temperature) + dd_source = openmc.stats.muir(e0=DDmean*1e6, m_rat=4, kt=ion_temperature) + dt_source = openmc.stats.muir(e0=DTmean*1e6, m_rat=5, kt=ion_temperature) # todo look into combining distributions openmc.data.combine_distributions() return [tt_source, dd_source, dt_source], [ strength_TT, @@ -76,7 +78,7 @@ def get_neutron_energy_distribution( elif ["D"] == sorted(set(fuel.keys())): strength_DD = 1.0 - dd_source = openmc.stats.muir(e0=2.5e6, m_rat=4, kt=ion_temperature) + dd_source = openmc.stats.muir(e0=DDmean*1e6, m_rat=4, kt=ion_temperature) return [dd_source], [strength_DD] elif ["T"] == sorted(set(fuel.keys())): diff --git a/tests/test_fuel_types.py b/tests/test_fuel_types.py index 110ff02..17dc1c5 100644 --- a/tests/test_fuel_types.py +++ b/tests/test_fuel_types.py @@ -1,33 +1,38 @@ import pytest -from openmc_plasma_source.fuel_types import Fuel, fuel_types +from openmc_plasma_source import get_neutron_energy_distribution -@pytest.mark.parametrize("energy,mass", [(2.5e7, 5), (15, 30)]) -def test_fuel_with_correct_inputs(energy, mass): +@pytest.mark.parametrize("temperature, fuel", [ + (2e3, {'D': 1.}), + (2e3, {'T': 1.}), + (2e3, {'T': 0.5, 'D': 0.5}), + (2e3, {'T': 0.2, 'D': 0.8}), +]) +def test_fuel_with_correct_inputs(temperature, fuel): # Should accept any non-zero positive inputs to either variable - fuel = Fuel(energy, mass) - assert fuel.mean_energy == energy - assert fuel.mass_of_reactants == mass + get_neutron_energy_distribution(temperature, fuel) -@pytest.mark.parametrize( - "energy,mass", [(2.5e7, -5), (-12, 30), (1e7, 0), (0, 4), (-12, -12)] -) -def test_fuel_with_bad_inputs(energy, mass): +@pytest.mark.parametrize("temperature, fuel", [ + (2e3, {'D': 1.1}), + (2e3, {'T': 0.9}), + (2e3, {'T': -0.5, 'D': 0.5}), + (2e3, {'T': -0.2, 'D': -0.8}), +]) +def test_fuel_with_bad_inputs(temperature, fuel): # Should reject any negative numbers and zeros. with pytest.raises(ValueError): - fuel = Fuel(energy, mass) + get_neutron_energy_distribution(temperature, fuel) -@pytest.mark.parametrize("fuel_type", ["DT", "DD"]) -def test_fuel_types(fuel_type): - # Should accept 'DD' and 'DT' - assert isinstance(fuel_types[fuel_type], Fuel) - - -@pytest.mark.parametrize("fuel_type", ["dt", "dd", "Dt", "dD", "hello world", 5]) -def test_incorrect_fuel_types(fuel_type): - # Should reject everything except 'DT' and 'DD' - with pytest.raises(KeyError): - my_fuel = fuel_types[fuel_type] +@pytest.mark.parametrize("temperature, fuel", [ + (2e3, {'DD': 1.1}), + (2e3, {'DT': 0.9}), + (2e3, {'He3': -0.5, 'D': 0.5}), + (2e3, {1: -0.2, 'D': -0.8}), +]) +def test_fuel_with_incorrect_isotopese(temperature, fuel): + # Should reject anything which is not 'D' or 'T'. + with pytest.raises(ValueError): + get_neutron_energy_distribution(temperature, fuel) diff --git a/tests/test_plotting.py b/tests/test_plotting.py index 4f11a15..8027ea0 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -2,13 +2,13 @@ import numpy as np import pytest -from openmc_plasma_source import TokamakSource +from openmc_plasma_source import tokamak_source from openmc_plasma_source import plotting as ops_plt @pytest.fixture def tokamak_source(): - return TokamakSource( + return tokamak_source( elongation=1.557, ion_density_centre=1.09e20, ion_density_peaking_factor=1, @@ -128,13 +128,13 @@ def test_scatter_tokamak_source_kwargs(tokamak_source, kwargs): def test_scatter_tokamak_not_source(): - """Ensure failure when given non-TokamakSource to plot""" + """Ensure failure when given non-tokamak_source to plot""" with pytest.raises(ValueError) as excinfo: fig = plt.figure() ax = fig.gca() ops_plt.scatter_tokamak_source("hello world", ax=ax) plt.close() - assert "TokamakSource" in str(excinfo.value) + assert "tokamak_source" in str(excinfo.value) @pytest.mark.parametrize("quantity", ["coucou", "ion_density", 17]) @@ -241,13 +241,13 @@ def test_plot_tokamak_source_3D_kwargs(tokamak_source, kwargs): def test_plot_tokamak_source_3D_not_source(): - """Ensure failure when given non-TokamakSource to plot""" + """Ensure failure when given non-tokamak_source to plot""" with pytest.raises(ValueError) as excinfo: fig = plt.figure() ax = fig.add_subplot(1, 1, 1, projection="3d") ops_plt.plot_tokamak_source_3D("hello world", ax=ax) plt.close() - assert "TokamakSource" in str(excinfo.value) + assert "tokamak_source" in str(excinfo.value) @pytest.mark.parametrize("quantity", ["coucou", "ion_density", 17]) diff --git a/tests/test_point_source.py b/tests/test_point_source.py index 0c383d3..9533f66 100644 --- a/tests/test_point_source.py +++ b/tests/test_point_source.py @@ -2,11 +2,11 @@ import openmc import pytest -from openmc_plasma_source import FusionPointSource +from openmc_plasma_source import fusion_point_source def test_creation(): - my_source = FusionPointSource() + my_source = fusion_point_source() # Ensure it is of type openmc.IndependentSource assert isinstance(my_source, openmc.IndependentSource) @@ -22,7 +22,7 @@ def test_creation(): ) def test_coordinate(coordinate): # Should allow any tuple of length 3 containing numbers - my_source = FusionPointSource(coordinate=coordinate) + my_source = fusion_point_source(coordinate=coordinate) assert np.array_equal(my_source.coordinate, coordinate) @@ -33,13 +33,13 @@ def test_bad_coordinate(coordinate): # Should reject iterables of length != 3, anything non-tuple, and anything # that can't convert to float with pytest.raises(ValueError): - FusionPointSource(coordinate=coordinate) + fusion_point_source(coordinate=coordinate) @pytest.mark.parametrize("temperature", [20000, 1e4, 0.1, 25000.0]) def test_temperature(temperature): # Should accept any positive float - my_source = FusionPointSource(temperature=temperature) + my_source = fusion_point_source(temperature=temperature) assert my_source.temperature == temperature @@ -47,13 +47,13 @@ def test_temperature(temperature): def test_bad_temperature(temperature): # Should reject negative floats and anything that isn't convertible to float with pytest.raises(ValueError): - FusionPointSource(temperature=temperature) + fusion_point_source(temperature=temperature) @pytest.mark.parametrize("fuel", ["DT", "DD"]) def test_fuel(fuel): # Should accept either 'DD' or 'DT' - my_source = FusionPointSource(fuel=fuel) + my_source = fusion_point_source(fuel=fuel) assert my_source.fuel_type == fuel @@ -61,4 +61,4 @@ def test_fuel(fuel): def test_wrong_fuel(fuel): # Should reject fuel types besides those listed in fuel_types.py with pytest.raises((KeyError, TypeError)): - FusionPointSource(fuel=fuel) + fusion_point_source(fuel=fuel) diff --git a/tests/test_ring_source.py b/tests/test_ring_source.py index a8c8158..1c6d758 100644 --- a/tests/test_ring_source.py +++ b/tests/test_ring_source.py @@ -2,11 +2,11 @@ import openmc import pytest -from openmc_plasma_source import FusionRingSource +from openmc_plasma_source import fusion_ring_source def test_creation(): - my_source = FusionRingSource(radius=1.0, z_placement=1.0) + my_source = fusion_ring_source(radius=1.0, z_placement=1.0) # Ensure it is of type openmc.IndependentSource assert isinstance(my_source, openmc.IndependentSource) @@ -20,7 +20,7 @@ def test_creation(): @pytest.mark.parametrize("radius", [1, 5.6, 1e5, 7.0]) def test_radius(radius): # should allow any positive float - my_source = FusionRingSource(radius=radius) + my_source = fusion_ring_source(radius=radius) assert my_source.radius == radius @@ -28,13 +28,13 @@ def test_radius(radius): def test_bad_radius(radius): # should reject any negative float or anything not convertible to float with pytest.raises(ValueError): - my_source = FusionRingSource(radius=radius) + fusion_ring_source(radius=radius) @pytest.mark.parametrize("angles", [(1, 2), (0.0, np.pi), (np.pi, 0.0)]) def test_angles(angles): # Should allow any tuple of length 2 with contents convertible to float - my_source = FusionRingSource(radius=1.0, angles=angles) + my_source = fusion_ring_source(radius=1.0, angles=angles) assert np.array_equal(my_source.angles, angles) @@ -43,13 +43,13 @@ def test_bad_angles(angles): # Should reject iterables of length != 2, anything non tuple, and anything # that can't convert to float with pytest.raises(ValueError): - FusionRingSource(radius=1.0, angles=angles) + fusion_ring_source(radius=1.0, angles=angles) @pytest.mark.parametrize("temperature", [20000.0, 1e4, 0.1, 25000]) def test_temperature(temperature): # Should accept any positive float - my_source = FusionRingSource(radius=1.0, temperature=temperature) + my_source = fusion_ring_source(radius=1.0, temperature=temperature) assert my_source.temperature == temperature @@ -57,13 +57,13 @@ def test_temperature(temperature): def test_bad_temperature(temperature): # Should reject negative floats and anything that isn't convertible to float with pytest.raises(ValueError): - FusionRingSource(radius=1.0, temperature=temperature) + fusion_ring_source(radius=1.0, temperature=temperature) @pytest.mark.parametrize("fuel", ["DT", "DD"]) def test_fuel(fuel): # Should accept either 'DD' or 'DT' - my_source = FusionRingSource(radius=1.0, fuel=fuel) + my_source = fusion_ring_source(radius=1.0, fuel=fuel) assert my_source.fuel_type == fuel @@ -71,10 +71,10 @@ def test_fuel(fuel): def test_wrong_fuel(fuel): # Should reject fuel types besides those listed in fuel_types.py with pytest.raises((KeyError, TypeError)): - FusionRingSource(radius=1.0, fuel=fuel) + fusion_ring_source(radius=1.0, fuel=fuel) @pytest.mark.parametrize("z", ["coucou", [5, 2.0]]) def test_wrong_z_placement(z): with pytest.raises((TypeError)): - FusionRingSource(radius=1.0, z_placement=z) + fusion_ring_source(radius=1.0, z_placement=z) diff --git a/tests/test_tokamak_source.py b/tests/test_tokamak_source.py index 8e557db..28845c4 100644 --- a/tests/test_tokamak_source.py +++ b/tests/test_tokamak_source.py @@ -1,15 +1,14 @@ import numpy as np import pytest -from hypothesis import assume, given, settings +from hypothesis import given, settings from hypothesis import strategies as st -from openmc import IndependentSource - -from openmc_plasma_source import TokamakSource +import openmc +from openmc_plasma_source import tokamak_source @pytest.fixture def tokamak_args_dict(): - """Returns a dict of realistic inputs for TokamakSource""" + """Returns a dict of realistic inputs for tokamak_source""" args_dict = { "elongation": 1.557, "triangularity": 0.270, @@ -33,51 +32,49 @@ def tokamak_args_dict(): @pytest.fixture def tokamak_source_example(tokamak_args_dict): - """Returns a TokamakSource with realistic inputs""" - return TokamakSource(**tokamak_args_dict) + """Returns a tokamak_source with realistic inputs""" + return tokamak_source(**tokamak_args_dict) def test_creation(tokamak_source_example): - """Tests that the sources generated by TokamakSource are of + """Tests that the sources generated by tokamak_source are of type openmc.Source""" for source in tokamak_source_example.sources: - assert isinstance(source, IndependentSource) + assert isinstance(source, openmc.IndependentSource) @pytest.mark.parametrize( "minor_radius,major_radius", [(3.0, 10.0), (3.0, 100), (3.0, 3.00001)] ) def test_major_radius(tokamak_args_dict, minor_radius, major_radius): - """Checks that TokamakSource creation accepts valid major radius""" + """Checks that tokamak_source creation accepts valid major radius""" tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["major_radius"] = major_radius - tokamak_source = TokamakSource(**tokamak_args_dict) - assert tokamak_source.major_radius == major_radius + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( "minor_radius,major_radius", [(3, 3), (3, 1), (3, -5), (3, "hello world")] ) def test_bad_major_radius(tokamak_args_dict, minor_radius, major_radius): - """Checks that TokamakSource creation rejects invalid major radius""" + """Checks that tokamak_source creation rejects invalid major radius""" tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["major_radius"] = major_radius with pytest.raises(ValueError): - tokamak_source = TokamakSource(**tokamak_args_dict) + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( "major_radius,minor_radius", [(10.0, 3.0), (10.0, 9.9), (10.0, 0.1)] ) def test_minor_radius(tokamak_args_dict, major_radius, minor_radius): - """Checks that TokamakSource creation accepts valid minor radius""" + """Checks that tokamak_source creation accepts valid minor radius""" tokamak_args_dict["major_radius"] = major_radius tokamak_args_dict["minor_radius"] = minor_radius # Set shafranov factor to 0 and pedestal factor to 0.8*minor_radius for safety tokamak_args_dict["pedestal_radius"] = 0.8 * minor_radius tokamak_args_dict["shafranov_factor"] = 0.0 - tokamak_source = TokamakSource(**tokamak_args_dict) - assert tokamak_source.minor_radius == minor_radius + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -85,43 +82,41 @@ def test_minor_radius(tokamak_args_dict, major_radius, minor_radius): [(10.0, 10.0), (10.0, 20.0), (10.0, 0), (10.0, -6), (10.0, "hello world")], ) def test_bad_minor_radius(tokamak_args_dict, major_radius, minor_radius): - """Checks that TokamakSource creation rejects invalid minor radius""" + """Checks that tokamak_source creation rejects invalid minor radius""" tokamak_args_dict["major_radius"] = major_radius tokamak_args_dict["minor_radius"] = minor_radius with pytest.raises(ValueError): - tokamak_source = TokamakSource(**tokamak_args_dict) + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("elongation", [1.0, 1.667, 0.5, 20, 0.001]) def test_elongation(tokamak_args_dict, elongation): - """Checks that TokamakSource creation accepts valid elongation""" + """Checks that tokamak_source creation accepts valid elongation""" tokamak_args_dict["elongation"] = elongation - tokamak_source = TokamakSource(**tokamak_args_dict) - assert tokamak_source.elongation == elongation + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("elongation", [0, -5, "hello world"]) def test_bad_elongation(tokamak_args_dict, elongation): - """Checks that TokamakSource creation rejects invalid elongation""" + """Checks that tokamak_source creation rejects invalid elongation""" tokamak_args_dict["elongation"] = elongation with pytest.raises(ValueError): - tokamak_source = TokamakSource(**tokamak_args_dict) + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("triangularity", [0.0, 0.5, 0.9, 1.0, -0.5, -0.9, -1.0]) def test_triangularity(tokamak_args_dict, triangularity): - """Checks that TokamakSource creation accepts valid triangularity""" + """Checks that tokamak_source creation accepts valid triangularity""" tokamak_args_dict["triangularity"] = triangularity - tokamak_source = TokamakSource(**tokamak_args_dict) - assert tokamak_source.triangularity == triangularity + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("triangularity", [1.1, -1.1, 10, -10, "hello world"]) def test_bad_triangularity(tokamak_args_dict, triangularity): - """Checks that TokamakSource creation rejects invalid triangularity""" + """Checks that tokamak_source creation rejects invalid triangularity""" tokamak_args_dict["triangularity"] = triangularity with pytest.raises(ValueError): - tokamak_source = TokamakSource(**tokamak_args_dict) + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -137,13 +132,12 @@ def test_bad_triangularity(tokamak_args_dict, triangularity): ], ) def test_shafranov_factor(tokamak_args_dict, major_radius, minor_radius, shaf): - """Checks that TokamakSource creation accepts valid Shafranov factor""" + """Checks that tokamak_source creation accepts valid Shafranov factor""" tokamak_args_dict["major_radius"] = major_radius tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["pedestal_radius"] = 0.8 * minor_radius tokamak_args_dict["shafranov_factor"] = shaf - tokamak_source = TokamakSource(**tokamak_args_dict) - assert tokamak_source.shafranov_factor == shaf + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -158,13 +152,13 @@ def test_shafranov_factor(tokamak_args_dict, major_radius, minor_radius, shaf): ], ) def test_bad_shafranov_factor(tokamak_args_dict, major_radius, minor_radius, shaf): - """Checks that TokamakSource creation rejects invalid Shafranov factor""" + """Checks that tokamak_source creation rejects invalid Shafranov factor""" tokamak_args_dict["major_radius"] = major_radius tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["pedestal_radius"] = 0.8 * minor_radius tokamak_args_dict["shafranov_factor"] = shaf with pytest.raises((ValueError, TypeError)): - tokamak_source = TokamakSource(**tokamak_args_dict) + tokamak_source = tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -174,7 +168,7 @@ def test_angles(tokamak_args_dict, angles): """Checks that custom angles can be set""" # Note: should accept negative angles and angles in reverse order tokamak_args_dict["angles"] = angles - tokamak_source = TokamakSource(**tokamak_args_dict) + tokamak_source = tokamak_source(**tokamak_args_dict) assert np.array_equal(tokamak_source.angles, angles) for source in tokamak_source.sources: assert np.array_equal((source.space.phi.a, source.space.phi.b), angles) @@ -187,7 +181,7 @@ def test_bad_angles(tokamak_args_dict, angles): # Contents should convert to float tokamak_args_dict["angles"] = angles with pytest.raises(ValueError) as excinfo: - tokamak_source = TokamakSource(**tokamak_args_dict) + tokamak_source = tokamak_source(**tokamak_args_dict) def test_ion_density(tokamak_source_example): @@ -248,7 +242,7 @@ def test_bad_convert_a_alpha_to_R_Z(tokamak_source_example): @st.composite def tokamak_source_strategy(draw): - """Defines a hypothesis strategy that automatically generates a TokamakSource. + """Defines a hypothesis strategy that automatically generates a tokamak_source. Geometry attributes are varied, while plasma attributes are fixed. """ # Used to avoid generation of inappropriate float values @@ -293,7 +287,7 @@ def tokamak_source_strategy(draw): ) ) - return TokamakSource( + return tokamak_source( elongation=elongation, triangularity=triangularity, major_radius=major_radius, From 1188867b0c7d2b8512287efe951bf1981fab71d1 Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 6 Mar 2024 13:45:28 +0000 Subject: [PATCH 16/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/fuel_types.py | 8 ++--- tests/test_fuel_types.py | 45 +++++++++++++++----------- 2 files changed, 31 insertions(+), 22 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 45e6ac6..2b50426 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -52,7 +52,7 @@ def get_neutron_energy_distribution( # print("DTmean", DTmean) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) # print("DDmean", DDmean) - #todo make use of these DTvar and DDvar values in muir or gaussian distribution + # todo make use of these DTvar and DDvar values in muir or gaussian distribution if ["D", "T"] == sorted(set(fuel.keys())): @@ -66,8 +66,8 @@ def get_neutron_energy_distribution( dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) - dd_source = openmc.stats.muir(e0=DDmean*1e6, m_rat=4, kt=ion_temperature) - dt_source = openmc.stats.muir(e0=DTmean*1e6, m_rat=5, kt=ion_temperature) + dd_source = openmc.stats.muir(e0=DDmean * 1e6, m_rat=4, kt=ion_temperature) + dt_source = openmc.stats.muir(e0=DTmean * 1e6, m_rat=5, kt=ion_temperature) # todo look into combining distributions openmc.data.combine_distributions() return [tt_source, dd_source, dt_source], [ strength_TT, @@ -78,7 +78,7 @@ def get_neutron_energy_distribution( elif ["D"] == sorted(set(fuel.keys())): strength_DD = 1.0 - dd_source = openmc.stats.muir(e0=DDmean*1e6, m_rat=4, kt=ion_temperature) + dd_source = openmc.stats.muir(e0=DDmean * 1e6, m_rat=4, kt=ion_temperature) return [dd_source], [strength_DD] elif ["T"] == sorted(set(fuel.keys())): diff --git a/tests/test_fuel_types.py b/tests/test_fuel_types.py index 17dc1c5..b0a2c8a 100644 --- a/tests/test_fuel_types.py +++ b/tests/test_fuel_types.py @@ -3,35 +3,44 @@ from openmc_plasma_source import get_neutron_energy_distribution -@pytest.mark.parametrize("temperature, fuel", [ - (2e3, {'D': 1.}), - (2e3, {'T': 1.}), - (2e3, {'T': 0.5, 'D': 0.5}), - (2e3, {'T': 0.2, 'D': 0.8}), -]) +@pytest.mark.parametrize( + "temperature, fuel", + [ + (2e3, {"D": 1.0}), + (2e3, {"T": 1.0}), + (2e3, {"T": 0.5, "D": 0.5}), + (2e3, {"T": 0.2, "D": 0.8}), + ], +) def test_fuel_with_correct_inputs(temperature, fuel): # Should accept any non-zero positive inputs to either variable get_neutron_energy_distribution(temperature, fuel) -@pytest.mark.parametrize("temperature, fuel", [ - (2e3, {'D': 1.1}), - (2e3, {'T': 0.9}), - (2e3, {'T': -0.5, 'D': 0.5}), - (2e3, {'T': -0.2, 'D': -0.8}), -]) +@pytest.mark.parametrize( + "temperature, fuel", + [ + (2e3, {"D": 1.1}), + (2e3, {"T": 0.9}), + (2e3, {"T": -0.5, "D": 0.5}), + (2e3, {"T": -0.2, "D": -0.8}), + ], +) def test_fuel_with_bad_inputs(temperature, fuel): # Should reject any negative numbers and zeros. with pytest.raises(ValueError): get_neutron_energy_distribution(temperature, fuel) -@pytest.mark.parametrize("temperature, fuel", [ - (2e3, {'DD': 1.1}), - (2e3, {'DT': 0.9}), - (2e3, {'He3': -0.5, 'D': 0.5}), - (2e3, {1: -0.2, 'D': -0.8}), -]) +@pytest.mark.parametrize( + "temperature, fuel", + [ + (2e3, {"DD": 1.1}), + (2e3, {"DT": 0.9}), + (2e3, {"He3": -0.5, "D": 0.5}), + (2e3, {1: -0.2, "D": -0.8}), + ], +) def test_fuel_with_incorrect_isotopese(temperature, fuel): # Should reject anything which is not 'D' or 'T'. with pytest.raises(ValueError): From 01d326b1e88c107c0492f994f136272e50e0f39f Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 6 Mar 2024 13:54:59 +0000 Subject: [PATCH 17/44] [skip ci] ring source example working --- examples/ring_source_example.py | 4 +-- src/openmc_plasma_source/point_source.py | 2 +- src/openmc_plasma_source/ring_source.py | 34 +++++++++++++++--------- 3 files changed, 25 insertions(+), 15 deletions(-) diff --git a/examples/ring_source_example.py b/examples/ring_source_example.py index ee17cbf..18efd13 100644 --- a/examples/ring_source_example.py +++ b/examples/ring_source_example.py @@ -3,7 +3,7 @@ import openmc -from openmc_plasma_source import FusionRingSource +from openmc_plasma_source import fusion_ring_source # just making use of a local cross section xml file, replace with your own cross sections or comment out openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml" @@ -14,7 +14,7 @@ geometry = openmc.Geometry([cell]) # define the source -my_source = FusionRingSource( +my_source = fusion_ring_source( radius=700, angles=(0.0, 2 * math.pi), # 360deg source temperature=20000.0, diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 646f27f..eefa4ba 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -51,7 +51,7 @@ def fusion_point_source( # else: for energy_distribution, strength in zip(energy_distributions, strengths): - source = openmc.Source() + source = openmc.IndependentSource() source.energy = energy_distribution source.space = openmc.stats.Point(coordinate) source.angle = openmc.stats.Isotropic() diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index 996a633..07903b6 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -53,18 +53,28 @@ def fusion_ring_source( else: raise ValueError("Temperature must be a float strictly greater than 0.") - source = openmc.IndependentSource() - - # performed after the super init as these are Source attributes - source.space = openmc.stats.CylindricalIndependent( - r=openmc.stats.Discrete([radius], [1]), - phi=openmc.stats.Uniform(a=angles[0], b=angles[1]), - z=openmc.stats.Discrete([z_placement], [1]), - origin=(0.0, 0.0, 0.0), - ) - source.angle = openmc.stats.Isotropic() - source.energy = get_neutron_energy_distribution( + + sources = [] + + energy_distributions, strengths = get_neutron_energy_distribution( ion_temperature=temperature, fuel=fuel ) - return source + for energy_distribution, strength in zip(energy_distributions, strengths): + + source = openmc.IndependentSource() + + source.space = openmc.stats.CylindricalIndependent( + r=openmc.stats.Discrete([radius], [1]), + phi=openmc.stats.Uniform(a=angles[0], b=angles[1]), + z=openmc.stats.Discrete([z_placement], [1]), + origin=(0.0, 0.0, 0.0), + ) + + source.energy = energy_distribution + source.angle = openmc.stats.Isotropic() + source.strength = strength + sources.append(source) + + return sources + From 0a12cda2d66626ddf29679f2c6f691cad51f5f6d Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 6 Mar 2024 13:56:04 +0000 Subject: [PATCH 18/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/ring_source.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index 07903b6..780e77f 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -53,7 +53,6 @@ def fusion_ring_source( else: raise ValueError("Temperature must be a float strictly greater than 0.") - sources = [] energy_distributions, strengths = get_neutron_energy_distribution( @@ -77,4 +76,3 @@ def fusion_ring_source( sources.append(source) return sources - From 146310e8ba9faa07ffad6bbec5eca46dc6a38a6a Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 6 Mar 2024 20:25:38 +0000 Subject: [PATCH 19/44] [skip ci] progress toko source --- examples/tokamak_source_example.py | 4 +- src/openmc_plasma_source/fuel_types.py | 16 +- src/openmc_plasma_source/tokamak_source.py | 543 ++++++++++----------- 3 files changed, 263 insertions(+), 300 deletions(-) diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index af973a5..495958a 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -2,7 +2,7 @@ import openmc -from openmc_plasma_source import TokamakSource +from openmc_plasma_source import tokamak_source # just making use of a local cross section xml file, replace with your own cross sections or comment out openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml" @@ -14,7 +14,7 @@ geometry = openmc.Geometry([cell]) # create a plasma source -my_plasma = TokamakSource( +my_plasma = tokamak_source( elongation=1.557, ion_density_centre=1.09e20, ion_density_peaking_factor=1, diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 2b50426..22ec8ad 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -51,7 +51,8 @@ def get_neutron_energy_distribution( DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) # print("DTmean", DTmean) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) - # print("DDmean", DDmean) + print("DDmean", DDmean) + print("DDvar", DDvar) # todo make use of these DTvar and DDvar values in muir or gaussian distribution if ["D", "T"] == sorted(set(fuel.keys())): @@ -66,8 +67,17 @@ def get_neutron_energy_distribution( dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) - dd_source = openmc.stats.muir(e0=DDmean * 1e6, m_rat=4, kt=ion_temperature) - dt_source = openmc.stats.muir(e0=DTmean * 1e6, m_rat=5, kt=ion_temperature) + + DD_std_dev = np.sqrt(DDvar* 1e12) # power 12 as this is in MeV^2 and we need eV + dd_source = openmc.stats.Normal(mean_value=DDmean * 1e6, std_dev=DD_std_dev) + # normal could be done with Muir but in this case we have the mean and std dev from NeSST + # dd_source = openmc.stats.muir(e0=DDmean * 1e6, m_rat=4, kt=ion_temperature) + + DT_std_dev = np.sqrt(DTvar* 1e12) # power 12 as this is in MeV^2 and we need eV + dt_source = openmc.stats.Normal(mean_value=DTmean * 1e6, std_dev=DT_std_dev) + # normal could be done with Muir but in this case we have the mean and std dev from NeSST + # dt_source = openmc.stats.muir(e0=DTmean * 1e6, m_rat=5, kt=ion_temperature) + # todo look into combining distributions openmc.data.combine_distributions() return [tt_source, dd_source, dt_source], [ strength_TT, diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 05a26b5..21e3365 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -82,306 +82,259 @@ def tokamak_source( if abs(shafranov_factor) >= 0.5 * minor_radius: raise ValueError("Shafranov factor must be smaller than 0.5*minor radius") + if isinstance(major_radius, (int, float)) and major_radius > 0: + pass + else: + raise ValueError( + "Major radius must be a number within the specified bounds" + ) + + if isinstance(minor_radius, (int, float)) and minor_radius > 0: + self._minor_radius = minor_radius + else: + raise ValueError( + "Minor radius must be a number within the specified bounds" + ) + + if isinstance(elongation, (int, float)) and elongation > 0: + pass + else: + raise ValueError("Elongation must be a number within the specified bounds") + + if isinstance(triangularity, (int, float)) and -1.0 <= value <= 1.0: + pass + else: + raise ValueError( + "Triangularity must be a number within the specified bounds" + ) + + if mode not in ["H", "L", "A"]: + raise ValueError("Mode must be one of the following: ['H', 'L', 'A']") + + if isinstance(ion_density_centre, (int, float)) and value > 0: + pass + else: + raise ValueError("Ion density centre must be a number greater than 0") + + if isinstance(ion_density_peaking_factor, (int, float)): + self._ion_density_peaking_factor = value + else: + raise ValueError("Ion density peaking factor must be a number") + + if isinstance(ion_density_pedestal, (int, float)) and value > 0: + pass + else: + raise ValueError("Ion density pedestal must be a number greater than 0") + + if isinstance(ion_density_separatrix, (int, float)) and value > 0: + self._ion_density_separatrix = value + else: + raise ValueError("Ion density separatrix must be a number greater than 0") + + if ( + isinstance(angles, tuple) + and len(angles) == 2 + and all( + isinstance(angle, (int, float)) and -2 * np.pi <= angle <= 2 * np.pi + for angle in angles + ) + ): + pass + else: + raise ValueError( + "Angles must be a tuple of floats between zero and 2 * np.pi" + ) + # Create a list of souces - self.sample_sources() + sources = sample_sources() self.sources = self.make_openmc_sources() - @property - def major_radius(self): - return self._major_radius - - @major_radius.setter - def major_radius(self, value): - if isinstance(value, (int, float)) and value > 0: - self._major_radius = value - else: - raise ValueError( - "Major radius must be a number within the specified bounds" - ) - - @property - def minor_radius(self): - return self._minor_radius - - @minor_radius.setter - def minor_radius(self, value): - if isinstance(value, (int, float)) and value > 0: - self._minor_radius = value - else: - raise ValueError( - "Minor radius must be a number within the specified bounds" - ) - - @property - def elongation(self): - return self._elongation - - @elongation.setter - def elongation(self, value): - if isinstance(value, (int, float)) and value > 0: - self._elongation = value - else: - raise ValueError("Elongation must be a number within the specified bounds") - - @property - def triangularity(self): - return self._triangularity - - @triangularity.setter - def triangularity(self, value): - if isinstance(value, (int, float)) and -1.0 <= value <= 1.0: - self._triangularity = value - else: - raise ValueError( - "Triangularity must be a number within the specified bounds" - ) - - @property - def mode(self): - return self._mode - - @mode.setter - def mode(self, value): - if value in ["H", "L", "A"]: - self._mode = value - else: - raise ValueError("Mode must be one of the following: ['H', 'L', 'A']") - - @property - def ion_density_centre(self): - return self._ion_density_centre - - @ion_density_centre.setter - def ion_density_centre(self, value): - if isinstance(value, (int, float)) and value > 0: - self._ion_density_centre = value - else: - raise ValueError("Ion density centre must be a number greater than 0") - - @property - def ion_density_peaking_factor(self): - return self._ion_density_peaking_factor - - @ion_density_peaking_factor.setter - def ion_density_peaking_factor(self, value): - if isinstance(value, (int, float)): - self._ion_density_peaking_factor = value - else: - raise ValueError("Ion density peaking factor must be a number") - - @property - def ion_density_pedestal(self): - return self._ion_density_pedestal - - @ion_density_pedestal.setter - def ion_density_pedestal(self, value): - if isinstance(value, (int, float)) and value > 0: - self._ion_density_pedestal = value - else: - raise ValueError("Ion density pedestal must be a number greater than 0") - - @property - def ion_density_separatrix(self): - return self._ion_density_separatrix - - @ion_density_separatrix.setter - def ion_density_separatrix(self, value): - if isinstance(value, (int, float)) and value > 0: - self._ion_density_separatrix = value - else: - raise ValueError("Ion density separatrix must be a number greater than 0") - - @property - def angles(self): - return self._angles - - @angles.setter - def angles(self, value): - if ( - isinstance(value, tuple) - and len(value) == 2 - and all( - isinstance(angle, (int, float)) and -2 * np.pi <= angle <= 2 * np.pi - for angle in value - ) - ): - self._angles = value - else: - raise ValueError( - "Angles must be a tuple of floats between zero and 2 * np.pi" - ) - - # TODO setters and getters for the rest - - def _bounds_check(value, bounds): - return bounds[0] < value - - def ion_density(self, r): - """Computes the ion density at a given position. The ion density is - only dependent on the minor radius. - - Args: - r (float, ndarray): the minor radius (cm) - - Returns: - float, ndarray: ion density in m-3 - """ - - r = np.asarray(r) - if np.any(r < 0): - raise ValueError("Minor radius must not be negative") - - if self.mode == "L": - density = ( - self.ion_density_centre - * (1 - (r / self.major_radius) ** 2) ** self.ion_density_peaking_factor - ) - elif self.mode in ["H", "A"]: - density = np.where( - r < self.pedestal_radius, - ( - (self.ion_density_centre - self.ion_density_pedestal) - * (1 - (r / self.pedestal_radius) ** 2) - ** self.ion_density_peaking_factor - + self.ion_density_pedestal - ), - ( - (self.ion_density_pedestal - self.ion_density_separatrix) - * (self.major_radius - r) - / (self.major_radius - self.pedestal_radius) - + self.ion_density_separatrix - ), - ) - return density - - def ion_temperature(self, r): - """Computes the ion temperature at a given position. The ion - temperature is only dependent on the minor radius. - - Args: - r (float, ndarray): minor radius (cm) - - Returns: - float, ndarray: ion temperature (keV) - """ - - r = np.asarray(r) - if np.any(r < 0): - raise ValueError("Minor radius must not be negative") - - if self.mode == "L": - temperature = ( - self.ion_temperature_centre - * (1 - (r / self.major_radius) ** 2) + + +def _ion_density( + mode, + ion_density_centre, + ion_density_peaking_factor, + ion_density_pedestal, + major_radius, + pedestal_radius, + ion_density_separatrix, + r +): + """Computes the ion density at a given position. The ion density is + only dependent on the minor radius. + + Args: + r (float, ndarray): the minor radius (cm) + + Returns: + float, ndarray: ion density in m-3 + """ + + r = np.asarray(r) + if np.any(r < 0): + raise ValueError("Minor radius must not be negative") + + if mode == "L": + density = ( + ion_density_centre + * (1 - (r / major_radius) ** 2) ** ion_density_peaking_factor + ) + elif mode in ["H", "A"]: + density = np.where( + r < pedestal_radius, + ( + (ion_density_centre - ion_density_pedestal) + * (1 - (r / pedestal_radius) ** 2) + ** ion_density_peaking_factor + + ion_density_pedestal + ), + ( + (ion_density_pedestal - ion_density_separatrix) + * (major_radius - r) + / (major_radius - pedestal_radius) + + ion_density_separatrix + ), + ) + return density + +def _ion_temperature(self, r): + """Computes the ion temperature at a given position. The ion + temperature is only dependent on the minor radius. + + Args: + r (float, ndarray): minor radius (cm) + + Returns: + float, ndarray: ion temperature (keV) + """ + + r = np.asarray(r) + if np.any(r < 0): + raise ValueError("Minor radius must not be negative") + + if self.mode == "L": + temperature = ( + self.ion_temperature_centre + * (1 - (r / self.major_radius) ** 2) + ** self.ion_temperature_peaking_factor + ) + elif self.mode in ["H", "A"]: + temperature = np.where( + r < self.pedestal_radius, + ( + self.ion_temperature_pedestal + + (self.ion_temperature_centre - self.ion_temperature_pedestal) + * (1 - (r / self.pedestal_radius) ** self.ion_temperature_beta) ** self.ion_temperature_peaking_factor - ) - elif self.mode in ["H", "A"]: - temperature = np.where( - r < self.pedestal_radius, - ( - self.ion_temperature_pedestal - + (self.ion_temperature_centre - self.ion_temperature_pedestal) - * (1 - (r / self.pedestal_radius) ** self.ion_temperature_beta) - ** self.ion_temperature_peaking_factor - ), - ( - self.ion_temperature_separatrix - + (self.ion_temperature_pedestal - self.ion_temperature_separatrix) - * (self.major_radius - r) - / (self.major_radius - self.pedestal_radius) - ), - ) - return temperature - - def convert_a_alpha_to_R_Z(self, a, alpha): - """Converts (r, alpha) cylindrical coordinates to (R, Z) cartesian - coordinates. - - Args: - a (float, ndarray): minor radius (cm) - alpha (float, ndarray): angle (rad) - - Returns: - ((float, ndarray), (float, ndarray)): (R, Z) coordinates - """ - a = np.asarray(a) - alpha = np.asarray(alpha) - if np.any(a < 0): - raise ValueError("Radius 'a' must not be negative") - - shafranov_shift = self.shafranov_factor * (1.0 - (a / self.minor_radius) ** 2) - R = ( - self.major_radius - + a * np.cos(alpha + (self.triangularity * np.sin(alpha))) - + shafranov_shift + ), + ( + self.ion_temperature_separatrix + + (self.ion_temperature_pedestal - self.ion_temperature_separatrix) + * (self.major_radius - r) + / (self.major_radius - self.pedestal_radius) + ), + ) + return temperature + +def _convert_a_alpha_to_R_Z(self, a, alpha): + """Converts (r, alpha) cylindrical coordinates to (R, Z) cartesian + coordinates. + + Args: + a (float, ndarray): minor radius (cm) + alpha (float, ndarray): angle (rad) + + Returns: + ((float, ndarray), (float, ndarray)): (R, Z) coordinates + """ + a = np.asarray(a) + alpha = np.asarray(alpha) + if np.any(a < 0): + raise ValueError("Radius 'a' must not be negative") + + shafranov_shift = self.shafranov_factor * (1.0 - (a / self.minor_radius) ** 2) + R = ( + self.major_radius + + a * np.cos(alpha + (self.triangularity * np.sin(alpha))) + + shafranov_shift + ) + Z = self.elongation * a * np.sin(alpha) + return (R, Z) + +def _sample_sources(sample_size, minor_radius): + """Samples self.sample_size neutrons and creates attributes .densities + (ion density), .temperatures (ion temperature), .strengths + (neutron source density) and .RZ (coordinates) + """ + # create a sample of (a, alpha) coordinates + a = np.random.random(self.sample_size) * self.minor_radius + alpha = np.random.random(self.sample_size) * 2 * np.pi + + # compute densities, temperatures, neutron source densities and + # convert coordinates + self.densities = self.ion_density( + mode=mode, + ion_density_centre=ion_density_centre, + ion_density_peaking_factor=ion_density_peaking_factor, + ion_density_pedestal=ion_density_pedestal, + major_radius=major_radius, + pedestal_radius=pedestal_radius, + ion_density_separatrix=ion_density_separatrix, + r=a + ) + self.temperatures = self.ion_temperature(a) + self.neutron_source_density = neutron_source_density( + self.densities, self.temperatures + ) + self.strengths = self.neutron_source_density / sum(self.neutron_source_density) + self.RZ = self.convert_a_alpha_to_R_Z(a, alpha) + +def _make_openmc_sources(self): + """Creates a list of OpenMC Sources() objects. The created sources are + ring sources based on the .RZ coordinates between two angles. The + energy of the sources are Muir energy spectra with ion temperatures + based on .temperatures. The strength of the sources (their probability) + is based on .strengths. + + Args: + angles ((float, float), optional): rotation of the ring source. + Defaults to (0, 2*np.pi). + + Returns: + list: list of openmc.IndependentSource() + """ + + sources = [] + # create a ring source for each sample in the plasma source + for i in range(len(self.strengths)): + my_source = openmc.IndependentSource() + + # extract the RZ values accordingly + radius = openmc.stats.Discrete([self.RZ[0][i]], [1]) + z_values = openmc.stats.Discrete([self.RZ[1][i]], [1]) + angle = openmc.stats.Uniform(a=self.angles[0], b=self.angles[1]) + + # create a ring source + my_source.space = openmc.stats.CylindricalIndependent( + r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0) ) - Z = self.elongation * a * np.sin(alpha) - return (R, Z) - - def sample_sources(self): - """Samples self.sample_size neutrons and creates attributes .densities - (ion density), .temperatures (ion temperature), .strengths - (neutron source density) and .RZ (coordinates) - """ - # create a sample of (a, alpha) coordinates - a = np.random.random(self.sample_size) * self.minor_radius - alpha = np.random.random(self.sample_size) * 2 * np.pi - - # compute densities, temperatures, neutron source densities and - # convert coordinates - self.densities = self.ion_density(a) - self.temperatures = self.ion_temperature(a) - self.neutron_source_density = neutron_source_density( - self.densities, self.temperatures + + my_source.angle = openmc.stats.Isotropic() + my_source.energy = get_neutron_energy_distribution( + ion_temperature=self.temperature[i], fuel=self.fuel ) - self.strengths = self.neutron_source_density / sum(self.neutron_source_density) - self.RZ = self.convert_a_alpha_to_R_Z(a, alpha) - - def make_openmc_sources(self): - """Creates a list of OpenMC Sources() objects. The created sources are - ring sources based on the .RZ coordinates between two angles. The - energy of the sources are Muir energy spectra with ion temperatures - based on .temperatures. The strength of the sources (their probability) - is based on .strengths. - - Args: - angles ((float, float), optional): rotation of the ring source. - Defaults to (0, 2*np.pi). - - Returns: - list: list of openmc.IndependentSource() - """ - - sources = [] - # create a ring source for each sample in the plasma source - for i in range(len(self.strengths)): - my_source = openmc.IndependentSource() - - # extract the RZ values accordingly - radius = openmc.stats.Discrete([self.RZ[0][i]], [1]) - z_values = openmc.stats.Discrete([self.RZ[1][i]], [1]) - angle = openmc.stats.Uniform(a=self.angles[0], b=self.angles[1]) - - # create a ring source - my_source.space = openmc.stats.CylindricalIndependent( - r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0) - ) - - my_source.angle = openmc.stats.Isotropic() - my_source.energy = get_neutron_energy_distribution( - ion_temperature=self.temperature[i], fuel=self.fuel - ) - - # the strength of the source (its probability) is given by - # self.strengths - my_source.strength = self.strengths[i] - - # append to the list of sources - sources.append(my_source) - return sources - - -def neutron_source_density(ion_density, ion_temperature): + + # the strength of the source (its probability) is given by + # self.strengths + my_source.strength = self.strengths[i] + + # append to the list of sources + sources.append(my_source) + return sources + + +def _neutron_source_density(ion_density, ion_temperature): """Computes the neutron source density given ion density and ion temperature. @@ -399,7 +352,7 @@ def neutron_source_density(ion_density, ion_temperature): return ion_density**2 * DT_xs(ion_temperature) -def DT_xs(ion_temperature): +def _DT_xs(ion_temperature): """Sadler–Van Belle formula Ref : https://doi.org/10.1016/j.fusengdes.2012.02.025 From 05f984831bca6753a528ff15badb9bcffbb44d00 Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 6 Mar 2024 20:26:13 +0000 Subject: [PATCH 20/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/fuel_types.py | 10 ++++--- src/openmc_plasma_source/tokamak_source.py | 31 +++++++++------------- 2 files changed, 19 insertions(+), 22 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 22ec8ad..defc8fc 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -68,16 +68,20 @@ def get_neutron_energy_distribution( dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) - DD_std_dev = np.sqrt(DDvar* 1e12) # power 12 as this is in MeV^2 and we need eV + DD_std_dev = np.sqrt( + DDvar * 1e12 + ) # power 12 as this is in MeV^2 and we need eV dd_source = openmc.stats.Normal(mean_value=DDmean * 1e6, std_dev=DD_std_dev) # normal could be done with Muir but in this case we have the mean and std dev from NeSST # dd_source = openmc.stats.muir(e0=DDmean * 1e6, m_rat=4, kt=ion_temperature) - DT_std_dev = np.sqrt(DTvar* 1e12) # power 12 as this is in MeV^2 and we need eV + DT_std_dev = np.sqrt( + DTvar * 1e12 + ) # power 12 as this is in MeV^2 and we need eV dt_source = openmc.stats.Normal(mean_value=DTmean * 1e6, std_dev=DT_std_dev) # normal could be done with Muir but in this case we have the mean and std dev from NeSST # dt_source = openmc.stats.muir(e0=DTmean * 1e6, m_rat=5, kt=ion_temperature) - + # todo look into combining distributions openmc.data.combine_distributions() return [tt_source, dd_source, dt_source], [ strength_TT, diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 21e3365..96e3ea8 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -85,16 +85,12 @@ def tokamak_source( if isinstance(major_radius, (int, float)) and major_radius > 0: pass else: - raise ValueError( - "Major radius must be a number within the specified bounds" - ) + raise ValueError("Major radius must be a number within the specified bounds") if isinstance(minor_radius, (int, float)) and minor_radius > 0: self._minor_radius = minor_radius else: - raise ValueError( - "Minor radius must be a number within the specified bounds" - ) + raise ValueError("Minor radius must be a number within the specified bounds") if isinstance(elongation, (int, float)) and elongation > 0: pass @@ -104,9 +100,7 @@ def tokamak_source( if isinstance(triangularity, (int, float)) and -1.0 <= value <= 1.0: pass else: - raise ValueError( - "Triangularity must be a number within the specified bounds" - ) + raise ValueError("Triangularity must be a number within the specified bounds") if mode not in ["H", "L", "A"]: raise ValueError("Mode must be one of the following: ['H', 'L', 'A']") @@ -141,16 +135,13 @@ def tokamak_source( ): pass else: - raise ValueError( - "Angles must be a tuple of floats between zero and 2 * np.pi" - ) + raise ValueError("Angles must be a tuple of floats between zero and 2 * np.pi") # Create a list of souces sources = sample_sources() self.sources = self.make_openmc_sources() - def _ion_density( mode, ion_density_centre, @@ -159,7 +150,7 @@ def _ion_density( major_radius, pedestal_radius, ion_density_separatrix, - r + r, ): """Computes the ion density at a given position. The ion density is only dependent on the minor radius. @@ -185,8 +176,7 @@ def _ion_density( r < pedestal_radius, ( (ion_density_centre - ion_density_pedestal) - * (1 - (r / pedestal_radius) ** 2) - ** ion_density_peaking_factor + * (1 - (r / pedestal_radius) ** 2) ** ion_density_peaking_factor + ion_density_pedestal ), ( @@ -198,6 +188,7 @@ def _ion_density( ) return density + def _ion_temperature(self, r): """Computes the ion temperature at a given position. The ion temperature is only dependent on the minor radius. @@ -216,8 +207,7 @@ def _ion_temperature(self, r): if self.mode == "L": temperature = ( self.ion_temperature_centre - * (1 - (r / self.major_radius) ** 2) - ** self.ion_temperature_peaking_factor + * (1 - (r / self.major_radius) ** 2) ** self.ion_temperature_peaking_factor ) elif self.mode in ["H", "A"]: temperature = np.where( @@ -237,6 +227,7 @@ def _ion_temperature(self, r): ) return temperature + def _convert_a_alpha_to_R_Z(self, a, alpha): """Converts (r, alpha) cylindrical coordinates to (R, Z) cartesian coordinates. @@ -262,6 +253,7 @@ def _convert_a_alpha_to_R_Z(self, a, alpha): Z = self.elongation * a * np.sin(alpha) return (R, Z) + def _sample_sources(sample_size, minor_radius): """Samples self.sample_size neutrons and creates attributes .densities (ion density), .temperatures (ion temperature), .strengths @@ -281,7 +273,7 @@ def _sample_sources(sample_size, minor_radius): major_radius=major_radius, pedestal_radius=pedestal_radius, ion_density_separatrix=ion_density_separatrix, - r=a + r=a, ) self.temperatures = self.ion_temperature(a) self.neutron_source_density = neutron_source_density( @@ -290,6 +282,7 @@ def _sample_sources(sample_size, minor_radius): self.strengths = self.neutron_source_density / sum(self.neutron_source_density) self.RZ = self.convert_a_alpha_to_R_Z(a, alpha) + def _make_openmc_sources(self): """Creates a list of OpenMC Sources() objects. The created sources are ring sources based on the .RZ coordinates between two angles. The From da838c810d88a49d096fa16384c11bb4b5f44ae9 Mon Sep 17 00:00:00 2001 From: shimwell Date: Thu, 7 Mar 2024 14:03:02 +0000 Subject: [PATCH 21/44] tokamak source example working --- examples/tokamak_source_example.py | 4 +- src/openmc_plasma_source/fuel_types.py | 14 +- src/openmc_plasma_source/tokamak_source.py | 265 ++++++++++++++------- 3 files changed, 181 insertions(+), 102 deletions(-) diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index 495958a..d8c9666 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -14,7 +14,7 @@ geometry = openmc.Geometry([cell]) # create a plasma source -my_plasma = tokamak_source( +my_sources = tokamak_source( elongation=1.557, ion_density_centre=1.09e20, ion_density_peaking_factor=1, @@ -39,7 +39,7 @@ settings.run_mode = "fixed source" settings.batches = 10 settings.particles = 1000 -settings.source = my_plasma.sources +settings.source = my_sources model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 22ec8ad..6f62201 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -7,7 +7,7 @@ def get_neutron_energy_distribution( ion_temperature: float, fuel: dict, ) -> openmc.stats.Discrete: - """Finds the energy distribution + """Finds the energy distribution and their relative strengths. Parameters ---------- @@ -49,11 +49,7 @@ def get_neutron_energy_distribution( E_pspec = np.linspace(0, 12, num_of_vals) # E_pspec is exspected in MeV units DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) - # print("DTmean", DTmean) DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) - print("DDmean", DDmean) - print("DDvar", DDvar) - # todo make use of these DTvar and DDvar values in muir or gaussian distribution if ["D", "T"] == sorted(set(fuel.keys())): @@ -65,6 +61,8 @@ def get_neutron_energy_distribution( "tt", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"] ) + total_strength = sum([strength_DT, strength_DD, strength_TT]) + dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) @@ -80,9 +78,9 @@ def get_neutron_energy_distribution( # todo look into combining distributions openmc.data.combine_distributions() return [tt_source, dd_source, dt_source], [ - strength_TT, - strength_DD, - strength_DT, + strength_TT/total_strength, + strength_DD/total_strength, + strength_DT/total_strength, ] elif ["D"] == sorted(set(fuel.keys())): diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 21e3365..c5d3443 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -82,54 +82,49 @@ def tokamak_source( if abs(shafranov_factor) >= 0.5 * minor_radius: raise ValueError("Shafranov factor must be smaller than 0.5*minor radius") - if isinstance(major_radius, (int, float)) and major_radius > 0: - pass - else: - raise ValueError( - "Major radius must be a number within the specified bounds" - ) + if not isinstance(major_radius, (int, float)): + raise ValueError("Major radius must be a number") - if isinstance(minor_radius, (int, float)) and minor_radius > 0: - self._minor_radius = minor_radius - else: - raise ValueError( - "Minor radius must be a number within the specified bounds" - ) + if major_radius < 0: + raise ValueError("Major radius must greater than 0") - if isinstance(elongation, (int, float)) and elongation > 0: - pass - else: - raise ValueError("Elongation must be a number within the specified bounds") + if not isinstance(minor_radius, (int, float)): + raise ValueError("Minor radius must be a number") + if minor_radius < 0: + raise ValueError("Minor radius must greater than 0") - if isinstance(triangularity, (int, float)) and -1.0 <= value <= 1.0: - pass - else: - raise ValueError( - "Triangularity must be a number within the specified bounds" - ) + if not isinstance(elongation, (int, float)): + raise ValueError("Elongation must be a number") + if elongation < 0: + raise ValueError("Elongation must greater than 0") + + if not isinstance(triangularity, (int, float)): + raise ValueError("Triangularity must be a number") + if not triangularity <= 1.0: + raise ValueError("Triangularity must less than or equal to 1.") + if not triangularity >= -1.0: + raise ValueError("Triangularity must greater than or equal to -1.") if mode not in ["H", "L", "A"]: raise ValueError("Mode must be one of the following: ['H', 'L', 'A']") - if isinstance(ion_density_centre, (int, float)) and value > 0: - pass - else: - raise ValueError("Ion density centre must be a number greater than 0") + if not isinstance(ion_density_centre, (int, float)): + raise ValueError("ion_density_centre must be a number") + if ion_density_centre < 0: + raise ValueError("ion_density_centre must greater than 0") - if isinstance(ion_density_peaking_factor, (int, float)): - self._ion_density_peaking_factor = value - else: - raise ValueError("Ion density peaking factor must be a number") + if not isinstance(ion_density_peaking_factor, (int, float)): + raise ValueError("ion_density_peaking_factor must be a number") - if isinstance(ion_density_pedestal, (int, float)) and value > 0: - pass - else: - raise ValueError("Ion density pedestal must be a number greater than 0") + if not isinstance(ion_density_pedestal, (int, float)): + raise ValueError("ion_density_pedestal must be a number") + if ion_density_pedestal < 0: + raise ValueError("ion_density_pedestal must greater than 0") - if isinstance(ion_density_separatrix, (int, float)) and value > 0: - self._ion_density_separatrix = value - else: - raise ValueError("Ion density separatrix must be a number greater than 0") + if not isinstance(ion_density_separatrix, (int, float)): + raise ValueError("ion_density_separatrix must be a number") + if ion_density_separatrix < 0: + raise ValueError("ion_density_separatrix must greater than 0") if ( isinstance(angles, tuple) @@ -145,9 +140,61 @@ def tokamak_source( "Angles must be a tuple of floats between zero and 2 * np.pi" ) - # Create a list of souces - sources = sample_sources() - self.sources = self.make_openmc_sources() + # Create a list of sources + """Samples sample_size neutrons and creates attributes .densities + (ion density), .temperatures (ion temperature), .strengths + (neutron source density) and .RZ (coordinates) + """ + # create a sample of (a, alpha) coordinates + 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 + densities = _ion_density( + mode=mode, + ion_density_centre=ion_density_centre, + ion_density_peaking_factor=ion_density_peaking_factor, + ion_density_pedestal=ion_density_pedestal, + major_radius=major_radius, + pedestal_radius=pedestal_radius, + ion_density_separatrix=ion_density_separatrix, + r=a + ) + temperatures = _ion_temperature( + r=a, + mode=mode, + pedestal_radius=pedestal_radius, + ion_temperature_pedestal=ion_temperature_pedestal, + ion_temperature_centre=ion_temperature_centre, + ion_temperature_beta=ion_temperature_beta, + ion_temperature_peaking_factor=ion_temperature_peaking_factor, + ion_temperature_separatrix=ion_temperature_separatrix, + major_radius=major_radius, + ) + + neutron_source_density = _neutron_source_density(densities, temperatures) + + strengths = neutron_source_density / sum(neutron_source_density) + + RZ = _convert_a_alpha_to_R_Z( + a=a, + alpha=alpha, + shafranov_factor=shafranov_factor, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=triangularity, + elongation=elongation, + ) + + sources = _make_openmc_sources( + strengths=strengths, + angles=angles, + temperatures=temperatures, + fuel=fuel, + RZ=RZ, + ) + return sources @@ -198,7 +245,17 @@ def _ion_density( ) return density -def _ion_temperature(self, r): +def _ion_temperature( + r, + mode, + pedestal_radius, + ion_temperature_pedestal, + ion_temperature_centre, + ion_temperature_beta, + ion_temperature_peaking_factor, + ion_temperature_separatrix, + major_radius, + ): """Computes the ion temperature at a given position. The ion temperature is only dependent on the minor radius. @@ -213,37 +270,48 @@ def _ion_temperature(self, r): if np.any(r < 0): raise ValueError("Minor radius must not be negative") - if self.mode == "L": + if mode == "L": temperature = ( - self.ion_temperature_centre - * (1 - (r / self.major_radius) ** 2) - ** self.ion_temperature_peaking_factor + ion_temperature_centre + * (1 - (r / major_radius) ** 2) + ** ion_temperature_peaking_factor ) - elif self.mode in ["H", "A"]: + elif mode in ["H", "A"]: temperature = np.where( - r < self.pedestal_radius, + r < pedestal_radius, ( - self.ion_temperature_pedestal - + (self.ion_temperature_centre - self.ion_temperature_pedestal) - * (1 - (r / self.pedestal_radius) ** self.ion_temperature_beta) - ** self.ion_temperature_peaking_factor + ion_temperature_pedestal + + (ion_temperature_centre - ion_temperature_pedestal) + * (1 - (r / pedestal_radius) ** ion_temperature_beta) + ** ion_temperature_peaking_factor ), ( - self.ion_temperature_separatrix - + (self.ion_temperature_pedestal - self.ion_temperature_separatrix) - * (self.major_radius - r) - / (self.major_radius - self.pedestal_radius) + ion_temperature_separatrix + + (ion_temperature_pedestal - ion_temperature_separatrix) + * (major_radius - r) + / (major_radius - pedestal_radius) ), ) return temperature -def _convert_a_alpha_to_R_Z(self, a, alpha): +def _convert_a_alpha_to_R_Z( + a, + alpha, + shafranov_factor, + minor_radius, + major_radius, + triangularity, + elongation, +): """Converts (r, alpha) cylindrical coordinates to (R, Z) cartesian coordinates. Args: a (float, ndarray): minor radius (cm) alpha (float, ndarray): angle (rad) + shafranov_factor: + minor_radius: + major_radius: Returns: ((float, ndarray), (float, ndarray)): (R, Z) coordinates @@ -253,27 +321,27 @@ def _convert_a_alpha_to_R_Z(self, a, alpha): if np.any(a < 0): raise ValueError("Radius 'a' must not be negative") - shafranov_shift = self.shafranov_factor * (1.0 - (a / self.minor_radius) ** 2) + shafranov_shift = shafranov_factor * (1.0 - (a / minor_radius) ** 2) R = ( - self.major_radius - + a * np.cos(alpha + (self.triangularity * np.sin(alpha))) + major_radius + + a * np.cos(alpha + (triangularity * np.sin(alpha))) + shafranov_shift ) - Z = self.elongation * a * np.sin(alpha) + Z = elongation * a * np.sin(alpha) return (R, Z) def _sample_sources(sample_size, minor_radius): - """Samples self.sample_size neutrons and creates attributes .densities + """Samples sample_size neutrons and creates attributes .densities (ion density), .temperatures (ion temperature), .strengths (neutron source density) and .RZ (coordinates) """ # create a sample of (a, alpha) coordinates - a = np.random.random(self.sample_size) * self.minor_radius - alpha = np.random.random(self.sample_size) * 2 * np.pi + 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 - self.densities = self.ion_density( + densities = _ion_density( mode=mode, ion_density_centre=ion_density_centre, ion_density_peaking_factor=ion_density_peaking_factor, @@ -283,14 +351,20 @@ def _sample_sources(sample_size, minor_radius): ion_density_separatrix=ion_density_separatrix, r=a ) - self.temperatures = self.ion_temperature(a) - self.neutron_source_density = neutron_source_density( - self.densities, self.temperatures + temperatures = ion_temperature(a) + neutron_source_density = neutron_source_density( + densities, temperatures ) - self.strengths = self.neutron_source_density / sum(self.neutron_source_density) - self.RZ = self.convert_a_alpha_to_R_Z(a, alpha) - -def _make_openmc_sources(self): + strengths = neutron_source_density / sum(neutron_source_density) + RZ = _convert_a_alpha_to_R_Z(a, alpha) + +def _make_openmc_sources( + strengths, + angles, + temperatures, + fuel, + RZ, +): """Creates a list of OpenMC Sources() objects. The created sources are ring sources based on the .RZ coordinates between two angles. The energy of the sources are Muir energy spectra with ion temperatures @@ -307,30 +381,37 @@ def _make_openmc_sources(self): sources = [] # create a ring source for each sample in the plasma source - for i in range(len(self.strengths)): - my_source = openmc.IndependentSource() - + for i in range(len(strengths)): # extract the RZ values accordingly - radius = openmc.stats.Discrete([self.RZ[0][i]], [1]) - z_values = openmc.stats.Discrete([self.RZ[1][i]], [1]) - angle = openmc.stats.Uniform(a=self.angles[0], b=self.angles[1]) - - # create a ring source - my_source.space = openmc.stats.CylindricalIndependent( - r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0) + radius = openmc.stats.Discrete([RZ[0][i]], [1]) + z_values = openmc.stats.Discrete([RZ[1][i]], [1]) + angle = openmc.stats.Uniform(a=angles[0], b=angles[1]) + strength = strengths[i] + + energy_distributions, dist_strengths = get_neutron_energy_distribution( + ion_temperature=temperatures[i], + fuel=fuel, ) - my_source.angle = openmc.stats.Isotropic() - my_source.energy = get_neutron_energy_distribution( - ion_temperature=self.temperature[i], fuel=self.fuel - ) + # now we have potentially 3 distributions (DT, DD, TT) + for energy_distribution, dist_strength in zip(energy_distributions, dist_strengths): + + my_source = openmc.IndependentSource() + + # create a ring source + my_source.space = openmc.stats.CylindricalIndependent( + r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0) + ) + my_source.angle = openmc.stats.Isotropic() + + my_source.energy = energy_distribution - # the strength of the source (its probability) is given by - # self.strengths - my_source.strength = self.strengths[i] + # the strength of the source (its probability) is given by the + # strength of the energy distribution and the location distribution + my_source.strength = dist_strength * strength - # append to the list of sources - sources.append(my_source) + # append to the list of sources + sources.append(my_source) return sources @@ -349,7 +430,7 @@ def _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) + return ion_density**2 * _DT_xs(ion_temperature) def _DT_xs(ion_temperature): From 7e44d24e10ac69c412ff2216b2e938a687824e4b Mon Sep 17 00:00:00 2001 From: shimwell Date: Thu, 7 Mar 2024 14:29:25 +0000 Subject: [PATCH 22/44] fuel testing working --- src/openmc_plasma_source/fuel_types.py | 5 ++++- tests/test_point_source.py | 12 ++++-------- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 6f62201..7565d43 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -31,7 +31,10 @@ def get_neutron_energy_distribution( ) if sum_fuel_isotopes < 0.0: - raise ValueError(f"isotope must sum to be above 0. Not {sum_fuel_isotopes}") + raise ValueError(f"isotope fractions must sum to be above 0. Not {sum_fuel_isotopes}") + + if sum_fuel_isotopes != 1.: + raise ValueError(f"isotope fractions must sum to 1. Not {sum_fuel_isotopes}") for k, v in fuel.items(): if k not in ["D", "T"]: diff --git a/tests/test_point_source.py b/tests/test_point_source.py index 9533f66..d250b79 100644 --- a/tests/test_point_source.py +++ b/tests/test_point_source.py @@ -22,8 +22,7 @@ def test_creation(): ) def test_coordinate(coordinate): # Should allow any tuple of length 3 containing numbers - my_source = fusion_point_source(coordinate=coordinate) - assert np.array_equal(my_source.coordinate, coordinate) + fusion_point_source(coordinate=coordinate) @pytest.mark.parametrize( @@ -39,9 +38,7 @@ def test_bad_coordinate(coordinate): @pytest.mark.parametrize("temperature", [20000, 1e4, 0.1, 25000.0]) def test_temperature(temperature): # Should accept any positive float - my_source = fusion_point_source(temperature=temperature) - assert my_source.temperature == temperature - + fusion_point_source(temperature=temperature) @pytest.mark.parametrize("temperature", [-20000.0, "hello world", [10000]]) def test_bad_temperature(temperature): @@ -50,11 +47,10 @@ def test_bad_temperature(temperature): fusion_point_source(temperature=temperature) -@pytest.mark.parametrize("fuel", ["DT", "DD"]) +@pytest.mark.parametrize("fuel", [{"D": 0.5, "T": 0.5}, {"D": 1.}, {'T': 1.}]) def test_fuel(fuel): # Should accept either 'DD' or 'DT' - my_source = fusion_point_source(fuel=fuel) - assert my_source.fuel_type == fuel + fusion_point_source(fuel=fuel) @pytest.mark.parametrize("fuel", ["топливо", 5]) From 68efbe84752cfb7301b078947b69f11f442ffec9 Mon Sep 17 00:00:00 2001 From: shimwell Date: Thu, 7 Mar 2024 14:52:19 +0000 Subject: [PATCH 23/44] point source testing working --- tests/test_point_source.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/tests/test_point_source.py b/tests/test_point_source.py index d250b79..10702f9 100644 --- a/tests/test_point_source.py +++ b/tests/test_point_source.py @@ -9,12 +9,13 @@ def test_creation(): my_source = fusion_point_source() # Ensure it is of type openmc.IndependentSource - assert isinstance(my_source, openmc.IndependentSource) + for source in my_source: + assert isinstance(source, openmc.IndependentSource) - # Ensure it has space, angle, and energy set - assert isinstance(my_source.space, openmc.stats.Point) - assert isinstance(my_source.angle, openmc.stats.Isotropic) - assert isinstance(my_source.energy, openmc.stats.univariate.Normal) + # Ensure it has space, angle, and energy set + assert isinstance(source.space, openmc.stats.Point) + assert isinstance(source.angle, openmc.stats.Isotropic) + assert isinstance(source.energy, openmc.stats.univariate.Normal) or isinstance(source.energy, openmc.stats.univariate.Tabular) @pytest.mark.parametrize( @@ -53,8 +54,8 @@ def test_fuel(fuel): fusion_point_source(fuel=fuel) -@pytest.mark.parametrize("fuel", ["топливо", 5]) +@pytest.mark.parametrize("fuel", [{"топливо": 1.}]) def test_wrong_fuel(fuel): # Should reject fuel types besides those listed in fuel_types.py - with pytest.raises((KeyError, TypeError)): + with pytest.raises(ValueError): fusion_point_source(fuel=fuel) From cb18a5f85ee3aa59ef7c3d41117e6dc3ef7b9420 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 09:59:55 +0000 Subject: [PATCH 24/44] added plotting example --- .gitignore | 3 +- examples/plot_tokamak_ion_temperature.py | 42 +++++++++++++++++++ src/openmc_plasma_source/__init__.py | 2 +- src/openmc_plasma_source/tokamak_source.py | 47 +++++----------------- 4 files changed, 55 insertions(+), 39 deletions(-) create mode 100644 examples/plot_tokamak_ion_temperature.py diff --git a/.gitignore b/.gitignore index e16264b..4a64908 100644 --- a/.gitignore +++ b/.gitignore @@ -135,8 +135,9 @@ dmypy.json # vim swap files *.swp -# image files created by tests +# image files created by tests and examples tests/*.png +*.png # openmc output files *.h5 diff --git a/examples/plot_tokamak_ion_temperature.py b/examples/plot_tokamak_ion_temperature.py new file mode 100644 index 0000000..031d2dc --- /dev/null +++ b/examples/plot_tokamak_ion_temperature.py @@ -0,0 +1,42 @@ +import matplotlib.pyplot as plt +import numpy as np +from openmc_plasma_source import tokamak_convert_a_alpha_to_R_Z, tokamak_ion_temperature + +sample_size=1000 +minor_radius=292.258 +major_radius=906 + +# create a sample of (a, alpha) coordinates +a = np.random.random(sample_size) * minor_radius +alpha = np.random.random(sample_size) * 2 * np.pi + +temperatures = tokamak_ion_temperature( + r=a, + mode='L', + pedestal_radius=0.8 * minor_radius, + ion_temperature_pedestal=6.09, + ion_temperature_centre=45.9, + ion_temperature_beta=2, + ion_temperature_peaking_factor=8.06, + ion_temperature_separatrix=0.1, + major_radius=major_radius, + ) + +RZ = tokamak_convert_a_alpha_to_R_Z( + a=a, + alpha=alpha, + shafranov_factor=0.44789, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=0.270, + elongation=1.557, +) + +plt.scatter(RZ[0], RZ[1], c=temperatures) +plt.gca().set_aspect('equal') +plt.xlabel("R [cm]") +plt.ylabel("Z [cm]") +plt.colorbar(label='Ion temperature [eV]') + +plt.savefig('tokamak_source_ion_temperature.png') +print('written tokamak_source_ion_temperature.png') \ No newline at end of file diff --git a/src/openmc_plasma_source/__init__.py b/src/openmc_plasma_source/__init__.py index 8a1ddda..0dc1642 100644 --- a/src/openmc_plasma_source/__init__.py +++ b/src/openmc_plasma_source/__init__.py @@ -14,4 +14,4 @@ from .fuel_types import get_neutron_energy_distribution from .point_source import fusion_point_source from .ring_source import fusion_ring_source -from .tokamak_source import tokamak_source +from .tokamak_source import * diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 4e09956..8bca1fb 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -149,7 +149,7 @@ def tokamak_source( # compute densities, temperatures, neutron source densities and # convert coordinates - densities = _ion_density( + densities = tokamak_ion_density( mode=mode, ion_density_centre=ion_density_centre, ion_density_peaking_factor=ion_density_peaking_factor, @@ -159,7 +159,7 @@ def tokamak_source( ion_density_separatrix=ion_density_separatrix, r=a, ) - temperatures = _ion_temperature( + temperatures = tokamak_ion_temperature( r=a, mode=mode, pedestal_radius=pedestal_radius, @@ -171,11 +171,11 @@ def tokamak_source( major_radius=major_radius, ) - neutron_source_density = _neutron_source_density(densities, temperatures) + neutron_source_density = tokamak_neutron_source_density(densities, temperatures) strengths = neutron_source_density / sum(neutron_source_density) - RZ = _convert_a_alpha_to_R_Z( + RZ = tokamak_convert_a_alpha_to_R_Z( a=a, alpha=alpha, shafranov_factor=shafranov_factor, @@ -185,7 +185,7 @@ def tokamak_source( elongation=elongation, ) - sources = _make_openmc_sources( + sources = tokamak_make_openmc_sources( strengths=strengths, angles=angles, temperatures=temperatures, @@ -195,7 +195,7 @@ def tokamak_source( return sources -def _ion_density( +def tokamak_ion_density( mode, ion_density_centre, ion_density_peaking_factor, @@ -242,7 +242,7 @@ def _ion_density( return density -def _ion_temperature( +def tokamak_ion_temperature( r, mode, pedestal_radius, @@ -291,7 +291,7 @@ def _ion_temperature( return temperature -def _convert_a_alpha_to_R_Z( +def tokamak_convert_a_alpha_to_R_Z( a, alpha, shafranov_factor, @@ -328,34 +328,7 @@ def _convert_a_alpha_to_R_Z( return (R, Z) -def _sample_sources(sample_size, minor_radius): - """Samples sample_size neutrons and creates attributes .densities - (ion density), .temperatures (ion temperature), .strengths - (neutron source density) and .RZ (coordinates) - """ - # create a sample of (a, alpha) coordinates - 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 - densities = _ion_density( - mode=mode, - ion_density_centre=ion_density_centre, - ion_density_peaking_factor=ion_density_peaking_factor, - ion_density_pedestal=ion_density_pedestal, - major_radius=major_radius, - pedestal_radius=pedestal_radius, - ion_density_separatrix=ion_density_separatrix, - r=a, - ) - temperatures = ion_temperature(a) - neutron_source_density = neutron_source_density(densities, temperatures) - strengths = neutron_source_density / sum(neutron_source_density) - RZ = _convert_a_alpha_to_R_Z(a, alpha) - - -def _make_openmc_sources( +def tokamak_make_openmc_sources( strengths, angles, temperatures, @@ -413,7 +386,7 @@ def _make_openmc_sources( return sources -def _neutron_source_density(ion_density, ion_temperature): +def tokamak_neutron_source_density(ion_density, ion_temperature): """Computes the neutron source density given ion density and ion temperature. From 88427dcb9c4f8d7a4aecbcca2ecbaa29840e8d9c Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 10:00:13 +0000 Subject: [PATCH 25/44] [skip ci] Apply formatting changes --- examples/plot_tokamak_ion_temperature.py | 34 ++++++++++++------------ tests/test_ring_source.py | 7 +++-- 2 files changed, 22 insertions(+), 19 deletions(-) diff --git a/examples/plot_tokamak_ion_temperature.py b/examples/plot_tokamak_ion_temperature.py index 031d2dc..ade5d8d 100644 --- a/examples/plot_tokamak_ion_temperature.py +++ b/examples/plot_tokamak_ion_temperature.py @@ -2,25 +2,25 @@ import numpy as np from openmc_plasma_source import tokamak_convert_a_alpha_to_R_Z, tokamak_ion_temperature -sample_size=1000 -minor_radius=292.258 -major_radius=906 +sample_size = 1000 +minor_radius = 292.258 +major_radius = 906 # create a sample of (a, alpha) coordinates a = np.random.random(sample_size) * minor_radius alpha = np.random.random(sample_size) * 2 * np.pi temperatures = tokamak_ion_temperature( - r=a, - mode='L', - pedestal_radius=0.8 * minor_radius, - ion_temperature_pedestal=6.09, - ion_temperature_centre=45.9, - ion_temperature_beta=2, - ion_temperature_peaking_factor=8.06, - ion_temperature_separatrix=0.1, - major_radius=major_radius, - ) + r=a, + mode="L", + pedestal_radius=0.8 * minor_radius, + ion_temperature_pedestal=6.09, + ion_temperature_centre=45.9, + ion_temperature_beta=2, + ion_temperature_peaking_factor=8.06, + ion_temperature_separatrix=0.1, + major_radius=major_radius, +) RZ = tokamak_convert_a_alpha_to_R_Z( a=a, @@ -33,10 +33,10 @@ ) plt.scatter(RZ[0], RZ[1], c=temperatures) -plt.gca().set_aspect('equal') +plt.gca().set_aspect("equal") plt.xlabel("R [cm]") plt.ylabel("Z [cm]") -plt.colorbar(label='Ion temperature [eV]') +plt.colorbar(label="Ion temperature [eV]") -plt.savefig('tokamak_source_ion_temperature.png') -print('written tokamak_source_ion_temperature.png') \ No newline at end of file +plt.savefig("tokamak_source_ion_temperature.png") +print("written tokamak_source_ion_temperature.png") diff --git a/tests/test_ring_source.py b/tests/test_ring_source.py index 00dae82..e7a0008 100644 --- a/tests/test_ring_source.py +++ b/tests/test_ring_source.py @@ -25,6 +25,7 @@ def test_radius(radius): # should allow any positive float fusion_ring_source(radius=radius) + @pytest.mark.parametrize("radius", [-1.0, "hello world", [1e5]]) def test_bad_radius(radius): # should reject any negative float or anything not convertible to float @@ -51,6 +52,7 @@ def test_temperature(temperature): # Should accept any positive float fusion_ring_source(radius=1.0, temperature=temperature) + @pytest.mark.parametrize("temperature", [-20000.0, "hello world", [10000]]) def test_bad_temperature(temperature): # Should reject negative floats and anything that isn't convertible to float @@ -58,12 +60,13 @@ def test_bad_temperature(temperature): fusion_ring_source(radius=1.0, temperature=temperature) -@pytest.mark.parametrize("fuel", [{'D':0.5,'T':0.5}, {"D":1.}]) +@pytest.mark.parametrize("fuel", [{"D": 0.5, "T": 0.5}, {"D": 1.0}]) def test_fuel(fuel): # Should accept either 'DD' or 'DT' fusion_ring_source(radius=1.0, fuel=fuel) -@pytest.mark.parametrize("fuel", [{"топливо": 1.}]) + +@pytest.mark.parametrize("fuel", [{"топливо": 1.0}]) def test_wrong_fuel(fuel): # Should reject fuel types besides those listed in fuel_types.py with pytest.raises(ValueError): From de02672d6afdee135a3b58236aac645dc78b13c1 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 10:10:13 +0000 Subject: [PATCH 26/44] added more tokamak plotting examples --- .../plot_tokamak_neutron_source_density.py | 59 ++++++++++++++++++ .../plot_tokamak_neutron_source_strengths.py | 61 +++++++++++++++++++ 2 files changed, 120 insertions(+) create mode 100644 examples/plot_tokamak_neutron_source_density.py create mode 100644 examples/plot_tokamak_neutron_source_strengths.py diff --git a/examples/plot_tokamak_neutron_source_density.py b/examples/plot_tokamak_neutron_source_density.py new file mode 100644 index 0000000..066adad --- /dev/null +++ b/examples/plot_tokamak_neutron_source_density.py @@ -0,0 +1,59 @@ +import matplotlib.pyplot as plt +import numpy as np +from openmc_plasma_source import tokamak_ion_temperature, tokamak_convert_a_alpha_to_R_Z, tokamak_neutron_source_density, tokamak_ion_density + +sample_size=20000 +minor_radius=292.258 +major_radius=906 +mode = 'L' +ion_density_centre=45.9 + +# create a sample of (a, alpha) coordinates +a = np.random.random(sample_size) * minor_radius +alpha = np.random.random(sample_size) * 2 * np.pi + +temperatures = tokamak_ion_temperature( + r=a, + mode=mode, + pedestal_radius=0.8 * minor_radius, + ion_temperature_pedestal=6.09, + ion_temperature_centre=ion_density_centre, + ion_temperature_beta=2, + ion_temperature_peaking_factor=8.06, + ion_temperature_separatrix=0.1, + major_radius=major_radius, + ) + +densities = tokamak_ion_density( + mode=mode, + ion_density_centre=ion_density_centre, + ion_density_peaking_factor=1, + ion_density_pedestal=1.09e20, + major_radius=major_radius, + pedestal_radius=0.8 * minor_radius, + ion_density_separatrix=3e19, + r=a, + ) + +neutron_source_density = tokamak_neutron_source_density( + densities, temperatures +) + +RZ = tokamak_convert_a_alpha_to_R_Z( + a=a, + alpha=alpha, + shafranov_factor=0.44789, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=0.270, + elongation=1.557, +) + +plt.scatter(RZ[0], RZ[1], c=neutron_source_density) +plt.gca().set_aspect('equal') +plt.xlabel("R [cm]") +plt.ylabel("Z [cm]") +plt.colorbar(label='neutron source density') + +plt.savefig('tokamak_source_neutron_source_density.png') +print('written tokamak_source_neutron_source_density.png') \ No newline at end of file diff --git a/examples/plot_tokamak_neutron_source_strengths.py b/examples/plot_tokamak_neutron_source_strengths.py new file mode 100644 index 0000000..a95802f --- /dev/null +++ b/examples/plot_tokamak_neutron_source_strengths.py @@ -0,0 +1,61 @@ +import matplotlib.pyplot as plt +import numpy as np +from openmc_plasma_source import tokamak_ion_temperature, tokamak_convert_a_alpha_to_R_Z, tokamak_neutron_source_density, tokamak_ion_density + +sample_size=2000 +minor_radius=292.258 +major_radius=906 +mode = 'L' +ion_density_centre=45.9 + +# create a sample of (a, alpha) coordinates +a = np.random.random(sample_size) * minor_radius +alpha = np.random.random(sample_size) * 2 * np.pi + +temperatures = tokamak_ion_temperature( + r=a, + mode=mode, + pedestal_radius=0.8 * minor_radius, + ion_temperature_pedestal=6.09, + ion_temperature_centre=ion_density_centre, + ion_temperature_beta=2, + ion_temperature_peaking_factor=8.06, + ion_temperature_separatrix=0.1, + major_radius=major_radius, + ) + +densities = tokamak_ion_density( + mode=mode, + ion_density_centre=ion_density_centre, + ion_density_peaking_factor=1, + ion_density_pedestal=1.09e20, + major_radius=major_radius, + pedestal_radius=0.8 * minor_radius, + ion_density_separatrix=3e19, + r=a, + ) + +neutron_source_density = tokamak_neutron_source_density( + densities, temperatures +) + +strengths = neutron_source_density / sum(neutron_source_density) + +RZ = tokamak_convert_a_alpha_to_R_Z( + a=a, + alpha=alpha, + shafranov_factor=0.44789, + minor_radius=minor_radius, + major_radius=major_radius, + triangularity=0.270, + elongation=1.557, +) + +plt.scatter(RZ[0], RZ[1], c=strengths) +plt.gca().set_aspect('equal') +plt.xlabel("R [cm]") +plt.ylabel("Z [cm]") +plt.colorbar(label='neutron emission strength') + +plt.savefig('tokamak_source_neutron_emission_strength.png') +print('written tokamak_source_neutron_emission_strength.png') From 4645fc1c40a610107461f1046c661b70c4c1225d Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 10:10:52 +0000 Subject: [PATCH 27/44] removed internal plotting methods --- src/openmc_plasma_source/plotting/__init__.py | 1 - .../plotting/plot_tokamak_source.py | 153 ---------- tests/test_plotting.py | 261 ------------------ 3 files changed, 415 deletions(-) delete mode 100644 src/openmc_plasma_source/plotting/__init__.py delete mode 100644 src/openmc_plasma_source/plotting/plot_tokamak_source.py delete mode 100644 tests/test_plotting.py diff --git a/src/openmc_plasma_source/plotting/__init__.py b/src/openmc_plasma_source/plotting/__init__.py deleted file mode 100644 index 9ba21fa..0000000 --- a/src/openmc_plasma_source/plotting/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from .plot_tokamak_source import plot_tokamak_source_3D, scatter_tokamak_source diff --git a/src/openmc_plasma_source/plotting/plot_tokamak_source.py b/src/openmc_plasma_source/plotting/plot_tokamak_source.py deleted file mode 100644 index 107a2e1..0000000 --- a/src/openmc_plasma_source/plotting/plot_tokamak_source.py +++ /dev/null @@ -1,153 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from matplotlib import cm - - -def scatter_tokamak_source(source, ax=None, quantity=None, aspect="equal", **kwargs): - """Create a 2D scatter plot of the tokamak source. - See https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.scatter.html - for more arguments. - - Args: - source (ops.TokamakSource): the plasma source - ax (maplotlib.Axes, optional): Axes object on which to plot. If not - provided by the user, the current working Axes is retrieved using - matplotlib.pyplot.gca(). - aspect (str, optional): Set aspect ratio of the plot. May be set to 'auto', - 'equal', or a number denoting the ratio of the plot height to the plot - width. - quantity (str, optional): value by which the lines should be - coloured. Defaults to None. - **kwargs: Keyword arguments compatible with matplotlib.pyplot.scatter - - Raises: - ValueError: If the quantity is unknown - """ - - # Define possible quantities, and link to arrays within tokamak_source - # If given a non-TokamakSource, this step will fail with an AttributeError - try: - quantity_to_attribute = { - "ion_temperature": source.temperatures, - "neutron_source_density": source.neutron_source_density, - "strength": source.strengths, - } - except AttributeError as e: - raise ValueError( - f"openmc_plasma_source.scatter_tokamak_source: argument 'source' " - f"must be of type TokamakSource" - ) from e - - # For a given quantity, determine colours to plot for each point - # If given an incorrect quantity name, this step will fail with a KeyError - colours = None - if quantity is not None: - try: - colours = quantity_to_attribute[quantity] - except KeyError as e: - raise ValueError( - f"openmc_plasma_source.scatter_tokamak_source: Unknown 'quantity' " - f"provided, options are {quantity_to_attribute.keys()}" - ) from e - - # If not provided with an Axes instance, retrieve the current Axes in focus - if ax is None: - ax = plt.gca() - - # Scatter the source R and Z positions, optionally colouring using the chosen - # quantity. - ax.scatter(source.RZ[0], source.RZ[1], c=colours, **kwargs) - - # Set the aspect ratio on the axes. - # Defaults to 'equal', so 1m on the x-axis has the same width as 1m on the y-axis - ax.set_aspect(aspect) - - return ax - - -def plot_tokamak_source_3D( - tokamak_source, - ax=None, - quantity=None, - angles=[0, 0.5 * np.pi], - colorbar=None, - **kwargs, -): - """Creates a 3D plot of the tokamak source. - See https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html#matplotlib.pyplot.plot - for more arguments. - - Args: - tokamak_source (ops.TokamakSource): the plasma source - ax (maplotlib.Axes, optional): Axes object on which to plot. If not - provided by the user, a new Axes is created by calling - matplotlib.pyplot.axes(projection="3d"). - quantity ("str", optional): value by which the lines should be - coloured. Defaults to None. - angles (list, optional): iterable of two floats defining the coverage. - Defaults to [0, 1/2*np.pi]. - colorbar (str, optional): colorbar used if quantity is not None. - When None, matplotlib currently defaults to "viridis". - - Raises: - ValueError: If the quantity is unknown - """ - - # Define possible quantities, and link to arrays within tokamak_source - # If given a non-TokamakSource, this step will fail with an AttributeError - try: - quantity_to_attribute = { - "ion_temperature": tokamak_source.temperatures, - "neutron_source_density": tokamak_source.neutron_source_density, - "strength": tokamak_source.strengths, - } - except AttributeError as e: - raise ValueError( - f"openmc_plasma_source.plot_tokamak_source_3D: argument 'source' " - f"must be of type TokamakSource" - ) from e - - # For a given quantity, get the corresponding array from tokamak_source - # If given an incorrect quantity name, this step will fail with a KeyError - if quantity is not None: - try: - quantity_values = quantity_to_attribute[quantity] - except KeyError as e: - raise ValueError( - f"openmc_plasma_source.plot_tokamak_source_3D: Unknown 'quantity' " - f"provided, options are {quantity_to_attribute.keys()}" - ) from e - else: - quantity_values = np.ones(tokamak_source.sample_size) - - # Get the colours used to plot each curve - # If 'quantity' == None, all have the same colour, selected from the middle - # of the colormap. - cmap = cm.get_cmap(colorbar) - if quantity is not None: - colors = cmap(quantity_values / np.max(quantity_values)) - else: - colors = cmap(np.full(tokamak_source.sample_size, 0.5)) - - # If not provided with an Axes object, create a new one - if ax is None: - ax = plt.axes(projection="3d") - - # Get curves to plot - n_theta = 100 - theta = np.linspace(*angles, n_theta) - xs = np.outer(tokamak_source.RZ[0], np.sin(theta)) - ys = np.outer(tokamak_source.RZ[0], np.cos(theta)) - zs = tokamak_source.RZ[1] - - # Plot each curve in turn - for x, y, z, color in zip(xs, ys, zs, colors): - ax.plot(x, y, z, color=color, **kwargs) - - # Set plot bounds - major_radius = tokamak_source.major_radius - ax.set_xlim(-major_radius, major_radius) - ax.set_ylim(-major_radius, major_radius) - ax.set_zlim(-major_radius, major_radius) - - return ax diff --git a/tests/test_plotting.py b/tests/test_plotting.py deleted file mode 100644 index 8027ea0..0000000 --- a/tests/test_plotting.py +++ /dev/null @@ -1,261 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -import pytest - -from openmc_plasma_source import tokamak_source -from openmc_plasma_source import plotting as ops_plt - - -@pytest.fixture -def tokamak_source(): - return tokamak_source( - elongation=1.557, - ion_density_centre=1.09e20, - ion_density_peaking_factor=1, - ion_density_pedestal=1.09e20, - ion_density_separatrix=3e19, - ion_temperature_centre=45.9, - ion_temperature_peaking_factor=8.06, - ion_temperature_pedestal=6.09, - ion_temperature_separatrix=0.1, - major_radius=9.06, - minor_radius=2.92258, - pedestal_radius=0.8 * 2.92258, - mode="H", - shafranov_factor=0.44789, - triangularity=0.270, - ion_temperature_beta=6, - sample_size=500, - ) - - -def test_scatter_tokamak_source_defaults(tokamak_source): - """Ensure plotting is successful without providing additional args""" - plt.figure() - assert not plt.gca().collections # Check current ax is empty - ops_plt.scatter_tokamak_source(tokamak_source) - assert plt.gca().collections # Check current ax is not empty - # Save for viewing, clean up - plt.xlabel("R") - plt.ylabel("Z", rotation=0) - plt.savefig("tests/test_scatter_tokamak_source_defaults.png") - plt.close() - - -def test_scatter_tokamak_source_with_ax(tokamak_source): - """Ensure plotting is successful for user-provided ax""" - fig = plt.figure() - ax = fig.gca() - assert not ax.collections # Check ax is empty - ops_plt.scatter_tokamak_source(tokamak_source, ax=ax) - assert ax.collections # Check ax is not empty - # Save for viewing, clean up - ax.set_xlabel("R") - ax.set_ylabel("Z", rotation=0) - fig.savefig("tests/test_scatter_tokamak_source_with_ax.png") - plt.close(fig) - - -def test_scatter_tokamak_source_with_subplots(tokamak_source): - """Ensure plotting is successful for multiple user-provided ax""" - fig, (ax1, ax2) = plt.subplots(1, 2) - # Plot on the first axes - assert not ax1.collections # Check ax is empty - ops_plt.scatter_tokamak_source(tokamak_source, ax=ax1) - assert ax1.collections # Check ax is not empty - # Generate new data - tokamak_source.sample_sources() - tokamak_source.sources = tokamak_source.make_openmc_sources() - # Plot on the other axes - assert not ax2.collections # Check ax is empty - ops_plt.scatter_tokamak_source(tokamak_source, ax=ax2) - assert ax2.collections # Check ax is not empty - # Save for viewing, clean up - ax1.set_xlabel("R") - ax1.set_ylabel("Z", rotation=0) - ax2.set_xlabel("R") - fig.savefig("tests/test_scatter_tokamak_source_subplots.png") - plt.close(fig) - - -@pytest.mark.parametrize( - "quantity", ["ion_temperature", "neutron_source_density", "strength"] -) -def test_scatter_tokamak_source_quantities(tokamak_source, quantity): - """Plot with colours set by 'quantity'""" - fig = plt.figure() - ax = fig.gca() - assert not ax.collections # Check ax is empty - ops_plt.scatter_tokamak_source(tokamak_source, ax=ax, quantity=quantity) - assert ax.collections # Check ax is not empty - # Save for viewing, clean up - ax.set_xlabel("R") - ax.set_ylabel("Z", rotation=0) - fig.savefig(f"tests/test_scatter_tokamak_source_quantities_{quantity}.png") - plt.close(fig) - - -@pytest.mark.parametrize("aspect", ["equal", "auto", 2]) -def test_scatter_tokamak_source_aspect(tokamak_source, aspect): - """Plot with various aspect ratios""" - fig = plt.figure() - ax = fig.gca() - assert not ax.collections # Check ax is empty - ops_plt.scatter_tokamak_source(tokamak_source, ax=ax, aspect=aspect) - assert ax.collections # Check ax is not empty - # Save for viewing, clean up - ax.set_xlabel("R") - ax.set_ylabel("Z", rotation=0) - fig.savefig(f"tests/test_scatter_tokamak_source_aspect_{aspect}.png") - plt.close(fig) - - -@pytest.mark.parametrize("kwargs", [{"alpha": 0.2}, {"marker": "x"}]) -def test_scatter_tokamak_source_kwargs(tokamak_source, kwargs): - """Plot with a kwarg compatible with 'scatter'""" - fig = plt.figure() - ax = fig.gca() - assert not ax.collections # Check ax is empty - ops_plt.scatter_tokamak_source(tokamak_source, ax=ax, **kwargs) - assert ax.collections # Check ax is not empty - # Save for viewing, clean up - ax.set_xlabel("R") - ax.set_ylabel("Z", rotation=0) - fig.savefig( - f"tests/test_scatter_tokamak_source_kwargs_{list(kwargs.keys())[0]}.png" - ) - plt.close(fig) - - -def test_scatter_tokamak_not_source(): - """Ensure failure when given non-tokamak_source to plot""" - with pytest.raises(ValueError) as excinfo: - fig = plt.figure() - ax = fig.gca() - ops_plt.scatter_tokamak_source("hello world", ax=ax) - plt.close() - assert "tokamak_source" in str(excinfo.value) - - -@pytest.mark.parametrize("quantity", ["coucou", "ion_density", 17]) -def test_scatter_tokamak_wrong_quantity(tokamak_source, quantity): - """Ensure failure when incorrect quantity specified""" - with pytest.raises(ValueError) as excinfo: - fig = plt.figure() - ax = fig.gca() - ops_plt.scatter_tokamak_source(tokamak_source, ax=ax, quantity=quantity) - plt.close() - assert "quantity" in str(excinfo.value) - - -def test_plot_tokamak_source_3D_default(tokamak_source): - """Ensure plots correctly with default inputs""" - plt.figure() - ops_plt.plot_tokamak_source_3D(tokamak_source) - assert plt.gca().lines # Check current ax is not empty - # Save for viewing, clean up - plt.savefig("tests/test_plot_tokamak_source_3D_defaults.png") - plt.close() - - -def test_plot_tokamak_source_3D_with_ax(tokamak_source): - """Ensure plots correctly given ax instance""" - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - assert not ax.lines # Check ax is empty - ops_plt.plot_tokamak_source_3D(tokamak_source, ax=ax) - assert ax.lines # Check ax is not empty - # Save for viewing, clean up - fig.savefig("tests/test_plot_tokamak_source_3D_with_ax.png") - plt.close(fig) - - -@pytest.mark.parametrize( - "quantity", ["ion_temperature", "neutron_source_density", "strength"] -) -def test_plot_tokamak_source_3D_quantities(tokamak_source, quantity): - """Ensure plots correctly for each quantity""" - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - assert not ax.lines # Check ax is empty - ops_plt.plot_tokamak_source_3D(tokamak_source, ax=ax, quantity=quantity) - assert ax.lines # Check ax is not empty - # Save for viewing, clean up - fig.savefig(f"tests/test_plot_tokamak_source_3D_quantities_{quantity}.png") - plt.close(fig) - - -@pytest.mark.parametrize("colorbar", ["plasma", "rainbow"]) -def test_plot_tokamak_source_3D_colorbars(tokamak_source, colorbar): - """Ensure plots correctly given colorbar choice""" - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - assert not ax.lines # Check ax is empty - ops_plt.plot_tokamak_source_3D( - tokamak_source, ax=ax, quantity="ion_temperature", colorbar=colorbar - ) - assert ax.lines # Check ax is not empty - # Save for viewing, clean up - fig.savefig(f"tests/test_plot_tokamak_source_3D_colorbar_{colorbar}.png") - plt.close(fig) - - -@pytest.mark.parametrize( - "angles,name", - [ - [(0, np.pi / 4), "eighth"], - [(0, np.pi), "half"], - [(np.pi, 2 * np.pi), "half_offset"], - [(0, 2 * np.pi), "full"], - ], -) -def test_plot_tokamak_source_3D_angles(tokamak_source, angles, name): - """Ensure plots correctly given angles range""" - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - assert not ax.lines # Check ax is empty - ops_plt.plot_tokamak_source_3D( - tokamak_source, ax=ax, quantity="ion_temperature", angles=angles - ) - assert ax.lines # Check ax is not empty - # Save for viewing, clean up - fig.savefig(f"tests/test_plot_tokamak_source_3D_angles_{name}.png") - plt.close(fig) - - -@pytest.mark.parametrize("kwargs", [{"alpha": 0.2}, {"linestyle": "dotted"}]) -def test_plot_tokamak_source_3D_kwargs(tokamak_source, kwargs): - """Ensure plots correctly given additonal keyword arguments to pass on to plot""" - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - assert not ax.lines # Check ax is empty - ops_plt.plot_tokamak_source_3D( - tokamak_source, ax=ax, quantity="ion_temperature", **kwargs - ) - assert ax.lines # Check ax is not empty - # Save for viewing, clean up - fig.savefig( - f"tests/test_plot_tokamak_source_3D_kwargs_{list(kwargs.keys())[0]}.png" - ) - plt.close(fig) - - -def test_plot_tokamak_source_3D_not_source(): - """Ensure failure when given non-tokamak_source to plot""" - with pytest.raises(ValueError) as excinfo: - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - ops_plt.plot_tokamak_source_3D("hello world", ax=ax) - plt.close() - assert "tokamak_source" in str(excinfo.value) - - -@pytest.mark.parametrize("quantity", ["coucou", "ion_density", 17]) -def test_plot_tokamak_source_3D_wrong_quantity(tokamak_source, quantity): - """Ensure failure when incorrect quantity specified""" - with pytest.raises(ValueError) as excinfo: - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection="3d") - ops_plt.plot_tokamak_source_3D(tokamak_source, ax=ax, quantity=quantity) - plt.close() - assert "quantity" in str(excinfo.value) From bf025739e626edd907f96ff60d80b71bd9c5882b Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 10:21:10 +0000 Subject: [PATCH 28/44] [skip ci] Apply formatting changes --- .../plot_tokamak_neutron_source_density.py | 67 ++++++++++--------- .../plot_tokamak_neutron_source_strengths.py | 67 ++++++++++--------- 2 files changed, 70 insertions(+), 64 deletions(-) diff --git a/examples/plot_tokamak_neutron_source_density.py b/examples/plot_tokamak_neutron_source_density.py index 066adad..f942fec 100644 --- a/examples/plot_tokamak_neutron_source_density.py +++ b/examples/plot_tokamak_neutron_source_density.py @@ -1,44 +1,47 @@ import matplotlib.pyplot as plt import numpy as np -from openmc_plasma_source import tokamak_ion_temperature, tokamak_convert_a_alpha_to_R_Z, tokamak_neutron_source_density, tokamak_ion_density +from openmc_plasma_source import ( + tokamak_ion_temperature, + tokamak_convert_a_alpha_to_R_Z, + tokamak_neutron_source_density, + tokamak_ion_density, +) -sample_size=20000 -minor_radius=292.258 -major_radius=906 -mode = 'L' -ion_density_centre=45.9 +sample_size = 20000 +minor_radius = 292.258 +major_radius = 906 +mode = "L" +ion_density_centre = 45.9 # create a sample of (a, alpha) coordinates a = np.random.random(sample_size) * minor_radius alpha = np.random.random(sample_size) * 2 * np.pi temperatures = tokamak_ion_temperature( - r=a, - mode=mode, - pedestal_radius=0.8 * minor_radius, - ion_temperature_pedestal=6.09, - ion_temperature_centre=ion_density_centre, - ion_temperature_beta=2, - ion_temperature_peaking_factor=8.06, - ion_temperature_separatrix=0.1, - major_radius=major_radius, - ) + r=a, + mode=mode, + pedestal_radius=0.8 * minor_radius, + ion_temperature_pedestal=6.09, + ion_temperature_centre=ion_density_centre, + ion_temperature_beta=2, + ion_temperature_peaking_factor=8.06, + ion_temperature_separatrix=0.1, + major_radius=major_radius, +) densities = tokamak_ion_density( - mode=mode, - ion_density_centre=ion_density_centre, - ion_density_peaking_factor=1, - ion_density_pedestal=1.09e20, - major_radius=major_radius, - pedestal_radius=0.8 * minor_radius, - ion_density_separatrix=3e19, - r=a, - ) - -neutron_source_density = tokamak_neutron_source_density( - densities, temperatures + mode=mode, + ion_density_centre=ion_density_centre, + ion_density_peaking_factor=1, + ion_density_pedestal=1.09e20, + major_radius=major_radius, + pedestal_radius=0.8 * minor_radius, + ion_density_separatrix=3e19, + r=a, ) +neutron_source_density = tokamak_neutron_source_density(densities, temperatures) + RZ = tokamak_convert_a_alpha_to_R_Z( a=a, alpha=alpha, @@ -50,10 +53,10 @@ ) plt.scatter(RZ[0], RZ[1], c=neutron_source_density) -plt.gca().set_aspect('equal') +plt.gca().set_aspect("equal") plt.xlabel("R [cm]") plt.ylabel("Z [cm]") -plt.colorbar(label='neutron source density') +plt.colorbar(label="neutron source density") -plt.savefig('tokamak_source_neutron_source_density.png') -print('written tokamak_source_neutron_source_density.png') \ No newline at end of file +plt.savefig("tokamak_source_neutron_source_density.png") +print("written tokamak_source_neutron_source_density.png") diff --git a/examples/plot_tokamak_neutron_source_strengths.py b/examples/plot_tokamak_neutron_source_strengths.py index a95802f..21671a9 100644 --- a/examples/plot_tokamak_neutron_source_strengths.py +++ b/examples/plot_tokamak_neutron_source_strengths.py @@ -1,44 +1,47 @@ import matplotlib.pyplot as plt import numpy as np -from openmc_plasma_source import tokamak_ion_temperature, tokamak_convert_a_alpha_to_R_Z, tokamak_neutron_source_density, tokamak_ion_density +from openmc_plasma_source import ( + tokamak_ion_temperature, + tokamak_convert_a_alpha_to_R_Z, + tokamak_neutron_source_density, + tokamak_ion_density, +) -sample_size=2000 -minor_radius=292.258 -major_radius=906 -mode = 'L' -ion_density_centre=45.9 +sample_size = 2000 +minor_radius = 292.258 +major_radius = 906 +mode = "L" +ion_density_centre = 45.9 # create a sample of (a, alpha) coordinates a = np.random.random(sample_size) * minor_radius alpha = np.random.random(sample_size) * 2 * np.pi temperatures = tokamak_ion_temperature( - r=a, - mode=mode, - pedestal_radius=0.8 * minor_radius, - ion_temperature_pedestal=6.09, - ion_temperature_centre=ion_density_centre, - ion_temperature_beta=2, - ion_temperature_peaking_factor=8.06, - ion_temperature_separatrix=0.1, - major_radius=major_radius, - ) + r=a, + mode=mode, + pedestal_radius=0.8 * minor_radius, + ion_temperature_pedestal=6.09, + ion_temperature_centre=ion_density_centre, + ion_temperature_beta=2, + ion_temperature_peaking_factor=8.06, + ion_temperature_separatrix=0.1, + major_radius=major_radius, +) densities = tokamak_ion_density( - mode=mode, - ion_density_centre=ion_density_centre, - ion_density_peaking_factor=1, - ion_density_pedestal=1.09e20, - major_radius=major_radius, - pedestal_radius=0.8 * minor_radius, - ion_density_separatrix=3e19, - r=a, - ) - -neutron_source_density = tokamak_neutron_source_density( - densities, temperatures + mode=mode, + ion_density_centre=ion_density_centre, + ion_density_peaking_factor=1, + ion_density_pedestal=1.09e20, + major_radius=major_radius, + pedestal_radius=0.8 * minor_radius, + ion_density_separatrix=3e19, + r=a, ) +neutron_source_density = tokamak_neutron_source_density(densities, temperatures) + strengths = neutron_source_density / sum(neutron_source_density) RZ = tokamak_convert_a_alpha_to_R_Z( @@ -52,10 +55,10 @@ ) plt.scatter(RZ[0], RZ[1], c=strengths) -plt.gca().set_aspect('equal') +plt.gca().set_aspect("equal") plt.xlabel("R [cm]") plt.ylabel("Z [cm]") -plt.colorbar(label='neutron emission strength') +plt.colorbar(label="neutron emission strength") -plt.savefig('tokamak_source_neutron_emission_strength.png') -print('written tokamak_source_neutron_emission_strength.png') +plt.savefig("tokamak_source_neutron_emission_strength.png") +print("written tokamak_source_neutron_emission_strength.png") From 2710cb4dc3eeb12913e3fecc88e2b9be260b0ac0 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 10:39:09 +0000 Subject: [PATCH 29/44] [skip ci] fixed some tok tests --- tests/test_tokamak_source.py | 223 +++++++++++++++++------------------ 1 file changed, 111 insertions(+), 112 deletions(-) diff --git a/tests/test_tokamak_source.py b/tests/test_tokamak_source.py index 28845c4..ee17022 100644 --- a/tests/test_tokamak_source.py +++ b/tests/test_tokamak_source.py @@ -3,7 +3,7 @@ from hypothesis import given, settings from hypothesis import strategies as st import openmc -from openmc_plasma_source import tokamak_source +from openmc_plasma_source import tokamak_source, tokamak_ion_density @pytest.fixture @@ -39,7 +39,7 @@ def tokamak_source_example(tokamak_args_dict): def test_creation(tokamak_source_example): """Tests that the sources generated by tokamak_source are of type openmc.Source""" - for source in tokamak_source_example.sources: + for source in tokamak_source_example: assert isinstance(source, openmc.IndependentSource) @@ -50,7 +50,7 @@ def test_major_radius(tokamak_args_dict, minor_radius, major_radius): """Checks that tokamak_source creation accepts valid major radius""" tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["major_radius"] = major_radius - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -61,7 +61,7 @@ def test_bad_major_radius(tokamak_args_dict, minor_radius, major_radius): tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["major_radius"] = major_radius with pytest.raises(ValueError): - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -74,7 +74,7 @@ def test_minor_radius(tokamak_args_dict, major_radius, minor_radius): # Set shafranov factor to 0 and pedestal factor to 0.8*minor_radius for safety tokamak_args_dict["pedestal_radius"] = 0.8 * minor_radius tokamak_args_dict["shafranov_factor"] = 0.0 - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -86,14 +86,14 @@ def test_bad_minor_radius(tokamak_args_dict, major_radius, minor_radius): tokamak_args_dict["major_radius"] = major_radius tokamak_args_dict["minor_radius"] = minor_radius with pytest.raises(ValueError): - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("elongation", [1.0, 1.667, 0.5, 20, 0.001]) def test_elongation(tokamak_args_dict, elongation): """Checks that tokamak_source creation accepts valid elongation""" tokamak_args_dict["elongation"] = elongation - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("elongation", [0, -5, "hello world"]) @@ -101,14 +101,14 @@ def test_bad_elongation(tokamak_args_dict, elongation): """Checks that tokamak_source creation rejects invalid elongation""" tokamak_args_dict["elongation"] = elongation with pytest.raises(ValueError): - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("triangularity", [0.0, 0.5, 0.9, 1.0, -0.5, -0.9, -1.0]) def test_triangularity(tokamak_args_dict, triangularity): """Checks that tokamak_source creation accepts valid triangularity""" tokamak_args_dict["triangularity"] = triangularity - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize("triangularity", [1.1, -1.1, 10, -10, "hello world"]) @@ -116,7 +116,7 @@ def test_bad_triangularity(tokamak_args_dict, triangularity): """Checks that tokamak_source creation rejects invalid triangularity""" tokamak_args_dict["triangularity"] = triangularity with pytest.raises(ValueError): - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -137,7 +137,7 @@ def test_shafranov_factor(tokamak_args_dict, major_radius, minor_radius, shaf): tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["pedestal_radius"] = 0.8 * minor_radius tokamak_args_dict["shafranov_factor"] = shaf - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -158,7 +158,7 @@ def test_bad_shafranov_factor(tokamak_args_dict, major_radius, minor_radius, sha tokamak_args_dict["pedestal_radius"] = 0.8 * minor_radius tokamak_args_dict["shafranov_factor"] = shaf with pytest.raises((ValueError, TypeError)): - tokamak_source = tokamak_source(**tokamak_args_dict) + tokamak_source(**tokamak_args_dict) @pytest.mark.parametrize( @@ -167,10 +167,9 @@ def test_bad_shafranov_factor(tokamak_args_dict, major_radius, minor_radius, sha def test_angles(tokamak_args_dict, angles): """Checks that custom angles can be set""" # Note: should accept negative angles and angles in reverse order - tokamak_args_dict["angles"] = angles - tokamak_source = tokamak_source(**tokamak_args_dict) - assert np.array_equal(tokamak_source.angles, angles) - for source in tokamak_source.sources: + angles = tokamak_args_dict["angles"] + sources = tokamak_source(**tokamak_args_dict) + for source in sources: assert np.array_equal((source.space.phi.a, source.space.phi.b), angles) @@ -184,60 +183,60 @@ def test_bad_angles(tokamak_args_dict, angles): tokamak_source = tokamak_source(**tokamak_args_dict) -def test_ion_density(tokamak_source_example): - # test with values of r that are within acceptable ranges. - r = np.linspace(0.0, tokamak_source_example.minor_radius, 100) - density = tokamak_source_example.ion_density(r) - assert isinstance(r, np.ndarray) - assert len(density) == len(r) - assert np.all(np.isfinite(density)) +# def test_ion_density(tokamak_source_example): +# # test with values of r that are within acceptable ranges. +# r = np.linspace(0.0, tokamak_source_example.minor_radius, 100) +# density = tokamak_source_example.ion_density(r) +# assert isinstance(r, np.ndarray) +# assert len(density) == len(r) +# assert np.all(np.isfinite(density)) -def test_bad_ion_density(tokamak_source_example): - # It should fail if given a negative r - with pytest.raises(ValueError) as excinfo: - density = tokamak_source_example.ion_density([0, 5, -6]) - assert "must not be negative" in str(excinfo.value) +# def test_bad_ion_density(tokamak_source_example): +# # It should fail if given a negative r +# with pytest.raises(ValueError) as excinfo: +# density = tokamak_ion_density([0, 5, -6]) +# assert "must not be negative" in str(excinfo.value) -def test_ion_temperature(tokamak_source_example): - # test with values of r that are within acceptable ranges. - r = np.linspace(0.0, tokamak_source_example.minor_radius, 100) - temperature = tokamak_source_example.ion_temperature(r) - assert isinstance(temperature, np.ndarray) - assert len(temperature) == len(r) - assert np.all(np.isfinite(temperature)) +# def test_ion_temperature(tokamak_source_example): +# # test with values of r that are within acceptable ranges. +# r = np.linspace(0.0, tokamak_source_example.minor_radius, 100) +# temperature = tokamak_source_example.ion_temperature(r) +# assert isinstance(temperature, np.ndarray) +# assert len(temperature) == len(r) +# assert np.all(np.isfinite(temperature)) -def test_bad_ion_temperature(tokamak_source_example): - # It should fail if given a negative r - with pytest.raises(ValueError) as excinfo: - temperature = tokamak_source_example.ion_temperature([0, 5, -6]) - assert "must not be negative" in str(excinfo.value) - - -def test_convert_a_alpha_to_R_Z(tokamak_source_example): - # Similar to test_source_locations_are_within_correct_range - # Rather than going in detail, simply tests validity of inputs and outputs - # Test with suitable values for a and alpha - a = np.linspace(0.0, tokamak_source_example.minor_radius, 100) - alpha = np.linspace(0.0, 2 * np.pi, 100) - R, Z = tokamak_source_example.convert_a_alpha_to_R_Z(a, alpha) - assert isinstance(R, np.ndarray) - assert isinstance(Z, np.ndarray) - assert len(R) == len(a) - assert len(Z) == len(a) - assert np.all(np.isfinite(R)) - assert np.all(np.isfinite(Z)) - - -def test_bad_convert_a_alpha_to_R_Z(tokamak_source_example): - # Repeat test_convert_a_alpha_to_R_Z, but show that negative a breaks it - a = np.linspace(0.0, tokamak_source_example.minor_radius, 100) - alpha = np.linspace(0.0, 2 * np.pi, 100) - with pytest.raises(ValueError) as excinfo: - R, Z = tokamak_source_example.convert_a_alpha_to_R_Z(-a, alpha) - assert "must not be negative" in str(excinfo.value) +# def test_bad_ion_temperature(tokamak_source_example): +# # It should fail if given a negative r +# with pytest.raises(ValueError) as excinfo: +# temperature = tokamak_source_example.ion_temperature([0, 5, -6]) +# assert "must not be negative" in str(excinfo.value) + + +# def test_convert_a_alpha_to_R_Z(tokamak_source_example): +# # Similar to test_source_locations_are_within_correct_range +# # Rather than going in detail, simply tests validity of inputs and outputs +# # Test with suitable values for a and alpha +# a = np.linspace(0.0, tokamak_source_example.minor_radius, 100) +# alpha = np.linspace(0.0, 2 * np.pi, 100) +# R, Z = tokamak_source_example.convert_a_alpha_to_R_Z(a, alpha) +# assert isinstance(R, np.ndarray) +# assert isinstance(Z, np.ndarray) +# assert len(R) == len(a) +# assert len(Z) == len(a) +# assert np.all(np.isfinite(R)) +# assert np.all(np.isfinite(Z)) + + +# def test_bad_convert_a_alpha_to_R_Z(tokamak_source_example): +# # Repeat test_convert_a_alpha_to_R_Z, but show that negative a breaks it +# a = np.linspace(0.0, tokamak_source_example.minor_radius, 100) +# alpha = np.linspace(0.0, 2 * np.pi, 100) +# with pytest.raises(ValueError) as excinfo: +# R, Z = tokamak_source_example.convert_a_alpha_to_R_Z(-a, alpha) +# assert "must not be negative" in str(excinfo.value) @st.composite @@ -307,51 +306,51 @@ def tokamak_source_strategy(draw): ) -@given(tokamak_source=tokamak_source_strategy()) -@settings(max_examples=50) -def test_strengths_are_normalised(tokamak_source): - """Tests that the sum of the strengths attribute is equal to""" - assert pytest.approx(sum(tokamak_source.strengths)) == 1 - - -@given(tokamak_source=tokamak_source_strategy()) -@settings(max_examples=50) -def test_source_locations_are_within_correct_range(tokamak_source): - """Tests that each source has RZ locations within the expected range. - - As the function converting (a,alpha) coordinates to (R,Z) is not bijective, - we cannot convert back to validate each individual point. However, we can - determine whether the generated points are contained within the shell of - the last closed magnetic surface. See "Tokamak D-T neutron source models - for different plasma physics confinement modes", C. Fausser et al., Fusion - Engineering and Design, 2012 for more info. - """ - R_0 = tokamak_source.major_radius - A = tokamak_source.minor_radius - El = tokamak_source.elongation - delta = tokamak_source.triangularity - - def get_R_on_LCMS(alpha): - """Gets R on the last closed magnetic surface for a given alpha""" - return R_0 + A * np.cos(alpha + delta * np.sin(alpha)) - - approx_lt = lambda x, y: x < y or np.isclose(x, y) - approx_gt = lambda x, y: x > y or np.isclose(x, y) - - for source in tokamak_source.sources: - R, Z = source.space.r.x[0], source.space.z.x[0] - # First test that the point is contained with a simple box with - # lower left (r_min,-z_max) and upper right (r_max,z_max) - assert approx_gt(R, R_0 - A) - assert approx_lt(R, R_0 + A) - assert approx_lt(abs(Z), A * El) - # For a given Z, we can determine the two values of alpha where - # where a = minor_radius, and from there determine the upper and - # lower bounds for R. - alpha_1 = np.arcsin(abs(Z) / (El * A)) - alpha_2 = np.pi - alpha_1 - R_max, R_min = get_R_on_LCMS(alpha_1), get_R_on_LCMS(alpha_2) - assert approx_lt(R_max, R_0 + A) - assert approx_gt(R_min, R_0 - A) - assert approx_lt(R, R_max) - assert approx_gt(R, R_min) +# @given(tokamak_source=tokamak_source_strategy()) +# @settings(max_examples=50) +# def test_strengths_are_normalised(tokamak_source): +# """Tests that the sum of the strengths attribute is equal to""" +# assert pytest.approx(sum(tokamak_source.strengths)) == 1 + + +# @given(tokamak_source=tokamak_source_strategy()) +# @settings(max_examples=50) +# def test_source_locations_are_within_correct_range(tokamak_source): +# """Tests that each source has RZ locations within the expected range. + +# As the function converting (a,alpha) coordinates to (R,Z) is not bijective, +# we cannot convert back to validate each individual point. However, we can +# determine whether the generated points are contained within the shell of +# the last closed magnetic surface. See "Tokamak D-T neutron source models +# for different plasma physics confinement modes", C. Fausser et al., Fusion +# Engineering and Design, 2012 for more info. +# """ +# R_0 = tokamak_source.major_radius +# A = tokamak_source.minor_radius +# El = tokamak_source.elongation +# delta = tokamak_source.triangularity + +# def get_R_on_LCMS(alpha): +# """Gets R on the last closed magnetic surface for a given alpha""" +# return R_0 + A * np.cos(alpha + delta * np.sin(alpha)) + +# approx_lt = lambda x, y: x < y or np.isclose(x, y) +# approx_gt = lambda x, y: x > y or np.isclose(x, y) + +# for source in tokamak_source.sources: +# R, Z = source.space.r.x[0], source.space.z.x[0] +# # First test that the point is contained with a simple box with +# # lower left (r_min,-z_max) and upper right (r_max,z_max) +# assert approx_gt(R, R_0 - A) +# assert approx_lt(R, R_0 + A) +# assert approx_lt(abs(Z), A * El) +# # For a given Z, we can determine the two values of alpha where +# # where a = minor_radius, and from there determine the upper and +# # lower bounds for R. +# alpha_1 = np.arcsin(abs(Z) / (El * A)) +# alpha_2 = np.pi - alpha_1 +# R_max, R_min = get_R_on_LCMS(alpha_1), get_R_on_LCMS(alpha_2) +# assert approx_lt(R_max, R_0 + A) +# assert approx_gt(R_min, R_0 - A) +# assert approx_lt(R, R_max) +# assert approx_gt(R, R_min) From 4ecbda158f0cf1da3ba46e7d32c2fa0c8e60ca3b Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 8 Mar 2024 13:59:05 +0000 Subject: [PATCH 30/44] toko tests working --- src/openmc_plasma_source/tokamak_source.py | 46 ++-- tests/test_tokamak_source.py | 271 +++++++++++++-------- 2 files changed, 195 insertions(+), 122 deletions(-) diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 8bca1fb..a0f046d 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -73,6 +73,31 @@ def tokamak_source( """ # Perform sanity checks for inputs not caught by properties + + if not isinstance(major_radius, (int, float)): + raise ValueError("Major radius must be a number") + + if not isinstance(minor_radius, (int, float)): + raise ValueError("Minor radius must be a number") + + if not isinstance(elongation, (int, float)): + raise ValueError("Elongation must be a number") + + if not isinstance(triangularity, (int, float)): + raise ValueError("Triangularity must be a number") + + if not isinstance(ion_density_centre, (int, float)): + raise ValueError("ion_density_centre must be a number") + + if not isinstance(ion_density_peaking_factor, (int, float)): + raise ValueError("ion_density_peaking_factor must be a number") + + if not isinstance(ion_density_pedestal, (int, float)): + raise ValueError("ion_density_pedestal must be a number") + + if not isinstance(ion_density_separatrix, (int, float)): + raise ValueError("ion_density_separatrix must be a number") + if minor_radius >= major_radius: raise ValueError("Minor radius must be smaller than major radius") @@ -82,24 +107,15 @@ def tokamak_source( if abs(shafranov_factor) >= 0.5 * minor_radius: raise ValueError("Shafranov factor must be smaller than 0.5*minor radius") - if not isinstance(major_radius, (int, float)): - raise ValueError("Major radius must be a number") - if major_radius < 0: raise ValueError("Major radius must greater than 0") - if not isinstance(minor_radius, (int, float)): - raise ValueError("Minor radius must be a number") if minor_radius < 0: raise ValueError("Minor radius must greater than 0") - if not isinstance(elongation, (int, float)): - raise ValueError("Elongation must be a number") - if elongation < 0: + if elongation <= 0: raise ValueError("Elongation must greater than 0") - if not isinstance(triangularity, (int, float)): - raise ValueError("Triangularity must be a number") if not triangularity <= 1.0: raise ValueError("Triangularity must less than or equal to 1.") if not triangularity >= -1.0: @@ -108,21 +124,13 @@ def tokamak_source( if mode not in ["H", "L", "A"]: raise ValueError("Mode must be one of the following: ['H', 'L', 'A']") - if not isinstance(ion_density_centre, (int, float)): - raise ValueError("ion_density_centre must be a number") if ion_density_centre < 0: raise ValueError("ion_density_centre must greater than 0") - if not isinstance(ion_density_peaking_factor, (int, float)): - raise ValueError("ion_density_peaking_factor must be a number") - if not isinstance(ion_density_pedestal, (int, float)): - raise ValueError("ion_density_pedestal must be a number") if ion_density_pedestal < 0: raise ValueError("ion_density_pedestal must greater than 0") - if not isinstance(ion_density_separatrix, (int, float)): - raise ValueError("ion_density_separatrix must be a number") if ion_density_separatrix < 0: raise ValueError("ion_density_separatrix must greater than 0") @@ -278,7 +286,7 @@ def tokamak_ion_temperature( ( ion_temperature_pedestal + (ion_temperature_centre - ion_temperature_pedestal) - * (1 - (r / pedestal_radius) ** ion_temperature_beta) + * (1 - (np.abs(r / pedestal_radius)) ** ion_temperature_beta) ** ion_temperature_peaking_factor ), ( diff --git a/tests/test_tokamak_source.py b/tests/test_tokamak_source.py index ee17022..4b5deed 100644 --- a/tests/test_tokamak_source.py +++ b/tests/test_tokamak_source.py @@ -1,9 +1,9 @@ import numpy as np import pytest -from hypothesis import given, settings +from hypothesis import given, settings, HealthCheck from hypothesis import strategies as st import openmc -from openmc_plasma_source import tokamak_source, tokamak_ion_density +from openmc_plasma_source import tokamak_source, tokamak_ion_density, tokamak_ion_temperature, tokamak_convert_a_alpha_to_R_Z @pytest.fixture @@ -167,7 +167,7 @@ def test_bad_shafranov_factor(tokamak_args_dict, major_radius, minor_radius, sha def test_angles(tokamak_args_dict, angles): """Checks that custom angles can be set""" # Note: should accept negative angles and angles in reverse order - angles = tokamak_args_dict["angles"] + tokamak_args_dict["angles"] = angles sources = tokamak_source(**tokamak_args_dict) for source in sources: assert np.array_equal((source.space.phi.a, source.space.phi.b), angles) @@ -180,63 +180,119 @@ def test_bad_angles(tokamak_args_dict, angles): # Contents should convert to float tokamak_args_dict["angles"] = angles with pytest.raises(ValueError) as excinfo: - tokamak_source = tokamak_source(**tokamak_args_dict) - - -# def test_ion_density(tokamak_source_example): -# # test with values of r that are within acceptable ranges. -# r = np.linspace(0.0, tokamak_source_example.minor_radius, 100) -# density = tokamak_source_example.ion_density(r) -# assert isinstance(r, np.ndarray) -# assert len(density) == len(r) -# assert np.all(np.isfinite(density)) - - -# def test_bad_ion_density(tokamak_source_example): -# # It should fail if given a negative r -# with pytest.raises(ValueError) as excinfo: -# density = tokamak_ion_density([0, 5, -6]) -# assert "must not be negative" in str(excinfo.value) - - -# def test_ion_temperature(tokamak_source_example): -# # test with values of r that are within acceptable ranges. -# r = np.linspace(0.0, tokamak_source_example.minor_radius, 100) -# temperature = tokamak_source_example.ion_temperature(r) -# assert isinstance(temperature, np.ndarray) -# assert len(temperature) == len(r) -# assert np.all(np.isfinite(temperature)) + tokamak_source(**tokamak_args_dict) -# def test_bad_ion_temperature(tokamak_source_example): -# # It should fail if given a negative r -# with pytest.raises(ValueError) as excinfo: -# temperature = tokamak_source_example.ion_temperature([0, 5, -6]) -# assert "must not be negative" in str(excinfo.value) +def test_ion_density(tokamak_args_dict, tokamak_source_example): + # test with values of r that are within acceptable ranges. + r = np.linspace(0.0, tokamak_args_dict['minor_radius'], 100) + density = tokamak_ion_density( + r=r, + mode='L', + ion_density_centre=tokamak_args_dict['ion_density_centre'], + ion_density_peaking_factor=tokamak_args_dict['ion_density_peaking_factor'], + ion_density_pedestal=tokamak_args_dict['ion_density_pedestal'], + major_radius=tokamak_args_dict['major_radius'], + pedestal_radius=tokamak_args_dict['pedestal_radius'], + ion_density_separatrix=tokamak_args_dict['ion_density_separatrix'], + ) + assert isinstance(r, np.ndarray) + assert len(density) == len(r) + assert np.all(np.isfinite(density)) -# def test_convert_a_alpha_to_R_Z(tokamak_source_example): -# # Similar to test_source_locations_are_within_correct_range -# # Rather than going in detail, simply tests validity of inputs and outputs -# # Test with suitable values for a and alpha -# a = np.linspace(0.0, tokamak_source_example.minor_radius, 100) -# alpha = np.linspace(0.0, 2 * np.pi, 100) -# R, Z = tokamak_source_example.convert_a_alpha_to_R_Z(a, alpha) -# assert isinstance(R, np.ndarray) -# assert isinstance(Z, np.ndarray) -# assert len(R) == len(a) -# assert len(Z) == len(a) -# assert np.all(np.isfinite(R)) -# assert np.all(np.isfinite(Z)) +def test_bad_ion_density(tokamak_args_dict, tokamak_source_example): + # It should fail if given a negative r + with pytest.raises(ValueError) as excinfo: + r=[0, 5, -6] + tokamak_ion_density( + r=r, + mode='L', + ion_density_centre=tokamak_args_dict['ion_density_centre'], + ion_density_peaking_factor=tokamak_args_dict['ion_density_peaking_factor'], + ion_density_pedestal=tokamak_args_dict['ion_density_pedestal'], + major_radius=tokamak_args_dict['major_radius'], + pedestal_radius=tokamak_args_dict['pedestal_radius'], + ion_density_separatrix=tokamak_args_dict['ion_density_separatrix'], + ) + assert "must not be negative" in str(excinfo.value) + + +def test_ion_temperature(tokamak_args_dict, tokamak_source_example): + # test with values of r that are within acceptable ranges. + r = np.linspace(0.0, 2.9, 100) + temperature = tokamak_ion_temperature( + r=r, + mode=tokamak_args_dict['mode'], + pedestal_radius=tokamak_args_dict['pedestal_radius'], + ion_temperature_pedestal=tokamak_args_dict['ion_temperature_pedestal'], + ion_temperature_centre=tokamak_args_dict['ion_temperature_centre'], + ion_temperature_beta=tokamak_args_dict['ion_temperature_beta'], + ion_temperature_peaking_factor=tokamak_args_dict['ion_temperature_peaking_factor'], + ion_temperature_separatrix=tokamak_args_dict['ion_temperature_separatrix'], + major_radius=tokamak_args_dict['major_radius'], + ) + assert isinstance(temperature, np.ndarray) + assert len(temperature) == len(r) + assert np.all(np.isfinite(temperature)) -# def test_bad_convert_a_alpha_to_R_Z(tokamak_source_example): -# # Repeat test_convert_a_alpha_to_R_Z, but show that negative a breaks it -# a = np.linspace(0.0, tokamak_source_example.minor_radius, 100) -# alpha = np.linspace(0.0, 2 * np.pi, 100) -# with pytest.raises(ValueError) as excinfo: -# R, Z = tokamak_source_example.convert_a_alpha_to_R_Z(-a, alpha) -# assert "must not be negative" in str(excinfo.value) +def test_bad_ion_temperature(tokamak_args_dict): + # It should fail if given a negative r + with pytest.raises(ValueError) as excinfo: + r=[0, 5, -6] + tokamak_ion_temperature( + r=r, + mode=tokamak_args_dict['mode'], + pedestal_radius=tokamak_args_dict['pedestal_radius'], + ion_temperature_pedestal=tokamak_args_dict['ion_temperature_pedestal'], + ion_temperature_centre=tokamak_args_dict['ion_temperature_centre'], + ion_temperature_beta=tokamak_args_dict['ion_temperature_beta'], + ion_temperature_peaking_factor=tokamak_args_dict['ion_temperature_peaking_factor'], + ion_temperature_separatrix=tokamak_args_dict['ion_temperature_separatrix'], + major_radius=tokamak_args_dict['major_radius'], + ) + assert "must not be negative" in str(excinfo.value) + + +def test_convert_a_alpha_to_R_Z(tokamak_args_dict): + # Similar to test_source_locations_are_within_correct_range + # Rather than going in detail, simply tests validity of inputs and outputs + # Test with suitable values for a and alpha + a = np.linspace(0.0, 2.9, 100) + alpha = np.linspace(0.0, 2 * np.pi, 100) + R, Z = tokamak_convert_a_alpha_to_R_Z( + a=a, + alpha=alpha, + shafranov_factor=tokamak_args_dict['shafranov_factor'], + minor_radius=tokamak_args_dict['minor_radius'], + major_radius=tokamak_args_dict['major_radius'], + triangularity=tokamak_args_dict['triangularity'], + elongation=tokamak_args_dict['elongation'], + ) + assert isinstance(R, np.ndarray) + assert isinstance(Z, np.ndarray) + assert len(R) == len(a) + assert len(Z) == len(a) + assert np.all(np.isfinite(R)) + assert np.all(np.isfinite(Z)) + + +def test_bad_convert_a_alpha_to_R_Z(tokamak_args_dict): + # Repeat test_convert_a_alpha_to_R_Z, but show that negative a breaks it + a = np.linspace(0.0, 2.9, 100) + alpha = np.linspace(0.0, 2 * np.pi, 100) + with pytest.raises(ValueError) as excinfo: + tokamak_convert_a_alpha_to_R_Z( + a=-a, + alpha=alpha, + shafranov_factor=tokamak_args_dict['shafranov_factor'], + minor_radius=tokamak_args_dict['minor_radius'], + major_radius=tokamak_args_dict['major_radius'], + triangularity=tokamak_args_dict['triangularity'], + elongation=tokamak_args_dict['elongation'], + ) + assert "must not be negative" in str(excinfo.value) @st.composite @@ -303,54 +359,63 @@ def tokamak_source_strategy(draw): ion_temperature_separatrix=0.1, mode="H", ion_temperature_beta=6, - ) + ), { + 'major_radius':major_radius, + 'minor_radius':minor_radius, + 'elongation':elongation, + 'triangularity':triangularity, + } -# @given(tokamak_source=tokamak_source_strategy()) -# @settings(max_examples=50) -# def test_strengths_are_normalised(tokamak_source): -# """Tests that the sum of the strengths attribute is equal to""" -# assert pytest.approx(sum(tokamak_source.strengths)) == 1 - - -# @given(tokamak_source=tokamak_source_strategy()) -# @settings(max_examples=50) -# def test_source_locations_are_within_correct_range(tokamak_source): -# """Tests that each source has RZ locations within the expected range. - -# As the function converting (a,alpha) coordinates to (R,Z) is not bijective, -# we cannot convert back to validate each individual point. However, we can -# determine whether the generated points are contained within the shell of -# the last closed magnetic surface. See "Tokamak D-T neutron source models -# for different plasma physics confinement modes", C. Fausser et al., Fusion -# Engineering and Design, 2012 for more info. -# """ -# R_0 = tokamak_source.major_radius -# A = tokamak_source.minor_radius -# El = tokamak_source.elongation -# delta = tokamak_source.triangularity - -# def get_R_on_LCMS(alpha): -# """Gets R on the last closed magnetic surface for a given alpha""" -# return R_0 + A * np.cos(alpha + delta * np.sin(alpha)) - -# approx_lt = lambda x, y: x < y or np.isclose(x, y) -# approx_gt = lambda x, y: x > y or np.isclose(x, y) - -# for source in tokamak_source.sources: -# R, Z = source.space.r.x[0], source.space.z.x[0] -# # First test that the point is contained with a simple box with -# # lower left (r_min,-z_max) and upper right (r_max,z_max) -# assert approx_gt(R, R_0 - A) -# assert approx_lt(R, R_0 + A) -# assert approx_lt(abs(Z), A * El) -# # For a given Z, we can determine the two values of alpha where -# # where a = minor_radius, and from there determine the upper and -# # lower bounds for R. -# alpha_1 = np.arcsin(abs(Z) / (El * A)) -# alpha_2 = np.pi - alpha_1 -# R_max, R_min = get_R_on_LCMS(alpha_1), get_R_on_LCMS(alpha_2) -# assert approx_lt(R_max, R_0 + A) -# assert approx_gt(R_min, R_0 - A) -# assert approx_lt(R, R_max) -# assert approx_gt(R, R_min) +@given(tokamak_source=tokamak_source_strategy()) +@settings(max_examples=30, suppress_health_check=(HealthCheck.too_slow,)) +def test_strengths_are_normalised(tokamak_source): + """Tests that the sum of the strengths attribute is equal to""" + local_strength = 0 + all_sources = tokamak_source[0] + for source in all_sources: + local_strength = local_strength + source.strength + assert pytest.approx(local_strength) == 1 + + +@given(tokamak_source=tokamak_source_strategy()) +@settings(max_examples=2) +def test_source_locations_are_within_correct_range(tokamak_source): + """Tests that each source has RZ locations within the expected range. + + As the function converting (a,alpha) coordinates to (R,Z) is not bijective, + we cannot convert back to validate each individual point. However, we can + determine whether the generated points are contained within the shell of + the last closed magnetic surface. See "Tokamak D-T neutron source models + for different plasma physics confinement modes", C. Fausser et al., Fusion + Engineering and Design, 2012 for more info. + """ + R_0 = tokamak_source[1]['major_radius'] + A = tokamak_source[1]['minor_radius'] + El = tokamak_source[1]['elongation'] + delta = tokamak_source[1]['triangularity'] + + def get_R_on_LCMS(alpha): + """Gets R on the last closed magnetic surface for a given alpha""" + return R_0 + A * np.cos(alpha + delta * np.sin(alpha)) + + approx_lt = lambda x, y: x < y or np.isclose(x, y) + approx_gt = lambda x, y: x > y or np.isclose(x, y) + + for source in tokamak_source[0]: + R, Z = source.space.r.x[0], source.space.z.x[0] + # First test that the point is contained with a simple box with + # lower left (r_min,-z_max) and upper right (r_max,z_max) + assert approx_gt(R, R_0 - A) + assert approx_lt(R, R_0 + A) + assert approx_lt(abs(Z), A * El) + # For a given Z, we can determine the two values of alpha where + # where a = minor_radius, and from there determine the upper and + # lower bounds for R. + alpha_1 = np.arcsin(abs(Z) / (El * A)) + alpha_2 = np.pi - alpha_1 + R_max, R_min = get_R_on_LCMS(alpha_1), get_R_on_LCMS(alpha_2) + assert approx_lt(R_max, R_0 + A) + assert approx_gt(R_min, R_0 - A) + assert approx_lt(R, R_max) + assert approx_gt(R, R_min) From df2b7defeb7ab15d9145af51b0d912c7b8e3eb38 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 13:59:22 +0000 Subject: [PATCH 31/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/tokamak_source.py | 1 - tests/test_tokamak_source.py | 113 +++++++++++---------- 2 files changed, 61 insertions(+), 53 deletions(-) diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index a0f046d..8b763ae 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -127,7 +127,6 @@ def tokamak_source( if ion_density_centre < 0: raise ValueError("ion_density_centre must greater than 0") - if ion_density_pedestal < 0: raise ValueError("ion_density_pedestal must greater than 0") diff --git a/tests/test_tokamak_source.py b/tests/test_tokamak_source.py index 4b5deed..e969597 100644 --- a/tests/test_tokamak_source.py +++ b/tests/test_tokamak_source.py @@ -3,7 +3,12 @@ from hypothesis import given, settings, HealthCheck from hypothesis import strategies as st import openmc -from openmc_plasma_source import tokamak_source, tokamak_ion_density, tokamak_ion_temperature, tokamak_convert_a_alpha_to_R_Z +from openmc_plasma_source import ( + tokamak_source, + tokamak_ion_density, + tokamak_ion_temperature, + tokamak_convert_a_alpha_to_R_Z, +) @pytest.fixture @@ -185,16 +190,16 @@ def test_bad_angles(tokamak_args_dict, angles): def test_ion_density(tokamak_args_dict, tokamak_source_example): # test with values of r that are within acceptable ranges. - r = np.linspace(0.0, tokamak_args_dict['minor_radius'], 100) + r = np.linspace(0.0, tokamak_args_dict["minor_radius"], 100) density = tokamak_ion_density( r=r, - mode='L', - ion_density_centre=tokamak_args_dict['ion_density_centre'], - ion_density_peaking_factor=tokamak_args_dict['ion_density_peaking_factor'], - ion_density_pedestal=tokamak_args_dict['ion_density_pedestal'], - major_radius=tokamak_args_dict['major_radius'], - pedestal_radius=tokamak_args_dict['pedestal_radius'], - ion_density_separatrix=tokamak_args_dict['ion_density_separatrix'], + mode="L", + ion_density_centre=tokamak_args_dict["ion_density_centre"], + ion_density_peaking_factor=tokamak_args_dict["ion_density_peaking_factor"], + ion_density_pedestal=tokamak_args_dict["ion_density_pedestal"], + major_radius=tokamak_args_dict["major_radius"], + pedestal_radius=tokamak_args_dict["pedestal_radius"], + ion_density_separatrix=tokamak_args_dict["ion_density_separatrix"], ) assert isinstance(r, np.ndarray) assert len(density) == len(r) @@ -204,16 +209,16 @@ def test_ion_density(tokamak_args_dict, tokamak_source_example): def test_bad_ion_density(tokamak_args_dict, tokamak_source_example): # It should fail if given a negative r with pytest.raises(ValueError) as excinfo: - r=[0, 5, -6] + r = [0, 5, -6] tokamak_ion_density( r=r, - mode='L', - ion_density_centre=tokamak_args_dict['ion_density_centre'], - ion_density_peaking_factor=tokamak_args_dict['ion_density_peaking_factor'], - ion_density_pedestal=tokamak_args_dict['ion_density_pedestal'], - major_radius=tokamak_args_dict['major_radius'], - pedestal_radius=tokamak_args_dict['pedestal_radius'], - ion_density_separatrix=tokamak_args_dict['ion_density_separatrix'], + mode="L", + ion_density_centre=tokamak_args_dict["ion_density_centre"], + ion_density_peaking_factor=tokamak_args_dict["ion_density_peaking_factor"], + ion_density_pedestal=tokamak_args_dict["ion_density_pedestal"], + major_radius=tokamak_args_dict["major_radius"], + pedestal_radius=tokamak_args_dict["pedestal_radius"], + ion_density_separatrix=tokamak_args_dict["ion_density_separatrix"], ) assert "must not be negative" in str(excinfo.value) @@ -223,14 +228,16 @@ def test_ion_temperature(tokamak_args_dict, tokamak_source_example): r = np.linspace(0.0, 2.9, 100) temperature = tokamak_ion_temperature( r=r, - mode=tokamak_args_dict['mode'], - pedestal_radius=tokamak_args_dict['pedestal_radius'], - ion_temperature_pedestal=tokamak_args_dict['ion_temperature_pedestal'], - ion_temperature_centre=tokamak_args_dict['ion_temperature_centre'], - ion_temperature_beta=tokamak_args_dict['ion_temperature_beta'], - ion_temperature_peaking_factor=tokamak_args_dict['ion_temperature_peaking_factor'], - ion_temperature_separatrix=tokamak_args_dict['ion_temperature_separatrix'], - major_radius=tokamak_args_dict['major_radius'], + mode=tokamak_args_dict["mode"], + pedestal_radius=tokamak_args_dict["pedestal_radius"], + ion_temperature_pedestal=tokamak_args_dict["ion_temperature_pedestal"], + ion_temperature_centre=tokamak_args_dict["ion_temperature_centre"], + ion_temperature_beta=tokamak_args_dict["ion_temperature_beta"], + ion_temperature_peaking_factor=tokamak_args_dict[ + "ion_temperature_peaking_factor" + ], + ion_temperature_separatrix=tokamak_args_dict["ion_temperature_separatrix"], + major_radius=tokamak_args_dict["major_radius"], ) assert isinstance(temperature, np.ndarray) assert len(temperature) == len(r) @@ -240,17 +247,19 @@ def test_ion_temperature(tokamak_args_dict, tokamak_source_example): def test_bad_ion_temperature(tokamak_args_dict): # It should fail if given a negative r with pytest.raises(ValueError) as excinfo: - r=[0, 5, -6] + r = [0, 5, -6] tokamak_ion_temperature( r=r, - mode=tokamak_args_dict['mode'], - pedestal_radius=tokamak_args_dict['pedestal_radius'], - ion_temperature_pedestal=tokamak_args_dict['ion_temperature_pedestal'], - ion_temperature_centre=tokamak_args_dict['ion_temperature_centre'], - ion_temperature_beta=tokamak_args_dict['ion_temperature_beta'], - ion_temperature_peaking_factor=tokamak_args_dict['ion_temperature_peaking_factor'], - ion_temperature_separatrix=tokamak_args_dict['ion_temperature_separatrix'], - major_radius=tokamak_args_dict['major_radius'], + mode=tokamak_args_dict["mode"], + pedestal_radius=tokamak_args_dict["pedestal_radius"], + ion_temperature_pedestal=tokamak_args_dict["ion_temperature_pedestal"], + ion_temperature_centre=tokamak_args_dict["ion_temperature_centre"], + ion_temperature_beta=tokamak_args_dict["ion_temperature_beta"], + ion_temperature_peaking_factor=tokamak_args_dict[ + "ion_temperature_peaking_factor" + ], + ion_temperature_separatrix=tokamak_args_dict["ion_temperature_separatrix"], + major_radius=tokamak_args_dict["major_radius"], ) assert "must not be negative" in str(excinfo.value) @@ -264,11 +273,11 @@ def test_convert_a_alpha_to_R_Z(tokamak_args_dict): R, Z = tokamak_convert_a_alpha_to_R_Z( a=a, alpha=alpha, - shafranov_factor=tokamak_args_dict['shafranov_factor'], - minor_radius=tokamak_args_dict['minor_radius'], - major_radius=tokamak_args_dict['major_radius'], - triangularity=tokamak_args_dict['triangularity'], - elongation=tokamak_args_dict['elongation'], + shafranov_factor=tokamak_args_dict["shafranov_factor"], + minor_radius=tokamak_args_dict["minor_radius"], + major_radius=tokamak_args_dict["major_radius"], + triangularity=tokamak_args_dict["triangularity"], + elongation=tokamak_args_dict["elongation"], ) assert isinstance(R, np.ndarray) assert isinstance(Z, np.ndarray) @@ -286,11 +295,11 @@ def test_bad_convert_a_alpha_to_R_Z(tokamak_args_dict): tokamak_convert_a_alpha_to_R_Z( a=-a, alpha=alpha, - shafranov_factor=tokamak_args_dict['shafranov_factor'], - minor_radius=tokamak_args_dict['minor_radius'], - major_radius=tokamak_args_dict['major_radius'], - triangularity=tokamak_args_dict['triangularity'], - elongation=tokamak_args_dict['elongation'], + shafranov_factor=tokamak_args_dict["shafranov_factor"], + minor_radius=tokamak_args_dict["minor_radius"], + major_radius=tokamak_args_dict["major_radius"], + triangularity=tokamak_args_dict["triangularity"], + elongation=tokamak_args_dict["elongation"], ) assert "must not be negative" in str(excinfo.value) @@ -360,10 +369,10 @@ def tokamak_source_strategy(draw): mode="H", ion_temperature_beta=6, ), { - 'major_radius':major_radius, - 'minor_radius':minor_radius, - 'elongation':elongation, - 'triangularity':triangularity, + "major_radius": major_radius, + "minor_radius": minor_radius, + "elongation": elongation, + "triangularity": triangularity, } @@ -390,10 +399,10 @@ def test_source_locations_are_within_correct_range(tokamak_source): for different plasma physics confinement modes", C. Fausser et al., Fusion Engineering and Design, 2012 for more info. """ - R_0 = tokamak_source[1]['major_radius'] - A = tokamak_source[1]['minor_radius'] - El = tokamak_source[1]['elongation'] - delta = tokamak_source[1]['triangularity'] + R_0 = tokamak_source[1]["major_radius"] + A = tokamak_source[1]["minor_radius"] + El = tokamak_source[1]["elongation"] + delta = tokamak_source[1]["triangularity"] def get_R_on_LCMS(alpha): """Gets R on the last closed magnetic surface for a given alpha""" From 6768bd7171d1b36235bc116d1279469d3d4eef90 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 8 Mar 2024 14:02:29 +0000 Subject: [PATCH 32/44] tests are slow but passing --- tests/test_tokamak_source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_tokamak_source.py b/tests/test_tokamak_source.py index 4b5deed..985ea8d 100644 --- a/tests/test_tokamak_source.py +++ b/tests/test_tokamak_source.py @@ -379,7 +379,7 @@ def test_strengths_are_normalised(tokamak_source): @given(tokamak_source=tokamak_source_strategy()) -@settings(max_examples=2) +@settings(max_examples=50, suppress_health_check=(HealthCheck.too_slow,)) def test_source_locations_are_within_correct_range(tokamak_source): """Tests that each source has RZ locations within the expected range. From da849c894445447bfcd3620193bc10e8b0ad5306 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 8 Mar 2024 14:19:08 +0000 Subject: [PATCH 33/44] commented out plotly part of examples --- .github/workflows/ci.yml | 3 +++ examples/point_source_example.py | 25 ++++++++++++++++--------- examples/ring_source_example.py | 14 +++++++++++++- examples/tokamak_source_example.py | 21 +++++++++++++++++---- 4 files changed, 49 insertions(+), 14 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a475db9..89b2bec 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -31,3 +31,6 @@ jobs: python examples/point_source_example.py python examples/ring_source_example.py python examples/tokamak_source_example.py + python examples/plot_tokamak_ion_temperature.py + python examples/plot_tokamak_neutron_source_density.py + python examples/plot_tokamak_neutron_source_strengths.py diff --git a/examples/point_source_example.py b/examples/point_source_example.py index d37d404..f9c83b5 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -2,8 +2,6 @@ import numpy as np import openmc -from openmc_source_plotter import plot_source_energy - from openmc_plasma_source import fusion_point_source # just making use of a local cross section xml file, replace with your own cross sections or comment out @@ -28,11 +26,20 @@ settings.particles = 1000 settings.source = my_source +model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) -plot = plot_source_energy( - this=settings, - n_samples=100000, - energy_bins=np.linspace(0, 16e6, 1000), - yaxis_type="log", -) -plot.show() +model.run() + + +# optionally if you would like to plot the energy of particles then another package can be used +# https://github.com/fusion-energy/openmc_source_plotter + +# from openmc_source_plotter import plot_source_energy + +# plot = plot_source_energy( +# this=settings, +# n_samples=1000000, # increase this value for a smoother plot +# energy_bins=np.linspace(0, 16e6, 1000), +# yaxis_type="log", +# ) +# plot.show() diff --git a/examples/ring_source_example.py b/examples/ring_source_example.py index 18efd13..2c959dc 100644 --- a/examples/ring_source_example.py +++ b/examples/ring_source_example.py @@ -2,7 +2,6 @@ from pathlib import Path import openmc - from openmc_plasma_source import fusion_ring_source # just making use of a local cross section xml file, replace with your own cross sections or comment out @@ -31,3 +30,16 @@ model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) model.run() + + +# optionally if you would like to plot the location of particles then another package can be used +# https://github.com/fusion-energy/openmc_source_plotter + +# from openmc_source_plotter import plot_source_position + +# plot = plot_source_position( +# this=settings, +# n_samples = 2000, +# ) + +# plot.show() diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index d8c9666..8f35a29 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -9,7 +9,7 @@ # minimal geometry -sphere_surface = openmc.Sphere(r=1000.0, boundary_type="vacuum") +sphere_surface = openmc.Sphere(r=100000.0, boundary_type="vacuum") cell = openmc.Cell(region=-sphere_surface) geometry = openmc.Geometry([cell]) @@ -24,9 +24,9 @@ ion_temperature_peaking_factor=8.06, ion_temperature_pedestal=6.09, ion_temperature_separatrix=0.1, - major_radius=9.06, - minor_radius=2.92258, - pedestal_radius=0.8 * 2.92258, + major_radius=906, + minor_radius=292.258, + pedestal_radius=0.8 * 292.258, mode="H", shafranov_factor=0.44789, triangularity=0.270, @@ -45,3 +45,16 @@ model = openmc.model.Model(materials=None, geometry=geometry, settings=settings) model.run() + + +# optionally if you would like to plot the direction of particles then another package can be used +# https://github.com/fusion-energy/openmc_source_plotter + +# from openmc_source_plotter import plot_source_direction + +# plot = plot_source_direction( +# this=settings, +# n_samples=200 +# ) + +# plot.show() From 08779d4a3dd62ffc82edad71113b86b32df65601 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 23:11:56 +0000 Subject: [PATCH 34/44] checking for incorrect values is more concise --- pyproject.toml | 2 +- src/openmc_plasma_source/point_source.py | 8 +- src/openmc_plasma_source/ring_source.py | 11 +-- src/openmc_plasma_source/tokamak_source.py | 96 +++++++--------------- tests/test_tokamak_source.py | 14 ++-- 5 files changed, 47 insertions(+), 84 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 7fc3d6f..9c318f1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,7 +18,7 @@ classifiers = [ authors = [ { name="Rémi Delaporte-Mathurin" }, ] -requires-python = ">=3.8" +requires-python = ">=3.9" keywords = ["python", "neutron", "fusion", "source", "openmc", "energy", "tokamak"] dependencies = [ "numpy>=1.9", diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index eefa4ba..9e896ab 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -10,10 +10,10 @@ def fusion_point_source( temperature: float = 20000.0, fuel: dict = {"D": 0.5, "T": 0.5}, ): - """An openmc.Source object with some presets to make it more convenient - for fusion simulations using a point source. All attributes can be changed - after initialization if required. Default isotropic point source at the - origin with a Muir energy distribution. + """Creates a list of openmc.IndependentSource objects representing an ICF source. + + Resulting ICF (Inertial Confinement Fusion) source will have an energy + distribution according to the fuel composition. Args: coordinate (tuple[float,float,float]): Location of the point source. diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index 751e4af..8ced493 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -12,11 +12,12 @@ def fusion_ring_source( z_placement: float = 0, temperature: float = 20000.0, fuel: dict = {"D": 0.5, "T": 0.5}, -): - """An openmc.Source object with some presets to make it more convenient - for fusion simulations using a ring source. All attributes can be changed - after initialization if required. Default isotropic ring source with a - realistic energy distribution. +)-> list[openmc.IndependentSource]: + """Creates a list of openmc.IndependentSource objects in a ring shape. + + Useful for simulations where all the plasma parameters are not known and + this simplified geometry will suffice. Resulting ring source will have an + energy distribution according to the fuel composition. Args: radius (float): the inner radius of the ring source, in metres diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 8b763ae..f24cfd3 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -2,6 +2,7 @@ import numpy as np import openmc +import openmc.checkvalue as cv from .fuel_types import get_neutron_energy_distribution @@ -26,20 +27,20 @@ def tokamak_source( angles: Tuple[float, float] = (0, 2 * np.pi), sample_size: int = 1000, fuel: dict = {"D": 0.5, "T": 0.5}, -): - """Plasma neutron source sampling. - This class greatly relies on models described in [1] +) -> list[openmc.IndependentSource]: + """Creates a list of openmc.IndependentSource objects representing a tokamak plasma. + + Resulting sources will have an energy distribution according to the fuel + composition.This function greatly relies on models described in [1] [1] : Fausser et al, 'Tokamak D-T neutron source models for different plasma physics confinement modes', Fus. Eng. and Design, https://doi.org/10.1016/j.fusengdes.2012.02.025 Usage: - my_plasma = Plasma(**plasma_prms) - my_plasma.sample_sources() - print(my_plasma.RZ) - print(my_plasma.temperatures) - openmc_sources = my_plasma.make_openmc_sources() + my_source = tokamak_source(**plasma_prms) + my_settings = openmc.Settings() + my_settings.source = my_source Args: major_radius (float): Plasma major radius (cm) @@ -73,65 +74,26 @@ def tokamak_source( """ # Perform sanity checks for inputs not caught by properties - - if not isinstance(major_radius, (int, float)): - raise ValueError("Major radius must be a number") - - if not isinstance(minor_radius, (int, float)): - raise ValueError("Minor radius must be a number") - - if not isinstance(elongation, (int, float)): - raise ValueError("Elongation must be a number") - - if not isinstance(triangularity, (int, float)): - raise ValueError("Triangularity must be a number") - - if not isinstance(ion_density_centre, (int, float)): - raise ValueError("ion_density_centre must be a number") - - if not isinstance(ion_density_peaking_factor, (int, float)): - raise ValueError("ion_density_peaking_factor must be a number") - - if not isinstance(ion_density_pedestal, (int, float)): - raise ValueError("ion_density_pedestal must be a number") - - if not isinstance(ion_density_separatrix, (int, float)): - raise ValueError("ion_density_separatrix must be a number") - - if minor_radius >= major_radius: - raise ValueError("Minor radius must be smaller than major radius") - - if pedestal_radius >= minor_radius: - raise ValueError("Pedestal radius must be smaller than minor radius") - - if abs(shafranov_factor) >= 0.5 * minor_radius: - raise ValueError("Shafranov factor must be smaller than 0.5*minor radius") - - if major_radius < 0: - raise ValueError("Major radius must greater than 0") - - if minor_radius < 0: - raise ValueError("Minor radius must greater than 0") - - if elongation <= 0: - raise ValueError("Elongation must greater than 0") - - if not triangularity <= 1.0: - raise ValueError("Triangularity must less than or equal to 1.") - if not triangularity >= -1.0: - raise ValueError("Triangularity must greater than or equal to -1.") - - if mode not in ["H", "L", "A"]: - raise ValueError("Mode must be one of the following: ['H', 'L', 'A']") - - if ion_density_centre < 0: - raise ValueError("ion_density_centre must greater than 0") - - if ion_density_pedestal < 0: - raise ValueError("ion_density_pedestal must greater than 0") - - if ion_density_separatrix < 0: - raise ValueError("ion_density_separatrix must greater than 0") + cv.check_type('major_radius', major_radius, (int, float)) + cv.check_type('minor_radius', minor_radius, (int, float)) + cv.check_type('elongation', elongation, (int, float)) + cv.check_type('triangularity', triangularity, (int, float)) + cv.check_type('ion_density_centre', ion_density_centre, (int, float)) + cv.check_type('ion_density_peaking_factor', ion_density_peaking_factor, (int, float)) + cv.check_type('ion_density_pedestal', ion_density_pedestal, (int, float)) + cv.check_type('ion_density_separatrix', ion_density_separatrix, (int, float)) + cv.check_less_than('minor_radius', minor_radius, major_radius) + cv.check_less_than('pedestal_radius', pedestal_radius, minor_radius) + cv.check_less_than('shafranov_factor', abs(shafranov_factor), 0.5 * minor_radius) + cv.check_greater_than('major_radius', major_radius, 0) + cv.check_greater_than('minor_radius', minor_radius, 0) + cv.check_greater_than('elongation', elongation, 0) + cv.check_less_than('triangularity', triangularity, 1., True) + cv.check_greater_than('triangularity', triangularity, -1., True) + cv.check_value('mode', mode, ["H", "L", "A"]) + cv.check_greater_than('ion_density_centre', ion_density_centre, 0) + cv.check_greater_than('ion_density_pedestal', ion_density_pedestal, 0) + cv.check_greater_than('ion_density_separatrix', ion_density_separatrix, 0) if ( isinstance(angles, tuple) diff --git a/tests/test_tokamak_source.py b/tests/test_tokamak_source.py index a459eed..9c43724 100644 --- a/tests/test_tokamak_source.py +++ b/tests/test_tokamak_source.py @@ -65,7 +65,7 @@ def test_bad_major_radius(tokamak_args_dict, minor_radius, major_radius): """Checks that tokamak_source creation rejects invalid major radius""" tokamak_args_dict["minor_radius"] = minor_radius tokamak_args_dict["major_radius"] = major_radius - with pytest.raises(ValueError): + with pytest.raises((ValueError, TypeError)): tokamak_source(**tokamak_args_dict) @@ -90,7 +90,7 @@ def test_bad_minor_radius(tokamak_args_dict, major_radius, minor_radius): """Checks that tokamak_source creation rejects invalid minor radius""" tokamak_args_dict["major_radius"] = major_radius tokamak_args_dict["minor_radius"] = minor_radius - with pytest.raises(ValueError): + with pytest.raises((ValueError, TypeError)): tokamak_source(**tokamak_args_dict) @@ -105,7 +105,7 @@ def test_elongation(tokamak_args_dict, elongation): def test_bad_elongation(tokamak_args_dict, elongation): """Checks that tokamak_source creation rejects invalid elongation""" tokamak_args_dict["elongation"] = elongation - with pytest.raises(ValueError): + with pytest.raises((ValueError, TypeError)): tokamak_source(**tokamak_args_dict) @@ -120,7 +120,7 @@ def test_triangularity(tokamak_args_dict, triangularity): def test_bad_triangularity(tokamak_args_dict, triangularity): """Checks that tokamak_source creation rejects invalid triangularity""" tokamak_args_dict["triangularity"] = triangularity - with pytest.raises(ValueError): + with pytest.raises((ValueError, TypeError)): tokamak_source(**tokamak_args_dict) @@ -184,11 +184,11 @@ def test_bad_angles(tokamak_args_dict, angles): # It should fail when given something that isn't a 2-tuple or similar # Contents should convert to float tokamak_args_dict["angles"] = angles - with pytest.raises(ValueError) as excinfo: + with pytest.raises((ValueError, TypeError)): tokamak_source(**tokamak_args_dict) -def test_ion_density(tokamak_args_dict, tokamak_source_example): +def test_ion_density(tokamak_args_dict): # test with values of r that are within acceptable ranges. r = np.linspace(0.0, tokamak_args_dict["minor_radius"], 100) density = tokamak_ion_density( @@ -206,7 +206,7 @@ def test_ion_density(tokamak_args_dict, tokamak_source_example): assert np.all(np.isfinite(density)) -def test_bad_ion_density(tokamak_args_dict, tokamak_source_example): +def test_bad_ion_density(tokamak_args_dict): # It should fail if given a negative r with pytest.raises(ValueError) as excinfo: r = [0, 5, -6] From 915337c241f68868dab61588ff2bc230b5bd912f Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 8 Mar 2024 23:12:14 +0000 Subject: [PATCH 35/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/ring_source.py | 2 +- src/openmc_plasma_source/tokamak_source.py | 42 +++++++++++----------- 2 files changed, 23 insertions(+), 21 deletions(-) diff --git a/src/openmc_plasma_source/ring_source.py b/src/openmc_plasma_source/ring_source.py index 8ced493..698d8c5 100644 --- a/src/openmc_plasma_source/ring_source.py +++ b/src/openmc_plasma_source/ring_source.py @@ -12,7 +12,7 @@ def fusion_ring_source( z_placement: float = 0, temperature: float = 20000.0, fuel: dict = {"D": 0.5, "T": 0.5}, -)-> list[openmc.IndependentSource]: +) -> list[openmc.IndependentSource]: """Creates a list of openmc.IndependentSource objects in a ring shape. Useful for simulations where all the plasma parameters are not known and diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index f24cfd3..cbf3e9d 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -74,26 +74,28 @@ def tokamak_source( """ # Perform sanity checks for inputs not caught by properties - cv.check_type('major_radius', major_radius, (int, float)) - cv.check_type('minor_radius', minor_radius, (int, float)) - cv.check_type('elongation', elongation, (int, float)) - cv.check_type('triangularity', triangularity, (int, float)) - cv.check_type('ion_density_centre', ion_density_centre, (int, float)) - cv.check_type('ion_density_peaking_factor', ion_density_peaking_factor, (int, float)) - cv.check_type('ion_density_pedestal', ion_density_pedestal, (int, float)) - cv.check_type('ion_density_separatrix', ion_density_separatrix, (int, float)) - cv.check_less_than('minor_radius', minor_radius, major_radius) - cv.check_less_than('pedestal_radius', pedestal_radius, minor_radius) - cv.check_less_than('shafranov_factor', abs(shafranov_factor), 0.5 * minor_radius) - cv.check_greater_than('major_radius', major_radius, 0) - cv.check_greater_than('minor_radius', minor_radius, 0) - cv.check_greater_than('elongation', elongation, 0) - cv.check_less_than('triangularity', triangularity, 1., True) - cv.check_greater_than('triangularity', triangularity, -1., True) - cv.check_value('mode', mode, ["H", "L", "A"]) - cv.check_greater_than('ion_density_centre', ion_density_centre, 0) - cv.check_greater_than('ion_density_pedestal', ion_density_pedestal, 0) - cv.check_greater_than('ion_density_separatrix', ion_density_separatrix, 0) + cv.check_type("major_radius", major_radius, (int, float)) + cv.check_type("minor_radius", minor_radius, (int, float)) + cv.check_type("elongation", elongation, (int, float)) + cv.check_type("triangularity", triangularity, (int, float)) + cv.check_type("ion_density_centre", ion_density_centre, (int, float)) + cv.check_type( + "ion_density_peaking_factor", ion_density_peaking_factor, (int, float) + ) + cv.check_type("ion_density_pedestal", ion_density_pedestal, (int, float)) + cv.check_type("ion_density_separatrix", ion_density_separatrix, (int, float)) + cv.check_less_than("minor_radius", minor_radius, major_radius) + cv.check_less_than("pedestal_radius", pedestal_radius, minor_radius) + cv.check_less_than("shafranov_factor", abs(shafranov_factor), 0.5 * minor_radius) + cv.check_greater_than("major_radius", major_radius, 0) + cv.check_greater_than("minor_radius", minor_radius, 0) + cv.check_greater_than("elongation", elongation, 0) + cv.check_less_than("triangularity", triangularity, 1.0, True) + cv.check_greater_than("triangularity", triangularity, -1.0, True) + cv.check_value("mode", mode, ["H", "L", "A"]) + cv.check_greater_than("ion_density_centre", ion_density_centre, 0) + cv.check_greater_than("ion_density_pedestal", ion_density_pedestal, 0) + cv.check_greater_than("ion_density_separatrix", ion_density_separatrix, 0) if ( isinstance(angles, tuple) From 22d89e9659a9a0b71e9f93bd23cde363a73b38fa Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 22 Mar 2024 18:02:46 +0000 Subject: [PATCH 36/44] added neutron mean and std functions; --- src/openmc_plasma_source/fuel_types.py | 108 ++++++++++++++++++++++--- 1 file changed, 97 insertions(+), 11 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 4670e5f..3a20e7a 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -3,6 +3,96 @@ import openmc +def neutron_energy_mean(ion_temperature: float, reaction: str) -> float: + """Calculates the mean energy of the neutron emitted during DD or DT + fusion accounting for temperature of the incident ions. Based on Ballabio + fits, see Table III of L. Ballabio et al 1998 Nucl. Fusion 38 1723 + + Args: + ion_temperature (float): the temperature of the ions in eV + reaction (str): the two isotope that fuse, can be either 'DD' or 'DT' + + Raises: + ValueError: if the reaction is not 'DD' or 'DT' then a ValueError is raised + + Returns: + float: the mean neutron energy in eV + """ + + # values from Ballabio paper + if reaction == "DD": + a_1 = 4.69515 + a_2 = -0.040729 + a_3 = 0.47 + a_4 = 0.81844 + mean = 2.4495e6 # in units of eV + elif reaction == "DT": + a_1 = 5.30509 + a_2 = 2.4736e-3 + a_3 = 1.84 + a_4 = 1.3818 + mean = 14.021e6 # in units of eV + else: + raise ValueError(f'reaction must be either "DD" or "DT" not {reaction}') + + ion_temperature_kev = ion_temperature / 1e3 # Ballabio equation accepts KeV units + mean_delta = ( + a_1 + * ion_temperature_kev ** (2./3.) + / (1. + a_2 * ion_temperature_kev**a_3) + + a_4 * ion_temperature_kev + ) + mean_delta *= 1e3 # converting back to eV + return mean + mean_delta + + +def neutron_energy_std_dev(ion_temperature: float, reaction: str) -> float: + """Calculates the standard deviation of the neutron energy emitted during DD + or DT fusion accounting for temperature of the incident ions. Based on + Ballabio fits, see Table III of L. Ballabio et al 1998 Nucl. Fusion 38 1723 + + Args: + ion_temperature (float): the temperature of the ions in eV + reaction (str): the two isotope that fuse, can be either 'DD' or 'DT' + + Raises: + ValueError: if the reaction is not 'DD' or 'DT' then a ValueError is raised + + Returns: + float: the mean neutron energy in eV + """ + + # values from Ballabio paper + if reaction == "DD": + w_0 = 82.542 + a_1 = 1.7013e-3 + a_2 = 0.16888 + a_3 = 0.49 + a_4 = 7.9460e-4 + elif reaction == "DT": + w_0 = 177.259 + a_1 = 5.1068e-4 + a_2 = 7.6223e-3 + a_3 = 1.78 + a_4 = 8.7691e-5 + else: + raise ValueError(f'reaction must be either "DD" or "DT" not {reaction}') + + ion_temperature_kev = ion_temperature / 1e3 # Ballabio equation accepts KeV units + delta = ( + a_1 + * ion_temperature_kev ** (2./3.) + / (1. + a_2 * ion_temperature_kev**a_3) + + a_4 * ion_temperature_kev + ) + + # 2.3548200450309493 on the line below comes from equation 2* math.sqrt(math.log(2)*2) + variance = ((w_0 * (1 + delta))**2 * ion_temperature_kev) / 2.3548200450309493**2 + variance *= 1e6 # converting keV^2 back to eV^2 + std_dev = np.sqrt(variance) + return std_dev + + def get_neutron_energy_distribution( ion_temperature: float, fuel: dict, @@ -53,8 +143,10 @@ def get_neutron_energy_distribution( # single grid for TT neutrons E_pspec = np.linspace(0, 12, num_of_vals) # E_pspec is exspected in MeV units - DTmean, DTvar = nst.DTprimspecmoments(ion_temperature_kev) - DDmean, DDvar = nst.DDprimspecmoments(ion_temperature_kev) + DDmean = neutron_energy_mean(ion_temperature=ion_temperature, reaction="DD") + DTmean = neutron_energy_mean(ion_temperature=ion_temperature, reaction="DT") + DD_std_dev = neutron_energy_std_dev(ion_temperature=ion_temperature, reaction="DD") + DT_std_dev = neutron_energy_std_dev(ion_temperature=ion_temperature, reaction="DT") if ["D", "T"] == sorted(set(fuel.keys())): strength_DT = 1.0 @@ -65,22 +157,16 @@ def get_neutron_energy_distribution( "tt", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"] ) - total_strength = sum([strength_DT, strength_DD, strength_TT]) + total_strength = sum([strength_TT, strength_DD, strength_DT]) dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev) tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT) - DD_std_dev = np.sqrt( - DDvar * 1e12 - ) # power 12 as this is in MeV^2 and we need eV - dd_source = openmc.stats.Normal(mean_value=DDmean * 1e6, std_dev=DD_std_dev) + dd_source = openmc.stats.Normal(mean_value=DDmean, std_dev=DD_std_dev) # normal could be done with Muir but in this case we have the mean and std dev from NeSST # dd_source = openmc.stats.muir(e0=DDmean * 1e6, m_rat=4, kt=ion_temperature) - DT_std_dev = np.sqrt( - DTvar * 1e12 - ) # power 12 as this is in MeV^2 and we need eV - dt_source = openmc.stats.Normal(mean_value=DTmean * 1e6, std_dev=DT_std_dev) + dt_source = openmc.stats.Normal(mean_value=DTmean, std_dev=DT_std_dev) # normal could be done with Muir but in this case we have the mean and std dev from NeSST # dt_source = openmc.stats.muir(e0=DTmean * 1e6, m_rat=5, kt=ion_temperature) From 6e5aa3b1a3e389cbf5ccc7076a2fe8d974017357 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 22 Mar 2024 18:03:06 +0000 Subject: [PATCH 37/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/fuel_types.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index 3a20e7a..3612e0c 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -38,8 +38,8 @@ def neutron_energy_mean(ion_temperature: float, reaction: str) -> float: ion_temperature_kev = ion_temperature / 1e3 # Ballabio equation accepts KeV units mean_delta = ( a_1 - * ion_temperature_kev ** (2./3.) - / (1. + a_2 * ion_temperature_kev**a_3) + * ion_temperature_kev ** (2.0 / 3.0) + / (1.0 + a_2 * ion_temperature_kev**a_3) + a_4 * ion_temperature_kev ) mean_delta *= 1e3 # converting back to eV @@ -81,13 +81,13 @@ def neutron_energy_std_dev(ion_temperature: float, reaction: str) -> float: ion_temperature_kev = ion_temperature / 1e3 # Ballabio equation accepts KeV units delta = ( a_1 - * ion_temperature_kev ** (2./3.) - / (1. + a_2 * ion_temperature_kev**a_3) + * ion_temperature_kev ** (2.0 / 3.0) + / (1.0 + a_2 * ion_temperature_kev**a_3) + a_4 * ion_temperature_kev ) # 2.3548200450309493 on the line below comes from equation 2* math.sqrt(math.log(2)*2) - variance = ((w_0 * (1 + delta))**2 * ion_temperature_kev) / 2.3548200450309493**2 + variance = ((w_0 * (1 + delta)) ** 2 * ion_temperature_kev) / 2.3548200450309493**2 variance *= 1e6 # converting keV^2 back to eV^2 std_dev = np.sqrt(variance) return std_dev From bbdb65c902aca33c273f768c0d4f4ca17c6e103e Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Sat, 30 Mar 2024 22:33:22 +0000 Subject: [PATCH 38/44] removed commented code --- src/openmc_plasma_source/point_source.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/openmc_plasma_source/point_source.py b/src/openmc_plasma_source/point_source.py index 9e896ab..3716f62 100644 --- a/src/openmc_plasma_source/point_source.py +++ b/src/openmc_plasma_source/point_source.py @@ -42,14 +42,6 @@ def fusion_point_source( ion_temperature=temperature, fuel=fuel ) - # if isinstance(energy_distributions, openmc.stats.Normal) or isinstance(energy_distributions, openmc.stats.Discrete) or isinstance(energy_distributions, openmc.stats.Tabular): - # source = openmc.Source() - # source.energy = energy_distributions - # source.space = openmc.stats.Point(coordinate) - # source.angle = openmc.stats.Isotropic() - # return source - - # else: for energy_distribution, strength in zip(energy_distributions, strengths): source = openmc.IndependentSource() source.energy = energy_distribution From 184aebca8726b9f2efb9e0be20ff01e6b8afb7c9 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Sat, 30 Mar 2024 23:07:28 +0000 Subject: [PATCH 39/44] moving to eV for nesst --- pyproject.toml | 2 +- src/openmc_plasma_source/fuel_types.py | 10 ++-- src/openmc_plasma_source/point_source.py | 2 +- src/openmc_plasma_source/tokamak_source.py | 55 +++++++++++++--------- 4 files changed, 39 insertions(+), 30 deletions(-) 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 From a308e6d8af0c9b704d3d42f1b8a274edd6d416f0 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Sat, 30 Mar 2024 23:20:42 +0000 Subject: [PATCH 40/44] planning isotope density --- src/openmc_plasma_source/tokamak_source.py | 35 ++++++++++++---------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 097494b..c869574 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -157,19 +157,23 @@ def tokamak_source( #TODO loop through the fuel reactions # compute neutron source densities - neutron_source_density = tokamak_neutron_source_density(densities, temperatures) + # loop through each reaction in the fuel + # could introduce a species density - strengths = neutron_source_density / sum(neutron_source_density) + for + neutron_source_density = tokamak_neutron_source_density(densities, temperatures, reaction=, fuel) + strengths = neutron_source_density / sum(neutron_source_density) - sources = tokamak_make_openmc_sources( - strengths=strengths, - angles=angles, - temperatures=temperatures, - fuel=fuel, - RZ=RZ, - ) - return sources + + sources = tokamak_make_openmc_sources( + strengths=strengths, + angles=angles, + temperatures=temperatures, + fuel=fuel, + RZ=RZ, + ) + return sources def tokamak_ion_density( @@ -383,12 +387,13 @@ def tokamak_neutron_source_density(ion_density, ion_temperature, reaction='DT'): ion_temperature = np.asarray(ion_temperature) return ion_density**2 * _DT_xs(ion_temperature) #TODO account for reaction + #TODO get a fuel fraction from fuel va if reaction == 'DT': - nst.reac_DT(ion_temperature) - if reaction == 'DD' - nst.reac_DD(ion_temperature) - if reaction == 'TT' - nst.reac_TT(ion_temperature) + ion_density**2 * nst.reac_DT(ion_temperature) # could use _DT_xs instead + if reaction == 'DD': + ion_density**2 * nst.reac_DD(ion_temperature) + if reaction == 'TT': + ion_density**2 * nst.reac_TT(ion_temperature) #TODO consider replace with NeSST or getting DD version as well From e219420621be50fde1362d20ba7e5a8479f2c8b8 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Apr 2024 12:52:40 +0100 Subject: [PATCH 41/44] adding nesst back --- pyproject.toml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 9c318f1..cbecd43 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" + ] [project.urls] @@ -36,7 +36,8 @@ write_to = "src/_version.py" [project.optional-dependencies] tests = [ "pytest>=5.4.3", - "hypothesis" + "hypothesis", + "NeSST" ] [tool.setuptools] From d1695efc3753e39e779d05fc307c19f4abbeddd8 Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 2 Apr 2024 14:06:19 +0000 Subject: [PATCH 42/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/fuel_types.py | 14 ++++++------- src/openmc_plasma_source/tokamak_source.py | 23 ++++++++++++---------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index e270446..1ca478a 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,7 +160,7 @@ 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) @@ -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) diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 338b121..033a2b8 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -132,7 +132,7 @@ 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( @@ -158,11 +158,12 @@ def tokamak_source( elongation=elongation, ) - neutron_source_density = tokamak_neutron_source_density(fuel_densities, temperatures, reaction) + neutron_source_density = tokamak_neutron_source_density( + fuel_densities, temperatures, reaction + ) strengths = neutron_source_density / sum(neutron_source_density) - sources = tokamak_make_openmc_sources( strengths=strengths, angles=angles, @@ -380,18 +381,20 @@ 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 == 'DT': - return ion_density * nst.reac_DT(ion_temperature) # could use _DT_xs instead - if reaction == 'DD': + if reaction == "DT": + return ion_density * nst.reac_DT(ion_temperature) # could use _DT_xs instead + if reaction == "DD": return ion_density * nst.reac_DD(ion_temperature) - if reaction == 'TT': + if reaction == "TT": return ion_density * nst.reac_TT(ion_temperature) -#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 @@ -400,7 +403,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, From 4c906fc6a74a8b08a31e87e7c92cfadc1a923d96 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Apr 2024 17:52:04 +0100 Subject: [PATCH 43/44] RZ values are nested? --- examples/point_source_example.py | 18 ++-- examples/tokamak_source_example.py | 22 ++--- src/openmc_plasma_source/fuel_types.py | 14 +-- src/openmc_plasma_source/tokamak_source.py | 104 +++++++++++---------- 4 files changed, 83 insertions(+), 75 deletions(-) diff --git a/examples/point_source_example.py b/examples/point_source_example.py index f9c83b5..46d15c2 100644 --- a/examples/point_source_example.py +++ b/examples/point_source_example.py @@ -34,12 +34,12 @@ # optionally if you would like to plot the energy of particles then another package can be used # https://github.com/fusion-energy/openmc_source_plotter -# from openmc_source_plotter import plot_source_energy - -# plot = plot_source_energy( -# this=settings, -# n_samples=1000000, # increase this value for a smoother plot -# energy_bins=np.linspace(0, 16e6, 1000), -# yaxis_type="log", -# ) -# plot.show() +from openmc_source_plotter import plot_source_energy + +plot = plot_source_energy( + this=settings, + n_samples=1000000, # increase this value for a smoother plot + energy_bins=np.linspace(0, 16e6, 1000), + yaxis_type="log", +) +plot.show() diff --git a/examples/tokamak_source_example.py b/examples/tokamak_source_example.py index 8f35a29..4711476 100644 --- a/examples/tokamak_source_example.py +++ b/examples/tokamak_source_example.py @@ -17,20 +17,20 @@ my_sources = tokamak_source( elongation=1.557, ion_density_centre=1.09e20, - ion_density_peaking_factor=1, ion_density_pedestal=1.09e20, + ion_density_peaking_factor=1, ion_density_separatrix=3e19, - ion_temperature_centre=45.9, + ion_temperature_centre=45.9e3, + ion_temperature_pedestal=6.09e3, + ion_temperature_separatrix=0.1e3, ion_temperature_peaking_factor=8.06, - ion_temperature_pedestal=6.09, - ion_temperature_separatrix=0.1, + ion_temperature_beta=6, major_radius=906, minor_radius=292.258, pedestal_radius=0.8 * 292.258, mode="H", shafranov_factor=0.44789, triangularity=0.270, - ion_temperature_beta=6, fuel={"D": 0.5, "T": 0.5}, ) @@ -50,11 +50,11 @@ # optionally if you would like to plot the direction of particles then another package can be used # https://github.com/fusion-energy/openmc_source_plotter -# from openmc_source_plotter import plot_source_direction +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() +plot.show() diff --git a/src/openmc_plasma_source/fuel_types.py b/src/openmc_plasma_source/fuel_types.py index e270446..ef3ed03 100644 --- a/src/openmc_plasma_source/fuel_types.py +++ b/src/openmc_plasma_source/fuel_types.py @@ -164,12 +164,12 @@ def get_neutron_energy_distribution( 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_source], [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_source], [strength_DD] + return {'DD': dd_source}, {'DD' :strength_DD} # DT, DD and TT reaction else: @@ -201,8 +201,8 @@ def get_neutron_energy_distribution( # dt_source = openmc.stats.muir(e0=DTmean * 1e6, m_rat=5, kt=ion_temperature) # todo look into combining distributions openmc.data.combine_distributions() - return [tt_source, dd_source, dt_source], [ - strength_TT / total_strength, - strength_DD / total_strength, - strength_DT / total_strength, - ] + return { + '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 338b121..3b30b94 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -4,8 +4,8 @@ import openmc import openmc.checkvalue as cv import NeSST as nst -from .fuel_types import get_neutron_energy_distribution - +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, @@ -54,14 +54,14 @@ def tokamak_source( ion_density_pedestal (float): Ion density at pedestal (m-3) ion_density_separatrix (float): Ion density at separatrix (m-3) ion_temperature_centre (float): Ion temperature at the plasma - centre (keV) + centre (eV) ion_temperature_peaking_factor (float): Ion temperature peaking factor (referred in [1] as ion temperature exponent alpha_T) ion_temperature_beta (float): Ion temperature beta exponent (referred in [1] as ion temperature exponent beta_T) - ion_temperature_pedestal (float): Ion temperature at pedestal (keV) + ion_temperature_pedestal (float): Ion temperature at pedestal (eV) ion_temperature_separatrix (float): Ion temperature at separatrix - (keV) + (eV) pedestal_radius (float): Minor radius at pedestal (cm) shafranov_factor (float): Shafranov factor (referred in [1] as esh) also known as outward radial displacement of magnetic surfaces @@ -139,11 +139,11 @@ def tokamak_source( r=a, mode=mode, pedestal_radius=pedestal_radius, - ion_temperature_pedestal=ion_temperature_pedestal, - ion_temperature_centre=ion_temperature_centre, + 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, + ion_temperature_separatrix=ion_temperature_separatrix/1e3, major_radius=major_radius, ) @@ -158,19 +158,27 @@ def tokamak_source( elongation=elongation, ) - neutron_source_density = tokamak_neutron_source_density(fuel_densities, temperatures, reaction) - - strengths = neutron_source_density / sum(neutron_source_density) - - - sources = tokamak_make_openmc_sources( - strengths=strengths, - angles=angles, - temperatures=temperatures, - fuel=fuel, - RZ=RZ, - ) - return sources + reactions = get_reactions_from_fuel(fuel) + + neutron_source_density={} + total_source_density = 0 + for reaction in reactions: + 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 + + sources = tokamak_make_openmc_sources( + strengths=neutron_source_density, + angles=angles, + temperatures=temperatures, + fuel=fuel, + RZ=RZ, + ) + all_sources = all_sources + sources + return all_sources def tokamak_ion_density( @@ -333,38 +341,37 @@ def tokamak_make_openmc_sources( sources = [] # create a ring source for each sample in the plasma source - for i in range(len(strengths)): + for i, (RZ_val, temperature) in enumerate(zip(RZ, temperatures)): # extract the RZ values accordingly - radius = openmc.stats.Discrete([RZ[0][i]], [1]) - z_values = openmc.stats.Discrete([RZ[1][i]], [1]) + radius = openmc.stats.Discrete([RZ_val[0]], [1]) + z_values = openmc.stats.Discrete([RZ_val[1]], [1]) angle = openmc.stats.Uniform(a=angles[0], b=angles[1]) - strength = strengths[i] - energy_distributions, dist_strengths = get_neutron_energy_distribution( - ion_temperature=temperatures[i], + energy_distributions_and_dist_strengths = get_neutron_energy_distribution( + ion_temperature=temperature, fuel=fuel, ) # now we have potentially 3 distributions (DT, DD, TT) - for energy_distribution, dist_strength in zip( - energy_distributions, dist_strengths - ): - my_source = openmc.IndependentSource() + for reaction, (energy_distribution, dist_strength) in energy_distributions_and_dist_strengths.items(): - # create a ring source - my_source.space = openmc.stats.CylindricalIndependent( - r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0) - ) - my_source.angle = openmc.stats.Isotropic() + if dist_strength * strengths[reaction][i] > 0.: + my_source = openmc.IndependentSource() - my_source.energy = energy_distribution + # create a ring source + my_source.space = openmc.stats.CylindricalIndependent( + r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0) + ) + my_source.angle = openmc.stats.Isotropic() - # the strength of the source (its probability) is given by the - # strength of the energy distribution and the location distribution - my_source.strength = dist_strength * strength + my_source.energy = energy_distribution - # append to the list of sources - sources.append(my_source) + # the strength of the source (its probability) is given by the + # strength of the energy distribution and the location distribution + my_source.strength = dist_strength * strengths[reaction][i] + + # append to the list of sources + sources.append(my_source) return sources @@ -383,12 +390,13 @@ def tokamak_neutron_source_density(ion_density, ion_temperature, reaction): ion_density = np.asarray(ion_density[reaction[0]]) * np.asarray(ion_density[reaction[1]]) ion_temperature = np.asarray(ion_temperature) - if reaction == 'DT': - return ion_density * nst.reac_DT(ion_temperature) # could use _DT_xs instead - if reaction == 'DD': - return ion_density * nst.reac_DD(ion_temperature) - if reaction == 'TT': - return ion_density * nst.reac_TT(ion_temperature) + if reaction == ['DD']: + return ion_density * reac_DD(ion_temperature) + 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 #TODO consider replace with NeSST or getting DD version as well From e860a74affcef7226f603e1013cb58a6b20faca4 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 19 Apr 2024 12:02:58 +0000 Subject: [PATCH 44/44] [skip ci] Apply formatting changes --- src/openmc_plasma_source/tokamak_source.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/openmc_plasma_source/tokamak_source.py b/src/openmc_plasma_source/tokamak_source.py index 9dfa843..d751918 100644 --- a/src/openmc_plasma_source/tokamak_source.py +++ b/src/openmc_plasma_source/tokamak_source.py @@ -123,7 +123,7 @@ def tokamak_source( rng = np.random.default_rng(self.sample_seed) a = rng.random(self.sample_size) * self.minor_radius alpha = rng.random(self.sample_size) * 2 * np.pi - + # compute densities, temperatures densities = tokamak_ion_density( mode=mode, @@ -282,7 +282,7 @@ def tokamak_ion_temperature( + (ion_temperature_pedestal - ion_temperature_separatrix) * (major_radius - r) / (major_radius - pedestal_radius) - ) + ), ) return temperature