Skip to content

reactivity

Reactivity calculations.

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

Represent a PEP 604 union type

E.g. for int | str

ReactivityMethod

Bases: enum.Enum

Reactivity calculation method.

Source code in tokamak_neutron_source/reactivity.py 
35
36
37
38
39
40
class ReactivityMethod(Enum):
    """Reactivity calculation method."""

    XS = auto()
    BOSCH_HALE = auto()
    AUTO = auto()

AUTO = <ReactivityMethod.AUTO: 3> class-attribute

Reactivity calculation method.

BOSCH_HALE = <ReactivityMethod.BOSCH_HALE: 2> class-attribute

Reactivity calculation method.

XS = <ReactivityMethod.XS: 1> class-attribute

Reactivity calculation method.

AneutronicReactions

Bases: tokamak_neutron_source.reactions.ReactionEnumMixin, enum.Enum

Aneutronic reaction channels.

Source code in tokamak_neutron_source/reactions.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
class AneutronicReactions(ReactionEnumMixin, Enum):
    """Aneutronic reaction channels."""

    D_D = ReactionData(
        label="D + D → T + p",
        total_energy=E_DD_TP_FUSION,
        num_neutrons=0,  # no neutrons in aneutronic branch
        cross_section=DD_TP_XS,
        bosch_hale_coefficients=BOSCH_HALE_DD_TP,
        ballabio_spectrum=None,
    )
    D_He3 = ReactionData(
        label="D + ³He → ⁴He + p",
        total_energy=E_DHE3_FUSION,
        num_neutrons=0,
        cross_section=DHE3_HEP_XS,
        bosch_hale_coefficients=None,
        ballabio_spectrum=None,
    )

D_D = <AneutronicReactions.D_D: ReactionData(label='D + D → T + p', total_energy=6.461016407480568e-13, num_neutrons=0, cross_section=<tokamak_neutron_source.reactivity_data.ReactionCrossSection object at 0x7fb504d05c50>, bosch_hale_coefficients=BoschHaleCoefficients(name='D + D --> T + p', t_min=0.2, t_max=100.0, bg=31.397, mrc2=937814.0, c=array([5.65718e-12, 3.41267e-03, 1.99167e-03, 0.00000e+00, 1.05060e-05,0.00000e+00, 0.00000e+00])), ballabio_spectrum=None)> class-attribute

Aneutronic reaction channels.

D_He3 = <AneutronicReactions.D_He3: ReactionData(label='D + ³He → ⁴He + p', total_energy=2.940668400408501e-12, num_neutrons=0, cross_section=<tokamak_neutron_source.reactivity_data.ReactionCrossSection object at 0x7fb504da04d0>, bosch_hale_coefficients=None, ballabio_spectrum=None)> class-attribute

Aneutronic reaction channels.

BoschHaleCoefficients

Bosch-Hale parameterisation dataclass.

H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611 DOI 10.1088/0029-5515/32/4/I07

Source code in tokamak_neutron_source/reactivity_data.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
@dataclass
class BoschHaleCoefficients:
    """
    Bosch-Hale parameterisation dataclass.

    H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611
    DOI 10.1088/0029-5515/32/4/I07
    """

    name: str
    t_min: float  # [keV]
    t_max: float  # [keV]
    bg: float  # [keV**0.5]
    mrc2: float  # [keV]
    c: npt.NDArray

ReactionCrossSection

Fusion reaction cross-section.

Source code in tokamak_neutron_source/reactivity_data.py 
18
19
20
21
22
23
24
25
26
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
class ReactionCrossSection:
    """Fusion reaction cross-section."""

    def __init__(self, file_name: str):
        """
        Parameters
        ----------
        file_name:
            Cross-sectional data file (ENDF format)

        Raises
        ------
        ReactivityError
            Data file path is not a file
        """
        path = get_tns_path("data/cross_sections")
        path = Path(path, file_name)
        if not path.is_file():
            raise ReactivityError(f"Cross-section data file {path} is not a file!")

        file = path.as_posix()
        # Read in the cross section (in barn) as a function of energy (MeV).
        energy, sigma = np.genfromtxt(file, comments="#", skip_footer=2, unpack=True)

        split = file_name.split(".", maxsplit=1)[0].split("_")
        collider, target = split[:2]
        self.name = f"{collider} + {target} -> {split[2]} + {split[3]}"

        mass_1, mass_2 = MOLAR_MASSES[collider], MOLAR_MASSES[target]

        self.reduced_mass = raw_uc(mass_1 * mass_2 / (mass_1 + mass_2), "amu", "kg")

        # Convert to center of mass frame
        # NOTE MC: Choice of target/collider thing makes Bosch-Hale line up...
        energy *= mass_2 / (mass_1 + mass_2)

        # Convert to kev / m^2
        self._cross_section = interp1d(energy * 1e3, sigma * 1e-28)

    def __call__(self, temp_kev: float | npt.NDArray) -> float | npt.NDArray:
        """Get cross section at a give temperature"""  # noqa: DOC201
        return self._cross_section(temp_kev)

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 0x7fb53dd00150>, 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 0x7fb504da0250>, 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 0x7fb504d603d0>, bosch_hale_coefficients=None, ballabio_spectrum=None)> class-attribute

