Skip to content

openmc_interface

OpenMC neutron source interface.

AllReactions = tokamak_neutron_source.reactions.Reactions | tokamak_neutron_source.reactions.AneutronicReactions module-attribute

Represent a PEP 604 union type

E.g. for int | str

EnergySpectrumMethod

Bases: enum.Enum

Energy spectrum calculation method.

Source code in tokamak_neutron_source/energy.py 
23
24
25
26
27
28
29
class EnergySpectrumMethod(Enum):
    """Energy spectrum calculation method."""

    DATA = auto()
    BALLABIO_GAUSSIAN = auto()
    BALLABIO_M_GAUSSIAN = auto()
    AUTO = auto()

AUTO = <EnergySpectrumMethod.AUTO: 4> class-attribute

Energy spectrum calculation method.

BALLABIO_GAUSSIAN = <EnergySpectrumMethod.BALLABIO_GAUSSIAN: 2> class-attribute

Energy spectrum calculation method.

BALLABIO_M_GAUSSIAN = <EnergySpectrumMethod.BALLABIO_M_GAUSSIAN: 3> class-attribute

Energy spectrum calculation method.

DATA = <EnergySpectrumMethod.DATA: 1> class-attribute

Energy spectrum calculation method.

QuietTTSpectrumWarnings

Bases: tokamak_neutron_source.tools.WarningFilter

Filter away all duplicate warnings from the energy and energy_data module.

Source code in tokamak_neutron_source/tools.py 
325
326
327
328
329
330
331
class QuietTTSpectrumWarnings(WarningFilter):
    """Filter away all duplicate warnings from the energy and energy_data module."""

    def __init__(self):
        super().__init__(
            "tokamak_neutron_source.energy", "tokamak_neutron_source.energy_data"
        )

Reactions

Bases: tokamak_neutron_source.reactions.ReactionEnumMixin, enum.Enum

Neutronic reaction channels.

Source code in tokamak_neutron_source/reactions.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
class Reactions(ReactionEnumMixin, Enum):
    """Neutronic reaction channels."""

    D_T = ReactionData(
        label="D + T → ⁴He + n",
        total_energy=E_DT_FUSION,
        num_neutrons=1,
        cross_section=DT_XS,
        bosch_hale_coefficients=BOSCH_HALE_DT_4HEN,
        ballabio_spectrum=BALLABIO_DT_NEUTRON,
    )
    D_D = ReactionData(
        label="D + D → ³He + n",
        total_energy=E_DD_HE3N_FUSION,
        num_neutrons=1,
        cross_section=DD_HE3N_XS,
        bosch_hale_coefficients=BOSCH_HALE_DD_3HEN,
        ballabio_spectrum=BALLABIO_DD_NEUTRON,
    )
    T_T = ReactionData(
        label="T + T → ⁴He + 2n",
        total_energy=E_TT_FUSION,
        num_neutrons=2,
        cross_section=TT_XS,
        bosch_hale_coefficients=None,
        ballabio_spectrum=None,
    )

D_D = <Reactions.D_D: ReactionData(label='D + D → ³He + n', total_energy=5.237367559215132e-13, num_neutrons=1, cross_section=<tokamak_neutron_source.reactivity_data.ReactionCrossSection object at 0x7ff2008dfb10>, bosch_hale_coefficients=BoschHaleCoefficients(name='D + D --> 3He + n', t_min=0.2, t_max=100.0, bg=31.397, mrc2=937814.0, c=array([ 5.43360e-12, 5.85778e-03, 7.68222e-03, 0.00000e+00,-2.96400e-06, 0.00000e+00, 0.00000e+00])), ballabio_spectrum=BallabioEnergySpectrum(energy_0=2449.5, omega_0=82.542, energy_shift_coeffs=BallabioCoefficients(a1=4.69515, a2=-0.040729, a3=0.47, a4=0.81844), width_correction_coeffs=BallabioCoefficients(a1=0.0017013, a2=0.16888, a3=0.49, a4=0.0007946)))> class-attribute

Neutronic reaction channels.

