main

Neutromak user-facing API.

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

Represent a PEP 604 union type

E.g. for int | str

CUSTOM_ORDER = [<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)))>, <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)))>, <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)>, <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)>, <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)>] module-attribute

Built-in mutable sequence.

If no argument is given, the constructor creates a new empty list. The argument must be an iterable if specified.

TYPE_CHECKING = False module-attribute

bool(x) -> bool

Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.

annotations = _Feature((3, 7, 0, 'beta', 1), None, 16777216) module-attribute

TokamakNeutronSource

Tokamak neutron source object.

Parameters:

Name	Type	Description	Default
transport	TransportInformation	Plasma profile information and species composition.	required
flux_map	FluxMap	Magneto-hydrodynamic equilibrium poloidal magnetic flux map containing LCFS geometry and the psi_norm inside of it.	required
source_type	AllReactions | list[AllReactions] | None	Which neutronic reaction(s) to include in the neutron source.	None
cell_side_length	float	Discretisation in [m] for sampling the 2-D poloidal plane within the LCFS. Square cells of side length cell_side_length are created.	0.1
total_fusion_power	float | None	If specified, will be used to normalise the source strength to a prescribed fusion power	None
reactivity_method	ReactivityMethod	Which method to use when calculating reactivities	<ReactivityMethod.AUTO: 3>

Source code in tokamak_neutron_source/main.py

