Merge pull request #8 from DanShort12/serialise_plasma_source

Serialise plasma source

.flake8

@@ -0,0 +1,3 @@
+[flake8]
+max-line-length = 88
+exclude = __init__.py

.github/workflows/python.yml

@@ -20,6 +20,16 @@ jobs:
       uses: actions/setup-python@v2
       with:
         python-version: ${{ matrix.python-version }}
+    - name: Install OpenMC
+      run: |
+        sudo apt-get install -y g++ cmake libhdf5-dev
+        git clone --recurse-submodules https://github.com/openmc-dev/openmc.git
+        cd openmc
+        git checkout
+        mkdir build && cd build
+        cmake ..
+        make
+        sudo make install
     - name: Install plasma source
       run: |
         pip install -r requirements-develop.txt

.gitignore

@@ -39,4 +39,11 @@ dist/
 # Python
 __pycache__/
 *.py[cod]
+.venv
 
+# Outputs
+*.xml
+*.h5
+
+# VS Code
+.vscode

README.md

@@ -8,42 +8,106 @@ The plasma source is based on a paper by [C. Fausser et al](https://www.scienced
 
 ## Installation
 
+### Installing from PyPI
+
 ```pip install parametric_plasma_source```
 
+### Installing from source
+
+Installation of the parametric plasma source from source requires cmake to build the underlying C++ code. This can be obtained from
+your OS's package manager by e.g. `sudo apt-get install cmake` or from cmake source.
+
+If you intend to develop the code then it is recommended to work in a virtual environment.
+
+The requirements for developing the code can be installed by running:
+
+```pip install -r requirements-develop.txt```
+
+The package can be built and installed in editable mode by:
+
+```pip install -e .```
+
 ## Usage
 
-The parametric plasma source can be imported an used in Python 3 in the following manner.
+The parametric plasma source can be sampled either directly in Python 3 or sampled in an OpenMC simulation.
+
+For a better understanding of the varibles take a look at the [C. Fausser et al](https://www.sciencedirect.com/science/article/pii/S0920379612000853) paper.
+
+### Sampling in Python
+
+The parametric plasma source can be imported an used in Python 3 in the following manner:
 
 ```[python]
-from parametric_plasma_source import Plasma
-my_plasma = Plasma(major_radius=6,
-                   minor_radius=1.5,
-                   elongation = 2.0
-                   triangularity = 0.55)
-my_plasma.export_plasma_source('custom_openmc_plasma_source.so')
+from parametric_plasma_source import PlasmaSource
+from random import random
+
+plasma_params = {
+    "elongation": 1.557,
+    "ion_density_origin": 1.09e20,
+    "ion_density_peaking_factor": 1,
+    "ion_density_pedestal": 1.09e20,
+    "ion_density_separatrix": 3e19,
+    "ion_temperature_origin": 45.9,
+    "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,
+    "plasma_id": 1,
+    "shafranov_shift": 0.44789,
+    "triangularity": 0.270,
+    "ion_temperature_beta": 6,
+}
+
+my_plasma = PlasmaSource(**plasma_params)
+sample = my_plasma.sample([random(), random(), random(), random(), random(), random(), random(), random()])
+particle_x, particle_y, particle_z = sample[0], sample[1], sample[2]
+particle_x_dir, particle_y_dir, particle_z_dir = sample[3], sample[4], sample[5]
+particle_energy_mev = sample[6]
 ```
 
-In the above example the major_radius, minor_radius, elongation and triangularity while the other varibles are kept as the default values.
+### Sampling in OpenMC
+
+The parametric plasma source also contains a plugin library for OpenMC to allow the source to be sampled in an OpenMC simulation.
 
-There are a number of additional arguments that can be passed to the Plasma class on construction. Units are in SI (e.g. meters not cm)
+When using the OpenMC sampling the inputs must be provided in meters where applicable (the sampling will convert to cm).
 
 ```[python]
-ion_density_pedistal = 1.09e+20
-ion_density_seperatrix = 3e+19
-ion_density_origin = 1.09e+20
-ion_temperature_pedistal = 6.09
-ion_temperature_seperatrix = 0.1
-ion_temperature_origin = 45.9
-pedistal_radius = 0.8
-ion_density_peaking_factor = 1
-ion_temperature_peaking_factor = 8.06
-minor_radius = 1.56
-major_radius = 2.5
-elongation = 2.0
-triangularity = 0.55
-shafranov_shift = 0.0
-number_of_bins = 100
-plasma_type = 1
+from parametric_plasma_source import PlasmaSource, SOURCE_SAMPLING_PATH
+import openmc
+
+plasma_params = {
+    "elongation": 1.557,
+    "ion_density_origin": 1.09e20,
+    "ion_density_peaking_factor": 1,
+    "ion_density_pedestal": 1.09e20,
+    "ion_density_separatrix": 3e19,
+    "ion_temperature_origin": 45.9,
+    "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,
+    "plasma_id": 1,
+    "shafranov_shift": 0.44789,
+    "triangularity": 0.270,
+    "ion_temperature_beta": 6,
+}
+
+my_plasma = PlasmaSource(**plasma_params)
+settings = openmc.Settings()
+settings.run_mode = "fixed source"
+settings.batches = 10
+settings.particles = 1000
+source = openmc.Source()
+source.library = SOURCE_SAMPLING_PATH
+source.parameters = str(my_plasma)
+settings.source = source
+settings.export_to_xml()
 ```
 
-For a better understanding of the varibles take a look at the [C. Fausser et al](https://www.sciencedirect.com/science/article/pii/S0920379612000853) paper.
+## Running Tests
+
+The tests are run by executing `pytest tests` from within your virtual environment.

data/example_source.txt

@@ -0,0 +1 @@
+major_radius=4.5, minor_radius=1.5, elongation=2, triangularity=0.55, shafranov_shift=0, pedestal_radius=1.2, ion_density_pedestal=1.09e+20, ion_density_separatrix=3e+19, ion_density_origin=1.09e+20, ion_density_alpha=1, ion_temperature_pedestal=6.09, ion_temperature_separatrix=0.1, ion_temperature_origin=45.9, ion_temperature_alpha=8.06, ion_temperature_beta=6, plasma_type=plasma, plasma_id=1, number_of_bins=100, minimum_toroidal_angle=0, maximum_toroidal_angle=360

examples/build_and_run_model.py

@@ -0,0 +1,74 @@
+"""
+OpenMC Plasma Source example.
+
+An example of how to generate a settings file containing the parameterised
+plasma source sampling library and how to generate the parameterisation in the
+OpenMC settings.xml file.
+"""
+
+from parametric_plasma_source import PlasmaSource, SOURCE_SAMPLING_PATH
+import openmc
+
+plasma_params = {
+    "elongation": 1.557,
+    "ion_density_origin": 1.09e20,
+    "ion_density_peaking_factor": 1,
+    "ion_density_pedestal": 1.09e20,
+    "ion_density_separatrix": 3e19,
+    "ion_temperature_origin": 45.9,
+    "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,
+    "plasma_id": 1,
+    "shafranov_shift": 0.44789,
+    "triangularity": 0.270,
+    "ion_temperature_beta": 6,
+}
+
+plasma = PlasmaSource(**plasma_params)
+
+# Create a single material
+iron = openmc.Material()
+iron.set_density("g/cm3", 5.0)
+iron.add_element("Fe", 1.0)
+mats = openmc.Materials([iron])
+
+# Create a 5 cm x 5 cm box filled with iron
+cells = []
+inner_box1 = openmc.ZCylinder(r=600.0)
+inner_box2 = openmc.ZCylinder(r=1400.0)
+outer_box = openmc.model.rectangular_prism(4000.0, 4000.0, boundary_type="vacuum")
+cells += [openmc.Cell(fill=iron, region=-inner_box1)]
+cells += [openmc.Cell(fill=None, region=+inner_box1 & -inner_box2)]
+cells += [openmc.Cell(fill=iron, region=+inner_box2 & outer_box)]
+geometry = openmc.Geometry(cells)
+
+# Tell OpenMC we're going to use our custom source
+settings = openmc.Settings()
+settings.run_mode = "fixed source"
+settings.batches = 10
+settings.particles = 1000
+source = openmc.Source()
+source.library = SOURCE_SAMPLING_PATH
+source.parameters = str(plasma)
+settings.source = source
+
+# Finally, define a mesh tally so that we can see the resulting flux
+mesh = openmc.RegularMesh()
+mesh.lower_left = (-2000.0, -2000.0)
+mesh.upper_right = (2000.0, 2000.0)
+mesh.dimension = (50, 50)
+
+tally = openmc.Tally()
+tally.filters = [openmc.MeshFilter(mesh)]
+tally.scores = ["flux"]
+tallies = openmc.Tallies([tally])
+
+model = openmc.model.Model(
+    materials=mats, geometry=geometry, settings=settings, tallies=tallies
+)
+
+model.run()

examples/plot_flux.py

@@ -0,0 +1,18 @@
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+import openmc
+
+mpl.use("TkAgg")
+
+# Get the flux from the statepoint
+with openmc.StatePoint("statepoint.10.h5") as sp:
+    flux = np.log(sp.get_tally(scores=["flux"]).mean)
+    flux.shape = (50, 50)
+
+# Plot the flux
+fig, ax = plt.subplots()
+ax.imshow(flux, origin="lower", extent=(-2000.0, 2000.0, -2000.0, 2000.0))
+ax.set_xlabel("x [cm]")
+ax.set_ylabel("y [cm]")
+plt.show()

parametric_plasma_source/__init__.py

@@ -0,0 +1,9 @@
+import os
+
+PLASMA_SOURCE_PATH = os.path.dirname(__file__)
+SOURCE_SAMPLING_PATH = os.sep.join([PLASMA_SOURCE_PATH, "source_sampling.so"])
+
+try:
+    from .plasma_source import *
+except ImportError:
+    print("The plasma_source module could not be found. Please compile before using.")