From 672a9f7013fafcd260b5ea95be5497187839fcb9 Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 30 Nov 2022 20:44:23 +0000 Subject: [PATCH 1/8] modern setup --- pyproject.toml | 55 +++++++++++++++++-- setup.cfg | 48 ---------------- setup.py | 4 -- .../regular_mesh_plotter}/__init__.py | 0 .../regular_mesh_plotter}/core.py | 0 .../regular_mesh_plotter}/utils.py | 0 6 files changed, 49 insertions(+), 58 deletions(-) delete mode 100644 setup.cfg delete mode 100644 setup.py rename {regular_mesh_plotter => src/regular_mesh_plotter}/__init__.py (100%) rename {regular_mesh_plotter => src/regular_mesh_plotter}/core.py (100%) rename {regular_mesh_plotter => src/regular_mesh_plotter}/utils.py (100%) diff --git a/pyproject.toml b/pyproject.toml index 1c00c55..24bdf41 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,10 +1,53 @@ [build-system] -requires = [ - "setuptools >= 45", - "wheel", - "setuptools_scm[toml] >= 6.2", -] +requires = ["setuptools >= 65.4.0", "setuptools_scm[toml]>=7.0.5"] build-backend = "setuptools.build_meta" +[project] +name = "regular_mesh_plotter" +authors = [ + { name="Jonathan Shimwell", email="mail@jshimwell.com" }, +] +license = {file = "LICENSE.txt"} +description = "A Python package for creating publication quality plots of regular mesh tallies with the underlying geometry" +readme = "README.md" +requires-python = ">=3.8" +keywords = ["regular", "mesh", "openmc", "tally", "plot", "slice"] +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", +] +dependencies = [ + "numpy>=1.21.1", + "matplotlib>=3.4.2", + "trimesh", + "shapely", + "scipy", + "pandas", + "dagmc_geometry_slice_plotter", + "openmc_tally_unit_converter", + "setuptools_scm", +] +dynamic = ["version"] + + [tool.setuptools_scm] -write_to = "regular_mesh_plotter/_version.py" +write_to = "src/regular_mesh_plotter/_version.py" + + +[project.optional-dependencies] +tests = [ + "pytest", + "dagmc_h5m_file_inspector", + "pytest-cov>=2.12.1", + "openmc-dagmc-wrapper", + "openmc_plasma_source", + "regular_mesh_plotter", +] + +[project.urls] +"Homepage" = "https://github.com/fusion-energy/regular_mesh_plotter" +"Bug Tracker" = "https://github.com/fusion-energy/regular_mesh_plotter/issues" + +[tool.setuptools] +package-dir = {"" = "src"} diff --git a/setup.cfg b/setup.cfg deleted file mode 100644 index 702be5a..0000000 --- a/setup.cfg +++ /dev/null @@ -1,48 +0,0 @@ -[metadata] -name = regular_mesh_plotter -author = The Fusion Energy Development Team -author_email = mail@jshimwell.com -description = A Python package for creating publication quality plots of regular mesh tallies with the underlying geometry -long_description = file: README.md -long_description_content_type = text/markdown -url = https://github.com/fusion-energy/regular_mesh_plotter -license = MIT -license_file = LICENSE.txt -classifiers = - Natural Language :: English - Topic :: Scientific/Engineering - Programming Language :: Python :: 3 - Programming Language :: Python :: 3.7 - Programming Language :: Python :: 3.8 - Programming Language :: Python :: 3.9 - Programming Language :: Python :: 3.10 - License :: OSI Approved :: MIT License - Operating System :: OS Independent -project_urls = - Source = https://github.com/fusion-energy/regular_mesh_plotter - Tracker = https://github.com/fusion-energy/regular_mesh_plotter/issues - -[options] -packages = find: -python_requires= >=3.6 -install_requires= - numpy >= 1.21.1 - matplotlib >= 3.4.2 - trimesh - shapely - scipy - pandas - dagmc_geometry_slice_plotter - openmc_tally_unit_converter - setuptools_scm - -[options.extras_require] -tests = - pytest >= 5.4.3 - pytest-cov>=2.12.1 - openmc-dagmc-wrapper - openmc_plasma_source - stl_to_h5m - -[flake8] -per-file-ignores = __init__.py:F401 diff --git a/setup.py b/setup.py deleted file mode 100644 index 1abbd06..0000000 --- a/setup.py +++ /dev/null @@ -1,4 +0,0 @@ -import setuptools - -if __name__ == "__main__": - setuptools.setup() diff --git a/regular_mesh_plotter/__init__.py b/src/regular_mesh_plotter/__init__.py similarity index 100% rename from regular_mesh_plotter/__init__.py rename to src/regular_mesh_plotter/__init__.py diff --git a/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py similarity index 100% rename from regular_mesh_plotter/core.py rename to src/regular_mesh_plotter/core.py diff --git a/regular_mesh_plotter/utils.py b/src/regular_mesh_plotter/utils.py similarity index 100% rename from regular_mesh_plotter/utils.py rename to src/regular_mesh_plotter/utils.py From 9164d533ef9db6ad7982ee06c4b7b44b4fe643b4 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 25 Jan 2023 18:46:21 +0000 Subject: [PATCH 2/8] using with openmc_geometry_plot to make hybrid plots --- .../plot_regular_mesh_tally_with_geometry.py | 103 +++-- src/regular_mesh_plotter/__init__.py | 14 +- src/regular_mesh_plotter/core.py | 404 +++++++----------- 3 files changed, 240 insertions(+), 281 deletions(-) diff --git a/examples/plot_regular_mesh_tally_with_geometry.py b/examples/plot_regular_mesh_tally_with_geometry.py index cd63873..bf7db81 100644 --- a/examples/plot_regular_mesh_tally_with_geometry.py +++ b/examples/plot_regular_mesh_tally_with_geometry.py @@ -1,29 +1,84 @@ -import regular_mesh_plotter as rmp +import regular_mesh_plotter import openmc -import matplotlib.pyplot as plt +import openmc_geometry_plot # extends openmc.Geometry class with plotting functions + + +# MATERIALS + +mat_1 = openmc.Material() +mat_1.add_element('Li', 1) +mat_1.set_density('g/cm3', 3.2720171e-2) # around 11 g/cm3 + +my_materials = openmc.Materials([mat_1]) + + +# GEOMETRY + +# surfaces +inner_surface = openmc.Sphere(r=500) +outer_surface = openmc.Sphere(r=1000, boundary_type='vacuum') + +# cells +inner_region = -inner_surface +inner_cell = openmc.Cell(region=inner_region) + +outer_region = -outer_surface & +inner_surface +outer_cell = openmc.Cell(region=outer_region) +outer_cell.fill = mat_1 + +my_geometry = openmc.Geometry([inner_cell, outer_cell]) + + +# SIMULATION SETTINGS -# loads in the statepoint file containing tallies -statepoint = openmc.StatePoint(filepath="statepoint.2.h5") - -# gets one tally from the available tallies -my_tally = statepoint.get_tally(name="neutron_effective_dose_on_2D_mesh_xy") - -# creates a plot of the mesh -rmp.plot_regular_mesh_tally_with_geometry( - tally=my_tally, - dagmc_file_or_trimesh_object="dagmc.h5m", - std_dev_or_tally_value="tally_value", - filename="plot_regular_mesh_tally_with_geometry.png", - label="", - x_label="X [cm]", - y_label="Y [cm]", - # plane_origin=[0, 0, 0], - plane_normal=[0, 0, 1], - rotate_mesh=0, - rotate_geometry=0, - required_units=None, - source_strength=None, +# Instantiate a Settings object +my_settings = openmc.Settings() +my_settings.batches = 10 +my_settings.inactive = 0 +my_settings.particles = 500 +my_settings.run_mode = 'fixed source' + +# Create a DT point source +source = openmc.Source() +source.space = openmc.stats.Point((0, 0, 0)) +source.angle = openmc.stats.Isotropic() +source.energy = openmc.stats.Discrete([14e6], [1]) +my_settings.source = source + + +my_tallies = openmc.Tallies() + +mesh = openmc.RegularMesh().from_domain( + my_geometry, # the corners of the mesh are being set automatically to surround the geometry + dimension=[10, 10, 10] +) +mesh_filter = openmc.MeshFilter(mesh) + +mesh_tally_1 = openmc.Tally(name='mesh_tally') +mesh_tally_1.filters = [mesh_filter] +mesh_tally_1.scores = ['flux'] +my_tallies.append(mesh_tally_1) + +model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies) +sp_filename = model.run() + +statepoint = openmc.StatePoint(sp_filename) + +# extracts the mesh tally by name +my_mesh_tally = statepoint.get_tally(name='mesh_tally') + +# converts the tally result into a VTK file +plot = mesh.plot_slice( + dataset=my_mesh_tally.mean, + slice_index=5, + view_direction='x' ) +plot.figure.savefig('test.png') + +# this section adds a contour line according the the material ids +material_ids = my_geometry.get_slice_of_material_ids(view_direction='x') +plot2 = my_geometry.get_outline_contour(outline_data=material_ids, view_direction='x') +import matplotlib.pyplot as plt -plt.show() +plt.savefig('test2.png') diff --git a/src/regular_mesh_plotter/__init__.py b/src/regular_mesh_plotter/__init__.py index 35d2821..12e2f90 100644 --- a/src/regular_mesh_plotter/__init__.py +++ b/src/regular_mesh_plotter/__init__.py @@ -11,16 +11,4 @@ __all__ = ["__version__"] -from .core import plot_regular_mesh_values -from .core import plot_regular_mesh_values_with_geometry - -from .core import plot_regular_mesh_tally -from .core import plot_regular_mesh_tally_with_geometry - -from .core import plot_regular_mesh_dose_tally -from .core import plot_regular_mesh_dose_tally_with_geometry - -from .utils import reshape_values_to_mesh_shape -from .utils import get_tally_extent -from .utils import get_values_from_tally -from .utils import get_std_dev_or_value_from_tally +from .core import * \ No newline at end of file diff --git a/src/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py index efda23c..0188b09 100644 --- a/src/regular_mesh_plotter/core.py +++ b/src/regular_mesh_plotter/core.py @@ -1,254 +1,170 @@ -from typing import List, Optional - -import dagmc_geometry_slice_plotter as dgsp -import matplotlib.pyplot as plt import numpy as np -import openmc_tally_unit_converter as otuc -from matplotlib import transforms - -from .utils import ( - get_std_dev_or_value_from_tally, - get_tally_extent, - get_values_from_tally, - reshape_values_to_mesh_shape, -) - - -def plot_regular_mesh_values( - values: np.ndarray, - filename: Optional[str] = None, - label="", - title=None, - extent=None, - x_label="X [cm]", - y_label="Y [cm]", - rotate_plot: float = 0, - **kwargs -): - - if rotate_plot != 0: - x_center = sum(extent[:2]) / 2 - y_center = sum(extent[2:]) / 2 - base = plt.gca().transData - rot = transforms.Affine2D().rotate_deg_around(x_center, y_center, rotate_plot) - - image_map = plt.imshow( - values, extent=extent, transform=rot + base, origin="lower", **kwargs - ) - else: - image_map = plt.imshow(values, extent=extent, origin="lower", **kwargs) - - plt.xlabel(x_label) - plt.ylabel(y_label) - if title: - plt.title(title) - - # image_map = fig.imshow(values, norm=scale, vmin=vmin) - plt.colorbar(image_map, label=label) - if filename: # TODO should we not let the users do that? - plt.savefig(filename, dpi=300) - - -def plot_regular_mesh_values_with_geometry( - values: np.ndarray, - dagmc_file_or_trimesh_object, - filename: Optional[str] = None, - label="", - title=None, - extent=None, - x_label="X [cm]", - y_label="Y [cm]", - plane_origin: List[float] = None, - plane_normal: List[float] = [0, 0, 1], - rotate_mesh: float = 0, - rotate_geometry: float = 0, - **kwargs -): - - slice = dgsp.plot_slice_of_dagmc_geometry( - dagmc_file_or_trimesh_object=dagmc_file_or_trimesh_object, - plane_origin=plane_origin, - plane_normal=plane_normal, - rotate_plot=rotate_geometry, - ) - - plot_regular_mesh_values( - values=values, - filename=filename, - label=label, - title=title, - extent=extent, - x_label=x_label, - y_label=y_label, - rotate_plot=rotate_mesh, +import typing +import openmc + +class RegularMesh(openmc.RegularMesh): + + def slice_of_data( + self, + dataset: np.ndarray, + slice_index: int = 0, + view_direction: str = 'z', + volume_normalization: bool = True + ): + """Obtains the dataset values on an axis aligned 2D slice through the + mesh. Useful for producing plots of slice data + + Parameters + ---------- + datasets : numpy.ndarray + 1-D array of data values. Will be reshaped to fill the mesh and + should therefore have the same number of entries to fill the mesh + slice_index : int + The index of the mesh slice to extract. + view_direction : str + The axis to view the slice from. Supports negative and positive axis + values. Acceptable values are 'x', 'y', 'z', '-x', '-y', '-z'. + volume_normalization : bool, optional + Whether or not to normalize the data by the volume of the mesh + elements. + + Returns + ------- + np.array() + the 2D array of dataset values + """ + + if volume_normalization: + dataset = dataset.flatten() / self.volumes.flatten() + + reshaped_ds = dataset.reshape(self.dimension, order="F") + + + if view_direction == "x": + # vertical axis is z, horizontal axis is -y + transposed_ds = reshaped_ds.transpose(0, 1, 2)[slice_index] + rotated_ds = np.rot90(transposed_ds, 1) + aligned_ds = np.fliplr(rotated_ds) + elif view_direction == "-x": + # vertical axis is z, horizontal axis is y + transposed_ds = reshaped_ds.transpose(0, 1, 2)[slice_index] + aligned_ds = np.rot90(transposed_ds, 1) + elif view_direction == "y": + # vertical axis is z, horizontal axis is x + transposed_ds = reshaped_ds.transpose(1, 2, 0)[slice_index] + aligned_ds = np.flipud(transposed_ds) + elif view_direction == "-y": + # vertical axis is z, horizontal axis is -x + transposed_ds = reshaped_ds.transpose(1, 2, 0)[slice_index] + aligned_ds = np.flipud(transposed_ds) + aligned_ds = np.fliplr(aligned_ds) + elif view_direction == "z": + # vertical axis is y, horizontal axis is -x + transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] + aligned_ds = np.rot90(transposed_ds, 1) + aligned_ds = np.fliplr(aligned_ds) + elif view_direction == "-z": + # vertical axis is y, horizontal axis is x + transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] + aligned_ds = np.rot90(transposed_ds, 1) + else: + msg = 'view_direction is not one of the acceptable options {supported_view_dirs}' + raise ValueError(msg) + + return aligned_ds + + + def plot_slice( + self, + dataset: np.ndarray, + slice_index: typing.Optional[int] = None, + view_direction: str = 'z', + axes: typing.Optional['matplotlib.Axes'] = None, + volume_normalization: bool = True, **kwargs - ) - - -def plot_regular_mesh_tally_with_geometry( - tally, - dagmc_file_or_trimesh_object, - std_dev_or_tally_value="tally_value", - filename: Optional[str] = None, - label="", - title=None, - x_label="X [cm]", - y_label="Y [cm]", - plane_origin: List[float] = None, - plane_normal: List[float] = [0, 0, 1], - rotate_mesh: float = 0, - rotate_geometry: float = 0, - required_units=None, - source_strength: float = None, - **kwargs -): - - if required_units is not None: - values = otuc.process_tally( - tally=tally, required_units=required_units, source_strength=source_strength - ) - else: - values = get_values_from_tally(tally) - - value = get_std_dev_or_value_from_tally(tally, values, std_dev_or_tally_value) - - extent = get_tally_extent(tally) - - dgsp.plot_slice_of_dagmc_geometry( - dagmc_file_or_trimesh_object=dagmc_file_or_trimesh_object, - plane_origin=plane_origin, - plane_normal=plane_normal, - rotate_plot=rotate_geometry, - ) - - plot_regular_mesh_values( - values=value, - filename=filename, - label=label, - title=title, - extent=extent, - x_label=x_label, - y_label=y_label, - rotate_plot=rotate_mesh, + ): + """Create a slice plot of the dataset on the RegularMesh. + + Parameters + ---------- + datasets : numpy.ndarray + 1-D array of data values. Will be reshaped to fill the mesh and + should therefore have the same number of entries to fill the mesh + slice_index : int + The index of the mesh slice to extract. If not set then the index + that slices through the middle of the mesh will be used + view_direction : str + The axis to view the slice from. Supports negative and positive axis + values. Acceptable values are 'x', 'y', 'z', '-x', '-y', '-z'. + axes : matplotlib.Axes, optional + Axes to draw to + volume_normalization : bool, optional + Whether or not to normalize the data by the volume of the mesh + elements. **kwargs - ) - - -def plot_regular_mesh_tally( - tally, - filename: Optional[str] = None, - label="", - title=None, - x_label="X [cm]", - y_label="Y [cm]", - rotate_plot: float = 0, - required_units: str = None, - source_strength: float = None, - std_dev_or_tally_value="tally_value", - **kwargs -): - - if required_units is not None: - values = otuc.process_tally( - tally=tally, required_units=required_units, source_strength=source_strength + Keyword arguments passed to :func:`matplotlib.pyplot.imshow` + + Returns + ------- + matplotlib.image.AxesImage + Resulting image + """ + + import matplotlib.pyplot as plt + + # gets the axis labels and bounding box index + if 'x' in view_direction: + x_label = 'Y [cm]' + y_label = 'Z [cm]' + bb_index = 0 + + if 'y' in view_direction: + x_label = 'X [cm]' + y_label = 'Z [cm]' + bb_index = 1 + + if 'z' in view_direction: + x_label = 'X [cm]' + y_label = 'Y [cm]' + bb_index = 2 + + # selecting mid index on the mesh for the slice + if slice_index is None: + slice_index = int(self.dimension[bb_index]/2) + + # slice_of_data also checks the view_direction is acceptable + image_slice = self.slice_of_data( + dataset=dataset, + slice_index=slice_index, + view_direction=view_direction, + volume_normalization=volume_normalization, ) - else: - values = get_values_from_tally(tally) - - value = get_std_dev_or_value_from_tally(tally, values, std_dev_or_tally_value) - - value = reshape_values_to_mesh_shape(tally, value) - extent = get_tally_extent(tally) - - plot_regular_mesh_values( - values=value, - filename=filename, - label=label, - title=title, - extent=extent, - x_label=x_label, - y_label=y_label, - rotate_plot=rotate_plot, - **kwargs - ) - - -def plot_regular_mesh_dose_tally( - tally, - filename: Optional[str] = None, - label="", - title=None, - x_label="X [cm]", - y_label="Y [cm]", - rotate_plot: float = 0, - required_units="picosievert / source_particle", - source_strength: float = None, - std_dev_or_tally_value: str = "tally_value", - **kwargs -): - - if required_units is not None: - values = otuc.process_dose_tally( - tally=tally, required_units=required_units, source_strength=source_strength + # gets the extent of the plot + x_min = self.lower_left[bb_index] + x_max = self.upper_right[bb_index] + y_min = self.lower_left[bb_index] + y_max = self.upper_right[bb_index] + + if axes is None: + fig, axes = plt.subplots() + axes.set_xlabel(x_label) + axes.set_ylabel(y_label) + axes.set_title(f'View direction {view_direction}') + + return axes.imshow( + X=image_slice, + extent=(x_min, x_max, y_min, y_max), + aspect='auto', + **kwargs ) - else: - values = get_values_from_tally(tally) - value = get_std_dev_or_value_from_tally(tally, values, std_dev_or_tally_value) - value = reshape_values_to_mesh_shape(tally, value) +# def get_cell_ids_for_regularmesh_slice(mesh, geometry): + # loop through the centroids of the mesh + # find the material at each point + # build up a pixel map of material ids - extent = get_tally_extent(tally) - plot_regular_mesh_values( - values=value, - filename=filename, - label=label, - title=title, - extent=extent, - x_label=x_label, - y_label=y_label, - rotate_plot=rotate_plot, - **kwargs - ) - - -def plot_regular_mesh_dose_tally_with_geometry( - tally, - dagmc_file_or_trimesh_object, - filename: Optional[str] = None, - label="", - title=None, - x_label="X [cm]", - y_label="Y [cm]", - plane_origin: List[float] = None, - plane_normal: List[float] = [0, 0, 1], - rotate_mesh: float = 0, - rotate_geometry: float = 0, - required_units="picosievert / source_particle", - source_strength: float = None, - std_dev_or_tally_value: str = "tally_value", - **kwargs -): - - dgsp.plot_slice_of_dagmc_geometry( - dagmc_file_or_trimesh_object=dagmc_file_or_trimesh_object, - plane_origin=plane_origin, - plane_normal=plane_normal, - rotate_plot=rotate_geometry, - ) - - plot_regular_mesh_dose_tally( - tally=tally, - filename=filename, - label=label, - title=title, - x_label=x_label, - y_label=y_label, - rotate_plot=rotate_mesh, - required_units=required_units, - source_strength=source_strength, - std_dev_or_tally_value=std_dev_or_tally_value, - **kwargs - ) +openmc.RegularMesh = RegularMesh + From 855cdfd424629d4cf4ec7dc01533d08b79763c39 Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 25 Jan 2023 18:47:49 +0000 Subject: [PATCH 3/8] [skip ci] Apply formatting changes --- .../plot_regular_mesh_tally_with_geometry.py | 32 ++++++------ src/regular_mesh_plotter/__init__.py | 2 +- src/regular_mesh_plotter/core.py | 50 ++++++++----------- 3 files changed, 37 insertions(+), 47 deletions(-) diff --git a/examples/plot_regular_mesh_tally_with_geometry.py b/examples/plot_regular_mesh_tally_with_geometry.py index bf7db81..dab528c 100644 --- a/examples/plot_regular_mesh_tally_with_geometry.py +++ b/examples/plot_regular_mesh_tally_with_geometry.py @@ -7,8 +7,8 @@ # MATERIALS mat_1 = openmc.Material() -mat_1.add_element('Li', 1) -mat_1.set_density('g/cm3', 3.2720171e-2) # around 11 g/cm3 +mat_1.add_element("Li", 1) +mat_1.set_density("g/cm3", 3.2720171e-2) # around 11 g/cm3 my_materials = openmc.Materials([mat_1]) @@ -17,7 +17,7 @@ # surfaces inner_surface = openmc.Sphere(r=500) -outer_surface = openmc.Sphere(r=1000, boundary_type='vacuum') +outer_surface = openmc.Sphere(r=1000, boundary_type="vacuum") # cells inner_region = -inner_surface @@ -37,7 +37,7 @@ my_settings.batches = 10 my_settings.inactive = 0 my_settings.particles = 500 -my_settings.run_mode = 'fixed source' +my_settings.run_mode = "fixed source" # Create a DT point source source = openmc.Source() @@ -50,14 +50,14 @@ my_tallies = openmc.Tallies() mesh = openmc.RegularMesh().from_domain( - my_geometry, # the corners of the mesh are being set automatically to surround the geometry - dimension=[10, 10, 10] + my_geometry, # the corners of the mesh are being set automatically to surround the geometry + dimension=[10, 10, 10], ) mesh_filter = openmc.MeshFilter(mesh) -mesh_tally_1 = openmc.Tally(name='mesh_tally') +mesh_tally_1 = openmc.Tally(name="mesh_tally") mesh_tally_1.filters = [mesh_filter] -mesh_tally_1.scores = ['flux'] +mesh_tally_1.scores = ["flux"] my_tallies.append(mesh_tally_1) model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies) @@ -66,19 +66,15 @@ statepoint = openmc.StatePoint(sp_filename) # extracts the mesh tally by name -my_mesh_tally = statepoint.get_tally(name='mesh_tally') +my_mesh_tally = statepoint.get_tally(name="mesh_tally") # converts the tally result into a VTK file -plot = mesh.plot_slice( - dataset=my_mesh_tally.mean, - slice_index=5, - view_direction='x' -) -plot.figure.savefig('test.png') +plot = mesh.plot_slice(dataset=my_mesh_tally.mean, slice_index=5, view_direction="x") +plot.figure.savefig("test.png") # this section adds a contour line according the the material ids -material_ids = my_geometry.get_slice_of_material_ids(view_direction='x') -plot2 = my_geometry.get_outline_contour(outline_data=material_ids, view_direction='x') +material_ids = my_geometry.get_slice_of_material_ids(view_direction="x") +plot2 = my_geometry.get_outline_contour(outline_data=material_ids, view_direction="x") import matplotlib.pyplot as plt -plt.savefig('test2.png') +plt.savefig("test2.png") diff --git a/src/regular_mesh_plotter/__init__.py b/src/regular_mesh_plotter/__init__.py index 12e2f90..801c6ac 100644 --- a/src/regular_mesh_plotter/__init__.py +++ b/src/regular_mesh_plotter/__init__.py @@ -11,4 +11,4 @@ __all__ = ["__version__"] -from .core import * \ No newline at end of file +from .core import * diff --git a/src/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py index 0188b09..59a324e 100644 --- a/src/regular_mesh_plotter/core.py +++ b/src/regular_mesh_plotter/core.py @@ -2,14 +2,14 @@ import typing import openmc -class RegularMesh(openmc.RegularMesh): +class RegularMesh(openmc.RegularMesh): def slice_of_data( self, dataset: np.ndarray, slice_index: int = 0, - view_direction: str = 'z', - volume_normalization: bool = True + view_direction: str = "z", + volume_normalization: bool = True, ): """Obtains the dataset values on an axis aligned 2D slice through the mesh. Useful for producing plots of slice data @@ -39,7 +39,6 @@ def slice_of_data( reshaped_ds = dataset.reshape(self.dimension, order="F") - if view_direction == "x": # vertical axis is z, horizontal axis is -y transposed_ds = reshaped_ds.transpose(0, 1, 2)[slice_index] @@ -68,20 +67,19 @@ def slice_of_data( transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] aligned_ds = np.rot90(transposed_ds, 1) else: - msg = 'view_direction is not one of the acceptable options {supported_view_dirs}' + msg = "view_direction is not one of the acceptable options {supported_view_dirs}" raise ValueError(msg) return aligned_ds - def plot_slice( self, dataset: np.ndarray, slice_index: typing.Optional[int] = None, - view_direction: str = 'z', - axes: typing.Optional['matplotlib.Axes'] = None, + view_direction: str = "z", + axes: typing.Optional["matplotlib.Axes"] = None, volume_normalization: bool = True, - **kwargs + **kwargs, ): """Create a slice plot of the dataset on the RegularMesh. @@ -113,24 +111,24 @@ def plot_slice( import matplotlib.pyplot as plt # gets the axis labels and bounding box index - if 'x' in view_direction: - x_label = 'Y [cm]' - y_label = 'Z [cm]' + if "x" in view_direction: + x_label = "Y [cm]" + y_label = "Z [cm]" bb_index = 0 - if 'y' in view_direction: - x_label = 'X [cm]' - y_label = 'Z [cm]' + if "y" in view_direction: + x_label = "X [cm]" + y_label = "Z [cm]" bb_index = 1 - if 'z' in view_direction: - x_label = 'X [cm]' - y_label = 'Y [cm]' + if "z" in view_direction: + x_label = "X [cm]" + y_label = "Y [cm]" bb_index = 2 # selecting mid index on the mesh for the slice if slice_index is None: - slice_index = int(self.dimension[bb_index]/2) + slice_index = int(self.dimension[bb_index] / 2) # slice_of_data also checks the view_direction is acceptable image_slice = self.slice_of_data( @@ -150,21 +148,17 @@ def plot_slice( fig, axes = plt.subplots() axes.set_xlabel(x_label) axes.set_ylabel(y_label) - axes.set_title(f'View direction {view_direction}') + axes.set_title(f"View direction {view_direction}") return axes.imshow( - X=image_slice, - extent=(x_min, x_max, y_min, y_max), - aspect='auto', - **kwargs + X=image_slice, extent=(x_min, x_max, y_min, y_max), aspect="auto", **kwargs ) # def get_cell_ids_for_regularmesh_slice(mesh, geometry): - # loop through the centroids of the mesh - # find the material at each point - # build up a pixel map of material ids +# loop through the centroids of the mesh +# find the material at each point +# build up a pixel map of material ids openmc.RegularMesh = RegularMesh - From abdc12e1361dab9e2f61f97f29922d60ca699ec1 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 26 Jan 2023 09:34:19 +0000 Subject: [PATCH 4/8] fixed import issue --- src/regular_mesh_plotter/core.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py index 59a324e..9bfd2a4 100644 --- a/src/regular_mesh_plotter/core.py +++ b/src/regular_mesh_plotter/core.py @@ -162,3 +162,4 @@ def plot_slice( openmc.RegularMesh = RegularMesh +openmc.mesh.RegularMesh = RegularMesh From 5521569d852af486e32192843432a50a295dd6d6 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 26 Jan 2023 18:32:17 +0000 Subject: [PATCH 5/8] reduce duplication --- src/regular_mesh_plotter/core.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/src/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py index 9bfd2a4..3423d33 100644 --- a/src/regular_mesh_plotter/core.py +++ b/src/regular_mesh_plotter/core.py @@ -7,8 +7,8 @@ class RegularMesh(openmc.RegularMesh): def slice_of_data( self, dataset: np.ndarray, - slice_index: int = 0, view_direction: str = "z", + slice_index: typing.Optional[int] = None, volume_normalization: bool = True, ): """Obtains the dataset values on an axis aligned 2D slice through the @@ -19,11 +19,11 @@ def slice_of_data( datasets : numpy.ndarray 1-D array of data values. Will be reshaped to fill the mesh and should therefore have the same number of entries to fill the mesh - slice_index : int - The index of the mesh slice to extract. view_direction : str The axis to view the slice from. Supports negative and positive axis values. Acceptable values are 'x', 'y', 'z', '-x', '-y', '-z'. + slice_index : int + The index of the mesh slice to extract. volume_normalization : bool, optional Whether or not to normalize the data by the volume of the mesh elements. @@ -34,6 +34,11 @@ def slice_of_data( the 2D array of dataset values """ + bb_index_to_view_direction = {'x': 0, 'y': 1, 'z': 2} + + if slice_index is None: + slice_index = int(self.dimension[bb_index_to_view_direction[view_direction]] / 2) + if volume_normalization: dataset = dataset.flatten() / self.volumes.flatten() @@ -67,7 +72,7 @@ def slice_of_data( transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] aligned_ds = np.rot90(transposed_ds, 1) else: - msg = "view_direction is not one of the acceptable options {supported_view_dirs}" + msg = "view_direction of {view_direction} is not one of the acceptable options ({supported_view_dirs})" raise ValueError(msg) return aligned_ds @@ -110,21 +115,24 @@ def plot_slice( import matplotlib.pyplot as plt + bb_index_to_view_direction = {'x': 0, 'y': 1, 'z': 2} + bb_index = bb_index_to_view_direction[view_direction] + + if slice_index is None: + slice_index = int(self.dimension[bb_index] / 2) + # gets the axis labels and bounding box index if "x" in view_direction: x_label = "Y [cm]" y_label = "Z [cm]" - bb_index = 0 if "y" in view_direction: x_label = "X [cm]" y_label = "Z [cm]" - bb_index = 1 if "z" in view_direction: x_label = "X [cm]" y_label = "Y [cm]" - bb_index = 2 # selecting mid index on the mesh for the slice if slice_index is None: From cbd2b83cac71577442610ec574eb50722cc3ee8f Mon Sep 17 00:00:00 2001 From: shimwell Date: Thu, 26 Jan 2023 18:33:39 +0000 Subject: [PATCH 6/8] [skip ci] Apply formatting changes --- src/regular_mesh_plotter/core.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py index 3423d33..0ebc051 100644 --- a/src/regular_mesh_plotter/core.py +++ b/src/regular_mesh_plotter/core.py @@ -34,10 +34,12 @@ def slice_of_data( the 2D array of dataset values """ - bb_index_to_view_direction = {'x': 0, 'y': 1, 'z': 2} + bb_index_to_view_direction = {"x": 0, "y": 1, "z": 2} if slice_index is None: - slice_index = int(self.dimension[bb_index_to_view_direction[view_direction]] / 2) + slice_index = int( + self.dimension[bb_index_to_view_direction[view_direction]] / 2 + ) if volume_normalization: dataset = dataset.flatten() / self.volumes.flatten() @@ -115,7 +117,7 @@ def plot_slice( import matplotlib.pyplot as plt - bb_index_to_view_direction = {'x': 0, 'y': 1, 'z': 2} + bb_index_to_view_direction = {"x": 0, "y": 1, "z": 2} bb_index = bb_index_to_view_direction[view_direction] if slice_index is None: From ae0bc5ae5ec60489cea341c639fc54dc605a34c7 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 1 Feb 2023 18:12:54 +0000 Subject: [PATCH 7/8] reduce maintence burden --- README.md | 122 ++------ ...eate_statepoint_file_with_meshes_openmc.py | 93 ------ ...tatepoint_file_with_meshes_openmc_dagmc.py | 107 ------- examples/example.stl | 86 ------ examples/examples_from_readme.py | 77 ----- examples/plot_regular_mesh_dose_tally.py | 24 -- ...t_regular_mesh_dose_tally_with_geometry.py | 27 -- examples/plot_regular_mesh_tally.py | 20 -- .../plot_regular_mesh_tally_with_geometry.py | 45 ++- examples/plot_regular_mesh_values.py | 76 ----- .../plot_regular_mesh_values_with_geometry.py | 79 ----- src/regular_mesh_plotter/core.py | 284 +++++++----------- src/regular_mesh_plotter/utils.py | 95 ------ 13 files changed, 167 insertions(+), 968 deletions(-) delete mode 100644 examples/create_statepoint_file_with_meshes_openmc.py delete mode 100644 examples/create_statepoint_file_with_meshes_openmc_dagmc.py delete mode 100644 examples/example.stl delete mode 100644 examples/examples_from_readme.py delete mode 100644 examples/plot_regular_mesh_dose_tally.py delete mode 100644 examples/plot_regular_mesh_dose_tally_with_geometry.py delete mode 100644 examples/plot_regular_mesh_tally.py delete mode 100644 examples/plot_regular_mesh_values.py delete mode 100644 examples/plot_regular_mesh_values_with_geometry.py delete mode 100644 src/regular_mesh_plotter/utils.py diff --git a/README.md b/README.md index 2501d0c..01208dc 100644 --- a/README.md +++ b/README.md @@ -6,122 +6,34 @@ [![codecov](https://codecov.io/gh/fusion-energy/regular_mesh_plotter/branch/main/graph/badge.svg)](https://codecov.io/gh/fusion-energy/regular_mesh_plotter) -## A minimal Python package that plots 2D mesh tally results with the underlying DAGMC geometry +## A minimal Python package that extracts 2D mesh tally results for plotting convenience. -# Installation +This package is deployed on [xsplot.com](https://www.xsplot.com) as part of the ```openmc_plot``` suite of plotting apps -```bash -pip install regular_mesh_plotter -``` - -Mesh results in the form of at OpenMC.tally objects can be plotted with a single API call. - -A Matplotlib.pyplot object is returned by all functions so one can make changes -to the legend, axis, colour map etc. However some key options are accessible -in the function call directly. - -There are additional options that allow - -- rotation of the mesh tally results -- rotation of the DAGMC geometry slice -- saving the plot as an image file -- specifying contour lines TODO -- changing axis and colour bar labels -- changing colour scale applied -- truncation of values -- The plane_normal of the DAGMC geometry - -The resulting plots can be used to show dose maps, activation, reaction rate -and other mesh tally results. - -The examples below require a mesh tally that can be read in with OpenMC in the following way. - -```python -import openmc - -# loads in the statepoint file containing tallies -statepoint = openmc.StatePoint(filepath="statepoint.2.h5") +# Local install -# gets one tally from the available tallies -my_tally = statepoint.get_tally(name="neutron_effective_dose_on_2D_mesh_xy") -``` - -Example 1 shows a OpenMC tally plotted -```python - -import regular_mesh_plotter as rmp -import matplotlib.pyplot as plt +First you will need openmc installed, then you can install this package with pip -rmp.plot_regular_mesh_tally( - tally=my_tally, - std_dev_or_tally_value="tally_value", - x_label="X [cm]", - y_label="Y [cm]", -) - -plt.savefig('openmc_mesh_tally_plot.png') +```bash +pip install regular_mesh_plotter ``` +# Usage -Example 4 shows a OpenMC tally plotted with an underlying DAGMC geometry -```python -plot_regular_mesh_tally_with_geometry( - tally=my_tally, - dagmc_file_or_trimesh_object='dagmc.h5m', - std_dev_or_tally_value="tally_value", - x_label="X [cm]", - y_label="Y [cm]", -) -``` +The package can be used from within your own python script to make plots or via a GUI that is also bundled into the package install. -Example 5 shows how to rotate the underlying DAGMC geometry and mesh tally. -This is sometimes necessary as the slice and mesh can get out of alignment -when changing the plane normal -```python -plot_regular_mesh_tally_with_geometry( - tally, - dagmc_file_or_trimesh_object, - std_dev_or_tally_value="tally_value", - filename: Optional[str] = None, - scale=None, # LogNorm(), - vmin=None, - label="", - x_label="X [cm]", - y_label="Y [cm]", - plane_origin: List[float] = None, - plane_normal: List[float] = [0, 0, 1], - rotate_mesh: float = 0, - rotate_geometry: float = 0, - required_units=None, - source_strength: float = None, -) -``` +## Python API script usage -Example 6 shows how to plot a dose tally with the underlying DAGMC geometry. -This also includes unit conversion from the base tally units to the requested -units. -```python -plot_regular_mesh_dose_tally_with_geometry( - tally, - dagmc_file_or_trimesh_object, - filename: Optional[str] = None, - scale=None, # LogNorm(), - vmin=None, - label="", - x_label="X [cm]", - y_label="Y [cm]", - plane_origin: List[float] = None, - plane_normal: List[float] = [0, 0, 1], - rotate_mesh: float = 0, - rotate_geometry: float = 0, - required_units="picosievert / source_particle", - source_strength: float = None, - std_dev_or_tally_value: str = "tally_value", -): -``` +See the [examples folder](https://github.com/fusion-energy/regular_mesh_plotter/tree/master/examples) for example scripts -Additional examples can be found in the [examples folder in the GitHub repository](https://github.com/fusion-energy/regular_mesh_plotter/tree/main/examples) +## Graphical User Interface (GUI) usage +After installing run ```openmc_mesh_plotter``` command from the terminal and the GUI should launch in a new browser window. # Related packages +[openmc_plot](https://github.com/fusion-energy/openmc_plot)A single package that includes all the various plotters. + If you want to plot the DAGMC geometry without a mesh tally then take a look at the [dagmc_geometry_slice_plotter](https://github.com/fusion-energy/dagmc_geometry_slice_plotter) package + +If you want to plot the Native CSG geometry without a mesh tally then take a look at +the [dagmc_geometry_slice_plotter](https://github.com/fusion-energy/openmc_geometry_plot) package diff --git a/examples/create_statepoint_file_with_meshes_openmc.py b/examples/create_statepoint_file_with_meshes_openmc.py deleted file mode 100644 index 508015b..0000000 --- a/examples/create_statepoint_file_with_meshes_openmc.py +++ /dev/null @@ -1,93 +0,0 @@ -# This minimal example makes a 3D volume and exports the shape to a stp file - -import openmc -import openmc_dagmc_wrapper as odw -import openmc_plasma_source as ops -import openmc_data_downloader as odd - - -# MATERIALS -breeder_material = openmc.Material(1, "PbLi") # Pb84.2Li15.8 -# breeder_material.add_element("Pb", 84.2, percent_type="ao") -breeder_material.add_element( - "Li", - 15.8, - percent_type="ao", - enrichment=50.0, - enrichment_target="Li6", - enrichment_type="ao", -) # 50% enriched -breeder_material.set_density("atom/b-cm", 3.2720171e-2) # around 11 g/cm3 - -iron = openmc.Material(name="iron") -iron.set_density("g/cm3", 7.75) -iron.add_element("Li", 0.95, percent_type="wo") - -materials = openmc.Materials([breeder_material, iron]) - -odd.just_in_time_library_generator(libraries="TENDL-2019", materials=materials) - -# GEOMETRY - -# surfaces -vessel_inner = openmc.Sphere(r=500) -first_wall_outer_surface = openmc.Sphere(r=510) -breeder_blanket_outer_surface = openmc.Sphere(r=610, boundary_type="vacuum") - -# cells -inner_vessel_region = -vessel_inner -inner_vessel_cell = openmc.Cell(region=inner_vessel_region) - -first_wall_region = -first_wall_outer_surface & +vessel_inner -first_wall_cell = openmc.Cell(region=first_wall_region) -first_wall_cell.fill = iron - -breeder_blanket_region = +first_wall_outer_surface & -breeder_blanket_outer_surface -breeder_blanket_cell = openmc.Cell(region=breeder_blanket_region) -breeder_blanket_cell.fill = breeder_material - -universe = openmc.Universe( - cells=[inner_vessel_cell, first_wall_cell, breeder_blanket_cell] -) -geometry = openmc.Geometry(universe) - -tally1 = odw.MeshTally2D( - tally_type="neutron_effective_dose", - plane="xy", - mesh_resolution=(10, 5), - bounding_box=[(-100, -100, 0), (100, 100, 1)], -) - -tally2 = odw.MeshTally3D( - mesh_resolution=(100, 100, 100), - bounding_box=[(-100, -100, 0), (100, 100, 1)], - tally_type="neutron_effective_dose", -) - -tally3 = odw.MeshTally2D( - tally_type="neutron_flux", - plane="xy", - mesh_resolution=(10, 5), - bounding_box=[(-100, -100, 0), (100, 100, 1)], -) - -tallies = openmc.Tallies( - [ - tally1, - tally2, - tally3, - ] -) - -settings = odw.FusionSettings() -settings.batches = 3 -settings.particles = 1000000 -# assigns a ring source of DT energy neutrons to the source using the -# openmc_plasma_source package -settings.source = ops.FusionPointSource() - - -my_model = openmc.model.Model( - materials=materials, geometry=geometry, settings=settings, tallies=tallies -) -statepoint_file = my_model.run() diff --git a/examples/create_statepoint_file_with_meshes_openmc_dagmc.py b/examples/create_statepoint_file_with_meshes_openmc_dagmc.py deleted file mode 100644 index fed8b24..0000000 --- a/examples/create_statepoint_file_with_meshes_openmc_dagmc.py +++ /dev/null @@ -1,107 +0,0 @@ -# This minimal example makes a 3D volume and exports the shape to a stp file -# A surrounding volume called a graveyard is needed for neutronics simulations - -import openmc -import openmc_dagmc_wrapper as odw -import openmc_plasma_source as ops -from stl_to_h5m import stl_to_h5m -from dagmc_bounding_box import DagmcBoundingBox - -# code used to create example.stl -# import paramak -# my_shape = paramak.ExtrudeStraightShape( -# points=[(1, 1), (1, 200), (600, 200), (600, 1)], -# distance=180, -# ) -# my_shape.export_stl("example.stl") - -# This script converts the CAD stl files generated into h5m files that can be -# used in DAGMC enabled codes. h5m files created in this way are imprinted, -# merged, faceted and ready for use in OpenMC. One of the key aspects of this -# is the assignment of materials to the volumes present in the CAD files. - -stl_to_h5m( - files_with_tags=[("example.stl", "mat1")], - h5m_filename="dagmc.h5m", -) - -my_corners = DagmcBoundingBox("dagmc.h5m").corners() - -# makes use of the previously created neutronics geometry (h5m file) and assigns -# actual materials to the material tags. Sets simulation intensity and specifies -# the neutronics results to record (know as tallies). - -geometry = odw.Geometry( - h5m_filename="dagmc.h5m", -) - -materials = odw.Materials( - h5m_filename="dagmc.h5m", correspondence_dict={"mat1": "FLiNaK"} -) - -tally1 = odw.MeshTally2D( - tally_type="neutron_effective_dose", - plane="xy", - mesh_resolution=(10, 5), - bounding_box=my_corners, -) -tally2 = odw.MeshTally2D( - tally_type="neutron_effective_dose", - plane="yz", - mesh_resolution=(10, 5), - bounding_box=my_corners, -) -tally3 = odw.MeshTally2D( - tally_type="neutron_effective_dose", - plane="xz", - mesh_resolution=(10, 5), - bounding_box=my_corners, -) -tally4 = odw.MeshTally2D( - tally_type="neutron_effective_dose", - plane="xy", - mesh_resolution=(10, 10), - bounding_box=my_corners, -) -tally5 = odw.MeshTally2D( - tally_type="neutron_effective_dose", - plane="yz", - mesh_resolution=(10, 5), - bounding_box=my_corners, -) -tally6 = odw.MeshTally2D( - tally_type="neutron_effective_dose", - plane="xz", - mesh_resolution=(10, 5), - bounding_box=my_corners, -) - -# tally2 = odw.MeshTally3D( -# mesh_resolution=(100, 100, 100), -# bounding_box=my_corners, -# tally_type="neutron_effective_dose", -# ) - -tallies = openmc.Tallies( - [ - tally1, - tally2, - tally3, - tally4, - tally5, - tally6, - ] -) - -settings = odw.FusionSettings() -settings.batches = 2 -settings.particles = 1000 -# assigns a ring source of DT energy neutrons to the source using the -# openmc_plasma_source package -settings.source = ops.FusionPointSource() - - -my_model = openmc.Model( - materials=materials, geometry=geometry, settings=settings, tallies=tallies -) -statepoint_file = my_model.run() diff --git a/examples/example.stl b/examples/example.stl deleted file mode 100644 index ddc4f41..0000000 --- a/examples/example.stl +++ /dev/null @@ -1,86 +0,0 @@ -solid - facet normal -1.000000e+00 0.000000e+00 -0.000000e+00 - outer loop - vertex 1.000000e+00 9.000000e+01 2.000000e+02 - vertex 1.000000e+00 9.000000e+01 1.000000e+00 - vertex 1.000000e+00 -9.000000e+01 1.000000e+00 - endloop - endfacet - facet normal -1.000000e+00 0.000000e+00 0.000000e+00 - outer loop - vertex 1.000000e+00 -9.000000e+01 2.000000e+02 - vertex 1.000000e+00 9.000000e+01 2.000000e+02 - vertex 1.000000e+00 -9.000000e+01 1.000000e+00 - endloop - endfacet - facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 - outer loop - vertex 1.000000e+00 -9.000000e+01 1.000000e+00 - vertex 6.000000e+02 9.000000e+01 1.000000e+00 - vertex 6.000000e+02 -9.000000e+01 1.000000e+00 - endloop - endfacet - facet normal 0.000000e+00 0.000000e+00 -1.000000e+00 - outer loop - vertex 1.000000e+00 -9.000000e+01 1.000000e+00 - vertex 1.000000e+00 9.000000e+01 1.000000e+00 - vertex 6.000000e+02 9.000000e+01 1.000000e+00 - endloop - endfacet - facet normal 0.000000e+00 -0.000000e+00 1.000000e+00 - outer loop - vertex 6.000000e+02 -9.000000e+01 2.000000e+02 - vertex 6.000000e+02 9.000000e+01 2.000000e+02 - vertex 1.000000e+00 9.000000e+01 2.000000e+02 - endloop - endfacet - facet normal 0.000000e+00 0.000000e+00 1.000000e+00 - outer loop - vertex 6.000000e+02 -9.000000e+01 2.000000e+02 - vertex 1.000000e+00 9.000000e+01 2.000000e+02 - vertex 1.000000e+00 -9.000000e+01 2.000000e+02 - endloop - endfacet - facet normal 0.000000e+00 1.000000e+00 0.000000e+00 - outer loop - vertex 6.000000e+02 9.000000e+01 2.000000e+02 - vertex 1.000000e+00 9.000000e+01 1.000000e+00 - vertex 1.000000e+00 9.000000e+01 2.000000e+02 - endloop - endfacet - facet normal 0.000000e+00 1.000000e+00 0.000000e+00 - outer loop - vertex 6.000000e+02 9.000000e+01 2.000000e+02 - vertex 6.000000e+02 9.000000e+01 1.000000e+00 - vertex 1.000000e+00 9.000000e+01 1.000000e+00 - endloop - endfacet - facet normal 1.000000e+00 0.000000e+00 -0.000000e+00 - outer loop - vertex 6.000000e+02 9.000000e+01 1.000000e+00 - vertex 6.000000e+02 9.000000e+01 2.000000e+02 - vertex 6.000000e+02 -9.000000e+01 2.000000e+02 - endloop - endfacet - facet normal 1.000000e+00 0.000000e+00 0.000000e+00 - outer loop - vertex 6.000000e+02 -9.000000e+01 1.000000e+00 - vertex 6.000000e+02 9.000000e+01 1.000000e+00 - vertex 6.000000e+02 -9.000000e+01 2.000000e+02 - endloop - endfacet - facet normal -0.000000e+00 -1.000000e+00 0.000000e+00 - outer loop - vertex 6.000000e+02 -9.000000e+01 2.000000e+02 - vertex 1.000000e+00 -9.000000e+01 2.000000e+02 - vertex 1.000000e+00 -9.000000e+01 1.000000e+00 - endloop - endfacet - facet normal 0.000000e+00 -1.000000e+00 -0.000000e+00 - outer loop - vertex 6.000000e+02 -9.000000e+01 2.000000e+02 - vertex 1.000000e+00 -9.000000e+01 1.000000e+00 - vertex 6.000000e+02 -9.000000e+01 1.000000e+00 - endloop - endfacet -endsolid diff --git a/examples/examples_from_readme.py b/examples/examples_from_readme.py deleted file mode 100644 index 1cf2073..0000000 --- a/examples/examples_from_readme.py +++ /dev/null @@ -1,77 +0,0 @@ -# this examples requires additional python packages - -from dagmc_geometry_slice_plotter import plot_slice_of_dagmc_geometry -from stl_to_h5m import stl_to_h5m -import matplotlib.pyplot as plt -import paramak - -my_reactor = paramak.SubmersionTokamak( - inner_bore_radial_thickness=30, - inboard_tf_leg_radial_thickness=30, - center_column_shield_radial_thickness=30, - divertor_radial_thickness=80, - inner_plasma_gap_radial_thickness=50, - plasma_radial_thickness=200, - outer_plasma_gap_radial_thickness=50, - firstwall_radial_thickness=30, - blanket_rear_wall_radial_thickness=30, - number_of_tf_coils=16, - rotation_angle=360, - support_radial_thickness=90, - inboard_blanket_radial_thickness=30, - outboard_blanket_radial_thickness=30, - elongation=2.00, - triangularity=0.50, - pf_coil_case_thicknesses=[10, 10, 10, 10], - pf_coil_radial_thicknesses=[20, 50, 50, 20], - pf_coil_vertical_thicknesses=[20, 50, 50, 20], - pf_coil_radial_position=[500, 550, 550, 500], - pf_coil_vertical_position=[270, 100, -100, -270], - rear_blanket_to_tf_gap=50, - outboard_tf_coil_radial_thickness=30, - outboard_tf_coil_poloidal_thickness=30, -) - - -stl_filenames = my_reactor.export_stl() - - -stl_to_h5m( - files_with_tags=[ - (stl_filename, name) - for name, stl_filename in zip(my_reactor.name, stl_filenames) - ], - h5m_filename="dagmc.h5m", -) - - -plot_slice_of_dagmc_geometry( - dagmc_file_or_trimesh_object="dagmc.h5m", - plane_normal=[0, 0, 1], - output_filename="my_plot1.png", -) - - -plot_slice_of_dagmc_geometry( - dagmc_file_or_trimesh_object="dagmc.h5m", - plane_origin=[0, 0, 300], - plane_normal=[0, 0, 1], - output_filename="my_plot2.png", -) - - -plot_slice_of_dagmc_geometry( - dagmc_file_or_trimesh_object="dagmc.h5m", - plane_normal=[0, 1, 0], - rotate_plot=45, - output_filename="my_plot3.png", -) - -plt.savefig("big.png", dpi=600) - -plot_slice_of_dagmc_geometry( - dagmc_file_or_trimesh_object="dagmc.h5m", - plane_normal=[1, 0, 0], - rotate_plot=270, - output_filename="my_plot4.png", -) diff --git a/examples/plot_regular_mesh_dose_tally.py b/examples/plot_regular_mesh_dose_tally.py deleted file mode 100644 index 9bcbf83..0000000 --- a/examples/plot_regular_mesh_dose_tally.py +++ /dev/null @@ -1,24 +0,0 @@ -import regular_mesh_plotter as rmp -import openmc -import matplotlib.pyplot as plt - -# loads in the statepoint file containing tallies -statepoint = openmc.StatePoint(filepath="statepoint.2.h5") - -# gets one tally from the available tallies -my_tally = statepoint.get_tally(name="neutron_effective_dose_on_2D_mesh_xy") - -# creates a plot of the mesh tally -rmp.plot_regular_mesh_dose_tally( - tally=my_tally, # the openmc tally object to plot, must be a 2d mesh tally - filename="plot_regular_mesh_dose_tally.png", # the filename of the picture file saved - x_label="X [cm]", - y_label="Y [cm]", - rotate_plot=0, - required_units="picosievert / source_particle", - source_strength=None, - label="Effective dose [picosievert / source_particle]", -) - -# displays the plot -plt.show() diff --git a/examples/plot_regular_mesh_dose_tally_with_geometry.py b/examples/plot_regular_mesh_dose_tally_with_geometry.py deleted file mode 100644 index 40b6f7d..0000000 --- a/examples/plot_regular_mesh_dose_tally_with_geometry.py +++ /dev/null @@ -1,27 +0,0 @@ -import regular_mesh_plotter as rmp -import openmc -import matplotlib.pyplot as plt - -# loads in the statepoint file containing tallies -statepoint = openmc.StatePoint(filepath="statepoint.2.h5") - -# gets one tally from the available tallies -my_tally = statepoint.get_tally(name="neutron_effective_dose_on_2D_mesh_xy") - -# creates a plot of the mesh -rmp.plot_regular_mesh_dose_tally_with_geometry( - tally=my_tally, - dagmc_file_or_trimesh_object="dagmc.h5m", - filename="plot_regular_mesh_dose_tally_with_geometry.png", - label="Effective dose [picosievert / second]", - x_label="X [cm]", - y_label="Y [cm]", - plane_origin=None, # this could be skipped as it defaults to None, which uses the center of the mesh - plane_normal=[0, 0, 1], - rotate_mesh=0, - rotate_geometry=0, - required_units="picosievert / second", - source_strength=1e20, -) - -plt.show() diff --git a/examples/plot_regular_mesh_tally.py b/examples/plot_regular_mesh_tally.py deleted file mode 100644 index 3f99d29..0000000 --- a/examples/plot_regular_mesh_tally.py +++ /dev/null @@ -1,20 +0,0 @@ -import regular_mesh_plotter as rmp -import openmc - -# loads in the statepoint file containing tallies -statepoint = openmc.StatePoint(filepath="statepoint.2.h5") - -# gets one tally from the available tallies -my_tally = statepoint.get_tally(name="neutron_effective_dose_on_2D_mesh_xy") - -# creates a plot of the mesh -rmp.plot_regular_mesh_tally( - tally=my_tally, - filename="neutron_effective_dose_on_2D_mesh_xy.png", - label="", - x_label="X [cm]", - y_label="Y [cm]", - # rotate_plot: float = 0, - # required_units: str = None, - # source_strength: float = None, -) diff --git a/examples/plot_regular_mesh_tally_with_geometry.py b/examples/plot_regular_mesh_tally_with_geometry.py index dab528c..1d59781 100644 --- a/examples/plot_regular_mesh_tally_with_geometry.py +++ b/examples/plot_regular_mesh_tally_with_geometry.py @@ -1,14 +1,14 @@ -import regular_mesh_plotter import openmc import openmc_geometry_plot # extends openmc.Geometry class with plotting functions +import regular_mesh_plotter # extends openmc.Mesh class with plotting functions # MATERIALS mat_1 = openmc.Material() mat_1.add_element("Li", 1) -mat_1.set_density("g/cm3", 3.2720171e-2) # around 11 g/cm3 +mat_1.set_density("g/cm3", 0.5) my_materials = openmc.Materials([mat_1]) @@ -41,7 +41,7 @@ # Create a DT point source source = openmc.Source() -source.space = openmc.stats.Point((0, 0, 0)) +source.space = openmc.stats.Point((100, 0, 0)) source.angle = openmc.stats.Isotropic() source.energy = openmc.stats.Discrete([14e6], [1]) my_settings.source = source @@ -51,13 +51,13 @@ mesh = openmc.RegularMesh().from_domain( my_geometry, # the corners of the mesh are being set automatically to surround the geometry - dimension=[10, 10, 10], + dimension=[20, 20, 20], ) mesh_filter = openmc.MeshFilter(mesh) mesh_tally_1 = openmc.Tally(name="mesh_tally") mesh_tally_1.filters = [mesh_filter] -mesh_tally_1.scores = ["flux"] +mesh_tally_1.scores = ["absorption"] my_tallies.append(mesh_tally_1) model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies) @@ -68,13 +68,36 @@ # extracts the mesh tally by name my_mesh_tally = statepoint.get_tally(name="mesh_tally") -# converts the tally result into a VTK file -plot = mesh.plot_slice(dataset=my_mesh_tally.mean, slice_index=5, view_direction="x") -plot.figure.savefig("test.png") +# gets a 2d slice of data to later plot +data_slice = mesh.slice_of_data( + dataset=my_mesh_tally.mean, view_direction="x" +) # this section adds a contour line according the the material ids -material_ids = my_geometry.get_slice_of_material_ids(view_direction="x") -plot2 = my_geometry.get_outline_contour(outline_data=material_ids, view_direction="x") +my_geometry.view_direction="x" +material_ids = my_geometry.get_slice_of_material_ids() + + import matplotlib.pyplot as plt +import numpy as np -plt.savefig("test2.png") +plt.imshow( + data_slice, + extent=mesh.get_mpl_plot_extent(), + interpolation="none", +) + +# gets unique levels for outlines contour plot and for the color scale +levels = np.unique([item for sublist in material_ids for item in sublist]) + +plt.contour( + material_ids, + origin="upper", + colors="k", + linestyles="solid", + levels=levels, + linewidths=0.5, + extent=my_geometry.get_mpl_plot_extent(), +) +plt.show() +# plt.savefig("mesh_with_geometry.png") diff --git a/examples/plot_regular_mesh_values.py b/examples/plot_regular_mesh_values.py deleted file mode 100644 index 1514f57..0000000 --- a/examples/plot_regular_mesh_values.py +++ /dev/null @@ -1,76 +0,0 @@ -import regular_mesh_plotter as rmp -import numpy as np - -values = np.array( - [ - [ - 1.64565869e09, - 2.59505642e09, - 3.06732422e09, - 2.64141597e09, - 1.76520128e09, - 2.14909584e09, - 3.70513916e09, - 4.99737560e09, - 3.86795536e09, - 2.21272467e09, - ], - [ - 2.92406594e09, - 5.91396360e09, - 9.39595883e09, - 5.94102629e09, - 2.78174231e09, - 3.47563407e09, - 1.02570496e10, - 2.50416310e10, - 1.01248003e10, - 3.34937674e09, - ], - [ - 3.99684099e09, - 1.72147289e10, - 3.68431465e11, - 1.64968908e10, - 3.74168705e09, - 3.90640820e09, - 1.65505774e10, - 3.66837062e11, - 1.65568070e10, - 3.65646431e09, - ], - [ - 3.27796129e09, - 1.01456714e10, - 2.45757058e10, - 1.00180127e10, - 3.45105436e09, - 2.63214911e09, - 5.95924816e09, - 9.43589769e09, - 5.71740137e09, - 2.75163850e09, - ], - [ - 2.37324680e09, - 3.77339226e09, - 5.01889988e09, - 3.58500172e09, - 2.16754228e09, - 1.81599509e09, - 2.57229036e09, - 3.09622197e09, - 2.50136006e09, - 1.72280196e09, - ], - ] -) - -# creates a plot of the array -rmp.plot_regular_mesh_values( - values=values, - filename="plot_regular_mesh_values.png", - label="legend label", - x_label="X [cm]", - y_label="Y [cm]", -) diff --git a/examples/plot_regular_mesh_values_with_geometry.py b/examples/plot_regular_mesh_values_with_geometry.py deleted file mode 100644 index 80ea5c4..0000000 --- a/examples/plot_regular_mesh_values_with_geometry.py +++ /dev/null @@ -1,79 +0,0 @@ -import regular_mesh_plotter as rmp -import numpy as np - -values = np.array( - [ - [ - 1.64565869e09, - 2.59505642e09, - 3.06732422e09, - 2.64141597e09, - 1.76520128e09, - 2.14909584e09, - 3.70513916e09, - 4.99737560e09, - 3.86795536e09, - 2.21272467e09, - ], - [ - 2.92406594e09, - 5.91396360e09, - 9.39595883e09, - 5.94102629e09, - 2.78174231e09, - 3.47563407e09, - 1.02570496e10, - 2.50416310e10, - 1.01248003e10, - 3.34937674e09, - ], - [ - 3.99684099e09, - 1.72147289e10, - 3.68431465e11, - 1.64968908e10, - 3.74168705e09, - 3.90640820e09, - 1.65505774e10, - 3.66837062e11, - 1.65568070e10, - 3.65646431e09, - ], - [ - 3.27796129e09, - 1.01456714e10, - 2.45757058e10, - 1.00180127e10, - 3.45105436e09, - 2.63214911e09, - 5.95924816e09, - 9.43589769e09, - 5.71740137e09, - 2.75163850e09, - ], - [ - 2.37324680e09, - 3.77339226e09, - 5.01889988e09, - 3.58500172e09, - 2.16754228e09, - 1.81599509e09, - 2.57229036e09, - 3.09622197e09, - 2.50136006e09, - 1.72280196e09, - ], - ] -) - - -rmp.plot_regular_mesh_values_with_geometry( - values=values, - dagmc_file_or_trimesh_object="example.stl", - filename="plot_regular_mesh_values_with_geometry.png", - label="legend label", - x_label="X [cm]", - y_label="Y [cm]", - plane_normal=[1, 0, 0], - extent=[0, 700, 0, 300], -) diff --git a/src/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py index 0ebc051..5c3129b 100644 --- a/src/regular_mesh_plotter/core.py +++ b/src/regular_mesh_plotter/core.py @@ -3,173 +3,121 @@ import openmc -class RegularMesh(openmc.RegularMesh): - def slice_of_data( - self, - dataset: np.ndarray, - view_direction: str = "z", - slice_index: typing.Optional[int] = None, - volume_normalization: bool = True, - ): - """Obtains the dataset values on an axis aligned 2D slice through the - mesh. Useful for producing plots of slice data - - Parameters - ---------- - datasets : numpy.ndarray - 1-D array of data values. Will be reshaped to fill the mesh and - should therefore have the same number of entries to fill the mesh - view_direction : str - The axis to view the slice from. Supports negative and positive axis - values. Acceptable values are 'x', 'y', 'z', '-x', '-y', '-z'. - slice_index : int - The index of the mesh slice to extract. - volume_normalization : bool, optional - Whether or not to normalize the data by the volume of the mesh - elements. - - Returns - ------- - np.array() - the 2D array of dataset values - """ - - bb_index_to_view_direction = {"x": 0, "y": 1, "z": 2} - - if slice_index is None: - slice_index = int( - self.dimension[bb_index_to_view_direction[view_direction]] / 2 - ) - - if volume_normalization: - dataset = dataset.flatten() / self.volumes.flatten() - - reshaped_ds = dataset.reshape(self.dimension, order="F") - - if view_direction == "x": - # vertical axis is z, horizontal axis is -y - transposed_ds = reshaped_ds.transpose(0, 1, 2)[slice_index] - rotated_ds = np.rot90(transposed_ds, 1) - aligned_ds = np.fliplr(rotated_ds) - elif view_direction == "-x": - # vertical axis is z, horizontal axis is y - transposed_ds = reshaped_ds.transpose(0, 1, 2)[slice_index] - aligned_ds = np.rot90(transposed_ds, 1) - elif view_direction == "y": - # vertical axis is z, horizontal axis is x - transposed_ds = reshaped_ds.transpose(1, 2, 0)[slice_index] - aligned_ds = np.flipud(transposed_ds) - elif view_direction == "-y": - # vertical axis is z, horizontal axis is -x - transposed_ds = reshaped_ds.transpose(1, 2, 0)[slice_index] - aligned_ds = np.flipud(transposed_ds) - aligned_ds = np.fliplr(aligned_ds) - elif view_direction == "z": - # vertical axis is y, horizontal axis is -x - transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] - aligned_ds = np.rot90(transposed_ds, 1) - aligned_ds = np.fliplr(aligned_ds) - elif view_direction == "-z": - # vertical axis is y, horizontal axis is x - transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] - aligned_ds = np.rot90(transposed_ds, 1) - else: - msg = "view_direction of {view_direction} is not one of the acceptable options ({supported_view_dirs})" - raise ValueError(msg) - - return aligned_ds - - def plot_slice( - self, - dataset: np.ndarray, - slice_index: typing.Optional[int] = None, - view_direction: str = "z", - axes: typing.Optional["matplotlib.Axes"] = None, - volume_normalization: bool = True, - **kwargs, - ): - """Create a slice plot of the dataset on the RegularMesh. - - Parameters - ---------- - datasets : numpy.ndarray - 1-D array of data values. Will be reshaped to fill the mesh and - should therefore have the same number of entries to fill the mesh - slice_index : int - The index of the mesh slice to extract. If not set then the index - that slices through the middle of the mesh will be used - view_direction : str - The axis to view the slice from. Supports negative and positive axis - values. Acceptable values are 'x', 'y', 'z', '-x', '-y', '-z'. - axes : matplotlib.Axes, optional - Axes to draw to - volume_normalization : bool, optional - Whether or not to normalize the data by the volume of the mesh - elements. - **kwargs - Keyword arguments passed to :func:`matplotlib.pyplot.imshow` - - Returns - ------- - matplotlib.image.AxesImage - Resulting image - """ - - import matplotlib.pyplot as plt - - bb_index_to_view_direction = {"x": 0, "y": 1, "z": 2} - bb_index = bb_index_to_view_direction[view_direction] - - if slice_index is None: - slice_index = int(self.dimension[bb_index] / 2) - - # gets the axis labels and bounding box index - if "x" in view_direction: - x_label = "Y [cm]" - y_label = "Z [cm]" - - if "y" in view_direction: - x_label = "X [cm]" - y_label = "Z [cm]" - - if "z" in view_direction: - x_label = "X [cm]" - y_label = "Y [cm]" - - # selecting mid index on the mesh for the slice - if slice_index is None: - slice_index = int(self.dimension[bb_index] / 2) - - # slice_of_data also checks the view_direction is acceptable - image_slice = self.slice_of_data( - dataset=dataset, - slice_index=slice_index, - view_direction=view_direction, - volume_normalization=volume_normalization, +def get_mpl_plot_extent(self, view_direction:str= 'x'): + """Returns the (x_min, x_max, y_min, y_max) of the bounding box. The + view_direction is taken into account and can be set using + openmc.Geometry.view_direction property is taken into account and can be + set to 'x', 'y' or 'z'.""" + + bb = (self.lower_left, self.upper_right) + + x_min = self.get_side_extent("left",view_direction, bb) + x_max = self.get_side_extent("right",view_direction, bb) + y_min = self.get_side_extent("bottom",view_direction, bb) + y_max = self.get_side_extent("top",view_direction, bb) + + return (x_min, x_max, y_min, y_max) + + +def get_side_extent(self, side: str, view_direction:str= 'x', bb=None): + + if bb is None: + bb = (self.lower_left, self.upper_right) + + avail_extents = {} + avail_extents[("left", "x")] = bb[0][1] + avail_extents[("right", "x")] = bb[1][1] + avail_extents[("top", "x")] = bb[1][2] + avail_extents[("bottom", "x")] = bb[0][2] + avail_extents[("left", "y")] = bb[0][0] + avail_extents[("right", "y")] = bb[1][0] + avail_extents[("top", "y")] = bb[1][2] + avail_extents[("bottom", "y")] = bb[0][2] + avail_extents[("left", "z")] = bb[0][0] + avail_extents[("right", "z")] = bb[1][0] + avail_extents[("top", "z")] = bb[1][1] + avail_extents[("bottom", "z")] = bb[0][1] + return avail_extents[(side, view_direction)] + +def slice_of_data( + self, + dataset: np.ndarray, + view_direction: str = "z", + slice_index: typing.Optional[int] = None, + volume_normalization: bool = True, +): + """Obtains the dataset values on an axis aligned 2D slice through the + mesh. Useful for producing plots of slice data + + Parameters + ---------- + datasets : numpy.ndarray + 1-D array of data values. Will be reshaped to fill the mesh and + should therefore have the same number of entries to fill the mesh + view_direction : str + The axis to view the slice from. Supports negative and positive axis + values. Acceptable values are 'x', 'y', 'z', '-x', '-y', '-z'. + slice_index : int + The index of the mesh slice to extract. + volume_normalization : bool, optional + Whether or not to normalize the data by the volume of the mesh + elements. + + Returns + ------- + np.array() + the 2D array of dataset values + """ + + bb_index_to_view_direction = {"x": 0, "y": 1, "z": 2} + + if slice_index is None: + slice_index = int( + self.dimension[bb_index_to_view_direction[view_direction]] / 2 ) - # gets the extent of the plot - x_min = self.lower_left[bb_index] - x_max = self.upper_right[bb_index] - y_min = self.lower_left[bb_index] - y_max = self.upper_right[bb_index] - - if axes is None: - fig, axes = plt.subplots() - axes.set_xlabel(x_label) - axes.set_ylabel(y_label) - axes.set_title(f"View direction {view_direction}") - - return axes.imshow( - X=image_slice, extent=(x_min, x_max, y_min, y_max), aspect="auto", **kwargs - ) - - -# def get_cell_ids_for_regularmesh_slice(mesh, geometry): -# loop through the centroids of the mesh -# find the material at each point -# build up a pixel map of material ids - - -openmc.RegularMesh = RegularMesh -openmc.mesh.RegularMesh = RegularMesh + if volume_normalization: + dataset = dataset.flatten() / self.volumes.flatten() + + reshaped_ds = dataset.reshape(self.dimension, order="F") + + if view_direction == "x": + # vertical axis is z, horizontal axis is -y + transposed_ds = reshaped_ds.transpose(0, 1, 2)[slice_index] + rotated_ds = np.rot90(transposed_ds, 1) + aligned_ds = np.fliplr(rotated_ds) + elif view_direction == "-x": + # vertical axis is z, horizontal axis is y + transposed_ds = reshaped_ds.transpose(0, 1, 2)[slice_index] + aligned_ds = np.rot90(transposed_ds, 1) + elif view_direction == "y": + # vertical axis is z, horizontal axis is x + transposed_ds = reshaped_ds.transpose(1, 2, 0)[slice_index] + aligned_ds = np.flipud(transposed_ds) + elif view_direction == "-y": + # vertical axis is z, horizontal axis is -x + transposed_ds = reshaped_ds.transpose(1, 2, 0)[slice_index] + aligned_ds = np.flipud(transposed_ds) + aligned_ds = np.fliplr(aligned_ds) + elif view_direction == "z": + # vertical axis is y, horizontal axis is -x + transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] + aligned_ds = np.rot90(transposed_ds, 1) + aligned_ds = np.fliplr(aligned_ds) + elif view_direction == "-z": + # vertical axis is y, horizontal axis is x + transposed_ds = reshaped_ds.transpose(2, 0, 1)[slice_index] + aligned_ds = np.rot90(transposed_ds, 1) + else: + msg = "view_direction of {view_direction} is not one of the acceptable options ({supported_view_dirs})" + raise ValueError(msg) + + return aligned_ds + + +openmc.RegularMesh.slice_of_data = slice_of_data +openmc.mesh.RegularMesh.slice_of_data = slice_of_data +openmc.RegularMesh.get_mpl_plot_extent = get_mpl_plot_extent +openmc.mesh.RegularMesh.get_mpl_plot_extent = get_mpl_plot_extent +openmc.RegularMesh.get_side_extent = get_side_extent +openmc.mesh.RegularMesh.get_side_extent = get_side_extent diff --git a/src/regular_mesh_plotter/utils.py b/src/regular_mesh_plotter/utils.py deleted file mode 100644 index a8c11fc..0000000 --- a/src/regular_mesh_plotter/utils.py +++ /dev/null @@ -1,95 +0,0 @@ -from pathlib import Path -from typing import List, Optional, Tuple - -import matplotlib.pyplot as plt -import numpy as np -import openmc -import pandas as pd - - -def reshape_values_to_mesh_shape(tally, values): - tally_filter = tally.find_filter(filter_type=openmc.MeshFilter) - shape = tally_filter.mesh.dimension.tolist() - # 2d mesh has a shape in the form [1, 400, 400] - if 1 in shape: - shape.remove(1) - return values.reshape(shape) - - -def get_tally_extent( - tally, -): - - for filter in tally.filters: - if isinstance(filter, openmc.MeshFilter): - mesh_filter = filter - - extent_x = ( - min(mesh_filter.mesh.lower_left[0], mesh_filter.mesh.upper_right[0]), - max(mesh_filter.mesh.lower_left[0], mesh_filter.mesh.upper_right[0]), - ) - extent_y = ( - min(mesh_filter.mesh.lower_left[1], mesh_filter.mesh.upper_right[1]), - max(mesh_filter.mesh.lower_left[1], mesh_filter.mesh.upper_right[1]), - ) - extent_z = ( - min(mesh_filter.mesh.lower_left[2], mesh_filter.mesh.upper_right[2]), - max(mesh_filter.mesh.lower_left[2], mesh_filter.mesh.upper_right[2]), - ) - - if 1 in mesh_filter.mesh.dimension.tolist(): - print("2d mesh tally") - index_of_1d = mesh_filter.mesh.dimension.tolist().index(1) - print("index", index_of_1d) - if index_of_1d == 0: - return extent_y + extent_z - if index_of_1d == 1: - return extent_x + extent_z - if index_of_1d == 2: - return extent_x + extent_y - return None - - -def get_values_from_tally(tally): - """Return a numpy array of the openmc tally values (mean entry in - dataframe) and if present the standard deviation (std. dev. in the - dateframe) is also returned""" - - data_frame = tally.get_pandas_dataframe() - if "std. dev." in get_data_frame_columns(data_frame): - values = (np.array(data_frame["mean"]), np.array(data_frame["std. dev."])) - else: - values = np.array(data_frame["mean"]) - return values - - -def get_data_frame_columns(data_frame): - if isinstance(data_frame.columns, pd.MultiIndex): - data_frame_columns = data_frame.columns.get_level_values(0).to_list() - else: - data_frame_columns = data_frame.columns.to_list() - return data_frame_columns - - -def get_std_dev_or_value_from_tally(tally, values, std_dev_or_tally_value): - - if std_dev_or_tally_value == "std_dev": - value_index = 1 - elif std_dev_or_tally_value == "tally_value": - value_index = 0 - else: - msg = f'Value of std_dev_or_tally_value should be either "std_dev" or "value", not {type(values)}' - raise ValueError(msg) - - if isinstance(values, tuple): - value = reshape_values_to_mesh_shape(tally, values[value_index]) - elif isinstance(values, np.ndarray): - value = reshape_values_to_mesh_shape(tally, values) - else: # isinstance(values, np.ndarray): - # a pint unit object - value = reshape_values_to_mesh_shape(tally, values.magnitude) - # else: - # msg = f'Values to plot should be a numpy ndarry or a tuple or numpy ndarrys not a {type(values)}' - # raise ValueError(msg) - - return value From 95f2f4d25ce2382c9fa3277dd0733bb79bfe88c3 Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 1 Feb 2023 18:14:36 +0000 Subject: [PATCH 8/8] [skip ci] Apply formatting changes --- examples/plot_regular_mesh_tally_with_geometry.py | 8 +++----- src/regular_mesh_plotter/core.py | 14 +++++++------- tests/test_plot_regular_mesh_values.py | 3 --- .../test_plot_regular_mesh_values_with_geometry.py | 2 -- 4 files changed, 10 insertions(+), 17 deletions(-) diff --git a/examples/plot_regular_mesh_tally_with_geometry.py b/examples/plot_regular_mesh_tally_with_geometry.py index 1d59781..bf9183a 100644 --- a/examples/plot_regular_mesh_tally_with_geometry.py +++ b/examples/plot_regular_mesh_tally_with_geometry.py @@ -1,7 +1,7 @@ import openmc import openmc_geometry_plot # extends openmc.Geometry class with plotting functions -import regular_mesh_plotter # extends openmc.Mesh class with plotting functions +import regular_mesh_plotter # extends openmc.Mesh class with plotting functions # MATERIALS @@ -69,12 +69,10 @@ my_mesh_tally = statepoint.get_tally(name="mesh_tally") # gets a 2d slice of data to later plot -data_slice = mesh.slice_of_data( - dataset=my_mesh_tally.mean, view_direction="x" -) +data_slice = mesh.slice_of_data(dataset=my_mesh_tally.mean, view_direction="x") # this section adds a contour line according the the material ids -my_geometry.view_direction="x" +my_geometry.view_direction = "x" material_ids = my_geometry.get_slice_of_material_ids() diff --git a/src/regular_mesh_plotter/core.py b/src/regular_mesh_plotter/core.py index 5c3129b..3703fbb 100644 --- a/src/regular_mesh_plotter/core.py +++ b/src/regular_mesh_plotter/core.py @@ -3,7 +3,7 @@ import openmc -def get_mpl_plot_extent(self, view_direction:str= 'x'): +def get_mpl_plot_extent(self, view_direction: str = "x"): """Returns the (x_min, x_max, y_min, y_max) of the bounding box. The view_direction is taken into account and can be set using openmc.Geometry.view_direction property is taken into account and can be @@ -11,16 +11,15 @@ def get_mpl_plot_extent(self, view_direction:str= 'x'): bb = (self.lower_left, self.upper_right) - x_min = self.get_side_extent("left",view_direction, bb) - x_max = self.get_side_extent("right",view_direction, bb) - y_min = self.get_side_extent("bottom",view_direction, bb) - y_max = self.get_side_extent("top",view_direction, bb) + x_min = self.get_side_extent("left", view_direction, bb) + x_max = self.get_side_extent("right", view_direction, bb) + y_min = self.get_side_extent("bottom", view_direction, bb) + y_max = self.get_side_extent("top", view_direction, bb) return (x_min, x_max, y_min, y_max) -def get_side_extent(self, side: str, view_direction:str= 'x', bb=None): - +def get_side_extent(self, side: str, view_direction: str = "x", bb=None): if bb is None: bb = (self.lower_left, self.upper_right) @@ -39,6 +38,7 @@ def get_side_extent(self, side: str, view_direction:str= 'x', bb=None): avail_extents[("bottom", "z")] = bb[0][1] return avail_extents[(side, view_direction)] + def slice_of_data( self, dataset: np.ndarray, diff --git a/tests/test_plot_regular_mesh_values.py b/tests/test_plot_regular_mesh_values.py index 91cf03e..dde60d1 100644 --- a/tests/test_plot_regular_mesh_values.py +++ b/tests/test_plot_regular_mesh_values.py @@ -8,7 +8,6 @@ class TestPlotRegularMeshValues(unittest.TestCase): def setUp(self): - self.values = np.array( [ [ @@ -75,11 +74,9 @@ def setUp(self): ) def test_plot_regular_mesh_values(self): - plot_regular_mesh_values(values=self.values) def test_plot_regular_mesh_values_with_output(self): - os.system("rm test.png") plot_regular_mesh_values(values=self.values, filename="test.png") diff --git a/tests/test_plot_regular_mesh_values_with_geometry.py b/tests/test_plot_regular_mesh_values_with_geometry.py index 28b8717..8bbf77a 100644 --- a/tests/test_plot_regular_mesh_values_with_geometry.py +++ b/tests/test_plot_regular_mesh_values_with_geometry.py @@ -9,7 +9,6 @@ class TestPlotRegularMeshValuesWithGeometry(unittest.TestCase): def setUp(self): - self.values = np.array( [ [ @@ -76,7 +75,6 @@ def setUp(self): ) def test_plot_regular_mesh_values_with_geometry(self): - plot_regular_mesh_values_with_geometry( values=self.values, dagmc_file_or_trimesh_object="tests/example.stl", # this could be a h5m file