class TokamakNeutronSource:
    """
    Tokamak neutron source object.

    Parameters
    ----------
    transport:
        Plasma profile information and species composition.
    flux_map:
        Magneto-hydrodynamic equilibrium poloidal magnetic flux map containing LCFS
        geometry and the psi_norm inside of it.
    source_type:
        Which neutronic reaction(s) to include in the neutron source.
    cell_side_length:
        Discretisation in [m] for sampling the 2-D poloidal plane within the LCFS.
        Square cells of side length `cell_side_length` are created.
    total_fusion_power:
        If specified, will be used to normalise the source strength to a prescribed
        fusion power
    reactivity_method:
        Which method to use when calculating reactivities
    """

    def __init__(
        self,
        transport: TransportInformation,
        flux_map: FluxMap,
        source_type: AllReactions | list[AllReactions] | None = None,
        cell_side_length: float = 0.1,
        total_fusion_power: float | None = None,
        reactivity_method: ReactivityMethod = ReactivityMethod.AUTO,
    ):
        self.source_type = _parse_source_type(source_type)

        self.x, self.z, self.d_volume = sample_space_2d(
            flux_map.lcfs, flux_map.o_point, cell_side_length
        )
        self.cell_side_length = cell_side_length
        psi_norm = flux_map.psi_norm(self.x, self.z)

        self.temperature = transport.temperature_profile.value(psi_norm)
        density_d = transport.deuterium_density_profile.value(psi_norm)
        density_t = transport.tritium_density_profile.value(psi_norm)
        density_he3 = transport.helium3_density_profile.value(psi_norm)

        # All reactions (neutronic and aneutronic)
        self.strength = {}
        self.num_reactions_per_second = {}
        self.num_neutrons_per_second = {}
        for reaction in self.source_type:
            n1_n2_sigma = density_weighted_reactivity(
                self.temperature,
                density_d,
                density_t,
                density_he3,
                reaction=reaction,
                method=reactivity_method,
            )
            self.strength[reaction] = n1_n2_sigma * self.d_volume
            self.num_reactions_per_second[reaction] = sum(self.strength[reaction])
            self.num_neutrons_per_second[reaction] = (
                self.num_reactions_per_second[reaction] * reaction.num_neutrons
            )

        self.flux_map = flux_map
        self.transport = transport

        if total_fusion_power is not None:
            self.normalise_fusion_power(total_fusion_power)

    @property
    def source_rate(self) -> float:
        """
        The total source rate in [neutrons / s].
        """
        return sum(self.num_neutrons_per_second.values())

    @property
    def source_T_rate(self) -> float:  # noqa: N802
        """
        The T consumption rate in [tritons / s].

        Notes
        -----
        If you are using a "(n,Xt)" tally to calculate TBR, note that the definition
        of TBR is relative to the number of D-T reactions, not the total number
        of fusion reactions.

        To correctly scale your "(n,Xt)" tally in [1/particles], you should scale by:
            tbr *= source_rate / source_T_rate
        """
        return sum([
            self.num_reactions_per_second.get(Reactions.D_T, 0.0),
            2.0 * self.num_reactions_per_second.get(Reactions.T_T, 0.0),
        ])

    def calculate_total_fusion_power(self) -> float:
        """
        Calculate the total fusion power from all reaction channels.

        Returns
        -------
        total_fusion_power
            The total fusion power

        Notes
        -----
        The aneutronic fusion reactions are included here.
        """
        return sum(
            np.sum(
                self.strength[reaction] * reaction.total_energy,
            )
            for reaction in self.source_type
        )

    def normalise_fusion_power(self, total_fusion_power: float):
        """
        Renormalise the source strength to match a total fusion power across all
        channels.

        Parameters
        ----------
        total_fusion_power:
            The total fusion power to normalise to [W]

        Notes
        -----
        This is done assuming the provided total fusion power is for the
        same channels as the source_type. The ratios of the strengths of
        each reaction is assumed to be same as modelled here.
        """
        actual_fusion_power = self.calculate_total_fusion_power()
        scaling_factor = total_fusion_power / actual_fusion_power

        for reaction in self.source_type:
            self.strength[reaction] *= scaling_factor
            self.num_reactions_per_second[reaction] *= scaling_factor
            self.num_neutrons_per_second[reaction] *= scaling_factor

    def to_openmc_source(
        self,
        energy_method: EnergySpectrumMethod = EnergySpectrumMethod.AUTO,
    ) -> list[IndependentSource]:
        """
        Create an OpenMC tokamak neutron source.

        Parameters
        ----------
        energy_method:
            Which method to use when calculating neutron spectra

        Returns
        -------
        :
            List of native OpenMC source objects
        """
        from tokamak_neutron_source.openmc_interface import (  # noqa: PLC0415
            make_openmc_full_combined_source,
        )

        return make_openmc_full_combined_source(
            self.x,
            self.z,
            self.cell_side_length,
            self.temperature,
            self.strength,
            self.source_rate,
            energy_method,
        )

    def to_sdef_card(self, filename) -> str:
        """
        Create an SDEF card which MCNP/openmc can use to make a tokamak neutron source.

        Notes
        -----
        The position-dependence of neutron energies be captured by SDEF. Therefore the
        energy distribution of neutrons is averaged, and the same (frozen) distribution
        is used everywhere in the reactor.
        """
        write_mcnp_sdef_source(
            filename,
            self.x,
            self.z,
            self.cell_side_length,
            self.temperature,
            self.strength,
        )

    def to_h5_source(self):
        """
        Create a source in the HDF5 format such that the full distribution of neutron
        energies and position

        Returns
        -------
        :
            H5 format
        """
        raise NotImplementedError

    def plot(
        self, reactions: list[AllReactions] | None = None
    ) -> tuple[plt.Figure, plt.Axes]:
        """
        Plot the tokamak neutron source.

        Returns
        -------
        f:
            Matplotlib figure object
        ax:
            Matplotlib axes object

        Raises
        ------
        TNSError
            If the requested reactions are not in the source
        """
        if reactions is None:
            plot_reactions = self.source_type
        else:
            plot_reactions = []
            for reaction in reactions:
                if reaction in self.source_type:
                    plot_reactions.append(reaction)
                else:
                    logger.warning(
                        f"Cannot plot reaction {reaction}; it was not specified upon "
                        "instantiation."
                    )
            if len(plot_reactions) == 0:
                raise TNSError("No valid reactions to plot.")

        f, ax = plt.subplots(1, len(plot_reactions), figsize=[15, 8])
        if len(plot_reactions) == 1:
            ax = [ax]

        for axis, reaction in zip(ax, plot_reactions, strict=False):
            self.flux_map.plot(f=f, ax=axis)
            axis.set_title(f"{reaction.label} reaction")
            cm = axis.scatter(
                self.x,
                self.z,
                c=self.strength[reaction],
                cmap="inferno",
            )

            # Make a colorbar axis that matches the plotting area's height
            divider = make_axes_locatable(axis)
            cax = divider.append_axes("right", size="9%", pad=0.08)
            cb = f.colorbar(cm, cax=cax)

            # Put label on top
            cb.ax.set_title("[1/s]")
            axis.plot(self.flux_map.lcfs.x, self.flux_map.lcfs.z, color="r")
            axis.plot(self.flux_map.o_point.x, self.flux_map.o_point.z, "o", color="b")
        f.tight_layout()
        plt.show()
        return f, ax

