diff --git a/examples/example_get_particle_data.py b/examples/example_get_particle_data.py index aafc09d..0f6a88e 100644 --- a/examples/example_get_particle_data.py +++ b/examples/example_get_particle_data.py @@ -13,14 +13,7 @@ # sets the energy distribution to 100% 14MeV neutrons my_source.energy = openmc.stats.Discrete([14e6], [1]) -# makes an initial_source.h5 file with details of the particles -initial_source_filename = osp.create_initial_particles( - source=my_source, - number_of_particles=10, - openmc_exec="/home/jshim/miniconda3/envs/openmc_0_11_0/bin/openmc", -) - # gets the particle corrdiantes, energy and direction -data = osp.get_particle_data(initial_source_filename) +data = osp.sample_initial_particles(my_source) print(data) diff --git a/examples/example_plot_plasma_source_position.py b/examples/example_plot_plasma_source_position.py index 2a118d8..c0f7947 100644 --- a/examples/example_plot_plasma_source_position.py +++ b/examples/example_plot_plasma_source_position.py @@ -22,10 +22,6 @@ # plots the particle energy distribution -plot = osp.plot_source_position( - source=my_source, - number_of_particles=10, - openmc_exec="/home/jshim/miniconda3/envs/openmc_0_11_0/bin/openmc", -) +plot = osp.plot_source_position(source=my_source, n_samples=100) plot.show() diff --git a/examples/example_plot_source_direction.py b/examples/example_plot_source_direction.py index d568f6f..596a91c 100644 --- a/examples/example_plot_source_direction.py +++ b/examples/example_plot_source_direction.py @@ -8,10 +8,6 @@ my_source.angle = openmc.stats.Isotropic() # plots the particle energy distribution -plot = osp.plot_source_direction( - source=my_source, - number_of_particles=100, - openmc_exec="/home/jshim/miniconda3/envs/openmc_0_11_0/bin/openmc", -) +plot = osp.plot_source_direction(source=my_source, n_samples=200) plot.show() diff --git a/examples/example_plot_source_energy.py b/examples/example_plot_source_energy.py index 74f0c9c..ea8aa4e 100644 --- a/examples/example_plot_source_energy.py +++ b/examples/example_plot_source_energy.py @@ -11,8 +11,7 @@ # plots the particle energy distribution plot = osp.plot_source_energy( source=my_source, - number_of_particles=10000, - energy_bins=np.linspace(0, 20e6, 100), + n_samples=10000, ) plot.show() diff --git a/examples/example_plot_source_position.py b/examples/example_plot_source_position.py index 381c763..0ffda4b 100644 --- a/examples/example_plot_source_position.py +++ b/examples/example_plot_source_position.py @@ -19,8 +19,6 @@ ) # plots the particle energy distribution -plot = osp.plot_source_position( - source=my_source, openmc_exec="/home/jshim/miniconda3/envs/openmc_0_11_0/bin/openmc" -) +plot = osp.plot_source_position(source=my_source) plot.show() diff --git a/examples/example_plot_two_source_energies.py b/examples/example_plot_two_source_energies.py index ae4dd75..c97c443 100644 --- a/examples/example_plot_two_source_energies.py +++ b/examples/example_plot_two_source_energies.py @@ -14,10 +14,6 @@ my_dd_source.energy = openmc.stats.Muir(e0=2080000.0, m_rat=2.0, kt=20000.0) # plots the particle energy distribution -plot = osp.plot_source_energy( - source=[my_dt_source, my_dd_source], - number_of_particles=10000, - energy_bins=np.linspace(0, 20e6, 100), -) +plot = osp.plot_source_energy(source=[my_dt_source, my_dd_source], n_samples=10000) plot.show() diff --git a/openmc_source_plotter/__init__.py b/openmc_source_plotter/__init__.py index 91d707e..4117953 100644 --- a/openmc_source_plotter/__init__.py +++ b/openmc_source_plotter/__init__.py @@ -13,8 +13,7 @@ __all__ = ["__version__"] -from .utils import get_particle_data -from .utils import create_initial_particles from .core import plot_source_direction from .core import plot_source_energy from .core import plot_source_position +from .core import sample_initial_particles diff --git a/openmc_source_plotter/_version.py b/openmc_source_plotter/_version.py index 5787b28..79de5f2 100644 --- a/openmc_source_plotter/_version.py +++ b/openmc_source_plotter/_version.py @@ -1,5 +1,5 @@ # coding: utf-8 # file generated by setuptools_scm # don't change, don't track in version control -__version__ = version = "0.2.1.dev1+g15fa018" -__version_tuple__ = version_tuple = (0, 2, 1, "dev1", "g15fa018") +__version__ = version = "0.2.1.dev0+g1eea6a4.d20220717" +__version_tuple__ = version_tuple = (0, 2, 1, "dev0", "g1eea6a4.d20220717") diff --git a/openmc_source_plotter/core.py b/openmc_source_plotter/core.py index 4267dd6..fce31d7 100644 --- a/openmc_source_plotter/core.py +++ b/openmc_source_plotter/core.py @@ -7,27 +7,48 @@ import numpy as np import openmc +import openmc.lib import plotly.graph_objects as go -from .utils import create_initial_particles, get_particle_data, save_plot + +def sample_initial_particles(source: openmc.source, n_samples=1000, prn_seed=None): + + settings = openmc.Settings() + settings.particles = 1 + settings.batches = 1 + settings.export_to_xml() + + materials = openmc.Materials() + materials.export_to_xml() + + sph = openmc.Sphere(r=1, boundary_type="vacuum") + cell = openmc.Cell(region=-sph) + geometry = openmc.Geometry([cell]) + + geometry.export_to_xml() + # model.geometry = openmc.Geometry([cell]) + + # model = openmc.Model() + + openmc.lib.init() + particles = openmc.lib.sample_external_source( + n_samples=n_samples, prn_seed=prn_seed + ) + openmc.lib.finalize() + return particles def plot_source_energy( source: Union[openmc.Source, List[openmc.Source]], - number_of_particles: int = 2000, - energy_bins: Union[np.array, str] = "auto", - filename: str = None, + n_samples: int = 2000, + prn_seed: int = 1, ): - """makes a plot of the energy distribution OpenMC source(s) + """makes a plot of the initial creation postions of an OpenMC source(s) Args: source: The openmc.Source object or list of openmc.Source objects to plot. - number_of_particles: The number of source samples to obtain, more will - take longer but produce a smoother plot. - energy_bins: A numpy array of energy bins to use as energy bins. 'Auto' - is also accepted and this picks the bins for you using numpy - filename: the filename to save the plot as should end with the correct - extension supported by matplotlib (e.g .png) or plotly (e.g .html) + n_smaples: The number of source samples to obtain. + prn_seed: The pseudorandom number seed """ figure = go.Figure() @@ -37,16 +58,16 @@ def plot_source_energy( for single_source in source: - e_values = single_source.energy.sample(n_samples=number_of_particles) + e_values = single_source.energy.sample(n_samples=n_samples) # Calculate pdf for source energies - probability, bin_edges = np.histogram(e_values, bins=energy_bins, density=True) + probability, bin_edges = np.histogram(e_values, bins="auto", density=True) # Plot source energy histogram figure.add_trace( go.Scatter( - x=energy_bins[:-1], - y=probability * np.diff(energy_bins), + x=bin_edges[:-1], + y=probability * np.diff(bin_edges), line={"shape": "hv"}, hoverinfo="text", ) @@ -59,24 +80,20 @@ def plot_source_energy( showlegend=True, ) - plotting_package = "plotly" # not matplotlib option for now - save_plot(plotting_package=plotting_package, filename=filename, figure=figure) - return figure def plot_source_position( source: Union[openmc.Source, List[openmc.Source]], - number_of_particles: int = 2000, - openmc_exec="openmc", - filename: str = None, + n_samples: int = 2000, + prn_seed: int = 1, ): """makes a plot of the initial creation postions of an OpenMC source(s) Args: source: The openmc.Source object or list of openmc.Source objects to plot. - number_of_particles: The number of source samples to obtain. - openmc_exec: The path of the openmc executable to use + n_smaples: The number of source samples to obtain. + prn_seed: The pseudorandom number seed """ figure = go.Figure() @@ -85,53 +102,42 @@ def plot_source_position( source = [source] for single_source in source: - tmp_filename = tempfile.mkstemp(suffix=".h5", prefix=f"openmc_source_")[1] - create_initial_particles( - source=single_source, - number_of_particles=number_of_particles, - openmc_exec=openmc_exec, - output_source_filename=tmp_filename, - ) - data = get_particle_data(tmp_filename) + data = sample_initial_particles(single_source, n_samples, prn_seed) - text = ["Energy = " + str(i) + " eV" for i in data["e_values"]] + text = ["Energy = " + str(particle.E) + " eV" for particle in data] figure.add_trace( go.Scatter3d( - x=data["x_values"], - y=data["y_values"], - z=data["z_values"], + x=[particle.r[0] for particle in data], + y=[particle.r[1] for particle in data], + z=[particle.r[2] for particle in data], hovertext=text, text=text, mode="markers", marker={ "size": 2, - "color": data["e_values"], + # "color": data.E, }, ) ) figure.update_layout(title="Particle production coordinates - coloured by energy") - plotting_package = "plotly" # not matplotlib option for now - save_plot(plotting_package=plotting_package, filename=filename, figure=figure) - return figure def plot_source_direction( source: Union[openmc.Source, List[openmc.Source]], - number_of_particles: int = 2000, - openmc_exec="openmc", - filename: str = None, + n_samples: int = 2000, + prn_seed: int = 1, ): - """makes a plot of the initial creation directions of the particle source + """makes a plot of the initial creation postions of an OpenMC source(s) Args: source: The openmc.Source object or list of openmc.Source objects to plot. - number_of_particles: The number of source samples to obtain. - openmc_exec: The path of the openmc executable to use + n_smaples: The number of source samples to obtain. + prn_seed: The pseudorandom number seed """ figure = go.Figure() @@ -139,25 +145,18 @@ def plot_source_direction( source = [source] for single_source in source: - tmp_filename = tempfile.mkstemp(suffix=".h5", prefix=f"openmc_source_")[1] - create_initial_particles( - source=single_source, - number_of_particles=number_of_particles, - openmc_exec=openmc_exec, - output_source_filename=tmp_filename, - ) - data = get_particle_data(tmp_filename) + data = sample_initial_particles(single_source, n_samples, prn_seed) figure.add_trace( { "type": "cone", "cauto": False, - "x": data["x_values"], - "y": data["y_values"], - "z": data["z_values"], - "u": data["x_dir"], - "v": data["y_dir"], - "w": data["z_dir"], + "x": [particle.r[0] for particle in data], + "y": [particle.r[1] for particle in data], + "z": [particle.r[2] for particle in data], + "u": [particle.u[0] for particle in data], + "v": [particle.u[1] for particle in data], + "w": [particle.u[2] for particle in data], "cmin": 0, "cmax": 1, "anchor": "tail", @@ -171,7 +170,4 @@ def plot_source_direction( figure.update_layout(title="Particle initial directions") - plotting_package = "plotly" # not matplotlib option for now - save_plot(plotting_package=plotting_package, filename=filename, figure=figure) - return figure diff --git a/openmc_source_plotter/utils.py b/openmc_source_plotter/utils.py index e7e53de..e69de29 100644 --- a/openmc_source_plotter/utils.py +++ b/openmc_source_plotter/utils.py @@ -1,114 +0,0 @@ -#!/usr/bin/env python3 - -"""Provides utilities for creating h5 files containing initial source -information""" - -import h5py -import openmc -import shutil - - -def get_particle_data(input_filename: str): - - f = h5py.File(input_filename, "r") - dset = f["source_bank"] - - x_values = [] - y_values = [] - z_values = [] - x_dir = [] - y_dir = [] - z_dir = [] - e_values = [] - - for particle in dset: - x_values.append(particle[0][0]) - y_values.append(particle[0][1]) - z_values.append(particle[0][2]) - x_dir.append(particle[1][0]) - y_dir.append(particle[1][1]) - z_dir.append(particle[1][2]) - e_values.append(particle[2]) - - return { - "x_values": x_values, - "y_values": y_values, - "z_values": z_values, - "x_dir": x_dir, - "y_dir": y_dir, - "z_dir": z_dir, - "e_values": e_values, - } - - -def create_initial_particles( - source: openmc.source, - number_of_particles: int = 2000, - openmc_exec="openmc", - output_source_filename="initial_source.h5", -) -> str: - """Accepts an openmc source and creates an initial_source.h5 file that can - be used to find propties of the source particles such as initial x,y,z - coordinates, direction and energy. - - Note, I've found it easiest to install the latest version of openmc as this - has the write_initial_source for fixed source problems and then in a - new empty conda environment I have installed openmc version 0.11.0 with the - 'conda install -c conda-forge openmc==0.11.0' command. I then go back to - the conda environment with the latest openmc version and run the python - script while setting the path to openmc_0.11.0 as the openmc_exec argument. - In my case it is '/home/jshimwell/miniconda3/envs/openmc_0_11_0/bin/openmc' - - Args: - source: the openmc source to use - number_of_particles: the number of particles to sample - openmc_exec: the path of openmc executable or executable name if it - appears in your system $PATH. Defaults to 'openmc' which will use - the default openmc in your system $PATH environmental variable. - output_source_filename: the filename of the initial source h5 file - produced. Note openmc will write 'initial_source.h5' but this will - be moved to the output_source_filename specified - Returns: - The filename of the initial source file created (initial_source.h5) - """ - - # no real materials are needed for finding the source - mats = openmc.Materials([]) - - # just a minimal geometry - outer_surface = openmc.Sphere(r=100000, boundary_type="vacuum") - cell = openmc.Cell(region=-outer_surface) - universe = openmc.Universe(cells=[cell]) - geom = openmc.Geometry(universe) - - # Instantiate a Settings object - settings = openmc.Settings() - settings.run_mode = "fixed source" - settings.particles = number_of_particles - settings.batches = 1 - settings.inactive = 0 - settings.write_initial_source = True - settings.source = source - - model = openmc.model.Model(geom, mats, settings) - - model.export_to_xml() - - openmc.run(openmc_exec=openmc_exec) - - if output_source_filename != "initial_source.h5": - shutil.move("initial_source.h5", output_source_filename) - - return output_source_filename - - -def save_plot(plotting_package: str, filename: str, figure): - """Saves the matplotlib or plotly graph object as a file.""" - if filename: - if plotting_package == "matplotlib": - figure.savefig(filename, bbox_inches="tight", dpi=400) - elif plotting_package == "plotly": - if Path(filename).suffix == ".html": - figure.write_html(filename) - else: - figure.write_image(filename)