diff --git a/.github/workflows/ci_with_install.yml b/.github/workflows/ci_with_install.yml index fbc3ce0..6fb06be 100644 --- a/.github/workflows/ci_with_install.yml +++ b/.github/workflows/ci_with_install.yml @@ -20,7 +20,8 @@ jobs: testing: runs-on: ubuntu-latest container: - image: continuumio/miniconda3:4.9.2 + # image: openmc/openmc:latest-dagmc is missing the dagmc universe + image: openmc/openmc:develop-dagmc steps: - name: Checkout repository uses: actions/checkout@v2 @@ -29,17 +30,27 @@ jobs: run: | python setup.py install - - name: install openmc - run: | - conda install -c conda-forge openmc - - name: install packages for tests run: | pip install -r requirements-test.txt - name: Run test_utils run: | + cp tests/example.stl . pytest tests/ -v --cov=regular_mesh_plotter --cov-append --cov-report term --cov-report xml - name: Upload to codecov uses: codecov/codecov-action@v2 + + - name: Run examples + run: | + cd examples + python create_statepoint_file_with_meshes_openmc_dagmc.py + python create_statepoint_file_with_meshes_openmc.py + # python examples_from_readme.py # commented out as it requires entire neutronics workflow + python plot_regular_mesh_dose_tally.py + python plot_regular_mesh_dose_tally_with_geometry.py + python plot_regular_mesh_tally.py + python plot_regular_mesh_tally_with_geometry.py + python plot_regular_mesh_values.py + python plot_regular_mesh_values_with_geometry.py diff --git a/README.md b/README.md index 59d7134..13c4428 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ Mesh results in the form of Numpy arrays or 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 accessable +to the legend, axis, colour map etc. However some key options are accessible in the function call directly. There are additional options that allow @@ -38,26 +38,129 @@ and other mesh tally results. Example 1 shows a Numpy array plotted ```python -TODO +plot_regular_mesh_values( + values: np.ndarray, + filename: Optional[str] = None, + scale=None, # LogNorm(), + vmin=None, + label="", + base_plt=None, + extent=None, + x_label="X [cm]", + y_label="Y [cm]", + rotate_plot: float = 0, +) ``` Example 2 shows a Numpy array plotted with an underlying DAGMC geometry ```python -TODO +plot_regular_mesh_values_with_geometry( + values: np.ndarray, + dagmc_file_or_trimesh_object, + filename: Optional[str] = None, + scale=None, # LogNorm(), + vmin=None, + label="", + 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, +) ``` -Example 3 shows a OpenMC tally plotted with an underlying DAGMC geometry +Example 3 shows a OpenMC tally plotted ```python -TODO +plot_regular_mesh_tally( + 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, +) ``` -Example 4 shows how to rotate the underlying DAGMC geometry and mesh tally. +Example 4 shows a OpenMC tally plotted with an underlying DAGMC geometry +```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, +) +``` + +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 +when changing the plane normal ```python -TODO +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, +) ``` +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 cm **2 / simulated_particle", + source_strength: float = None, + std_dev_or_tally_value: str = "tally_value", +): +``` + +Additional examples can be found in the [examples folder in the GitHub repository](https://github.com/fusion-energy/regular_mesh_plotter/tree/main/examples) + # Related packages If you want to plot the DAGMC geometry without a mesh tally then take a look at diff --git a/examples/create_statepoint_file_with_meshes_openmc.py b/examples/create_statepoint_file_with_meshes_openmc.py index 4fdcd8d..4f4eaf3 100644 --- a/examples/create_statepoint_file_with_meshes_openmc.py +++ b/examples/create_statepoint_file_with_meshes_openmc.py @@ -81,7 +81,7 @@ ) settings = odw.FusionSettings() -settings.batches = 2 +settings.batches = 3 settings.particles = 1000000 # assigns a ring source of DT energy neutrons to the source using the # openmc_plasma_source package diff --git a/examples/create_statepoint_file_with_meshes_openmc_dagmc.py b/examples/create_statepoint_file_with_meshes_openmc_dagmc.py index 7d4dd3e..228ee7a 100644 --- a/examples/create_statepoint_file_with_meshes_openmc_dagmc.py +++ b/examples/create_statepoint_file_with_meshes_openmc_dagmc.py @@ -4,15 +4,15 @@ import openmc import openmc_dagmc_wrapper as odw import openmc_plasma_source as ops -import paramak from stl_to_h5m import stl_to_h5m -my_shape = paramak.ExtrudeStraightShape( - points=[(1, 1), (1, 200), (600, 200), (600, 1)], - distance=180, -) - -my_shape.export_stl("example.stl") +# 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, @@ -33,13 +33,13 @@ ) materials = odw.Materials( - h5m_filename="dagmc.h5m", correspondence_dict={"mat1": "eurofer"} + h5m_filename="dagmc.h5m", correspondence_dict={"mat1": "FLiNaK"} ) tally1 = odw.MeshTally2D( tally_type="neutron_effective_dose", plane="xy", - mesh_resolution=(10, ), + mesh_resolution=(10, 5), bounding_box="dagmc.h5m", ) tally2 = odw.MeshTally2D( @@ -57,7 +57,7 @@ tally4 = odw.MeshTally2D( tally_type="neutron_effective_dose", plane="xy", - mesh_resolution=(10, ), + mesh_resolution=(10, 10), bounding_box="dagmc.h5m", ) tally5 = odw.MeshTally2D( diff --git a/examples/example.stl b/examples/example.stl new file mode 100644 index 0000000..ddc4f41 --- /dev/null +++ b/examples/example.stl @@ -0,0 +1,86 @@ +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/plot_regular_mesh_dose_tally.py b/examples/plot_regular_mesh_dose_tally.py index b8a1625..87a0cae 100644 --- a/examples/plot_regular_mesh_dose_tally.py +++ b/examples/plot_regular_mesh_dose_tally.py @@ -10,15 +10,15 @@ # creates a plot of the mesh tally my_plot = 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 + filename="plot_regular_mesh_dose_tally.png", # the filename of the picture file saved scale=None, # LogNorm(), vmin=None, label="", x_label="X [cm]", y_label="Y [cm]", - rotate_plot = 0, - required_units='picosievert cm **2 / simulated_particle', - source_strength = None, + rotate_plot=0, + required_units="picosievert cm **2 / simulated_particle", + source_strength=None, ) # displays the plot diff --git a/examples/plot_regular_mesh_dose_tally_with_geometry.py b/examples/plot_regular_mesh_dose_tally_with_geometry.py index 6651661..21c07e7 100644 --- a/examples/plot_regular_mesh_dose_tally_with_geometry.py +++ b/examples/plot_regular_mesh_dose_tally_with_geometry.py @@ -11,18 +11,18 @@ my_plot = 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', + filename="plot_regular_mesh_dose_tally_with_geometry.png", scale=None, # LogNorm(), vmin=None, label="", 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 cm **2 / simulated_particle', - source_strength = None, + 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 cm **2 / simulated_particle", + source_strength=None, ) my_plot.show() diff --git a/examples/plot_regular_mesh_tally.py b/examples/plot_regular_mesh_tally.py index 7e5b5b1..56b540f 100644 --- a/examples/plot_regular_mesh_tally.py +++ b/examples/plot_regular_mesh_tally.py @@ -1,4 +1,3 @@ - import regular_mesh_plotter as rmp import openmc @@ -11,7 +10,7 @@ # creates a plot of the mesh rmp.plot_regular_mesh_tally( tally=my_tally, - filename= 'neutron_effective_dose_on_2D_mesh_xy.png', + filename="neutron_effective_dose_on_2D_mesh_xy.png", scale=None, # LogNorm(), vmin=None, label="", @@ -21,4 +20,4 @@ # rotate_plot: float = 0, # required_units: str = None, # source_strength: float = None, -) \ No newline at end of file +) diff --git a/examples/plot_regular_mesh_tally_with_geometry.py b/examples/plot_regular_mesh_tally_with_geometry.py index 484ff30..787c58c 100644 --- a/examples/plot_regular_mesh_tally_with_geometry.py +++ b/examples/plot_regular_mesh_tally_with_geometry.py @@ -1,4 +1,3 @@ - import regular_mesh_plotter as rmp import openmc @@ -11,20 +10,20 @@ # creates a plot of the mesh my_plot = 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', + dagmc_file_or_trimesh_object="dagmc.h5m", + std_dev_or_tally_value="tally_value", + filename="plot_regular_mesh_tally_with_geometry.png", scale=None, # LogNorm(), vmin=None, 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, + # plane_origin=[0, 0, 0], + plane_normal=[0, 0, 1], + rotate_mesh=0, + rotate_geometry=0, required_units=None, - source_strength = None + source_strength=None, ) my_plot.show() diff --git a/examples/plot_regular_mesh_values.py b/examples/plot_regular_mesh_values.py index 95b4cb7..e7b76a3 100644 --- a/examples/plot_regular_mesh_values.py +++ b/examples/plot_regular_mesh_values.py @@ -1,23 +1,78 @@ - import regular_mesh_plotter as rmp -import openmc - -# loads in the statepoint file containing tallies -statepoint = openmc.StatePoint(filepath="statepoint.2.h5") +import numpy as np -# gets one tally from the available tallies -my_tally = statepoint.get_tally(name="2_neutron_effective_dose") +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 mesh -plot_regular_mesh_values( - values: np.ndarray, - filename: Optional[str] = None, +# creates a plot of the array +rmp.plot_regular_mesh_values( + values=values, + filename="plot_regular_mesh_values.png", scale=None, # LogNorm(), vmin=None, - label="", - base_plt=None, - extent=None, + label="legend label", x_label="X [cm]", y_label="Y [cm]", - rotate_plot: float = 0, -) \ No newline at end of file +) diff --git a/examples/plot_regular_mesh_values_with_geometry.py b/examples/plot_regular_mesh_values_with_geometry.py index 868d4d1..0053ed3 100644 --- a/examples/plot_regular_mesh_values_with_geometry.py +++ b/examples/plot_regular_mesh_values_with_geometry.py @@ -1,27 +1,81 @@ - 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="2_neutron_effective_dose") +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, + ], + ] +) -def plot_regular_mesh_values_with_geometry( - values: np.ndarray, - dagmc_file_or_trimesh_object, - filename: Optional[str] = None, +rmp.plot_regular_mesh_values_with_geometry( + values=values, + dagmc_file_or_trimesh_object="example.stl", + filename="plot_regular_mesh_values_with_geometry.png", scale=None, # LogNorm(), vmin=None, - label="", - extent=None, + label="legend 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, -): \ No newline at end of file + plane_normal=[1, 0, 0], + extent=[0, 700, 0, 300], +) diff --git a/regular_mesh_plotter/__init__.py b/regular_mesh_plotter/__init__.py index 7b797fe..f0cd08b 100644 --- a/regular_mesh_plotter/__init__.py +++ b/regular_mesh_plotter/__init__.py @@ -1,8 +1,13 @@ - 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 .core import get_tally_extent + +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 diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index b63537c..9585052 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -1,14 +1,17 @@ -import trimesh +from typing import List, Optional + +import dagmc_geometry_slice_plotter as dgsp import matplotlib.pyplot as plt -from matplotlib import transforms -from typing import List, Optional, Tuple import numpy as np -from pathlib import Path -import dagmc_geometry_slice_plotter as dgsp +import openmc_post_processor as opp from matplotlib import transforms -import openmc_post_processor as opp -import openmc +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( @@ -91,45 +94,11 @@ def plot_regular_mesh_values_with_geometry( return both -def get_values_from_tally(tally): - data_frame = tally.get_pandas_dataframe() - if "std. dev." in data_frame.columns.to_list(): - values = ( - np.array(data_frame["mean"]), - np.array(data_frame["std. dev."]) - ) - else: - values = ( - np.array(data_frame["mean"]) - ) - return values - -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): - 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 def plot_regular_mesh_tally_with_geometry( tally, dagmc_file_or_trimesh_object, - std_dev_or_tally_value='tally_value', + std_dev_or_tally_value="tally_value", filename: Optional[str] = None, scale=None, # LogNorm(), vmin=None, @@ -141,18 +110,16 @@ def plot_regular_mesh_tally_with_geometry( rotate_mesh: float = 0, rotate_geometry: float = 0, required_units=None, - source_strength: float = None + source_strength: float = None, ): if required_units is not None: values = opp.process_tally( - tally=tally, - required_units=required_units, - source_strength=source_strength + 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) @@ -180,7 +147,6 @@ def plot_regular_mesh_tally_with_geometry( return plot - def plot_regular_mesh_tally( tally, filename: Optional[str] = None, @@ -193,18 +159,16 @@ def plot_regular_mesh_tally( rotate_plot: float = 0, required_units: str = None, source_strength: float = None, - std_dev_or_tally_value='tally_value', + std_dev_or_tally_value="tally_value", ): if required_units is not None: values = opp.process_tally( - tally=tally, - required_units=required_units, - source_strength=source_strength + 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) value = reshape_values_to_mesh_shape(tally, value) @@ -237,16 +201,14 @@ def plot_regular_mesh_dose_tally( x_label="X [cm]", y_label="Y [cm]", rotate_plot: float = 0, - required_units='picosievert cm **2 / simulated_particle', + required_units="picosievert cm **2 / simulated_particle", source_strength: float = None, - std_dev_or_tally_value: str = 'tally_value', + std_dev_or_tally_value: str = "tally_value", ): if required_units is not None: values = opp.process_dose_tally( - tally=tally, - required_units=required_units, - source_strength=source_strength + tally=tally, required_units=required_units, source_strength=source_strength ) else: values = get_values_from_tally(tally) @@ -286,9 +248,9 @@ def plot_regular_mesh_dose_tally_with_geometry( plane_normal: List[float] = [0, 0, 1], rotate_mesh: float = 0, rotate_geometry: float = 0, - required_units='picosievert cm **2 / simulated_particle', + required_units="picosievert cm **2 / simulated_particle", source_strength: float = None, - std_dev_or_tally_value: str = 'tally_value', + std_dev_or_tally_value: str = "tally_value", ): slice = dgsp.plot_slice_of_dagmc_geometry( @@ -314,46 +276,3 @@ def plot_regular_mesh_dose_tally_with_geometry( ) return both - - -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 diff --git a/regular_mesh_plotter/utils.py b/regular_mesh_plotter/utils.py new file mode 100644 index 0000000..a1f364a --- /dev/null +++ b/regular_mesh_plotter/utils.py @@ -0,0 +1,86 @@ +from pathlib import Path +from typing import List, Optional, Tuple + +import matplotlib.pyplot as plt +import numpy as np +import openmc + + +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 data_frame.columns.to_list(): + values = (np.array(data_frame["mean"]), np.array(data_frame["std. dev."])) + else: + values = np.array(data_frame["mean"]) + return values + + +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 diff --git a/requirements-test.txt b/requirements-test.txt index a270299..aff697a 100644 --- a/requirements-test.txt +++ b/requirements-test.txt @@ -1,2 +1,4 @@ - +openmc-dagmc-wrapper pytest-cov>=2.12.1 +openmc_plasma_source +stl_to_h5m diff --git a/setup.py b/setup.py index cb86a3f..d386804 100644 --- a/setup.py +++ b/setup.py @@ -39,6 +39,6 @@ "shapely", "scipy", "dagmc_geometry_slice_plotter", - "openmc_post_processor" + "openmc_post_processor", ], ) diff --git a/tests/test_plot_regular_mesh_values.py b/tests/test_plot_regular_mesh_values.py index 834d98b..2a9b08d 100644 --- a/tests/test_plot_regular_mesh_values.py +++ b/tests/test_plot_regular_mesh_values.py @@ -83,8 +83,8 @@ def test_plot_regular_mesh_values(self): def test_plot_regular_mesh_values_with_output(self): - os.system('rm test.png') + os.system("rm test.png") - plot_regular_mesh_values(values=self.values, filename='test.png') + plot_regular_mesh_values(values=self.values, filename="test.png") - assert Path('test.png').is_file() + assert Path("test.png").is_file() diff --git a/tests/test_plot_regular_mesh_values_with_geometry.py b/tests/test_plot_regular_mesh_values_with_geometry.py index a8c934b..10af4d6 100644 --- a/tests/test_plot_regular_mesh_values_with_geometry.py +++ b/tests/test_plot_regular_mesh_values_with_geometry.py @@ -79,7 +79,7 @@ def test_plot_regular_mesh_values_with_geometry(self): test_plot = plot_regular_mesh_values_with_geometry( values=self.values, - dagmc_file_or_trimesh_object='example.stl' # this could be a h5m file + dagmc_file_or_trimesh_object="example.stl", # this could be a h5m file ) assert isinstance(test_plot, type(matplotlib.pyplot)) @@ -88,9 +88,9 @@ def test_plot_regular_mesh_values_with_geometry(self): test_plot = plot_regular_mesh_values_with_geometry( values=self.values, - dagmc_file_or_trimesh_object='example.stl', # this could be a h5m file - filename='test.png' + dagmc_file_or_trimesh_object="example.stl", # this could be a h5m file + filename="test.png", ) assert isinstance(test_plot, type(matplotlib.pyplot)) - assert Path('test.png').is_file() + assert Path("test.png").is_file()