Neutronic reaction channels.

ReactivityError

Bases: tokamak_neutron_source.error.TNSError

Reactivity error class

Source code in tokamak_neutron_source/error.py 
29
30
class ReactivityError(TNSError):
    """Reactivity error class"""

density_weighted_reactivity(temp_kev, density_d, density_t, density_he3, reaction=<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 0x7fb504da0250>, 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)))>, method=<ReactivityMethod.AUTO: 3>)

Calculate the density-weighted thermal reactivity of a fusion reaction in Maxwellian plasmas, \t:math:n_1 n_2 <\sigma v>.

Parameters:

Name Type Description Default
temp_kev float | ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Temperature [keV]

 required
density_d float | ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Deuterium density [m^-3]

 required
density_t float | ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Tritium density [m^-3]

 required
density_he3 float | ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Helium-3 density [m^-3]

 required
reaction str | Reactions | AneutronicReactions

The fusion reaction

 <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 0x7fb504da0250>, 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)))>
method ReactivityMethod

The method to use when calculating the reactivity

 <ReactivityMethod.AUTO: 3>

Returns:

Type Description
float | ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Density-weighted reactivity of the reaction at the specified temperature(s) [1/m^3/s]
Source code in tokamak_neutron_source/reactivity.py 
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
81
82
83
84
85
86
87
88
89
90
def density_weighted_reactivity(
    temp_kev: float | npt.NDArray,
    density_d: float | npt.NDArray,
    density_t: float | npt.NDArray,
    density_he3: float | npt.NDArray,
    reaction: str | AllReactions = Reactions.D_T,
    method: ReactivityMethod = ReactivityMethod.AUTO,
) -> float | npt.NDArray:
    """
    Calculate the density-weighted thermal reactivity of a fusion reaction in
    Maxwellian plasmas, \\t:math:`n_1 n_2 <\\sigma v>`.

    Parameters
    ----------
    temp_kev:
        Temperature [keV]
    density_d:
        Deuterium density [m^-3]
    density_t:
        Tritium density [m^-3]
    density_he3:
        Helium-3 density [m^-3]
    reaction:
        The fusion reaction
    method:
        The method to use when calculating the reactivity

    Returns
    -------
    :
        Density-weighted reactivity of the reaction at the specified temperature(s)
        [1/m^3/s]
    """
    reaction = _parse_reaction(reaction)

    match reaction:
        case Reactions.D_D | AneutronicReactions.D_D:
            n1_n2 = density_d * density_d / 2
        case Reactions.D_T:
            n1_n2 = density_d * density_t
        case Reactions.T_T:
            n1_n2 = density_t * density_t / 2
        case AneutronicReactions.D_He3:
            n1_n2 = density_d * density_he3
        case _:
            raise NotImplementedError(f"Reaction {reaction} not implemented.")

    return n1_n2 * reactivity(temp_kev, reaction, method)

reactivity(temp_kev, reaction=<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 0x7fb504da0250>, 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)))>, method=<ReactivityMethod.AUTO: 3>)

Calculate the thermal reactivity of a fusion reaction in Maxwellian plasmas, \t:math:<\sigma v>.

Parameters:

Name Type Description Default
temp_kev float | ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Temperature [keV]

 required
reaction str | Reactions | AneutronicReactions

The fusion reaction

 <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 0x7fb504da0250>, 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)))>
method ReactivityMethod

The method to use when calculating the reactivity

 <ReactivityMethod.AUTO: 3>

Returns:

Type Description
float | ndarray[tuple[Any, ...], dtype[~_ScalarT]]

Reactivity of the reaction at the specified temperature(s) [m^3/s]
Source code in tokamak_neutron_source/reactivity.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
def reactivity(
    temp_kev: float | npt.NDArray,
    reaction: str | AllReactions = Reactions.D_T,
    method: ReactivityMethod = ReactivityMethod.AUTO,
) -> float | npt.NDArray:
    """
    Calculate the thermal reactivity of a fusion reaction in Maxwellian plasmas,
    \\t:math:`<\\sigma v>`.

    Parameters
    ----------
    temp_kev:
        Temperature [keV]
    reaction:
        The fusion reaction
    method:
        The method to use when calculating the reactivity

    Returns
    -------
    :
        Reactivity of the reaction at the specified temperature(s) [m^3/s]
    """
    reaction = _parse_reaction(reaction)

    match method:
        case ReactivityMethod.AUTO:
            if reaction.bosch_hale_coefficients is not None:
                return _reactivity_bosch_hale(temp_kev, reaction)
            return _reactivity_from_xs(temp_kev, reaction.cross_section)
        case ReactivityMethod.BOSCH_HALE:
            if reaction.bosch_hale_coefficients is not None:
                return _reactivity_bosch_hale(temp_kev, reaction)

            logger.warning(
                f"There is no Bosch-Hale parameterisation for reaction {reaction.name}, "
                "returning reactivity calculated by cross-section."
            )
            return _reactivity_from_xs(temp_kev, reaction.cross_section)
        case ReactivityMethod.XS:
            return _reactivity_from_xs(temp_kev, reaction.cross_section)

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
        )