source_T_rate property

The T consumption rate in [tritons / s].

Notes

If you are using a "(n,Xt)" tally to calculate TBR, note that the definition of TBR is relative to the number of D-T reactions, not the total number of fusion reactions.

To correctly scale your "(n,Xt)" tally in [1/particles], you should scale by: tbr *= source_rate / source_T_rate

source_rate property

The total source rate in [neutrons / s].

calculate_total_fusion_power()

Calculate the total fusion power from all reaction channels.

Returns:

Type	Description
total_fusion_power	The total fusion power

Notes

The aneutronic fusion reactions are included here.

Source code in tokamak_neutron_source/main.py

def calculate_total_fusion_power(self) -> float:
    """
    Calculate the total fusion power from all reaction channels.

    Returns
    -------
    total_fusion_power
        The total fusion power

    Notes
    -----
    The aneutronic fusion reactions are included here.
    """
    return sum(
        np.sum(
            self.strength[reaction] * reaction.total_energy,
        )
        for reaction in self.source_type
    )

normalise_fusion_power(total_fusion_power)

Renormalise the source strength to match a total fusion power across all channels.

Parameters:

Name	Type	Description	Default
total_fusion_power	float	The total fusion power to normalise to [W]	required

Notes

This is done assuming the provided total fusion power is for the same channels as the source_type. The ratios of the strengths of each reaction is assumed to be same as modelled here.

Source code in tokamak_neutron_source/main.py

def normalise_fusion_power(self, total_fusion_power: float):
    """
    Renormalise the source strength to match a total fusion power across all
    channels.

    Parameters
    ----------
    total_fusion_power:
        The total fusion power to normalise to [W]

    Notes
    -----
    This is done assuming the provided total fusion power is for the
    same channels as the source_type. The ratios of the strengths of
    each reaction is assumed to be same as modelled here.
    """
    actual_fusion_power = self.calculate_total_fusion_power()
    scaling_factor = total_fusion_power / actual_fusion_power

    for reaction in self.source_type:
        self.strength[reaction] *= scaling_factor
        self.num_reactions_per_second[reaction] *= scaling_factor
        self.num_neutrons_per_second[reaction] *= scaling_factor

to_openmc_source(energy_method=<EnergySpectrumMethod.AUTO: 4>)

Create an OpenMC tokamak neutron source.

Parameters:

Name	Type	Description	Default
energy_method	EnergySpectrumMethod	Which method to use when calculating neutron spectra	<EnergySpectrumMethod.AUTO: 4>

Returns:

Type	Description
list[IndependentSource]	List of native OpenMC source objects

Source code in tokamak_neutron_source/main.py

