Integrating PyDress to NeSST

.github/workflows/ci.yml

@@ -25,7 +25,7 @@ jobs:
 
       - name: install package with testing dependencies
         run: |
-          pip install .[test]
+          pip install .[test,dress]
           pytest
 
       - name: run examples

README.md

@@ -40,7 +40,7 @@ pip install -e .
 ```
 
 ## Current model specifications:
-- Primary spectrum models for DT, DD and TT
+- Primary spectrum models for DT, DD and TT, including wrapper for Monte Carlo code [pyDRESS](https://github.com/jacob-eri/pydress)
 - Elastic and inelastic (n,2n) processes for D and T
 - ENDF file interface using [ENDF-python](https://github.com/paulromano/endf-python)
 - Relativistic corrections to elastic scattering kernels

docs/guide.rst

@@ -51,6 +51,15 @@ Two shapes are available for the primary spectrum:
    # Ballabio (modified Gaussian) shape
    I_E = nst.QBallabio(Ein, mean, variance)
 
+Additionally, the DRESS Monte Carlo code can be used to generate primary spectra without assuming a particular shape:
+
+.. code-block:: python
+
+   Ein = np.linspace(10e6, 18e6, 200)  # energy array in eV
+
+   # DRESS DT primary spectrum, bin centres provided by Ein
+   I_E = nst.QDRESS_DT(Ein, Tion, n_samples = int(1e6))
+
 Yields and Reactivities
 ~~~~~~~~~~~~~~~~~~~~~~~
 

pyproject.toml

@@ -40,6 +40,9 @@ docs = [
     "myst-parser>=2.0",
     "nbsphinx>=0.9",
 ]
+dress = [
+    "pydress"
+]
 
 [tool.setuptools.package-data]
 "NeSST.data" = ["*.dat", "*.txt","*.csv","*.json"]

src/NeSST/__init__.py

@@ -1,5 +1,10 @@
 ##################
 from .constants import *
 from .core import *
+
+try:
+    from .dress_interface import *
+except ImportError:
+    pass
 from .fitting import *
 from .time_of_flight import *

src/NeSST/core.py

@@ -99,6 +99,44 @@ def QBallabio(Ein: npt.NDArray, mean: float, variance: float) -> npt.NDArray:
     return spec
 
 
+def QDress_DT(Ein: npt.NDArray, T_D: float, T_T: float | None = None, n_samples: int = int(1e6)) -> npt.NDArray:
+    """Calculates the DT primary neutron spectrum using the DRESS Monte Carlo code.
+
+    Args:
+        Ein (numpy.array): array of energy bin centres on which to compute the spectrum (eV)
+        T_D (float): temperature of deuterons in eV
+        T_T (float, optional): temperature of tritons in eV. If None, uses T_D for both.
+        n_samples (int): number of Monte Carlo samples (default 1e6)
+
+    Returns:
+        numpy.array: normalised DT spectrum (1/eV) evaluated at the bin centres of Ein
+    """
+    from NeSST.dress_interface import DRESS_DT_spec, Ecentres_to_edges
+
+    if T_T is None:
+        T_T = T_D
+
+    Ebins, _ = Ecentres_to_edges(Ein)
+    return DRESS_DT_spec(T_D, T_T, n_samples, Ebins)
+
+
+def QDress_DD(Ein: npt.NDArray, Tion: float, n_samples: int = int(1e6)) -> npt.NDArray:
+    """Calculates the DD primary neutron spectrum using the DRESS Monte Carlo code.
+
+    Args:
+        Ein (numpy.array): array of energy bin centres on which to compute the spectrum (eV)
+        Tion (float): temperature of the ions in eV
+        n_samples (int): number of Monte Carlo samples (default 1e6)
+
+    Returns:
+        numpy.array: normalised DD spectrum (1/eV) evaluated at the bin centres of Ein
+    """
+    from NeSST.dress_interface import DRESS_DD_spec, Ecentres_to_edges
+
+    Ebins, _ = Ecentres_to_edges(Ein)
+    return DRESS_DD_spec(Tion, n_samples, Ebins)
+
+
 # TT spectral shape
 def dNdE_TT(E: npt.NDArray, Tion: float) -> npt.NDArray:
     """Calculates the TT primary spectrum with Doppler broadening effect as

src/NeSST/dress_interface.py

@@ -0,0 +1,113 @@
+import numpy as np
+import numpy.typing as npt
+
+try:
+    import dress
+
+    dress_available = True
+except ImportError:
+    dress_available = False
+
+
+def _check_dress():
+    if not dress_available:
+        raise ImportError("pydress is required for DRESS interface functions. Install it with: pip install pydress")
+
+
+def Ecentres_to_edges(Ecentres: npt.NDArray) -> npt.NDArray:
+    """Convert energy bin centres to edges.
+
+    Args:
+        Ecentres (numpy.array): energy bin centres in eV
+
+    Returns:
+        numpy.array: energy bin edges in eV
+    """
+    dE = np.diff(Ecentres)
+    Eedges = np.empty(len(Ecentres) + 1)
+    Eedges[1:-1] = 0.5 * (Ecentres[:-1] + Ecentres[1:])
+    Eedges[0] = Ecentres[0] - 0.5 * dE[0]
+    Eedges[-1] = Ecentres[-1] + 0.5 * dE[-1]
+    return Eedges, Eedges[1:] - Eedges[:-1]
+
+
+def DRESS_DT_spec(T_D: float, T_T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray:
+    """Computes DT primary neutron spectrum using DRESS with separate D and T temperatures.
+
+    Args:
+        T_D (float): temperature of deuterons in eV
+        T_T (float): temperature of tritons in eV
+        n_samples (int): number of Monte Carlo samples
+        bins (numpy.array): energy bin edges in eV
+
+    Returns:
+        numpy.array: normalised DT spectrum (1/eV) on the bin centres
+    """
+    _check_dress()
+    T_D_keV = T_D / 1e3
+    T_T_keV = T_T / 1e3
+    bins_keV = bins / 1e3
+
+    reaction = dress.reactions.DTNHe4Reaction()
+    spec_calc = dress.SpectrumCalculator(reaction, n_samples=n_samples)
+
+    dist_a = dress.dists.MaxwellianDistribution(T_D_keV, spec_calc.reactant_a.particle)
+    dist_b = dress.dists.MaxwellianDistribution(T_T_keV, spec_calc.reactant_b.particle)
+
+    spec_calc.reactant_a.v = dist_a.sample(spec_calc.n_samples)
+    spec_calc.reactant_b.v = dist_b.sample(spec_calc.n_samples)
+
+    spec = spec_calc(bins=bins_keV, normalize=True)
+
+    # Convert from unit-sum to unit-area (1/eV)
+    bins_w = bins[1:] - bins[:-1]
+
+    return spec / bins_w
+
+
+def DRESS_DT_spec_single_T(T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray:
+    """Computes DT primary neutron spectrum using DRESS with a single ion temperature.
+
+    Wrapper around DRESS_DT_spec for the case where D and T share the same temperature.
+
+    Args:
+        T (float): ion temperature in eV
+        n_samples (int): number of Monte Carlo samples
+        bins (numpy.array): energy bin edges in eV
+
+    Returns:
+        numpy.array: normalised DT spectrum (1/eV) on the bin centres
+    """
+    return DRESS_DT_spec(T, T, n_samples, bins)
+
+
+def DRESS_DD_spec(T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray:
+    """Computes DD primary neutron spectrum using DRESS.
+
+    Args:
+        T (float): temperature of deuterons in eV
+        n_samples (int): number of Monte Carlo samples
+        bins (numpy.array): energy bin edges in eV
+
+    Returns:
+        numpy.array: normalised DD spectrum (1/eV) on the bin centres
+    """
+    _check_dress()
+    T_keV = T / 1e3
+    bins_keV = bins / 1e3
+
+    reaction = dress.reactions.DDNHe3Reaction()
+    spec_calc = dress.SpectrumCalculator(reaction, n_samples=n_samples)
+
+    dist_a = dress.dists.MaxwellianDistribution(T_keV, spec_calc.reactant_a.particle)
+    dist_b = dress.dists.MaxwellianDistribution(T_keV, spec_calc.reactant_b.particle)
+
+    spec_calc.reactant_a.v = dist_a.sample(spec_calc.n_samples)
+    spec_calc.reactant_b.v = dist_b.sample(spec_calc.n_samples)
+
+    spec = spec_calc(bins=bins_keV, normalize=True)
+
+    # Convert from unit-sum to unit-area (1/eV)
+    bins_w = bins[1:] - bins[:-1]
+
+    return spec / bins_w

tests/test_dress.py

@@ -0,0 +1,93 @@
+import numpy as np
+import pytest
+
+pytest.importorskip("dress", reason="pydress is required for DRESS tests")
+
+import NeSST as nst
+
+# Number of Monte Carlo samples used in tests - large enough for accurate
+# statistics but small enough to keep test runtime reasonable.
+_N_SAMPLES = int(5e5)
+
+# Tolerance for relative comparison of mean energy against Ballabio (%)
+_MEAN_RTOL = 0.05  # 5 %
+# Tolerance for relative comparison of standard deviation against Ballabio (%)
+_STDDEV_RTOL = 0.05  # 5 %
+
+
+def _spectrum_moments(Ein, spec):
+    """Return the mean and standard deviation of a normalised spectrum."""
+    _, dE = nst.Ecentres_to_edges(Ein)
+    total = np.sum(spec * dE)
+    mean = np.sum(Ein * spec * dE) / total
+    stddev = np.sqrt(np.sum((Ein - mean) ** 2 * spec * dE) / total)
+    return mean, stddev
+
+
+def test_QDress_DT_mean_matches_Ballabio():
+    """DRESS DT mean neutron energy should match Ballabio within 1 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES)
+
+    mean_dress, _ = _spectrum_moments(Ein, spec)
+
+    assert mean_dress - nst.E0_DT == pytest.approx(mean_ball - nst.E0_DT, rel=_MEAN_RTOL)
+
+
+def test_QDress_DT_stddev_matches_Ballabio():
+    """DRESS DT spectral width should match Ballabio within 5 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES)
+
+    _, stddev_dress = _spectrum_moments(Ein, spec)
+
+    assert stddev_dress == pytest.approx(stddev_ball, rel=_STDDEV_RTOL)
+
+
+def test_QDress_DT_separate_temperatures():
+    """Using different D and T temperatures should produce a different spectrum
+    from a single-temperature run."""
+    Tion = 5e3  # 5 keV in eV
+    T_D = 5e3
+    T_T = 10e3  # tritons hotter than deuterons
+
+    mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion)
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+
+    spec_single = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES)
+    spec_diff = nst.QDress_DT(Ein, T_D, T_T, n_samples=_N_SAMPLES)
+
+    # The spectra should not be identical when temperatures differ
+    assert not np.allclose(spec_single, spec_diff)
+
+
+def test_QDress_DD_mean_matches_Ballabio():
+    """DRESS DD mean neutron energy should match Ballabio within 1 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DDprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DD(Ein, Tion, n_samples=_N_SAMPLES)
+
+    mean_dress, _ = _spectrum_moments(Ein, spec)
+
+    assert mean_dress - nst.E0_DD == pytest.approx(mean_ball - nst.E0_DD, rel=_MEAN_RTOL)
+
+
+def test_QDress_DD_stddev_matches_Ballabio():
+    """DRESS DD spectral width should match Ballabio within 5 %."""
+    Tion = 5e3  # 5 keV in eV
+    mean_ball, stddev_ball, _ = nst.DDprimspecmoments(Tion)
+
+    Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500)
+    spec = nst.QDress_DD(Ein, Tion, n_samples=_N_SAMPLES)
+
+    _, stddev_dress = _spectrum_moments(Ein, spec)
+
+    assert stddev_dress == pytest.approx(stddev_ball, rel=_STDDEV_RTOL)