diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..fd1c60c --- /dev/null +++ b/.flake8 @@ -0,0 +1,3 @@ +[flake8] +max-line-length = 88 +exclude = __init__.py \ No newline at end of file diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index c46ced0..2c3af7a 100644 --- a/.github/workflows/python.yml +++ b/.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 diff --git a/.gitignore b/.gitignore index e6e4fb4..1fadb9e 100644 --- a/.gitignore +++ b/.gitignore @@ -39,4 +39,11 @@ dist/ # Python __pycache__/ *.py[cod] +.venv +# Outputs +*.xml +*.h5 + +# VS Code +.vscode diff --git a/README.md b/README.md index 84d791e..4cc5517 100644 --- a/README.md +++ b/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. diff --git a/data/example_source.txt b/data/example_source.txt new file mode 100644 index 0000000..ee90fb2 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/examples/build_and_run_model.py b/examples/build_and_run_model.py new file mode 100644 index 0000000..8cf8516 --- /dev/null +++ b/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() diff --git a/examples/plot_flux.py b/examples/plot_flux.py new file mode 100644 index 0000000..13381b8 --- /dev/null +++ b/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() diff --git a/parametric_plasma_source/Makefile b/parametric_plasma_source/Makefile deleted file mode 100644 index a060d7b..0000000 --- a/parametric_plasma_source/Makefile +++ /dev/null @@ -1,52 +0,0 @@ - -# Makefile to build dynamic sources for OpenMC -# this assumes that your source sampling filename is -# source_sampling.cpp -# -# you can add fortran, c and cpp dependencies to this source -# adding them in FC_DEPS, C_DEPS, and CXX_DEPS accordingly - -ifeq ($(FC), f77) -FC = gfortran -endif - -default: all - -ALL = source_sampling -# add your fortran depencies here -FC_DEPS = -# add your c dependencies here -C_DEPS = -# add your cpp dependencies here -CPP_DEPS = plasma_source.cpp -DEPS = $(FC_DEPS) $(C_DEPS) $(CPP_DEPS) -OPT_LEVEL = -O3 -FFLAGS = $(OPT_LEVEL) -fPIC -C_FLAGS = -fPIC -CXX_FLAGS = -fPIC - - #this directory will need changing if your openmc path is different -OPENMC_DIR = /opt/openmc -OPENMC_INC_DIR = $(OPENMC_DIR)/include -OPENMC_LIB_DIR = $(OPENMC_DIR)/lib -# setting the so name is important -LINK_FLAGS = $(OPT_LEVEL) -Wall -Wl,-soname,source_sampling.so -# setting shared is important -LINK_FLAGS += -L$(OPENMC_LIB_DIR) -lopenmc -lgfortran -fPIC -shared -OPENMC_INCLUDES = -I$(OPENMC_INC_DIR) -I$(OPENMC_DIR)/vendor/pugixml - -all: $(ALL) - -source_sampling: $(DEPS) - $(CXX) source_sampling.cpp $(DEPS) $(OPENMC_INCLUDES) $(LINK_FLAGS) -o $@.so -# make any fortran objects -%.o : %.F90 - $(FC) -c $(FFLAGS) $*.F90 -o $@ -# make any c objects -%.o : %.c - $(CC) -c $(FFLAGS) $*.c -o $@ -#make any cpp objects -%.o : %.cpp - $(CXX) -c $(FFLAGS) $*.cpp -o $@ -clean: - rm -rf *.o *.mod diff --git a/parametric_plasma_source/__init__.py b/parametric_plasma_source/__init__.py index e69de29..0d71a4b 100644 --- a/parametric_plasma_source/__init__.py +++ b/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.") diff --git a/parametric_plasma_source/build_python.py b/parametric_plasma_source/build_python.py deleted file mode 100644 index 6ee7e39..0000000 --- a/parametric_plasma_source/build_python.py +++ /dev/null @@ -1,56 +0,0 @@ -import os - -plasma_source_cpp = "plasma_source.cpp" -plasma_source_hpp = "plasma_source.hpp" -source_sampling_cpp = "source_sampling.cpp" -make_file = "Makefile" - -py_file = "source.py" - -def write_disclaimer(py_file_path): - disclaimer = "\"\"\"\n" - disclaimer += "Auto-generated file. To edit, run `python build_python.py`.\n" - disclaimer += "This will need to be run whenever a new C++ version is made available.\n" - disclaimer += "\"\"\"\n\n" - with open(py_file_path, "a+") as f_py: - f_py.write(disclaimer) - -def generate_variable(cpp_file_path, variable_name, py_file_path, end_str="\n\n"): - with open(cpp_file_path, "r") as f_cpp: - lines = f_cpp.readlines() - py_lines = variable_name - py_lines += " = (\n\"\"\"" - py_lines += "".join(lines) - py_lines += "\"\"\"\n)" - py_lines += end_str - with open(py_file_path, "a+") as f_py: - f_py.write(py_lines) - -def build(source_dir="."): - output_path = os.path.join(source_dir, py_file) - - if os.path.exists(output_path): - os.remove(output_path) - - write_disclaimer(output_path) - generate_variable( - os.path.join(source_dir, source_sampling_cpp), - "source_sampling_cpp", - output_path - ) - generate_variable( - os.path.join(source_dir, plasma_source_cpp), - "plasma_source_cpp", - output_path - ) - generate_variable( - os.path.join(source_dir, plasma_source_hpp), - "plasma_source_hpp", - output_path - ) - generate_variable( - os.path.join(source_dir, make_file), - "make_file", - output_path, - "\n" - ) diff --git a/parametric_plasma_source/compile.sh b/parametric_plasma_source/compile.sh deleted file mode 100644 index cd74266..0000000 --- a/parametric_plasma_source/compile.sh +++ /dev/null @@ -1,5 +0,0 @@ -#!/bin/bash - -make clean - -make \ No newline at end of file diff --git a/parametric_plasma_source/plasma.py b/parametric_plasma_source/plasma.py deleted file mode 100644 index 7cab511..0000000 --- a/parametric_plasma_source/plasma.py +++ /dev/null @@ -1,380 +0,0 @@ - -import os -from pathlib import Path -import tempfile -import shutil -from .source import source_sampling_cpp, plasma_source_cpp, plasma_source_hpp, make_file - - -class Plasma(): - - def __init__(self, - elongation=2., - major_radius=450, - minor_radius=150, - single_null=True, - triangularity=0.55, - ion_density_pedistal=1.09e20, - ion_density_seperatrix=3e19, - ion_density_origin=1.09e20, - ion_temperature_pedistal=6.09, - ion_temperature_seperatrix=0.1, - ion_temperature_origin=45.9, - pedistal_radius=120, # 0.8 * minor_radius - ion_density_peaking_factor=1, - ion_temperature_peaking_factor=8.06, - ion_temperature_beta=6.0, - shafranov_shift=0.0, - number_of_bins=100, - plasma_type=1, - openmc_install_directory = '/opt/openmc/' - ): - - # properties needed for plasma shapes - self.elongation = elongation - self.major_radius = major_radius - self.minor_radius = minor_radius - self.single_null = single_null - self.triangularity = triangularity - self.ion_density_pedistal = ion_density_pedistal # ions per m^3 - self.ion_density_seperatrix = ion_density_seperatrix - self.ion_density_origin = ion_density_origin - self.ion_temperature_pedistal = ion_temperature_pedistal - self.ion_temperature_seperatrix = ion_temperature_seperatrix - self.ion_temperature_origin = ion_temperature_origin - self.pedistal_radius = pedistal_radius # pedistal major rad - self.ion_density_peaking_factor = ion_density_peaking_factor - self.ion_temperature_peaking_factor = ion_temperature_peaking_factor - self.ion_temperature_beta = ion_temperature_beta - self.shafranov_shift = shafranov_shift - self.number_of_bins = number_of_bins - self.plasma_type = plasma_type # 0 = L mode anything else H/A mode - self.openmc_install_directory = openmc_install_directory - - self.plasma_source_cpp_file = plasma_source_cpp - self.plasma_source_hpp_file = plasma_source_hpp - self.source_sampling_cpp_file = source_sampling_cpp - self.plasma_make_file = make_file - - @property - def openmc_install_directory(self): - return self._openmc_install_directory - - @openmc_install_directory.setter - def openmc_install_directory(self, value): - if Path(value).exists() == False: - raise ValueError('openmc_install_directory is out of range') - else: - self._openmc_install_directory = value - - @property - def plasma_type(self): - return self._plasma_type - - @plasma_type.setter - def plasma_type(self, plasma_type): - if plasma_type < 0: - raise ValueError('plasma_type is out of range') - else: - self._plasma_type = plasma_type - - @property - def number_of_bins(self): - return self._number_of_bins - - @number_of_bins.setter - def number_of_bins(self, number_of_bins): - if number_of_bins < 0: - raise ValueError('number_of_bins is out of range') - else: - self._number_of_bins = number_of_bins - - @property - def shafranov_shift(self): - return self._shafranov_shift - - @shafranov_shift.setter - def shafranov_shift(self, shafranov_shift): - if shafranov_shift < 0: - raise ValueError('shafranov_shift is out of range') - else: - self._shafranov_shift = shafranov_shift - - @property - def ion_temperature_peaking_factor(self): - return self._ion_temperature_peaking_factor - - @ion_temperature_peaking_factor.setter - def ion_temperature_peaking_factor(self, ion_temperature_peaking_factor): - if ion_temperature_peaking_factor < 0: - raise ValueError('ion_temperature_peaking_factor is out of range') - else: - self._ion_temperature_peaking_factor = ion_temperature_peaking_factor - - @property - def ion_temperature_beta(self): - return self._ion_temperature_beta - - @ion_temperature_beta.setter - def ion_temperature_beta(self, ion_temperature_beta): - if ion_temperature_beta < 0: - raise ValueError('ion_temperature_beta is out of range') - else: - self._ion_temperature_beta = ion_temperature_beta - - @property - def ion_density_peaking_factor(self): - return self._ion_density_peaking_factor - - @ion_density_peaking_factor.setter - def ion_density_peaking_factor(self, ion_density_peaking_factor): - if ion_density_peaking_factor < 0: - raise ValueError('ion_density_peaking_factor is out of range') - else: - self._ion_density_peaking_factor = ion_density_peaking_factor - - @property - def pedistal_radius(self): - return self._pedistal_radius - - @pedistal_radius.setter - def pedistal_radius(self, pedistal_radius): - if pedistal_radius < 0: - raise ValueError('pedistal_radius is out of range') - else: - self._pedistal_radius = pedistal_radius - - @property - def ion_temperature_origin(self): - return self._ion_temperature_origin - - @ion_temperature_origin.setter - def ion_temperature_origin(self, ion_temperature_origin): - if ion_temperature_origin < 0: - raise ValueError('ion_temperature_origin is out of range') - else: - self._ion_temperature_origin = ion_temperature_origin - - @property - def ion_temperature_seperatrix(self): - return self._ion_temperature_seperatrix - - @ion_temperature_seperatrix.setter - def ion_temperature_seperatrix(self, ion_temperature_seperatrix): - if ion_temperature_seperatrix < 0: - raise ValueError('ion_temperature_seperatrix is out of range') - else: - self._ion_temperature_seperatrix = ion_temperature_seperatrix - - @property - def ion_temperature_pedistal(self): - return self._ion_temperature_pedistal - - @ion_temperature_pedistal.setter - def ion_temperature_pedistal(self, ion_temperature_pedistal): - if ion_temperature_pedistal < 0: - raise ValueError('ion_temperature_pedistal is out of range') - else: - self._ion_temperature_pedistal = ion_temperature_pedistal - - @property - def ion_density_origin(self): - return self._ion_density_origin - - @ion_density_origin.setter - def ion_density_origin(self, ion_density_origin): - if ion_density_origin < 0: - raise ValueError('ion_density_origin is out of range') - else: - self._ion_density_origin = ion_density_origin - - @property - def ion_density_seperatrix(self): - return self._ion_density_seperatrix - - @ion_density_seperatrix.setter - def ion_density_seperatrix(self, ion_density_seperatrix): - if ion_density_seperatrix < 0: - raise ValueError('ion_density_seperatrix is out of range') - else: - self._ion_density_seperatrix = ion_density_seperatrix - - @property - def triangularity(self): - return self._triangularity - - @triangularity.setter - def triangularity(self, triangularity): - if triangularity > 2000 or triangularity < -2000: - raise ValueError('triangularity is out of range') - else: - self._triangularity = triangularity - - @property - def single_null(self): - return self._single_null - - @single_null.setter - def single_null(self, single_null): - if type(single_null) != bool : - raise ValueError('single_null must be True or False') - else: - self._single_null = single_null - - @property - def minor_radius(self): - return self._minor_radius - - @minor_radius.setter - def minor_radius(self, minor_radius): - self._minor_radius = minor_radius - - @property - def major_radius(self): - return self._major_radius - - @major_radius.setter - def major_radius(self, major_radius): - if major_radius > 2000 or major_radius < 1: - raise ValueError('major_radius is out of range') - else: - self._major_radius = major_radius - - @property - def elongation(self): - return self._elongation - - @elongation.setter - def elongation(self, elongation): - if elongation > 4 or elongation < 0: - raise ValueError('elongation is out of range') - else: - self._elongation = elongation - - @property - def ion_density_pedistal(self): - return self._ion_density_pedistal - - @ion_density_pedistal.setter - def ion_density_pedistal(self, ion_density_pedistal): - if ion_density_pedistal > 10e22 or ion_density_pedistal < 1e4: - raise ValueError('ion_density_pedistal is out of range') - else: - self._ion_density_pedistal = ion_density_pedistal - - def replace_variable_value( - self, input_strings, variable_name, new_value, is_cpp=True - ): - """Replaces the value assigned to a variable with the provided value - - param: input_strings: The list of strings to search for the variable. - type input_strings: List[str] - - param: variable_name: The name of the variable to search for. - type variable_name: str - - param: new_value: The value to set the variable to. - type new_value: Union[str, float, int] - - param is_cpp: Whether the variable is a C++ variable. - Optional, by default True. - type is_cpp: bool - - raises ValueError: variable_name was not found in input_strings. - """ - if is_cpp: - strs_to_find = ("const", variable_name, "=", ";") - else: - strs_to_find = (variable_name, "=") - for idx, string in enumerate(input_strings): - if all(str_to_find in string for str_to_find in strs_to_find): - equals_idx = string.find("=") - if equals_idx >= 0: - if is_cpp: - new_value = f"{new_value};" - input_strings[idx] = string.replace( - input_strings[idx][equals_idx:], f" = {new_value}" - ) - break - else: - file_content = "\n".join(input_strings) - raise ValueError( - f"{variable_name} string not found in {file_content}" - ) - - def export_plasma_source(self, output_filename): - """Writes and compiles custom plasma source for the reactor - :param output_folder: the output folder where the .so complied plasma source will be created - :type output_folder: str - ... - :return: filename of the compiled source - :rtype: str - """ - - if self.openmc_install_directory is None: - raise ValueError('directory must be set to create .so file') - - temp_folder = Path(tempfile.mkdtemp()) - - Path(output_filename).parent.mkdir(parents=True, exist_ok=True) - - plasma_make_file_lines = self.plasma_make_file.split("\n") - self.replace_variable_value( - plasma_make_file_lines, "OPENMC_DIR", self.openmc_install_directory, False - ) - - with open(temp_folder/'Makefile', "w") as text_file: - text_file.write("\n".join(plasma_make_file_lines)) - - with open(temp_folder/'plasma_source.cpp', "w") as text_file: - text_file.write(self.plasma_source_cpp_file) - - with open(temp_folder/'plasma_source.hpp', "w") as text_file: - text_file.write(self.plasma_source_hpp_file) - - plasma_variables = { - "ion_density_pedistal": self.ion_density_pedistal, - "ion_density_seperatrix": self.ion_density_seperatrix, - "ion_density_origin": self.ion_density_origin, - "ion_temperature_pedistal": self.ion_temperature_pedistal, - "ion_temperature_seperatrix": self.ion_temperature_seperatrix, - "ion_temperature_origin": self.ion_temperature_origin, - "pedistal_radius": self.pedistal_radius, - "ion_density_peaking_factor": self.ion_density_peaking_factor, - "ion_temperature_peaking_factor": self.ion_temperature_peaking_factor, - "ion_temperature_beta": self.ion_temperature_beta, - "minor_radius": self.minor_radius, - "major_radius": self.major_radius, - "elongation": self.elongation, - "triangularity": self.triangularity, - "shafranov_shift": self.shafranov_shift, - "number_of_bins": self.number_of_bins, - "plasma_type": self.plasma_type, - } - - source_sampling_cpp_lines = self.source_sampling_cpp_file.split("\n") - [ - self.replace_variable_value(source_sampling_cpp_lines, *variable) - for variable in plasma_variables.items() - ] - - with open(temp_folder/'source_sampling.cpp', "w") as text_file: - text_file.write("\n".join(source_sampling_cpp_lines)) - - cwd = os.getcwd() - os.chdir(Path(temp_folder)) - - os.system('make clean') - os.system('make') - - os.chdir(cwd) - shutil.move(temp_folder/'source_sampling.so', output_filename) - print('parametric plasma source compiled and saved to ', output_filename) - shutil.rmtree(temp_folder) - - return output_filename - - -if __name__ == "__main__": - p = Plasma() - p.export_plasma_source("source_sampling.so") diff --git a/parametric_plasma_source/plasma_source.cpp b/parametric_plasma_source/plasma_source.cpp index b0b61b0..4af4941 100644 --- a/parametric_plasma_source/plasma_source.cpp +++ b/parametric_plasma_source/plasma_source.cpp @@ -1,5 +1,7 @@ -#include +#include +#include #include +#include #include #include "plasma_source.hpp" #include @@ -11,25 +13,25 @@ PlasmaSource::PlasmaSource() {} // large constructor PlasmaSource::PlasmaSource(const double ion_density_ped, const double ion_density_sep, - const double ion_density_origin, const double ion_temp_ped, - const double ion_temp_sep, const double ion_temp_origin, - const double pedistal_rad, const double ion_density_peak, - const double ion_temp_peak, const double ion_temp_beta, - const double minor_radius, const double major_radius, - const double elongation, const double triangularity, - const double shafranov, const std::string plasma_type, - const int plasma_id, const int number_of_bins, - const double min_toroidal_angle, - const double max_toroidal_angle ) { + const double ion_density_origin, const double ion_temp_ped, + const double ion_temp_sep, const double ion_temp_origin, + const double pedestal_rad, const double ion_density_peak, + const double ion_temp_peak, const double ion_temp_beta, + const double minor_radius, const double major_radius, + const double elongation, const double triangularity, + const double shafranov, const std::string plasma_type, + const int plasma_id, const int number_of_bins, + const double min_toroidal_angle, + const double max_toroidal_angle ) { // set params - ionDensityPedistal = ion_density_ped; - ionDensitySeperatrix = ion_density_sep; + ionDensityPedestal = ion_density_ped; + ionDensitySeparatrix = ion_density_sep; ionDensityOrigin = ion_density_origin; - ionTemperaturePedistal = ion_temp_ped; - ionTemperatureSeperatrix = ion_temp_sep; + ionTemperaturePedestal = ion_temp_ped; + ionTemperatureSeparatrix = ion_temp_sep; ionTemperatureOrigin = ion_temp_origin; - pedistalRadius = pedistal_rad; + pedestalRadius = pedestal_rad; ionDensityPeaking = ion_density_peak; ionTemperaturePeaking = ion_temp_peak; ionTemperatureBeta = ion_temp_beta; @@ -41,8 +43,8 @@ PlasmaSource::PlasmaSource(const double ion_density_ped, const double ion_densit plasmaType = plasma_type; plasmaId = plasma_id; numberOfBins = number_of_bins; - minToroidalAngle = min_toroidal_angle/180.*M_PI; - maxToroidalAngle = max_toroidal_angle/180.*M_PI; + minToroidalAngle = (min_toroidal_angle / 180.) * M_PI; + maxToroidalAngle = (max_toroidal_angle / 180.) * M_PI; setup_plasma_source(); } @@ -51,17 +53,17 @@ PlasmaSource::PlasmaSource(const double ion_density_ped, const double ion_densit PlasmaSource::~PlasmaSource(){} // main master sample function -void PlasmaSource::SampleSource(std::array random_numbers, - double &x, - double &y, - double &z, - double &u, - double &v, - double &w, - double &E) { +void PlasmaSource::sample(std::array random_numbers, + double &x, + double &y, + double &z, + double &u, + double &v, + double &w, + double &E) { double radius = 0.; int bin = 0; - sample_source_radial(random_numbers[0],random_numbers[1],radius,bin); + sample_radial(random_numbers[0],random_numbers[1],radius,bin); double r = 0.; convert_rad_to_rz(radius,random_numbers[2],r,z); convert_r_to_xy(r,random_numbers[3],x,y); @@ -74,7 +76,7 @@ void PlasmaSource::SampleSource(std::array random_numbers, /* * sample the pdf src_profile, to generate the sampled minor radius */ -void PlasmaSource::sample_source_radial(double rn_store1, double rn_store2, +void PlasmaSource::sample_radial(double rn_store1, double rn_store2, double &sampled_radius, int &sampled_bin) { for ( int i = 0 ; i < numberOfBins ; i++ ) { @@ -101,7 +103,7 @@ void PlasmaSource::sample_source_radial(double rn_store1, double rn_store2, * sample the energy of the neutrons, updates energy neutron in mev */ void PlasmaSource::sample_energy(const int bin_number, double random_number1, double random_number2, - double &energy_neutron) { + double &energy_neutron) { // generate the normally distributed number const double twopi = 6.28318530718; double sample1 = std::sqrt(-2.0*std::log(random_number1)); @@ -194,18 +196,18 @@ double PlasmaSource::ion_density(const double sample_radius) ion_dens = ionDensityOrigin* (1.0-std::pow(sample_radius/minorRadius,2)); } else { - if(sample_radius <= pedistalRadius) { - ion_dens += ionDensityPedistal; + if(sample_radius <= pedestalRadius) { + ion_dens += ionDensityPedestal; double product; - product = 1.0-std::pow(sample_radius/pedistalRadius,2); + product = 1.0-std::pow(sample_radius/pedestalRadius,2); product = std::pow(product,ionDensityPeaking); - ion_dens += (ionDensityOrigin-ionDensityPedistal)* + ion_dens += (ionDensityOrigin-ionDensityPedestal)* (product); } else { - ion_dens += ionDensitySeperatrix; + ion_dens += ionDensitySeparatrix; double product; - product = ionDensityPedistal - ionDensitySeperatrix; - ion_dens += product*(minorRadius-sample_radius)/(minorRadius-pedistalRadius); + product = ionDensityPedestal - ionDensitySeparatrix; + ion_dens += product*(minorRadius-sample_radius)/(minorRadius-pedestalRadius); } } @@ -225,18 +227,18 @@ double PlasmaSource::ion_temperature(const double sample_radius) (1.0-std::pow(sample_radius/minorRadius, ionTemperaturePeaking)); } else { - if(sample_radius <= pedistalRadius) { - ion_temp += ionTemperaturePedistal; + if(sample_radius <= pedestalRadius) { + ion_temp += ionTemperaturePedestal; double product; - product = 1.0-std::pow(sample_radius/pedistalRadius,ionTemperatureBeta); + product = 1.0-std::pow(sample_radius/pedestalRadius,ionTemperatureBeta); product = std::pow(product,ionTemperaturePeaking); ion_temp += (ionTemperatureOrigin- - ionTemperaturePedistal)*(product); + ionTemperaturePedestal)*(product); } else { - ion_temp += ionTemperatureSeperatrix; + ion_temp += ionTemperatureSeparatrix; double product; - product = ionTemperaturePedistal - ionTemperatureSeperatrix; - ion_temp += product*(minorRadius-sample_radius)/(minorRadius-pedistalRadius); + product = ionTemperaturePedestal - ionTemperatureSeparatrix; + ion_temp += product*(minorRadius-sample_radius)/(minorRadius-pedestalRadius); } } @@ -250,8 +252,8 @@ double PlasmaSource::dt_xs(double ion_temp) { double dt; double c[7]={2.5663271e-18,19.983026,2.5077133e-2, - 2.5773408e-3,6.1880463e-5,6.6024089e-2, - 8.1215505e-3}; + 2.5773408e-3,6.1880463e-5,6.6024089e-2, + 8.1215505e-3}; double u = 1.0-ion_temp*(c[2]+ion_temp*(c[3]-c[4]*ion_temp)) /(1.0+ion_temp*(c[5]+c[6]*ion_temp)); @@ -279,4 +281,70 @@ void PlasmaSource::isotropic_direction(const double random1, return; } +std::string PlasmaSource::to_string() { + std::stringstream result; + result << "major_radius=" << majorRadius << ", "; + result << "minor_radius=" << minorRadius << ", "; + result << "elongation=" << elongation << ", "; + result << "triangularity=" << triangularity << ", "; + result << "shafranov_shift=" << shafranov << ", "; + result << "pedestal_radius=" << pedestalRadius << ", "; + result << "ion_density_pedestal=" << ionDensityPedestal << ", "; + result << "ion_density_separatrix=" << ionDensitySeparatrix << ", "; + result << "ion_density_origin=" << ionDensityOrigin << ", "; + result << "ion_density_alpha=" << ionDensityPeaking << ", "; + result << "ion_temperature_pedestal=" << ionTemperaturePedestal << ", "; + result << "ion_temperature_separatrix=" << ionTemperatureSeparatrix << ", "; + result << "ion_temperature_origin=" << ionTemperatureOrigin << ", "; + result << "ion_temperature_alpha=" << ionTemperaturePeaking << ", "; + result << "ion_temperature_beta=" << ionTemperatureBeta << ", "; + result << "plasma_type=" << plasmaType << ", "; + result << "plasma_id=" << plasmaId << ", "; + result << "number_of_bins=" << numberOfBins << ", "; + result << "minimum_toroidal_angle=" << minToroidalAngle * 180.0 / M_PI << ", "; + result << "maximum_toroidal_angle=" << maxToroidalAngle * 180.0 / M_PI; + return result.str(); +} + +PlasmaSource PlasmaSource::from_string(std::string parameters) { + std::unordered_map parameter_mapping; + + std::stringstream ss(parameters); + std::string parameter; + while (std::getline(ss, parameter, ',')) { + parameter.erase(0, parameter.find_first_not_of(' ')); + std::string key = parameter.substr(0, parameter.find_first_of('=')); + std::string value = parameter.substr(parameter.find_first_of('=') + 1, parameter.length()); + parameter_mapping[key] = value; + } + + double major_radius = std::stod(parameter_mapping["major_radius"]); + double minor_radius = std::stod(parameter_mapping["minor_radius"]); + double elongation = std::stod(parameter_mapping["elongation"]); + double triangularity = std::stod(parameter_mapping["triangularity"]); + double shafranov_shift = std::stod(parameter_mapping["shafranov_shift"]); + double pedestal_radius = std::stod(parameter_mapping["pedestal_radius"]); + double ion_density_pedestal = std::stod(parameter_mapping["ion_density_pedestal"]); + double ion_density_separatrix = std::stod(parameter_mapping["ion_density_separatrix"]); + double ion_density_origin = std::stod(parameter_mapping["ion_density_origin"]); + double ion_temperature_pedestal = std::stod(parameter_mapping["ion_temperature_pedestal"]); + double ion_temperature_separatrix = std::stod(parameter_mapping["ion_temperature_separatrix"]); + double ion_temperature_origin = std::stod(parameter_mapping["ion_temperature_origin"]); + double ion_density_alpha = std::stod(parameter_mapping["ion_density_alpha"]); + double ion_temperature_alpha = std::stod(parameter_mapping["ion_temperature_alpha"]); + double ion_temperature_beta = std::stod(parameter_mapping["ion_temperature_beta"]); + std::string plasma_type = parameter_mapping["plasma_type"]; + int plasma_id = std::stoi(parameter_mapping["plasma_id"]); + int number_of_bins = std::stoi(parameter_mapping["number_of_bins"]); + double min_toroidal_angle = std::stod(parameter_mapping["minimum_toroidal_angle"]); + double max_toroidal_angle = std::stod(parameter_mapping["maximum_toroidal_angle"]); + + PlasmaSource plasma_source = PlasmaSource( + ion_density_pedestal, ion_density_separatrix, ion_density_origin, ion_temperature_pedestal, ion_temperature_separatrix, + ion_temperature_origin, pedestal_radius, ion_density_alpha, ion_temperature_alpha, ion_temperature_beta, minor_radius, + major_radius, elongation, triangularity, shafranov_shift, plasma_type, plasma_id, number_of_bins, min_toroidal_angle, + max_toroidal_angle); + return plasma_source; +} + } // end of namespace diff --git a/parametric_plasma_source/plasma_source.hpp b/parametric_plasma_source/plasma_source.hpp index bf1dbf1..8f81148 100644 --- a/parametric_plasma_source/plasma_source.hpp +++ b/parametric_plasma_source/plasma_source.hpp @@ -1,5 +1,6 @@ #include #include + namespace plasma_source { struct xs_params { @@ -7,116 +8,242 @@ struct xs_params { }; class PlasmaSource { -public: -// constructor -PlasmaSource(); -// destructor -~PlasmaSource(); -// large constructor -PlasmaSource(const double ion_density_ped, const double ion_density_sep, - const double ion_density_origin, const double ion_temp_ped, - const double ion_temp_sep, const double ion_temp_origin, - const double pedistal_rad, const double ion_density_peak, - const double ion_temp_peak, const double ion_temp_beta, - const double minor_radius, const double major_radius, - const double elongation, const double triangularity, - const double shafranov, const std::string plasma_type, - const int plasma_id, const int number_of_bins, - const double min_toroidal_angle = 0.0, - const double max_toridal_angle = 360.); - -// main sample fucnction -void SampleSource(std::array randoms, - double &x, - double &y, - double &z, - double &u, - double &v, - double &w, - double &E); - -/* - * Function to setup the plasma source in the first case. - */ -void setup_plasma_source(); - -/* - * function to calculate the ion density at a specific minor - * radius - */ -double ion_density(const double sample_radius); - -/* - * function to calculate the ion temperature at a specific minor - * radius - */ -double ion_temperature(const double sample_radius); - -/* - * function to determine the value of the dt xs cross sections at - * a specific ion temp - */ -double dt_xs(double ion_temp); - -/* - * sample the source, returns the minor radius sampled - * expects new rn_store every call - */ -void sample_source_radial(double rn_store1, double rn_store2, - double &sampled_radius, - int &sampled_bin); - -/* - * sample the neutron energy in MeV - */ -void sample_energy(const int bin_number, double random_number1, double random_number2, - double &energy_neutron); - -/* - * take the sampled minor radius and convert to cylindrical coordinates - */ -void convert_rad_to_rz(const double minor_sampled, - const double rn_store, - double &radius, - double &height); - -/* - * convert partial cylindrical coords to xyz - */ -void convert_r_to_xy(const double r, const double rn_store, - double &x, double &y); -/* - * get an isotropically direction vector - */ -void isotropic_direction(const double random1, const double random2, - double &u, double &v, double &w); - -private: - std::vector source_profile; - std::vector ion_kt; - - double ionDensityPedistal; - double ionDensitySeperatrix; - double ionDensityOrigin; - double ionTemperaturePedistal; - double ionTemperatureSeperatrix; - double ionTemperatureOrigin; - double pedistalRadius; - double ionDensityPeaking; - double ionTemperaturePeaking; - double ionTemperatureBeta; - double minorRadius; - double majorRadius; - double elongation; - double triangularity; - double shafranov; - double minToroidalAngle; - double maxToroidalAngle; - - std::string plasmaType; - int plasmaId; - double binWidth; - int numberOfBins; + public: + // constructor + PlasmaSource(); + + // destructor + ~PlasmaSource(); + + // large constructor + PlasmaSource(const double ion_density_ped, + const double ion_density_sep, + const double ion_density_origin, + const double ion_temp_ped, + const double ion_temp_sep, + const double ion_temp_origin, + const double pedestal_rad, + const double ion_density_peak, + const double ion_temp_peak, + const double ion_temp_beta, + const double minor_radius, + const double major_radius, + const double elongation, + const double triangularity, + const double shafranov, + const std::string plasma_type, + const int plasma_id, + const int number_of_bins, + const double min_toroidal_angle = 0.0, + const double max_toridal_angle = 360.); + + // main sample fucnction + void sample(std::array randoms, + double &x, + double &y, + double &z, + double &u, + double &v, + double &w, + double &E); + + /* + * Function to setup the plasma source in the first case. + */ + void setup_plasma_source(); + + /* + * function to calculate the ion density at a specific minor + * radius + */ + double ion_density(const double sample_radius); + + /* + * function to calculate the ion temperature at a specific minor + * radius + */ + double ion_temperature(const double sample_radius); + + /* + * function to determine the value of the dt xs cross sections at + * a specific ion temp + */ + double dt_xs(double ion_temp); + + /* + * sample the source, returns the minor radius sampled + * expects new rn_store every call + */ + void sample_radial(double rn_store1, + double rn_store2, + double &sampled_radius, + int &sampled_bin); + + /* + * sample the neutron energy in MeV + */ + void sample_energy(const int bin_number, double random_number1, double random_number2, + double &energy_neutron); + + /* + * take the sampled minor radius and convert to cylindrical coordinates + */ + void convert_rad_to_rz(const double minor_sampled, + const double rn_store, + double &radius, + double &height); + + /* + * convert partial cylindrical coords to xyz + */ + void convert_r_to_xy(const double r, const double rn_store, double &x, double &y); + + /* + * get an isotropically direction vector + */ + void isotropic_direction(const double random1, + const double random2, + double &u, + double &v, + double &w); + + /* + * get a key-value pair string representation of this instance of the source + */ + std::string to_string(); + + /* + * create a new source from the provided key-value pair string representation + * of the source + */ + static PlasmaSource from_string(std::string parameters); + + /* + * Getter for pedestal ion density + */ + const double &getIonDensityPedestal() const { return ionDensityPedestal; } + + /* + * Getter for separatrix ion density + */ + const double &getIonDensitySeparatrix() const { return ionDensitySeparatrix; } + + /* + * Getter for origin ion density + */ + const double &getIonDensityOrigin() const { return ionDensityOrigin; } + + /* + * Getter for pedestal ion temperature + */ + const double &getIonTemperaturePedestal() const { return ionTemperaturePedestal; } + /* + * Getter for separatrix ion temperature + */ + const double &getIonTemperatureSeparatrix() const { return ionTemperatureSeparatrix; } + + /* + * Getter for origin ion temperature + */ + const double &getIonTemperatureOrigin() const { return ionTemperatureOrigin; } + + /* + * Getter for pedestal radius + */ + const double &getPedestalRadius() const { return pedestalRadius; } + + /* + * Getter for ion density peaking factor + */ + const double &getIonDensityPeaking() const { return ionDensityPeaking; } + + /* + * Getter for ion temperature peaking factor + */ + const double &getIonTemperaturePeaking() const { return ionTemperaturePeaking; } + + /* + * Getter for ion temperature beta factor + */ + const double &getIonTemperatureBeta() const { return ionTemperatureBeta; } + + /* + * Getter for minor radius + */ + const double &getMinorRadius() const { return minorRadius; } + + /* + * Getter for major radius + */ + const double &getMajorRadius() const { return majorRadius; } + + /* + * Getter for elongation + */ + const double &getElongation() const { return elongation; } + + /* + * Getter for triangularity + */ + const double &getTriangularity() const { return triangularity; } + + /* + * Getter for shafranov shift + */ + const double &getShafranov() const { return shafranov; } + + /* + * Getter for minimum toroidal angle + */ + const double &getMinToroidalAngle() const { return minToroidalAngle; } + + /* + * Getter for maximum toroidal angle + */ + const double &getMaxToroidalAngle() const { return maxToroidalAngle; } + + /* + * Getter for plasma type + */ + const std::string &getPlasmaType() const { return plasmaType; } + + /* + * Getter for plasma identifier + */ + const int &getPlasmaId() const { return plasmaId; } + + /* + * Getter for number of bins + */ + const int &getNumberOfBins() const { return numberOfBins; } + + private: + std::vector source_profile; + std::vector ion_kt; + + double ionDensityPedestal; + double ionDensitySeparatrix; + double ionDensityOrigin; + double ionTemperaturePedestal; + double ionTemperatureSeparatrix; + double ionTemperatureOrigin; + double pedestalRadius; + double ionDensityPeaking; + double ionTemperaturePeaking; + double ionTemperatureBeta; + double minorRadius; + double majorRadius; + double elongation; + double triangularity; + double shafranov; + double minToroidalAngle; + double maxToroidalAngle; + + std::string plasmaType; + int plasmaId; + double binWidth; + int numberOfBins; }; -}// end of namespace \ No newline at end of file + +}// end of namespace diff --git a/parametric_plasma_source/plasma_source_pybind.cpp b/parametric_plasma_source/plasma_source_pybind.cpp index c220d2a..218e68e 100644 --- a/parametric_plasma_source/plasma_source_pybind.cpp +++ b/parametric_plasma_source/plasma_source_pybind.cpp @@ -1,4 +1,5 @@ #include +#include #include "plasma_source.hpp" namespace py = pybind11; @@ -13,36 +14,97 @@ PYBIND11_MODULE(plasma_source, m) { const double &, const double &, const double &, const double &, const double &, const double &, const double &, const std::string &, const int &, const int &, const double &, const double &>(), - py::arg("ion_density_pedistal")=1.09e20, - py::arg("ion_density_seperatrix")=3e19, - py::arg("ion_density_origin")=1.09e20, - py::arg("ion_temperature_pedistal")=6.09, - py::arg("ion_temperature_seperatrix")=0.1, - py::arg("ion_temperature_origin")=45.9, - py::arg("pedistal_radius")=0.8, - py::arg("ion_density_peaking_factor")=1.0, - py::arg("ion_temperature_peaking_factor")=8.06, - py::arg("ion_temperature_beta")=6.0, - py::arg("minor_radius")=1.5, - py::arg("major_radius")=4.5, - py::arg("elongation")=2.0, - py::arg("triangularity")=0.55, - py::arg("shafranov_shift")=0.0, + py::arg("ion_density_pedestal"), + py::arg("ion_density_separatrix"), + py::arg("ion_density_origin"), + py::arg("ion_temperature_pedestal"), + py::arg("ion_temperature_separatrix"), + py::arg("ion_temperature_origin"), + py::arg("pedestal_radius"), + py::arg("ion_density_peaking_factor"), + py::arg("ion_temperature_peaking_factor"), + py::arg("ion_temperature_beta"), + py::arg("minor_radius"), + py::arg("major_radius"), + py::arg("elongation"), + py::arg("triangularity"), + py::arg("shafranov_shift"), py::arg("plasma_type")="plasma", py::arg("plasma_id")=1, py::arg("number_of_bins")=100, - py::arg("min_toroidal_angle") = 0.0, - py::arg("max_toridal_angle") = 360.0) + py::arg("min_toroidal_angle")=0.0, + py::arg("max_toridal_angle")=360.0) + .def_property_readonly("ion_density_pedestal", &ps::PlasmaSource::getIonDensityPedestal) + .def_property_readonly("ion_density_separatrix", &ps::PlasmaSource::getIonDensitySeparatrix) + .def_property_readonly("ion_density_origin", &ps::PlasmaSource::getIonDensityOrigin) + .def_property_readonly("ion_temperature_pedestal", &ps::PlasmaSource::getIonTemperaturePedestal) + .def_property_readonly("ion_temperature_separatrix", &ps::PlasmaSource::getIonTemperatureSeparatrix) + .def_property_readonly("ion_temperature_origin", &ps::PlasmaSource::getIonTemperatureOrigin) + .def_property_readonly("pedestal_radius", &ps::PlasmaSource::getPedestalRadius) + .def_property_readonly("ion_density_peaking_factor", &ps::PlasmaSource::getIonDensityPeaking) + .def_property_readonly("ion_temperature_peaking_factor", &ps::PlasmaSource::getIonTemperaturePeaking) + .def_property_readonly("ion_temperature_beta", &ps::PlasmaSource::getIonTemperatureBeta) + .def_property_readonly("minor_radius", &ps::PlasmaSource::getMinorRadius) + .def_property_readonly("major_radius", &ps::PlasmaSource::getMajorRadius) + .def_property_readonly("elongation", &ps::PlasmaSource::getElongation) + .def_property_readonly("triangularity", &ps::PlasmaSource::getTriangularity) + .def_property_readonly("shafranov_shift", &ps::PlasmaSource::getShafranov) + .def_property_readonly("plasma_type", &ps::PlasmaSource::getPlasmaType) + .def_property_readonly("plasma_id", &ps::PlasmaSource::getPlasmaId) + .def_property_readonly("number_of_bins", &ps::PlasmaSource::getNumberOfBins) + .def_property_readonly("min_toroidal_angle", &ps::PlasmaSource::getMinToroidalAngle) + .def_property_readonly("max_toridal_angle", &ps::PlasmaSource::getMaxToroidalAngle) .def("ion_density", &ps::PlasmaSource::ion_density, "Calculate the ion density at a specific minor radius", py::arg("minor_radius")) .def("ion_temperature", &ps::PlasmaSource::ion_temperature, - "calculate the ion temperature at a specific minor radius", + "Calculate the ion temperature at a specific minor radius", py::arg("minor_radius")) .def("dt_xs", &ps::PlasmaSource::dt_xs, - "determine the value of the dt xs cross sections at a specific ion temperature", - py::arg("ion_temperature")); + "Determine the value of the dt xs cross sections at a specific ion temperature", + py::arg("ion_temperature")) + .def("sample", + [](ps::PlasmaSource &source, std::array random_numbers) { + double x, y, z; + double u, v, w; + double e; + source.sample(random_numbers, x, y, z, u, v, w, e); + return py::make_tuple(x, y, z, u, v, w, e); + }, + R"( + Sample the source. + + Parameters + ---------- + random_numbers : List[float[8]] + The set of eight random numbers to use when sampling the source. + + Returns + ------- + x : float + The initial x-coordinate of the particle. + y : float + The initial y-coordinate of the particle. + z : float + The initial z-coordinate of the particle. + x : float + The initial x direction of the particle. + y : float + The initial y direction of the particle. + z : float + The initial z direction of the particle. + e : float + The initial energy of the particle. + )", + py::arg("random_numbers")) + .def("__str__", + (std::string (ps::PlasmaSource::*)()) &ps::PlasmaSource::to_string, + "Write out the PlasmaSource as a key-value string") + .def_static("from_string", + &ps::PlasmaSource::from_string, + "Load the PlasmaSource from a key-value string.", + py::arg("parameters")); } diff --git a/parametric_plasma_source/source.py b/parametric_plasma_source/source.py deleted file mode 100644 index f744dd4..0000000 --- a/parametric_plasma_source/source.py +++ /dev/null @@ -1,567 +0,0 @@ -""" -Auto-generated file. To edit, run `python build_python.py`. -This will need to be run whenever a new C++ version is made available. -""" - -source_sampling_cpp = ( -"""#include -#include "openmc/random_lcg.h" -#include "openmc/source.h" -#include "openmc/particle.h" -#include "plasma_source.hpp" - - -// Spherical tokamak SOURCE -// units are in SI units -const double ion_density_pedistal = 1.09e+20; // ions per m^3 -const double ion_density_seperatrix = 3e+19; -const double ion_density_origin = 1.09e+20; -const double ion_temperature_pedistal = 6.09; -const double ion_temperature_seperatrix = 0.1; -const double ion_temperature_origin = 45.9; -const double ion_density_peaking_factor = 1; -const double ion_temperature_peaking_factor = 8.06; // check alpha or beta value from paper -const double ion_temperature_beta = 6.0; -const double minor_radius = 1.56; // metres -const double major_radius = 2.5; // metres -const double pedistal_radius = 0.8 * minor_radius; // pedistal minor rad in metres -const double elongation = 2.0; -const double triangularity = 0.55; -const double shafranov_shift = 0.0; //metres -const std::string name = "parametric_plasma_source"; -const int number_of_bins = 100; -const int plasma_type = 1; // 1 is default; //0 = L mode anything else H/A mode - - -plasma_source::PlasmaSource source = plasma_source::PlasmaSource(ion_density_pedistal, - ion_density_seperatrix, - ion_density_origin, - ion_temperature_pedistal, - ion_temperature_seperatrix, - ion_temperature_origin, - pedistal_radius, - ion_density_peaking_factor, - ion_temperature_peaking_factor, - ion_temperature_beta, - minor_radius, - major_radius, - elongation, - triangularity, - shafranov_shift, - name, - plasma_type, - number_of_bins, - 0.0, - 360.0); - -// you must have external C linkage here otherwise -// dlopen will not find the file -extern "C" openmc::Particle::Bank sample_source(uint64_t* seed) { - openmc::Particle::Bank particle; - // wgt - particle.particle = openmc::Particle::Type::neutron; - particle.wgt = 1.0; - // position - - std::array randoms = {openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed)}; - - double u,v,w,E; - source.SampleSource(randoms,particle.r.x,particle.r.y,particle.r.z, - u,v,w,E); - - particle.r.x *= 100.; - particle.r.y *= 100.; - particle.r.z *= 100.; - - // particle.E = 14.08e6; - particle.E = E*1e6; // convert from MeV -> eV - - particle.u = {u, - v, - w}; - - - particle.delayed_group = 0; - return particle; -} - - -""" -) - -plasma_source_cpp = ( -"""#include -#include -#include -#include "plasma_source.hpp" -#include -#include "openmc/random_lcg.h" - -#define RANDOM openmc::prn() - -namespace plasma_source { - -// default constructor -PlasmaSource::PlasmaSource() {} - -// large constructor -PlasmaSource::PlasmaSource(const double ion_density_ped, const double ion_density_sep, - const double ion_density_origin, const double ion_temp_ped, - const double ion_temp_sep, const double ion_temp_origin, - const double pedistal_rad, const double ion_density_peak, - const double ion_temp_peak, const double ion_temp_beta, - const double minor_radius, const double major_radius, - const double elongation, const double triangularity, - const double shafranov, const std::string plasma_type, - const int plasma_id, const int number_of_bins, - const double min_toroidal_angle, - const double max_toroidal_angle ) { - - // set params - ionDensityPedistal = ion_density_ped; - ionDensitySeperatrix = ion_density_sep; - ionDensityOrigin = ion_density_origin; - ionTemperaturePedistal = ion_temp_ped; - ionTemperatureSeperatrix = ion_temp_sep; - ionTemperatureOrigin = ion_temp_origin; - pedistalRadius = pedistal_rad; - ionDensityPeaking = ion_density_peak; - ionTemperaturePeaking = ion_temp_peak; - ionTemperatureBeta = ion_temp_beta; - minorRadius = minor_radius; - majorRadius = major_radius; - this->elongation = elongation; - this->triangularity = triangularity; - this->shafranov = shafranov; - plasmaType = plasma_type; - plasmaId = plasma_id; - numberOfBins = number_of_bins; - minToroidalAngle = min_toroidal_angle/180.*M_PI; - maxToroidalAngle = max_toroidal_angle/180.*M_PI; - - setup_plasma_source(); -} - -// destructor -PlasmaSource::~PlasmaSource(){} - -// main master sample function -void PlasmaSource::SampleSource(std::array random_numbers, - double &x, - double &y, - double &z, - double &u, - double &v, - double &w, - double &E) { - double radius = 0.; - int bin = 0; - sample_source_radial(random_numbers[0],random_numbers[1],radius,bin); - double r = 0.; - convert_rad_to_rz(radius,random_numbers[2],r,z); - convert_r_to_xy(r,random_numbers[3],x,y); - sample_energy(bin,random_numbers[4],random_numbers[5],E); - isotropic_direction(random_numbers[6],random_numbers[7], - u,v,w); - -} - -/* - * sample the pdf src_profile, to generate the sampled minor radius - */ -void PlasmaSource::sample_source_radial(double rn_store1, double rn_store2, - double &sampled_radius, int &sampled_bin) { - - for ( int i = 0 ; i < numberOfBins ; i++ ) { - if ( rn_store1 <= source_profile[i] ) { - if ( i > 0 ) { - sampled_radius = (float(i-1)*(binWidth)) + (binWidth*(rn_store2)); - sampled_bin = i; - return; - } else { - sampled_radius = binWidth*(rn_store2); - sampled_bin = i; - return; - } - } - } - - std::cerr << "error" << std::endl; - std::cerr << "Sample position greater than plasma radius" << std::endl; - exit(1); - return; -} - -/* - * sample the energy of the neutrons, updates energy neutron in mev - */ -void PlasmaSource::sample_energy(const int bin_number, double random_number1, double random_number2, - double &energy_neutron) { - // generate the normally distributed number - const double twopi = 6.28318530718; - double sample1 = std::sqrt(-2.0*std::log(random_number1)); - double sample2 = cos(twopi*(random_number2)); - energy_neutron = (5.59/2.35)*(ion_kt[bin_number])*sample1*sample2; - energy_neutron += 14.08; - // test energy limit - // if (energy_neutron < 15.5){energy_neutron = 15.5} else {} - return; -} - -/* - * convert the sampled radius to an rz coordinate by using plasma parameters - */ -void PlasmaSource::convert_rad_to_rz( const double minor_sampled, - const double rn_store, - double &radius, double &height) -{ - const double twopi = 6.28318530718; - - double alpha = twopi*(rn_store); - - double shift = shafranov*(1.0-std::pow(minor_sampled/(minorRadius),2)); - - radius = majorRadius + minor_sampled*cos(alpha+(triangularity*sin(alpha))) + shift; - height = elongation*minor_sampled*sin(alpha); - - - return; -} - - -/* - * convert rz_to_xyz - */ -void PlasmaSource::convert_r_to_xy(const double r, const double rn_store, - double &x, double &y) - -{ - double toroidal_extent = maxToroidalAngle - minToroidalAngle; - double toroidal_angle = toroidal_extent*rn_store + minToroidalAngle; - x = r*sin(toroidal_angle); - y = r*cos(toroidal_angle); - return; -} - -/* - * sets up the cumulatitive probability profile - * on the basis of the ion temp and ion density - * this portion is deterministic - */ -void PlasmaSource::setup_plasma_source() -{ - double ion_d; // ion density - double ion_t; // ion temp - - std::vector src_strength; // the source strength, n/m3 - double r; - - binWidth = minorRadius/float(numberOfBins); - double total = 0.0; // total source strength - - for (int i = 0 ; i < numberOfBins ; i++) { - r = binWidth * float(i); - ion_d = ion_density(r); - ion_t = ion_temperature(r); - src_strength.push_back(std::pow(ion_d,2)*dt_xs(ion_t)); - ion_kt.push_back(sqrt(ion_t/1000.0)); // convert to sqrt(MeV) - total += src_strength[i]; - } - - // normalise the source profile - double sum = 0 ; - for ( int i = 0 ; i < numberOfBins ; i++) { - sum += src_strength[i]; - source_profile.push_back(sum/total); - } - return; -} - -/* - * function that returns the ion density given the - * given the critical plasma parameters - */ -double PlasmaSource::ion_density(const double sample_radius) -{ - double ion_dens = 0.0; - - if( plasmaId == 0 ) { - ion_dens = ionDensityOrigin* - (1.0-std::pow(sample_radius/minorRadius,2)); - } else { - if(sample_radius <= pedistalRadius) { - ion_dens += ionDensityPedistal; - double product; - product = 1.0-std::pow(sample_radius/pedistalRadius,2); - product = std::pow(product,ionDensityPeaking); - ion_dens += (ionDensityOrigin-ionDensityPedistal)* - (product); - } else { - ion_dens += ionDensitySeperatrix; - double product; - product = ionDensityPedistal - ionDensitySeperatrix; - ion_dens += product*(minorRadius-sample_radius)/(minorRadius-pedistalRadius); - } - } - - return ion_dens; -} - -/* - * function that returns the ion density given the - * given the critical plasma parameters - */ -double PlasmaSource::ion_temperature(const double sample_radius) -{ - double ion_temp = 0.0; - - if( plasmaId == 0 ) { - ion_temp = ionTemperatureOrigin* - (1.0-std::pow(sample_radius/minorRadius, - ionTemperaturePeaking)); - } else { - if(sample_radius <= pedistalRadius) { - ion_temp += ionTemperaturePedistal; - double product; - product = 1.0-std::pow(sample_radius/pedistalRadius,ionTemperatureBeta); - product = std::pow(product,ionTemperaturePeaking); - ion_temp += (ionTemperatureOrigin- - ionTemperaturePedistal)*(product); - } else { - ion_temp += ionTemperatureSeperatrix; - double product; - product = ionTemperaturePedistal - ionTemperatureSeperatrix; - ion_temp += product*(minorRadius-sample_radius)/(minorRadius-pedistalRadius); - } - } - - return ion_temp; -} - -/* - * returns the dt cross section for a given ion temp - */ -double PlasmaSource::dt_xs(double ion_temp) -{ - double dt; - double c[7]={2.5663271e-18,19.983026,2.5077133e-2, - 2.5773408e-3,6.1880463e-5,6.6024089e-2, - 8.1215505e-3}; - - double u = 1.0-ion_temp*(c[2]+ion_temp*(c[3]-c[4]*ion_temp)) - /(1.0+ion_temp*(c[5]+c[6]*ion_temp)); - - dt = c[0]/(std::pow(u,5./6.)*std::pow(ion_temp,2./3.0)); - dt *= exp(-1.*c[1]*std::pow(u/ion_temp,1./3.)); - - return dt; -} - -/* - * returns the dt cross section for a given ion temp - */ -void PlasmaSource::isotropic_direction(const double random1, - const double random2, - double &u, double &v, - double &w) { - double t = 2*M_PI*random1; - double p = acos(1. - 2.*random2); - - u = sin(p)*cos(t); - v = sin(p)*sin(t); - w = cos(p); - - return; -} - -} // end of namespace -""" -) - -plasma_source_hpp = ( -"""#include -#include -namespace plasma_source { - -struct xs_params { - double c[7]; -}; - -class PlasmaSource { -public: -// constructor -PlasmaSource(); -// destructor -~PlasmaSource(); -// large constructor -PlasmaSource(const double ion_density_ped, const double ion_density_sep, - const double ion_density_origin, const double ion_temp_ped, - const double ion_temp_sep, const double ion_temp_origin, - const double pedistal_rad, const double ion_density_peak, - const double ion_temp_peak, const double ion_temp_beta, - const double minor_radius, const double major_radius, - const double elongation, const double triangularity, - const double shafranov, const std::string plasma_type, - const int plasma_id, const int number_of_bins, - const double min_toroidal_angle = 0.0, - const double max_toridal_angle = 360.); - -// main sample fucnction -void SampleSource(std::array randoms, - double &x, - double &y, - double &z, - double &u, - double &v, - double &w, - double &E); - -/* - * Function to setup the plasma source in the first case. - */ -void setup_plasma_source(); - -/* - * function to calculate the ion density at a specific minor - * radius - */ -double ion_density(const double sample_radius); - -/* - * function to calculate the ion temperature at a specific minor - * radius - */ -double ion_temperature(const double sample_radius); - -/* - * function to determine the value of the dt xs cross sections at - * a specific ion temp - */ -double dt_xs(double ion_temp); - -/* - * sample the source, returns the minor radius sampled - * expects new rn_store every call - */ -void sample_source_radial(double rn_store1, double rn_store2, - double &sampled_radius, - int &sampled_bin); - -/* - * sample the neutron energy in MeV - */ -void sample_energy(const int bin_number, double random_number1, double random_number2, - double &energy_neutron); - -/* - * take the sampled minor radius and convert to cylindrical coordinates - */ -void convert_rad_to_rz(const double minor_sampled, - const double rn_store, - double &radius, - double &height); - -/* - * convert partial cylindrical coords to xyz - */ -void convert_r_to_xy(const double r, const double rn_store, - double &x, double &y); -/* - * get an isotropically direction vector - */ -void isotropic_direction(const double random1, const double random2, - double &u, double &v, double &w); - -private: - std::vector source_profile; - std::vector ion_kt; - - double ionDensityPedistal; - double ionDensitySeperatrix; - double ionDensityOrigin; - double ionTemperaturePedistal; - double ionTemperatureSeperatrix; - double ionTemperatureOrigin; - double pedistalRadius; - double ionDensityPeaking; - double ionTemperaturePeaking; - double ionTemperatureBeta; - double minorRadius; - double majorRadius; - double elongation; - double triangularity; - double shafranov; - double minToroidalAngle; - double maxToroidalAngle; - - std::string plasmaType; - int plasmaId; - double binWidth; - int numberOfBins; - -}; -}// end of namespace""" -) - -make_file = ( -""" -# Makefile to build dynamic sources for OpenMC -# this assumes that your source sampling filename is -# source_sampling.cpp -# -# you can add fortran, c and cpp dependencies to this source -# adding them in FC_DEPS, C_DEPS, and CXX_DEPS accordingly - -ifeq ($(FC), f77) -FC = gfortran -endif - -default: all - -ALL = source_sampling -# add your fortran depencies here -FC_DEPS = -# add your c dependencies here -C_DEPS = -# add your cpp dependencies here -CPP_DEPS = plasma_source.cpp -DEPS = $(FC_DEPS) $(C_DEPS) $(CPP_DEPS) -OPT_LEVEL = -O3 -FFLAGS = $(OPT_LEVEL) -fPIC -C_FLAGS = -fPIC -CXX_FLAGS = -fPIC - - #this directory will need changing if your openmc path is different -OPENMC_DIR = /opt/openmc -OPENMC_INC_DIR = $(OPENMC_DIR)/include -OPENMC_LIB_DIR = $(OPENMC_DIR)/lib -# setting the so name is important -LINK_FLAGS = $(OPT_LEVEL) -Wall -Wl,-soname,source_sampling.so -# setting shared is important -LINK_FLAGS += -L$(OPENMC_LIB_DIR) -lopenmc -lgfortran -fPIC -shared -OPENMC_INCLUDES = -I$(OPENMC_INC_DIR) -I$(OPENMC_DIR)/vendor/pugixml - -all: $(ALL) - -source_sampling: $(DEPS) - $(CXX) source_sampling.cpp $(DEPS) $(OPENMC_INCLUDES) $(LINK_FLAGS) -o $@.so -# make any fortran objects -%.o : %.F90 - $(FC) -c $(FFLAGS) $*.F90 -o $@ -# make any c objects -%.o : %.c - $(CC) -c $(FFLAGS) $*.c -o $@ -#make any cpp objects -%.o : %.cpp - $(CXX) -c $(FFLAGS) $*.cpp -o $@ -clean: - rm -rf *.o *.mod -""" -) diff --git a/parametric_plasma_source/source_sampling.cpp b/parametric_plasma_source/source_sampling.cpp index f763d4a..49b71a5 100644 --- a/parametric_plasma_source/source_sampling.cpp +++ b/parametric_plasma_source/source_sampling.cpp @@ -1,89 +1,66 @@ -#include +#include + #include "openmc/random_lcg.h" #include "openmc/source.h" #include "openmc/particle.h" #include "plasma_source.hpp" +// defines a class that wraps our PlasmaSource and exposes it to OpenMC. +class SampledSource : public openmc::CustomSource { + private: + // the source that we will sample from + plasma_source::PlasmaSource source; -// Spherical tokamak SOURCE -// units are in SI units -const double ion_density_pedistal = 1.09e+20; // ions per m^3 -const double ion_density_seperatrix = 3e+19; -const double ion_density_origin = 1.09e+20; -const double ion_temperature_pedistal = 6.09; -const double ion_temperature_seperatrix = 0.1; -const double ion_temperature_origin = 45.9; -const double ion_density_peaking_factor = 1; -const double ion_temperature_peaking_factor = 8.06; // check alpha or beta value from paper -const double ion_temperature_beta = 6.0; -const double minor_radius = 1.56; // metres -const double major_radius = 2.5; // metres -const double pedistal_radius = 0.8 * minor_radius; // pedistal minor rad in metres -const double elongation = 2.0; -const double triangularity = 0.55; -const double shafranov_shift = 0.0; //metres -const std::string name = "parametric_plasma_source"; -const int number_of_bins = 100; -const int plasma_type = 1; // 1 is default; //0 = L mode anything else H/A mode + public: + // create a SampledSource as a wrapper for the source that we will sample from + SampledSource(plasma_source::PlasmaSource source) : source(source) { } + // the function that will be exposed in the openmc::CustomSource parent class + // so that the source can be sampled from. + // essentially wraps the sample_source method on the source and populates the + // relevant values in the openmc::Particle::Bank. + openmc::Particle::Bank sample(uint64_t* seed) { + openmc::Particle::Bank particle; + + // random numbers sampled from openmc::prn + std::array randoms = {openmc::prn(seed), + openmc::prn(seed), + openmc::prn(seed), + openmc::prn(seed), + openmc::prn(seed), + openmc::prn(seed), + openmc::prn(seed), + openmc::prn(seed)}; -plasma_source::PlasmaSource source = plasma_source::PlasmaSource(ion_density_pedistal, - ion_density_seperatrix, - ion_density_origin, - ion_temperature_pedistal, - ion_temperature_seperatrix, - ion_temperature_origin, - pedistal_radius, - ion_density_peaking_factor, - ion_temperature_peaking_factor, - ion_temperature_beta, - minor_radius, - major_radius, - elongation, - triangularity, - shafranov_shift, - name, - plasma_type, - number_of_bins, - 0.0, - 360.0); - -// you must have external C linkage here otherwise -// dlopen will not find the file -extern "C" openmc::Particle::Bank sample_source(uint64_t* seed) { - openmc::Particle::Bank particle; - // wgt - particle.particle = openmc::Particle::Type::neutron; - particle.wgt = 1.0; - // position + double u, v, w, E; - std::array randoms = {openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed), - openmc::prn(seed)}; + // sample from the source + this->source.sample(randoms, particle.r.x, particle.r.y, particle.r.z, u, v, w, E); - double u,v,w,E; - source.SampleSource(randoms,particle.r.x,particle.r.y,particle.r.z, - u,v,w,E); + // wgt + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + particle.delayed_group = 0; - particle.r.x *= 100.; - particle.r.y *= 100.; - particle.r.z *= 100.; - - // particle.E = 14.08e6; - particle.E = E*1e6; // convert from MeV -> eV + // position (note units converted m -> cm) + particle.r.x *= 100.; + particle.r.y *= 100.; + particle.r.z *= 100.; - particle.u = {u, - v, - w}; + // energy + particle.E = E * 1e6; // convert from MeV -> eV - - particle.delayed_group = 0; - return particle; -} + // direction + particle.u = { u, v, w }; + return particle; + } +}; +// A function to create a unique pointer to an instance of this class when generated +// via a plugin call using dlopen/dlsym. +// You must have external C linkage here otherwise dlopen will not find the file +extern "C" std::unique_ptr openmc_create_source(std::string parameters) { + plasma_source::PlasmaSource source = plasma_source::PlasmaSource::from_string(parameters); + return std::unique_ptr (new SampledSource(source)); +} diff --git a/requirements-develop.txt b/requirements-develop.txt index 0d73fcd..763b67d 100644 --- a/requirements-develop.txt +++ b/requirements-develop.txt @@ -1,2 +1,4 @@ +black +flake8 pytest wheel diff --git a/setup.py b/setup.py index 65fb269..ede8bd3 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,6 @@ from setuptools import setup, find_packages, Extension from setuptools.command.build_ext import build_ext -from parametric_plasma_source.build_python import build as build_python class CMakeExtention(Extension): def __init__(self, name, sourcedir=""): @@ -18,54 +17,50 @@ def run(self): try: subprocess.check_output(["cmake", "--version"]) except OSError: - raise RuntimeError("CMake must be installed to build the " - "following extentions: " - ", ".join(e.name for e in self.extensions)) + raise RuntimeError( + "CMake must be installed to build the " + "following extentions: " + ", ".join(e.name for e in self.extensions) + ) for ext in self.extensions: self.build_extension(ext) def build_extension(self, ext): - extdir = os.path.abspath( - os.path.dirname(self.get_ext_fullpath(ext.name)) - ) + extdir = os.path.abspath(os.path.dirname(self.get_ext_fullpath(ext.name))) if not extdir.endswith(os.path.sep): extdir += os.path.sep - cmake_args = ["-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=" + extdir, - "-DPYTHON_EXECUTABLE=" + sys.executable] + cmake_args = [ + "-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=" + extdir, + "-DCMAKE_ARCHIVE_OUTPUT_DIRECTORY=" + extdir, + "-DPYTHON_EXECUTABLE=" + sys.executable, + ] cfg = "Debug" if self.debug else "Release" - build_args = ["--target", "plasma_source"] - build_args += ["--config", cfg] + build_args = ["--config", cfg] cmake_args += ["-DCMAKE_BUILD_TYPE=" + cfg] build_args += ["--", "-j2"] env = os.environ.copy() - env["CXXFLAGS"] = "{} -DVERSION_INFO=\\\"{}\\\"".format( - env.get("CXXFLAGS", ""), - self.distribution.get_version() + env["CXXFLAGS"] = '{} -DVERSION_INFO=\\"{}\\"'.format( + env.get("CXXFLAGS", ""), self.distribution.get_version() ) if not os.path.exists(self.build_temp): os.makedirs(self.build_temp) subprocess.check_call( - ["cmake", ext.sourcedir] + cmake_args, - cwd=self.build_temp, - env=env + ["cmake", ext.sourcedir] + cmake_args, cwd=self.build_temp, env=env ) subprocess.check_call( - ["cmake", "--build", "."] + build_args, - cwd=self.build_temp + ["cmake", "--build", "."] + build_args, cwd=self.build_temp ) with open("README.md", "r") as fh: long_description = fh.read() -build_python(source_dir="parametric_plasma_source") - setup( name="parametric_plasma_source", version="0.0.6", @@ -77,6 +72,7 @@ def build_extension(self, ext): url="https://github.com/makeclean/parametric-plasma-source/", packages=find_packages(), ext_modules=[CMakeExtention("parametric_plasma_source/plasma_source")], + package_data={"parametric_plasma_source": ["source_sampling.so"]}, cmdclass=dict(build_ext=CMakeBuild), classifiers=[ "Programming Language :: Python :: 3", diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..4b8a656 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,28 @@ +import pytest + +from parametric_plasma_source import PlasmaSource + +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, +} + + +@pytest.fixture(scope="session") +def plasma_source(): + """Make a plasma source to use as a test fixture.""" + return PlasmaSource(**plasma_params) diff --git a/tests/test_Compile.py b/tests/test_Compile.py deleted file mode 100644 index 405a275..0000000 --- a/tests/test_Compile.py +++ /dev/null @@ -1,26 +0,0 @@ - - -import pytest -import unittest - -import os -from pathlib import Path - -from parametric_plasma_source.plasma import Plasma - -pytest.importorskip("openmc") - -class test_object_properties(unittest.TestCase): - - def test_compile(self): - - os.system('test_plasma.so') - test_plasma = Plasma() - test_plasma.export_plasma_source('test_plasma.so') - - assert Path('test_plasma.so').exists() == True - os.system('rm test_plasma.so') - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_data/baseline_source.txt b/tests/test_data/baseline_source.txt new file mode 100644 index 0000000..e249ef4 --- /dev/null +++ b/tests/test_data/baseline_source.txt @@ -0,0 +1 @@ +major_radius=9.06, minor_radius=2.92258, elongation=1.557, triangularity=0.27, shafranov_shift=0.44789, pedestal_radius=2.33806, 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 \ No newline at end of file diff --git a/tests/test_plasma_source.py b/tests/test_plasma_source.py index 56a5437..88cb045 100644 --- a/tests/test_plasma_source.py +++ b/tests/test_plasma_source.py @@ -1,33 +1,11 @@ """Tests for the methods in plasma_source.""" +import math +import os import pytest +from random import random -from parametric_plasma_source.plasma_source import PlasmaSource - -plasma_params = { - "elongation": 1.557, - "ion_density_origin": 1.09e20, - "ion_density_peaking_factor": 1, - "ion_density_pedistal": 1.09e20, - "ion_density_seperatrix": 3e19, - "ion_temperature_origin": 45.9, - "ion_temperature_peaking_factor": 8.06, - "ion_temperature_pedistal": 6.09, - "ion_temperature_seperatrix": 0.1, - "major_radius": 906.0, - "minor_radius": 292.258, - "pedistal_radius": 0.8 * 292.258, - "plasma_id": 1, - "shafranov_shift": 44.789, - "triangularity": 0.270, - "ion_temperature_beta": 6, -} - - -@pytest.fixture(scope="session") -def plasma_source(): - """Make a plasma source to use as a test fixture.""" - return PlasmaSource(**plasma_params) +from parametric_plasma_source import PlasmaSource class TestPlasmaSource: @@ -37,60 +15,111 @@ def test_ion_density_magnetic_origin(self, plasma_source): """Test the ion density at the magnetic origin.""" ion_density = plasma_source.ion_density(0.0) - assert pytest.approx(ion_density, 1.09e20) + assert ion_density == pytest.approx(1.09e20) def test_ion_density_inside_pedestal(self, plasma_source): """Test the ion density inside the pedestal.""" ion_density = plasma_source.ion_density(0.2) - assert pytest.approx(ion_density, 1.09e20) + assert ion_density == pytest.approx(1.09e20) def test_ion_density_outside_pedestal(self, plasma_source): """Test the ion density outside the pedestal.""" ion_density = plasma_source.ion_density(2.4) - assert pytest.approx(ion_density, 1.00628584e20) + assert ion_density == pytest.approx(1.00629067e20) def test_ion_density_boundary(self, plasma_source): """Test the ion density at the boundary.""" - boundary = plasma_params["minor_radius"] / 100.0 - ion_density = plasma_source.ion_density(boundary) + ion_density = plasma_source.ion_density(plasma_source.minor_radius) - assert pytest.approx( - ion_density, plasma_params["ion_density_seperatrix"] - ) + assert ion_density == pytest.approx(plasma_source.ion_density_separatrix) def test_ion_temperature_magnetic_origin(self, plasma_source): """Test the ion temperature at the magnetic origin.""" ion_temperature = plasma_source.ion_temperature(0.0) - assert pytest.approx( - ion_temperature, plasma_params["ion_temperature_origin"] - ) + assert ion_temperature == pytest.approx(plasma_source.ion_temperature_origin) def test_ion_temperature_inside_pedestal(self, plasma_source): """Test the ion temperature inside the pedestal.""" ion_temperature = plasma_source.ion_temperature(0.2) - assert pytest.approx(ion_temperature, 45.89987429) + assert ion_temperature == pytest.approx(45.89987429) def test_ion_temperature_outside_pedestal(self, plasma_source): """Test the ion temperature outside the pedestal.""" ion_temperature = plasma_source.ion_temperature(2.4) - assert pytest.approx(ion_temperature, 5.45525594) + assert ion_temperature == pytest.approx(5.45529258) def test_ion_temperature_boundary(self, plasma_source): """Test the ion temperature at the boundary.""" - boundary = plasma_params["minor_radius"] / 100.0 - ion_temperature = plasma_source.ion_temperature(boundary) + ion_temperature = plasma_source.ion_temperature(plasma_source.minor_radius) - assert pytest.approx( - ion_temperature, plasma_params["ion_temperature_seperatrix"] + assert ion_temperature == pytest.approx( + plasma_source.ion_temperature_separatrix ) def test_dt_cross_section(self, plasma_source): """Test the dt cross section at a specific temperature.""" dt_cross_section = plasma_source.dt_xs(4.25e7) - assert pytest.approx(dt_cross_section, 0.0) + assert dt_cross_section == pytest.approx(0.0) + + def test_source_to_from_string(self, plasma_source): + """Test the source can be converted to and from a string.""" + the_str = str(plasma_source) + new_source = PlasmaSource.from_string(the_str) + + assert id(plasma_source) != id(new_source) + assert str(new_source) == the_str + + def test_source_string_regression(self, plasma_source): + """Test the source string representation matches the baseline.""" + baseline_file = os.sep.join( + [os.path.dirname(__file__), "test_data", "baseline_source.txt"] + ) + with open(baseline_file, "r") as f: + baseline_str = f.read().strip() + + the_str = str(plasma_source) + assert baseline_str == the_str + + def test_sampling(self, plasma_source): + """Test the sampling function.""" + randoms = [r() for r in [random] * 8] + sample = plasma_source.sample(randoms) + + assert len(sample) == 7 + + # The 4th, 5th and 6th elements together define a unit vector + direction_length = math.sqrt(sample[3] ** 2 + sample[4] ** 2 + sample[5] ** 2) + assert direction_length == pytest.approx(1.0) + + def test_sampling_regression(self, plasma_source): + """Test the sampling function returns matching results.""" + randoms = [ + 0.17213994440390412, + 0.9186868218670968, + 0.5789738834800362, + 0.08876642179434446, + 0.9556278780110383, + 0.8967227763309567, + 0.5187262083328932, + 0.09064281320718603, + ] + + expected = ( + 4.9015857574526, + 7.8576247486517, + -0.19321843360054, + -0.5702309715680232, + -0.06740484811110535, + 0.8187143735856279, + 14.202333312096737, + ) + + sample = plasma_source.sample(randoms) + + assert sample == pytest.approx(expected)