Merge pull request #23 from SteSeg/neutronics

added neutronics tools

.github/workflows/ci.yml

@@ -8,20 +8,31 @@ jobs:
     steps:
       - uses: actions/checkout@v4
 
-      - name: Set up Python
-        uses: actions/setup-python@v2
+      - name: Set up Conda
+        uses: conda-incubator/setup-miniconda@v2
         with:
-          python-version: "3.10"
+          activate-environment: myenv
+          miniforge-version: latest
+          # use-mamba: true
+          channels: conda-forge
+
+      - name: Create Conda environment
+        shell: bash -l {0}
+        run: |
+          conda install -c conda-forge openmc>=0.14.0
 
       - name: Install dependencies
+        shell: bash -l {0}
         run: |
-          pip install .[tests]
+          python -m pip install .[tests]
 
       - name: Run tests
+        shell: bash -l {0}
         run: |
           python -m pytest test/ --cov libra_toolbox --cov-report xml --cov-report term
 
       - name: Run notebook examples
+        shell: bash -l {0}
         run: |
           jupyter nbconvert --to python --execute docs/examples/*.ipynb
 

libra_toolbox/__init__.py

@@ -1,2 +1,3 @@
 from . import tritium
 from . import neutron_detection
+from . import neutronics

libra_toolbox/_version.py

@@ -0,0 +1,16 @@
+# file generated by setuptools_scm
+# don't change, don't track in version control
+TYPE_CHECKING = False
+if TYPE_CHECKING:
+    from typing import Tuple, Union
+    VERSION_TUPLE = Tuple[Union[int, str], ...]
+else:
+    VERSION_TUPLE = object
+
+version: str
+__version__: str
+__version_tuple__: VERSION_TUPLE
+version_tuple: VERSION_TUPLE
+
+__version__ = version = '0.1.dev87+geef7af6.d20241209'
+__version_tuple__ = version_tuple = (0, 1, 'dev87', 'geef7af6.d20241209')

libra_toolbox/neutronics/__init__.py

@@ -0,0 +1 @@
+from .neutron_source import *

libra_toolbox/neutronics/neutron_source.py

@@ -0,0 +1,77 @@
+# building the MIT-VaultLab neutron generator
+# angular and energy distribution
+
+from pathlib import Path
+from collections.abc import Iterable
+import pandas as pd
+import numpy as np
+import h5py
+import openmc
+
+
+def A325_generator_diamond(center=(0, 0, 0), reference_uvw=(0, 0, 1)) -> Iterable[openmc.IndependentSource]:
+    """
+    Builds the MIT-VaultLab A-325 neutron generator in OpenMC
+    with data tabulated from John Ball and Shon Mackie characterization
+    via diamond detectors
+
+    Parameters
+    ----------
+    center : tuple, optional
+        coordinate position of the source (it is a point source),
+        by default (0, 0, 0)
+    reference_uvw : tuple, optional
+        direction for the polar angle (tuple or list of versors)
+    it is the same for the openmc.PolarAzimuthal class
+    more specifically, polar angle = 0 is the direction of the D accelerator
+    towards the Zr-T target, by default (0, 0, 1)
+
+    Returns
+    -------
+        list of openmc neutron sources with angular and energy distribution 
+        and total strength of 1
+    """
+
+    filename = "A325_generator_diamond.h5"
+    filename = str(Path(__file__).parent) / Path(filename)
+
+    with h5py.File(filename, "r") as source:
+        df = pd.DataFrame(source["values/table"][()]).drop(columns='index')
+        # energy values
+        energies = np.array(df["Energy (MeV)"]) * 1e6
+        # angle column names
+        angles = df.columns[1:]
+        # angular bins in [0, pi)
+        pbins = np.cos([np.deg2rad(float(a)) for a in angles] + [np.pi])
+        spectra = [np.array(df[col]) for col in angles]
+
+    # yield values for strengths
+    yields = np.sum(spectra, axis=-1) * np.diff(pbins)
+    yields /= np.sum(yields)
+
+    # azimuthal values
+    phi = openmc.stats.Uniform(a=0, b=2 * np.pi)
+
+    all_sources = []
+    for i, angle in enumerate(pbins[:-1]):
+
+        mu = openmc.stats.Uniform(a=pbins[i + 1], b=pbins[i])
+
+        space = openmc.stats.Point(center)
+        angle = openmc.stats.PolarAzimuthal(mu=mu, phi=phi, reference_uvw=reference_uvw)
+        energy = openmc.stats.Tabular(
+            energies, spectra[i], interpolation="linear-linear"
+        )
+        strength = yields[i]
+
+        my_source = openmc.Source(
+            space=space,
+            angle=angle,
+            energy=energy,
+            strength=strength,
+            particle="neutron",
+        )
+
+        all_sources.append(my_source)
+
+    return all_sources

pyproject.toml

@@ -15,10 +15,13 @@ description = "Design and analysis tools for LIBRA project"
 license = {file = "LICENSE"}
 requires-python = ">=3.6"
 dynamic = ["version"]
-dependencies = ["numpy", "pint", "scipy", "matplotlib", "sympy", "pandas"]
+dependencies = ["numpy", "pint", "scipy", "matplotlib", "sympy", "pandas", "openmc>=0.14.0", "h5py"]
 
 [project.optional-dependencies]
 tests = ["pytest>=5.4.3", "pytest-cov", "nbconvert", "ipykernel"]
 
 [tool.setuptools_scm]
-write_to = "libra_toolbox/_version.py"
+write_to = "libra_toolbox/_version.py"
+
+[tool.setuptools.package-data]
+mvng_source = ["*.h5"]

test/neutronics/test_neutron_source.py

@@ -0,0 +1,10 @@
+from libra_toolbox.neutronics.neutron_source import *
+from collections.abc import Iterable
+
+def test_get_avg_neutron_rate():
+
+    source = A325_generator_diamond((0, 0, 0), (0, 0, 1))
+
+    assert isinstance(source, Iterable)
+    for s in source:
+        assert isinstance(s, openmc.IndependentSource)