Improved energy distribution using NeSST

.github/workflows/ci.yml

@@ -31,3 +31,6 @@ jobs:
           python examples/point_source_example.py
           python examples/ring_source_example.py
           python examples/tokamak_source_example.py
+          python examples/plot_tokamak_ion_temperature.py
+          python examples/plot_tokamak_neutron_source_density.py
+          python examples/plot_tokamak_neutron_source_strengths.py

.gitignore

@@ -135,8 +135,9 @@ dmypy.json
 # vim swap files
 *.swp
 
-# image files created by tests
+# image files created by tests and examples
 tests/*.png
+*.png
 
 # openmc output files
 *.h5

examples/plot_tokamak_ion_temperature.py

@@ -0,0 +1,42 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from openmc_plasma_source import tokamak_convert_a_alpha_to_R_Z, tokamak_ion_temperature
+
+sample_size = 1000
+minor_radius = 292.258
+major_radius = 906
+
+# create a sample of (a, alpha) coordinates
+a = np.random.random(sample_size) * minor_radius
+alpha = np.random.random(sample_size) * 2 * np.pi
+
+temperatures = tokamak_ion_temperature(
+    r=a,
+    mode="L",
+    pedestal_radius=0.8 * minor_radius,
+    ion_temperature_pedestal=6.09,
+    ion_temperature_centre=45.9,
+    ion_temperature_beta=2,
+    ion_temperature_peaking_factor=8.06,
+    ion_temperature_separatrix=0.1,
+    major_radius=major_radius,
+)
+
+RZ = tokamak_convert_a_alpha_to_R_Z(
+    a=a,
+    alpha=alpha,
+    shafranov_factor=0.44789,
+    minor_radius=minor_radius,
+    major_radius=major_radius,
+    triangularity=0.270,
+    elongation=1.557,
+)
+
+plt.scatter(RZ[0], RZ[1], c=temperatures)
+plt.gca().set_aspect("equal")
+plt.xlabel("R [cm]")
+plt.ylabel("Z [cm]")
+plt.colorbar(label="Ion temperature [eV]")
+
+plt.savefig("tokamak_source_ion_temperature.png")
+print("written tokamak_source_ion_temperature.png")

examples/plot_tokamak_neutron_source_density.py

@@ -0,0 +1,62 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from openmc_plasma_source import (
+    tokamak_ion_temperature,
+    tokamak_convert_a_alpha_to_R_Z,
+    tokamak_neutron_source_density,
+    tokamak_ion_density,
+)
+
+sample_size = 20000
+minor_radius = 292.258
+major_radius = 906
+mode = "L"
+ion_density_centre = 45.9
+
+# create a sample of (a, alpha) coordinates
+a = np.random.random(sample_size) * minor_radius
+alpha = np.random.random(sample_size) * 2 * np.pi
+
+temperatures = tokamak_ion_temperature(
+    r=a,
+    mode=mode,
+    pedestal_radius=0.8 * minor_radius,
+    ion_temperature_pedestal=6.09,
+    ion_temperature_centre=ion_density_centre,
+    ion_temperature_beta=2,
+    ion_temperature_peaking_factor=8.06,
+    ion_temperature_separatrix=0.1,
+    major_radius=major_radius,
+)
+
+densities = tokamak_ion_density(
+    mode=mode,
+    ion_density_centre=ion_density_centre,
+    ion_density_peaking_factor=1,
+    ion_density_pedestal=1.09e20,
+    major_radius=major_radius,
+    pedestal_radius=0.8 * minor_radius,
+    ion_density_separatrix=3e19,
+    r=a,
+)
+
+neutron_source_density = tokamak_neutron_source_density(densities, temperatures)
+
+RZ = tokamak_convert_a_alpha_to_R_Z(
+    a=a,
+    alpha=alpha,
+    shafranov_factor=0.44789,
+    minor_radius=minor_radius,
+    major_radius=major_radius,
+    triangularity=0.270,
+    elongation=1.557,
+)
+
+plt.scatter(RZ[0], RZ[1], c=neutron_source_density)
+plt.gca().set_aspect("equal")
+plt.xlabel("R [cm]")
+plt.ylabel("Z [cm]")
+plt.colorbar(label="neutron source density")
+
+plt.savefig("tokamak_source_neutron_source_density.png")
+print("written tokamak_source_neutron_source_density.png")

examples/plot_tokamak_neutron_source_strengths.py

@@ -0,0 +1,64 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from openmc_plasma_source import (
+    tokamak_ion_temperature,
+    tokamak_convert_a_alpha_to_R_Z,
+    tokamak_neutron_source_density,
+    tokamak_ion_density,
+)
+
+sample_size = 2000
+minor_radius = 292.258
+major_radius = 906
+mode = "L"
+ion_density_centre = 45.9
+
+# create a sample of (a, alpha) coordinates
+a = np.random.random(sample_size) * minor_radius
+alpha = np.random.random(sample_size) * 2 * np.pi
+
+temperatures = tokamak_ion_temperature(
+    r=a,
+    mode=mode,
+    pedestal_radius=0.8 * minor_radius,
+    ion_temperature_pedestal=6.09,
+    ion_temperature_centre=ion_density_centre,
+    ion_temperature_beta=2,
+    ion_temperature_peaking_factor=8.06,
+    ion_temperature_separatrix=0.1,
+    major_radius=major_radius,
+)
+
+densities = tokamak_ion_density(
+    mode=mode,
+    ion_density_centre=ion_density_centre,
+    ion_density_peaking_factor=1,
+    ion_density_pedestal=1.09e20,
+    major_radius=major_radius,
+    pedestal_radius=0.8 * minor_radius,
+    ion_density_separatrix=3e19,
+    r=a,
+)
+
+neutron_source_density = tokamak_neutron_source_density(densities, temperatures)
+
+strengths = neutron_source_density / sum(neutron_source_density)
+
+RZ = tokamak_convert_a_alpha_to_R_Z(
+    a=a,
+    alpha=alpha,
+    shafranov_factor=0.44789,
+    minor_radius=minor_radius,
+    major_radius=major_radius,
+    triangularity=0.270,
+    elongation=1.557,
+)
+
+plt.scatter(RZ[0], RZ[1], c=strengths)
+plt.gca().set_aspect("equal")
+plt.xlabel("R [cm]")
+plt.ylabel("Z [cm]")
+plt.colorbar(label="neutron emission strength")
+
+plt.savefig("tokamak_source_neutron_emission_strength.png")
+print("written tokamak_source_neutron_emission_strength.png")

examples/point_source_example.py

@@ -1,6 +1,8 @@
-import openmc
-from openmc_plasma_source import FusionPointSource
 from pathlib import Path
+import numpy as np
+
+import openmc
+from openmc_plasma_source import fusion_point_source
 
 # just making use of a local cross section xml file, replace with your own cross sections or comment out
 openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml"
@@ -11,7 +13,11 @@
 geometry = openmc.Geometry([cell])
 
 # define the source
-my_source = FusionPointSource(coordinate=(0, 0, 0), temperature=20000.0, fuel="DT")
+my_source = fusion_point_source(
+    coordinate=(0, 0, 0),
+    temperature=20000.0,
+    fuel={"D": 0.5, "T": 0.5},
+)
 
 # Tell OpenMC we're going to use our custom source
 settings = openmc.Settings()
@@ -20,7 +26,20 @@
 settings.particles = 1000
 settings.source = my_source
 
-
 model = openmc.model.Model(materials=None, geometry=geometry, settings=settings)
 
 model.run()
+
+
+# optionally if you would like to plot the energy of particles then another package can be used
+# https://github.com/fusion-energy/openmc_source_plotter
+
+from openmc_source_plotter import plot_source_energy
+
+plot = plot_source_energy(
+    this=settings,
+    n_samples=1000000,  # increase this value for a smoother plot
+    energy_bins=np.linspace(0, 16e6, 1000),
+    yaxis_type="log",
+)
+plot.show()

examples/ring_source_example.py

@@ -1,8 +1,9 @@
-import openmc
-from openmc_plasma_source import FusionRingSource
 import math
 from pathlib import Path
 
+import openmc
+from openmc_plasma_source import fusion_ring_source
+
 # just making use of a local cross section xml file, replace with your own cross sections or comment out
 openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml"
 
@@ -12,10 +13,11 @@
 geometry = openmc.Geometry([cell])
 
 # define the source
-my_source = FusionRingSource(
+my_source = fusion_ring_source(
     radius=700,
     angles=(0.0, 2 * math.pi),  # 360deg source
     temperature=20000.0,
+    fuel={"D": 0.5, "T": 0.5},
 )
 
 settings = openmc.Settings()
@@ -28,3 +30,16 @@
 model = openmc.model.Model(materials=None, geometry=geometry, settings=settings)
 
 model.run()
+
+
+# optionally if you would like to plot the location of particles then another package can be used
+# https://github.com/fusion-energy/openmc_source_plotter
+
+# from openmc_source_plotter import plot_source_position
+
+# plot = plot_source_position(
+#     this=settings,
+#     n_samples = 2000,
+# )
+
+# plot.show()

examples/tokamak_source_example.py

@@ -1,44 +1,57 @@
 from pathlib import Path
+
 import openmc
-from openmc_plasma_source import TokamakSource
+
+from openmc_plasma_source import tokamak_source
 
 # just making use of a local cross section xml file, replace with your own cross sections or comment out
 openmc.config["cross_sections"] = Path(__file__).parent.resolve() / "cross_sections.xml"
 
 
 # minimal geometry
-sphere_surface = openmc.Sphere(r=1000.0, boundary_type="vacuum")
+sphere_surface = openmc.Sphere(r=100000.0, boundary_type="vacuum")
 cell = openmc.Cell(region=-sphere_surface)
 geometry = openmc.Geometry([cell])
 
 # create a plasma source
-my_plasma = TokamakSource(
+my_sources = tokamak_source(
     elongation=1.557,
     ion_density_centre=1.09e20,
-    ion_density_peaking_factor=1,
     ion_density_pedestal=1.09e20,
+    ion_density_peaking_factor=1,
     ion_density_separatrix=3e19,
-    ion_temperature_centre=45.9,
+    ion_temperature_centre=45.9e3,
+    ion_temperature_pedestal=6.09e3,
+    ion_temperature_separatrix=0.1e3,
     ion_temperature_peaking_factor=8.06,
-    ion_temperature_pedestal=6.09,
-    ion_temperature_separatrix=0.1,
-    major_radius=9.06,
-    minor_radius=2.92258,
-    pedestal_radius=0.8 * 2.92258,
+    ion_temperature_beta=6,
+    major_radius=906,
+    minor_radius=292.258,
+    pedestal_radius=0.8 * 292.258,
     mode="H",
     shafranov_factor=0.44789,
     triangularity=0.270,
-    ion_temperature_beta=6,
+    fuel={"D": 0.5, "T": 0.5},
 )
 
 # Tell OpenMC we're going to use our custom source
 settings = openmc.Settings()
 settings.run_mode = "fixed source"
 settings.batches = 10
 settings.particles = 1000
-settings.source = my_plasma.sources
+settings.source = my_sources
 
 
 model = openmc.model.Model(materials=None, geometry=geometry, settings=settings)
 
 model.run()
+
+
+# optionally if you would like to plot the direction of particles then another package can be used
+# https://github.com/fusion-energy/openmc_source_plotter
+
+from openmc_source_plotter import plot_source_direction
+
+plot = plot_source_direction(this=settings, n_samples=200)
+
+plot.show()

pyproject.toml

@@ -18,11 +18,12 @@ classifiers = [
 authors = [
   { name="Rémi Delaporte-Mathurin" },
 ]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 keywords = ["python", "neutron", "fusion", "source", "openmc", "energy", "tokamak"]
 dependencies = [
     "numpy>=1.9",
     "matplotlib>=3.2.2",
+    "NeSST>=1.0.0"
 ]
 
 [project.urls]
@@ -35,7 +36,8 @@ write_to = "src/_version.py"
 [project.optional-dependencies]
 tests = [
     "pytest>=5.4.3",
-    "hypothesis"
+    "hypothesis",
+    "NeSST>=1.0.0"
 ]
 
 [tool.setuptools]

src/openmc_plasma_source/__init__.py

@@ -1,7 +1,7 @@
 try:
-    from importlib.metadata import version, PackageNotFoundError
+    from importlib.metadata import PackageNotFoundError, version
 except (ModuleNotFoundError, ImportError):
-    from importlib_metadata import version, PackageNotFoundError
+    from importlib_metadata import PackageNotFoundError, version
 try:
     __version__ = version("openmc_plasma_source")
 except PackageNotFoundError:
@@ -11,6 +11,7 @@
 
 __all__ = ["__version__"]
 
-from .tokamak_source import TokamakSource
-from .ring_source import FusionRingSource
-from .point_source import FusionPointSource
+from .fuel_types import get_neutron_energy_distribution
+from .point_source import fusion_point_source
+from .ring_source import fusion_ring_source
+from .tokamak_source import *