D_T = <Reactions.D_T: ReactionData(label='D + T → ⁴He + n', total_energy=2.8183035155819573e-12, num_neutrons=1, cross_section=<tokamak_neutron_source.reactivity_data.ReactionCrossSection object at 0x7ff2008d4c50>, bosch_hale_coefficients=BoschHaleCoefficients(name='D + T --> 4He + n', t_min=0.2, t_max=100.0, bg=34.3827, mrc2=1124656.0, c=array([ 1.17302e-09, 1.51361e-02, 7.51886e-02, 4.60643e-03,1.35000e-02, -1.06750e-04, 1.36600e-05])), ballabio_spectrum=BallabioEnergySpectrum(energy_0=14021.0, omega_0=177.259, energy_shift_coeffs=BallabioCoefficients(a1=5.30509, a2=0.0024736, a3=1.84, a4=1.3818), width_correction_coeffs=BallabioCoefficients(a1=0.00051068, a2=0.0076223, a3=1.78, a4=8.7691e-05)))> class-attribute

Neutronic reaction channels.

T_T = <Reactions.T_T: ReactionData(label='T + T → ⁴He + 2n', total_energy=1.8157845541890245e-12, num_neutrons=2, cross_section=<tokamak_neutron_source.reactivity_data.ReactionCrossSection object at 0x7ff20ad03850>, bosch_hale_coefficients=None, ballabio_spectrum=None)> class-attribute

Neutronic reaction channels.

get_neutron_energy_spectrum(reaction, temp_kev, method)

Get a native OpenMC neutron energy spectrum.

Parameters:

Name Type Description Default
reaction Reactions

The neutronic reaction for which to retrieve the neutron spectrum

 required
temp_kev float

The ion temperature of the reactants

 required
method EnergySpectrumMethod

Which method to use when calculating the energy spectrum

 required

Returns:

Type Description
Tabular | Discrete

OpenMC tabular neutron energy distribution for the given reaction.
Notes

Log-linear interpolation is used within OpenMC.

Source code in tokamak_neutron_source/openmc_interface.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def get_neutron_energy_spectrum(
    reaction: Reactions, temp_kev: float, method: EnergySpectrumMethod
) -> Tabular | Discrete:
    """
    Get a native OpenMC neutron energy spectrum.

    Parameters
    ----------
    reaction:
        The neutronic reaction for which to retrieve the neutron spectrum
    temp_kev: float
        The ion temperature of the reactants
    method:
        Which method to use when calculating the energy spectrum

    Returns
    -------
    :
        OpenMC tabular neutron energy distribution for the given reaction.

    Notes
    -----
    Log-linear interpolation is used within OpenMC.
    """
    if (
        method is EnergySpectrumMethod.BALLABIO_GAUSSIAN
        and reaction.ballabio_spectrum is not None
    ):
        mean = reaction.ballabio_spectrum.mean_energy(temp_kev)
        std = reaction.ballabio_spectrum.std_deviation(temp_kev)
        return Normal(raw_uc(mean, "keV", "eV"), raw_uc(std, "keV", "eV"))
    energy, probability = energy_spectrum(temp_kev, reaction, method)
    energy_ev = raw_uc(energy, "keV", "eV")
    # Log-linear interpolation is not supported in OpenMC at present
    # see: https://github.com/openmc-dev/openmc/issues/2409
    return Tabular(energy_ev, probability, interpolation="linear-linear")

make_openmc_ring_source(r, z, half_cell_length, energy_distribution, strength)

Make a single OpenMC ring source with a square cross-section.

Parameters:

Name Type Description Default
r float

Radial position of the centroid of the 3-D ring [m]

 required
z float

Vertical position of the centroid of the 3-D ring [m]

 required
half_cell_length float

Half the square cell length [m]

 required
energy_distribution Univariate

Neutron energy distribution

 required
strength float

Strength of the source [number of neutrons]

 required

Returns:

Type Description
IndependentSource

An OpenMC IndependentSource object, or None if strength is zero.
Notes

The z values within the square cell are uniform, and the r values vary linearly with increasing radius.