def to_openmc_source(
    self,
    energy_method: EnergySpectrumMethod = EnergySpectrumMethod.AUTO,
) -> list[IndependentSource]:
    """
    Create an OpenMC tokamak neutron source.

    Parameters
    ----------
    energy_method:
        Which method to use when calculating neutron spectra

    Returns
    -------
    :
        List of native OpenMC source objects
    """
    from tokamak_neutron_source.openmc_interface import (  # noqa: PLC0415
        make_openmc_full_combined_source,
    )

    return make_openmc_full_combined_source(
        self.x,
        self.z,
        self.cell_side_length,
        self.temperature,
        self.strength,
        self.source_rate,
        energy_method,
    )

to_sdef_card(filename)

Create an SDEF card which MCNP/openmc can use to make a tokamak neutron source.

Notes

The position-dependence of neutron energies be captured by SDEF. Therefore the energy distribution of neutrons is averaged, and the same (frozen) distribution is used everywhere in the reactor.

Source code in tokamak_neutron_source/main.py

def to_sdef_card(self, filename) -> str:
    """
    Create an SDEF card which MCNP/openmc can use to make a tokamak neutron source.

    Notes
    -----
    The position-dependence of neutron energies be captured by SDEF. Therefore the
    energy distribution of neutrons is averaged, and the same (frozen) distribution
    is used everywhere in the reactor.
    """
    write_mcnp_sdef_source(
        filename,
        self.x,
        self.z,
        self.cell_side_length,
        self.temperature,
        self.strength,
    )

to_h5_source()

Create a source in the HDF5 format such that the full distribution of neutron energies and position

Returns:

Type	Description
	H5 format

Source code in tokamak_neutron_source/main.py

def to_h5_source(self):
    """
    Create a source in the HDF5 format such that the full distribution of neutron
    energies and position

    Returns
    -------
    :
        H5 format
    """
    raise NotImplementedError

plot(reactions=None)

Plot the tokamak neutron source.

Returns:

Name	Type	Description
f	Figure	Matplotlib figure object
ax	Axes	Matplotlib axes object

Raises:

Type	Description
TNSError	If the requested reactions are not in the source

Source code in tokamak_neutron_source/main.py

def plot(
    self, reactions: list[AllReactions] | None = None
) -> tuple[plt.Figure, plt.Axes]:
    """
    Plot the tokamak neutron source.

    Returns
    -------
    f:
        Matplotlib figure object
    ax:
        Matplotlib axes object

    Raises
    ------
    TNSError
        If the requested reactions are not in the source
    """
    if reactions is None:
        plot_reactions = self.source_type
    else:
        plot_reactions = []
        for reaction in reactions:
            if reaction in self.source_type:
                plot_reactions.append(reaction)
            else:
                logger.warning(
                    f"Cannot plot reaction {reaction}; it was not specified upon "
                    "instantiation."
                )
        if len(plot_reactions) == 0:
            raise TNSError("No valid reactions to plot.")

    f, ax = plt.subplots(1, len(plot_reactions), figsize=[15, 8])
    if len(plot_reactions) == 1:
        ax = [ax]

    for axis, reaction in zip(ax, plot_reactions, strict=False):
        self.flux_map.plot(f=f, ax=axis)
        axis.set_title(f"{reaction.label} reaction")
        cm = axis.scatter(
            self.x,
            self.z,
            c=self.strength[reaction],
            cmap="inferno",
        )

        # Make a colorbar axis that matches the plotting area's height
        divider = make_axes_locatable(axis)
        cax = divider.append_axes("right", size="9%", pad=0.08)
        cb = f.colorbar(cm, cax=cax)

        # Put label on top
        cb.ax.set_title("[1/s]")
        axis.plot(self.flux_map.lcfs.x, self.flux_map.lcfs.z, color="r")
        axis.plot(self.flux_map.o_point.x, self.flux_map.o_point.z, "o", color="b")
    f.tight_layout()
    plt.show()
    return f, ax

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.

ReactivityError

Bases: tokamak_neutron_source.error.TNSError

Reactivity error class

