mcnp_interface

MCNP neutron source (SDEF) 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

AneutronicReactions

Bases: tokamak_neutron_source.reactions.ReactionEnumMixin, enum.Enum

Aneutronic reaction channels.

Source code in tokamak_neutron_source/reactions.py

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 0x7ff2008b80d0>, 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 0x7ff200901310>, bosch_hale_coefficients=None, ballabio_spectrum=None)> class-attribute

Aneutronic reaction channels.

EnergySpectrumMethod

Bases: enum.Enum

Energy spectrum calculation method.

Source code in tokamak_neutron_source/energy.py

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.

Reactions

Bases: tokamak_neutron_source.reactions.ReactionEnumMixin, enum.Enum

Neutronic reaction channels.

Source code in tokamak_neutron_source/reactions.py

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.

write_mcnp_sdef_source(file, r, z, cell_side_length, temperature, strength)

Write an MCNP SDEF source for a ring source at (r,z).

Parameters:

Name	Type	Description	Default
file	str | Path	The file name stub to which to write the SDEF source	required
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	side length of square source cell	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 [arbitrary units]	required

Notes

Only Neutronic reactions are written to SDEF file Aneutronic reactions are ignored. The radial distribution bouldaries and probabilities are set to the SI3 and SP3 cards The DS4 card is used as the dependent distribution numbers for the vertical distributions

Source code in tokamak_neutron_source/mcnp_interface.py

def write_mcnp_sdef_source(
    file: str | Path,
    r: npt.NDArray,
    z: npt.NDArray,
    cell_side_length: float,
    temperature: npt.NDArray,
    strength: dict[AllReactions, npt.NDArray],
):
    """
    Write an MCNP SDEF source for a ring source at (r,z).

    Parameters
    ----------
    file:
        The file name stub to which to write the SDEF source
    r:
        Radial positions of the rings [m]
    z:
        Vertical positions of the rings [m]
    cell_side_length:
        side length of square source cell
    temperature:
        Ion temperatures at the rings [keV]
    strength:
        Dictionary of strengths for each reaction at the rings [arbitrary units]

    Notes
    -----
    Only Neutronic reactions are written to SDEF file Aneutronic reactions are ignored.
    The radial distribution bouldaries and probabilities are set to the SI3 and SP3 cards
    The DS4 card is used as the dependent distribution numbers
    for the vertical distributions

    """
    r = raw_uc(r, "m", "cm")
    z = raw_uc(z, "m", "cm")
    dr = raw_uc(cell_side_length, "m", "cm")

    # Half widths of 'cells'
    drad = dzed = dr / 2

    offset = 5  # First 5 distribution are reserved

    for reaction, r_data in strength.items():
        if reaction not in AneutronicReactions:
            short_react = re.findall(r"[DT]", reaction.label)
            file_name = f"{file}.{short_react[0]}{short_react[1]}"

            header = sdef_header(reaction, r_data, temperature)

            # Calculate the radial boundaries based on the ring centres
            # and 'cell width' (dr)
            r_bounds = np.unique(r) - drad
            r_bounds = np.append(r_bounds, r_bounds[-1] + dr)  # Add the last boundary

            # Identify where radial position changes
            # and therefore the range of each vertical distribution
            z_ints = np.concatenate([[-1], np.nonzero(r[1:] != r[:-1])[0], [len(r) - 1]])

            si_card = _text_wrap(
                "SI3 H " + " ".join(f"{ri:.5e}" for ri in r_bounds), new_lines=1
            )
            sp_card = _text_wrap(
                f"SP3 D {0.0:.5e} "
                + " ".join(
                    f"{np.sum(r_data[z_ints[i] + 1 : z_ints[i + 1] + 1]):.5e}"
                    for i in range(len(z_ints) - 1)
                ),
                new_lines=1,
            )

            ds_card = _text_wrap(
                "DS4 S "
                + " ".join(f"{i:d}" for i in range(offset, offset + len(r_bounds) - 1)),
                new_lines=1,
            )

            with open(file_name, "w") as sdef_file:
                sdef_file.write(
                    f"{header}{si_card}{sp_card}{ds_card}"
                    "C\nC 3. Neutron Emission Probability - Vertical Distribution\nC\n"
                )
                for si_card, sp_card in _si_sp_vertical_dist_cards(
                    offset, z, z_ints, dzed, r_data
                ):
                    sdef_file.write(f"{si_card}{sp_card}")

        else:
            logger.info(f"Skipping reaction {reaction.label} for MCNP SDEF source.")