Source code in tokamak_neutron_source/openmc_interface.py 
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def make_openmc_ring_source(
    r: float,
    z: float,
    half_cell_length: float,
    energy_distribution: Univariate,
    strength: float,
) -> IndependentSource:
    """
    Make a single OpenMC ring source with a square cross-section.

    Parameters
    ----------
    r:
        Radial position of the centroid of the  3-D ring [m]
    z:
        Vertical position of the centroid of the 3-D ring [m]
    half_cell_length:
        Half the square cell length [m]
    energy_distribution:
        Neutron energy distribution
    strength:
        Strength of the source [number of neutrons]

    Returns
    -------
    :
        An OpenMC IndependentSource object, or None if strength is zero.

    Notes
    -----
    The z values within the square cell are uniform, and the r values vary
    linearly with increasing radius.
    """
    if strength > 0:
        r_in, r_out = r - half_cell_length, r + half_cell_length
        z_down, z_up = z - half_cell_length, z + half_cell_length
        r_lim_cm = raw_uc([r_in, r_out], "m", "cm")
        z_lim_cm = raw_uc([z_down, z_up], "m", "cm")
        r_lim_prob = np.array(r_lim_cm) / sum(r_lim_cm)
        return IndependentSource(
            energy=energy_distribution,
            space=CylindricalIndependent(
                r=Tabular(r_lim_cm, r_lim_prob, interpolation="linear-linear"),
                phi=Uniform(0, 2 * np.pi),
                z=Uniform(*z_lim_cm),
                origin=(0.0, 0.0, 0.0),
            ),
            angle=Isotropic(),
            strength=strength,
        )
    return None

make_openmc_full_combined_source(r, z, cell_side_length, temperature, strength, source_rate, energy_spectrum_method=<EnergySpectrumMethod.AUTO: 4>)

Make an OpenMC source combining multiple reactions across the whole plasma.

Parameters:

Name Type Description Default
r ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Radial positions of the rings [m]

 required
z ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Vertical positions of the rings [m]

 required
cell_side_length float

Radial and vertical spacings of the rings [m]

 required
temperature ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Ion temperatures at the rings [keV]

 required
strength dict[Reactions | AneutronicReactions, ndarray[tuple[Any, ...], dtype[~_ScalarT]]]

Dictionary of strengths for each reaction at the rings [neutrons]

 required
source_rate float

Total source rate [neutrons/s]

 required
energy_spectrum_method EnergySpectrumMethod

Which method to use when calculating neutron spectra

 <EnergySpectrumMethod.AUTO: 4>

Returns:

Type Description
IndependentSource

A list of OpenMC IndependentSource objects, one per ring.
Source code in tokamak_neutron_source/openmc_interface.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
def make_openmc_full_combined_source(  # noqa: PLR0913, PLR0917
    r: npt.NDArray,
    z: npt.NDArray,
    cell_side_length: float,
    temperature: npt.NDArray,
    strength: dict[AllReactions, npt.NDArray],
    source_rate: float,
    energy_spectrum_method: EnergySpectrumMethod = EnergySpectrumMethod.AUTO,
) -> IndependentSource:
    """
    Make an OpenMC source combining multiple reactions across the whole plasma.

    Parameters
    ----------
    r:
        Radial positions of the rings [m]
    z:
        Vertical positions of the rings [m]
    cell_side_length:
        Radial and vertical spacings of the rings [m]
    temperature:
        Ion temperatures at the rings [keV]
    strength:
        Dictionary of strengths for each reaction at the rings [neutrons]
    source_rate:
        Total source rate [neutrons/s]
    energy_spectrum_method:
        Which method to use when calculating neutron spectra

    Returns
    -------
    :
        A list of OpenMC IndependentSource objects, one per ring.
    """
    sources = []
    # Neutronic reaction channels only
    # We multiply the T-T channel by 2 because it is 2n
    n_strength = {
        reaction: rate * reaction.num_neutrons
        for reaction, rate in strength.items()
        if isinstance(reaction, Reactions)
    }

    l_2 = cell_side_length / 2
    with QuietTTSpectrumWarnings():
        for i, (ri, zi, ti) in enumerate(zip(r, z, temperature, strict=False)):
            distributions = []
            weights = []

            for reaction, s in n_strength.items():
                if s[i] > 0.0:
                    distributions.append(
                        get_neutron_energy_spectrum(reaction, ti, energy_spectrum_method)
                    )
                    weights.append(s[i])

            local_strength = sum(weights)

            distribution = Mixture(np.array(weights) / local_strength, distributions)

            source = make_openmc_ring_source(
                ri,
                zi,
                l_2,
                distribution,
                local_strength / source_rate,
            )
            if source is not None:
                sources.append(source)

    return sources

