Improved energy distributions including DD, TT and DT reactions

.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

README.md

@@ -8,9 +8,9 @@ This python-based package offers a collection of pre-built [OpenMC](https://gith
 
 OpenMC is required to use this package.
 
-To install openmc-plasma-source, simply run:
+To install openmc_plasma_source, simply run:
 ```
-pip install openmc-plasma-source
+pip install openmc_plasma_source
 ```
 
 ## Usage
@@ -24,26 +24,27 @@ Each source has its own strength (or probability that a neutron spawns in this l
 The equations implemented here are taken from [this paper](https://doi.org/10.1016/j.fusengdes.2012.02.025).
 
 ```python
-from openmc_plasma_source import TokamakSource
+from openmc_plasma_source import tokamak_source
 
-my_source = 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
-  ).make_openmc_sources()
+    fuel={"D": 0.5, "T": 0.5},
+)
 ```
 
 For a more complete example check out the [example script](https://github.com/fusion-energy/openmc-plasma-source/blob/main/examples/tokamak_source_example.py).
@@ -57,13 +58,13 @@ For a more complete example check out the [example script](https://github.com/fu
 Create a ring source with temperature distribution of a 2000 eV plasma.
 
 ```python
-from openmc_plasma_source import FusionRingSource
+from openmc_plasma_source import fusion_ring_source
 
-my_plasma = FusionRingSource(
-    angles = (0., 6.28318530718),  # input is in radians
-    radius = 400,  # units in cm
-    temperature = 20000.,  # ion temperature in eV
-    fuel='DT'  # or 'DD'
+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},
 )
 ```
 ### Point Source
@@ -72,12 +73,12 @@ Create a point source with temperature distribution of a 2000 eV plasma.
 
 
 ```python
-from openmc_plasma_source import FusionPointSource
+from openmc_plasma_source import fusion_point_source
 
-my_plasma = FusionPointSource(
-    coordinate = (0, 0, 0),
-    temperature = 20000.,  # ion temperature in eV
-    fuel = 'DT'  # or 'DD'
+my_source = fusion_point_source(
+    coordinate=(0, 0, 0),
+    temperature=20000.0,
+    fuel={"D": 0.09, "T": 0.91},  # note this is mainly tritium fuel so that TT reactions are more likely
 )
 ```
 
@@ -86,5 +87,5 @@ my_plasma = FusionPointSource(
 To run the tests, simply run
 
 ```
-pytest tests/
+pytest tests
 ```

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 = 20000
+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, "DD")
+
+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 = 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, "DD")
+
+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,14 @@
 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.09,
+        "T": 0.91,
+    },  # note this is mainly tritium fuel so that TT reactions are more likely
+)
 
 # Tell OpenMC we're going to use our custom source
 settings = openmc.Settings()
@@ -20,7 +29,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=2000000,  # 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()
@@ -25,6 +27,19 @@
 # tell OpenMC we're going to use our custom source
 settings.source = my_source
 
-model = openmc.model.Model(materials=None, geometry=geometry, settings=settings)
+model = openmc.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()