Source code in tokamak_neutron_source/error.py

class ReactivityError(TNSError):
    """Reactivity error class"""

ReactivityMethod

Bases: enum.Enum

Reactivity calculation method.

Source code in tokamak_neutron_source/reactivity.py

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.

TNSError

Bases: builtins.Exception

Base exception class. Sub-class from this for module level Errors.

Source code in tokamak_neutron_source/error.py

class TNSError(Exception):
    """
    Base exception class. Sub-class from this for module level Errors.
    """

    def __str__(self) -> str:
        """
        Prettier handling of the Exception strings.

        Returns
        -------
        :
            The formatted exception string.
        """
        return fill(dedent(self.args[0]))

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 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)))>, 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 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)))>
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

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)

sample_space_2d(lcfs, o_point, cell_side_length)

Sample the 2-D poloidal plane within the LCFS.

Parameters:

Name	Type	Description	Default
lcfs	ClosedFluxSurface	Last closed flux surface	required
o_point	FluxPoint	O-point location	required
cell_side_length	float	Side length of square cells [m]	required

Returns:

Name	Type	Description
x	ndarray[tuple[Any, ...], dtype[~_ScalarT]]	Radial coordinates of sampled points [m]
z	ndarray[tuple[Any, ...], dtype[~_ScalarT]]	Vertical coordinates of sampled points [m]
d_volume	ndarray[tuple[Any, ...], dtype[~_ScalarT]]	Volumes of cells centred at points [m^3]

Notes

Creates points at the centres of square cells of fixed size (cell_side_length by cell_side_length). Only cells whose centres fall inside the LCFS polygon are kept. Cells are positioned such that the centre of one cell lies on the O-point.

Source code in tokamak_neutron_source/space.py

def sample_space_2d(
    lcfs: ClosedFluxSurface,
    o_point: FluxPoint,
    cell_side_length: float,
) -> tuple[npt.NDArray, npt.NDArray, npt.NDArray]:
    """
    Sample the 2-D poloidal plane within the LCFS.

    Parameters
    ----------
    lcfs:
        Last closed flux surface
    o_point:
        O-point location
    cell_side_length:
        Side length of square cells [m]

    Returns
    -------
    x:
        Radial coordinates of sampled points [m]
    z:
        Vertical coordinates of sampled points [m]
    d_volume:
        Volumes of cells centred at points [m^3]

    Notes
    -----
    Creates points at the centres of square cells of fixed size
    (cell_side_length by cell_side_length). Only cells whose centres fall
    inside the LCFS polygon are kept.
    Cells are positioned such that the  centre of one cell lies on
    the O-point.
    """
    # Get bounding box around LCFS (+ offset)
    polygon_path = Path(np.c_[lcfs.x, lcfs.z])
    off = 2.0 * cell_side_length  # Just to be sure
    x_min, x_max = np.min(lcfs.x) - off, np.max(lcfs.x) + off
    z_min, z_max = np.min(lcfs.z) - off, np.max(lcfs.z) + off

    # Shift grid so O-point lies on a cell center
    dx = (o_point.x - x_min) % cell_side_length - 0.5 * cell_side_length
    dz = (o_point.z - z_min) % cell_side_length - 0.5 * cell_side_length

    # Construct grid
    x_centers = np.arange(x_min + dx, x_max, cell_side_length) + 0.5 * cell_side_length
    z_centers = np.arange(z_min + dz, z_max, cell_side_length) + 0.5 * cell_side_length
    x, z = np.meshgrid(x_centers, z_centers, indexing="ij")
    points = np.c_[x.ravel(), z.ravel()]

    # Mask by LCFS polygon
    inside = polygon_path.contains_points(points)
    points = points[inside]

    # Volumes: toroidal rotation of each square cell
    d_volume = 2 * np.pi * points[:, 0] * cell_side_length**2

    return points[:, 0], points[:, 1], d_volume

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.")