energy_spectrum(temp_kev, reaction, method=<EnergySpectrumMethod.BALLABIO_M_GAUSSIAN: 3>)

Calculate the tabulated energy spectrum of a reaction at a given ion temperature.

Parameters:

Name Type Description Default
temp_kev float

Ion temperature

 required
reaction Reactions

Neutronic fusion reaction

 required
method EnergySpectrumMethod

Method to use to calculate the energy spectrum

 <EnergySpectrumMethod.BALLABIO_M_GAUSSIAN: 3>

Returns:

Name Type Description
energies ndarray[tuple[Any, ...], dtype[~_ScalarT]]

The energy bins of the probability distribution function
pdf ndarray[tuple[Any, ...], dtype[~_ScalarT]]

The PDF values
Source code in tokamak_neutron_source/energy.py 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def energy_spectrum(
    temp_kev: float,
    reaction: Reactions,
    method: EnergySpectrumMethod = EnergySpectrumMethod.BALLABIO_M_GAUSSIAN,
) -> tuple[npt.NDArray, npt.NDArray]:
    """
    Calculate the tabulated energy spectrum of a reaction at a given ion temperature.

    Parameters
    ----------
    temp_kev:
        Ion temperature
    reaction:
        Neutronic fusion reaction
    method:
        Method to use to calculate the energy spectrum

    Returns
    -------
    energies:
        The energy bins of the probability distribution function
    pdf:
        The PDF values
    """
    match method:
        case (
            EnergySpectrumMethod.BALLABIO_GAUSSIAN
            | EnergySpectrumMethod.BALLABIO_M_GAUSSIAN
        ):
            if reaction.ballabio_spectrum is not None:
                return _ballabio_spectrum(reaction.ballabio_spectrum, temp_kev, method)

            logger.warning(
                f"There is no Ballabio parameterisation for reaction {reaction.name}, "
                "returning energy spectrum calculated by data.",
                stacklevel=5,
            )
            return _data_spectrum(reaction, temp_kev)

        case EnergySpectrumMethod.DATA:
            return _data_spectrum(reaction, temp_kev)
        case EnergySpectrumMethod.AUTO:
            match reaction:
                case Reactions.D_D | Reactions.D_T:
                    return energy_spectrum(
                        temp_kev, reaction, EnergySpectrumMethod.BALLABIO_M_GAUSSIAN
                    )
                case Reactions.T_T:
                    return energy_spectrum(temp_kev, reaction, EnergySpectrumMethod.DATA)

raw_uc(value, unit_from, unit_to)

Raw unit converter Converts a value from one unit to another

Parameters:

Name Type Description Default
value ValueLikeT

value to convert

 required
unit_from str | Unit

unit to convert from

 required
unit_to str | Unit

unit to convert to

 required

Returns:

Type Description
ValueLikeT

converted value
Source code in tokamak_neutron_source/constants.py 
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
def raw_uc(
    value: ValueLikeT,
    unit_from: str | ureg.Unit,
    unit_to: str | ureg.Unit,
) -> ValueLikeT:
    """
    Raw unit converter
    Converts a value from one unit to another

    Parameters
    ----------
    value:
        value to convert
    unit_from:
        unit to convert from
    unit_to:
        unit to convert to

    Returns
    -------
    :
        converted value
    """
    try:
        return (
            ureg.Quantity(value, ureg.Unit(unit_from)).to(ureg.Unit(unit_to)).magnitude
        )
    except ValueError:
        # Catch scales on units eg the ridculousness of this unit: 10^19/m^3
        unit_from_q = ureg.Quantity(unit_from)
        unit_to_q = ureg.Quantity(unit_to)
        return (
            ureg.Quantity(value * unit_from_q).to(unit_to_q.units).magnitude
            / unit_to_q.magnitude
        )