sdef_header(reaction, reaction_data, temperature)

Create SDEF file header

Parameters:

Name	Type	Description	Default
reaction	Reactions	Reaction to be created	required
reaction_data	ndarray[tuple[Any, ...], dtype[~_ScalarT]]	strength of source	required
temperature	float	Ion temperature	required

Notes

For DT and DD reactions MCNP's built-in gaussian spectrums are used For TT reactions the tabulated data is used

Source code in tokamak_neutron_source/mcnp_interface.py

def sdef_header(
    reaction: Reactions, reaction_data: npt.NDArray, temperature: float
) -> str:
    """Create SDEF file header

    Parameters
    ----------
    reaction:
        Reaction to be created
    reaction_data:
        strength of source
    temperature:
        Ion temperature

    Notes
    -----
    For DT and DD reactions MCNP's built-in gaussian spectrums are used
    For TT reactions the tabulated data is used

    """  # noqa: DOC201
    strength = sum(reaction_data)
    ion_temp = mean_ion_temp(reaction_data, temperature)

    # Read authors and git address from CITATION.cff
    citation = load_citation()
    authors = "\n".join(
        (
            f"C    {author.get('given-names', '')}"
            f" {author.get('family-names', '')},"
            f" {author.get('affiliation', '')}"
        )
        for author in citation.get("authors", [])
    )
    gitaddr = citation.get("repository-code", "")

    if reaction in {Reactions.D_T, Reactions.D_D}:
        # -1 for D-T and -2 for D-D
        reaction_data = (
            f"SP2 -4 {raw_uc(ion_temp, 'keV', 'MeV'):5e} "
            f"{-1 if reaction == Reactions.D_T else -2}\n"
        )

    elif reaction == Reactions.T_T:
        energies, probabilities = energy_spectrum(
            ion_temp, reaction, EnergySpectrumMethod.DATA
        )
        reaction_data = "SI2 H " + _text_wrap(
            f"{0.0:.5e} " + " ".join(f"{e:.5e}" for e in raw_uc(energies, "keV", "MeV"))
        )
        reaction_data += "SP2 D " + _text_wrap(
            f"{0.0:.5e} " + " ".join(f"{p:.5e}" for p in probabilities)
        )

    return f"""C ============================
C SDEF Card for Tokamak Neutron Source generated by:
C {gitaddr}
C ============================
C
C ============================
C Authors:
{authors}
C ============================
C
C ============================
C Method:
C 1. Create a cylinder that encloses the entire torus.
C 2. Then slice the cylinder along the R-axis.
C 3. Finally, define the vertical distribution, assuming rotational symmetry.
C ============================
C
C ============================
C Reaction channel: {reaction.label}
C Total source neutrons: {strength:5e} n/s
C ============================
C
C 1. Neutron Emission Probability - Set up cylindrical source
C
sdef erg=d2 par=1 wgt=1
      pos = 0 0 0    $ Center = origin
      axs = 0 0 1    $ Cylinder points along the Z axis
      rad = d3       $ radial distribution defined by distribution 3
      ext = frad d4  $ extent distribution defined by distribution 4 which is dependent on distribution rad
{reaction_data}C
C 2. Neutron Emission Probability - Radial Distribution
C
"""  # noqa: E501

mean_ion_temp(strength, temperature)

Calculate the strength-weighted mean ion temperature.

Source code in tokamak_neutron_source/mcnp_interface.py

def mean_ion_temp(strength: npt.NDArray, temperature: npt.NDArray) -> float:
    """Calculate the strength-weighted mean ion temperature."""  # noqa: DOC201
    return np.sum(strength * temperature) / np.sum(strength)

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

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)

load_citation()

Load the CITATION.cff file.

Returns:

Type	Description
dict	The contents of the CITATION.cff file as a dictionary.

Source code in tokamak_neutron_source/tools.py

def load_citation() -> dict:
    """
    Load the CITATION.cff file.

    Returns
    -------
    :
        The contents of the CITATION.cff file as a dictionary.
    """
    with open(get_tns_path("data") / "CITATION.cff") as citation_file:
        return yaml.safe_load(citation_file)

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

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
        )