diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000..48da971 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,10 @@ +cff-version: 1.2.0 +message: "If you use this software, please cite it as below." +authors: +- family-names: "Shimwell" + given-names: "Jonathan" + orcid: "https://orcid.org/0000-0001-6909-0946" +title: "A Python package for plotting regular mesh tally results from neutronics simulations." +version: 0.1.0 +date-released: 2021-10-30 +url: "https://github.com/fusion-energy/regular_mesh_plotter" diff --git a/examples/create_statepoint_file_with_meshes_openmc.py b/examples/create_statepoint_file_with_meshes_openmc.py new file mode 100644 index 0000000..4fdcd8d --- /dev/null +++ b/examples/create_statepoint_file_with_meshes_openmc.py @@ -0,0 +1,94 @@ +# 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 +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 = 2 +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 new file mode 100644 index 0000000..7d4dd3e --- /dev/null +++ b/examples/create_statepoint_file_with_meshes_openmc_dagmc.py @@ -0,0 +1,104 @@ +# 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 +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") + +# 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", +) + +# 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": "eurofer"} +) + +tally1 = odw.MeshTally2D( + tally_type="neutron_effective_dose", + plane="xy", + mesh_resolution=(10, ), + bounding_box="dagmc.h5m", +) +tally2 = odw.MeshTally2D( + tally_type="neutron_effective_dose", + plane="yz", + mesh_resolution=(10, 5), + bounding_box="dagmc.h5m", +) +tally3 = odw.MeshTally2D( + tally_type="neutron_effective_dose", + plane="xz", + mesh_resolution=(10, 5), + bounding_box="dagmc.h5m", +) +tally4 = odw.MeshTally2D( + tally_type="neutron_effective_dose", + plane="xy", + mesh_resolution=(10, ), + bounding_box="dagmc.h5m", +) +tally5 = odw.MeshTally2D( + tally_type="neutron_effective_dose", + plane="yz", + mesh_resolution=(10, 5), + bounding_box="dagmc.h5m", +) +tally6 = odw.MeshTally2D( + tally_type="neutron_effective_dose", + plane="xz", + mesh_resolution=(10, 5), + bounding_box="dagmc.h5m", +) + +# tally2 = odw.MeshTally3D( +# mesh_resolution=(100, 100, 100), +# bounding_box="dagmc.h5m", +# 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/plot_mesh_with_geometry_slice.py b/examples/plot_mesh_with_geometry_slice.py deleted file mode 100644 index fb2622c..0000000 --- a/examples/plot_mesh_with_geometry_slice.py +++ /dev/null @@ -1,17 +0,0 @@ - -from regular_mesh_plotter import plot_regular_mesh_tally_with_geometry -import openmc - -# loads in the statepoint file containing tallies -statepoint = openmc.StatePoint(filepath="statepoint.3.h5") - -# gets one tally from the available tallies -my_tally = statepoint.get_tally(name="neutron_effective_dose_on_2D_mesh_xy") - -# creates the matplotlib mesh plot with geometry -plot = plot_regular_mesh_tally_with_geometry( - tally=my_tally, - mesh_file_or_trimesh_object="dagmc.h5m", -) - -plot.show() diff --git a/examples/plot_regular_mesh_dose_tally.py b/examples/plot_regular_mesh_dose_tally.py new file mode 100644 index 0000000..44d941d --- /dev/null +++ b/examples/plot_regular_mesh_dose_tally.py @@ -0,0 +1,25 @@ +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") + +# 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 + 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, +) + +# displays the plot +my_plot.show() diff --git a/examples/plot_regular_mesh_dose_tally_with_geometry.py b/examples/plot_regular_mesh_dose_tally_with_geometry.py new file mode 100644 index 0000000..2f683ec --- /dev/null +++ b/examples/plot_regular_mesh_dose_tally_with_geometry.py @@ -0,0 +1,28 @@ +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") + +# creates a plot of the mesh +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', + 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, +): + +my_plot.show() diff --git a/examples/plot_regular_mesh_tally.py b/examples/plot_regular_mesh_tally.py new file mode 100644 index 0000000..83f966f --- /dev/null +++ b/examples/plot_regular_mesh_tally.py @@ -0,0 +1,23 @@ + +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") + +# creates a plot of the mesh +plot_regular_mesh_tally( + tally, + filename: Optional[str] = None, + scale=None, # LogNorm(), + vmin=None, + label="", + base_plt=None, + x_label="X [cm]", + y_label="Y [cm]", + 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 new file mode 100644 index 0000000..9afd170 --- /dev/null +++ b/examples/plot_regular_mesh_tally_with_geometry.py @@ -0,0 +1,30 @@ + +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") + +# 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', + 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, + required_units=None, + source_strength = None +) + +my_plot.show() diff --git a/examples/plot_regular_mesh_values.py b/examples/plot_regular_mesh_values.py new file mode 100644 index 0000000..95b4cb7 --- /dev/null +++ b/examples/plot_regular_mesh_values.py @@ -0,0 +1,23 @@ + +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") + +# creates a plot of the mesh +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, +) \ 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 new file mode 100644 index 0000000..868d4d1 --- /dev/null +++ b/examples/plot_regular_mesh_values_with_geometry.py @@ -0,0 +1,27 @@ + +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") + + + +def 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, +): \ No newline at end of file diff --git a/regular_mesh_plotter/__init__.py b/regular_mesh_plotter/__init__.py index 49becfc..7b797fe 100644 --- a/regular_mesh_plotter/__init__.py +++ b/regular_mesh_plotter/__init__.py @@ -3,5 +3,6 @@ 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 - diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index 2ff1f32..287456e 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -1,7 +1,7 @@ import trimesh import matplotlib.pyplot as plt from matplotlib import transforms -from typing import List, Optional +from typing import List, Optional, Tuple import numpy as np from pathlib import Path import dagmc_geometry_slice_plotter as dgsp @@ -91,10 +91,41 @@ 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 = reshape_values_to_mesh_shape(tally, values[1]) + elif std_dev_or_tally_value == 'tally_value': + if isinstance(values, np.ndarray): + value = reshape_values_to_mesh_shape(tally, values) + elif isinstance(values, tuple): + value = reshape_values_to_mesh_shape(tally, values[0]) + else: + msg = f'Values to plot should be a numpy ndarry or a tuple or numpy ndarrys not a {type(values)}' + raise ValueError(msg) + else: + msg = f'Value of std_dev_or_tally_value should be either "std_dev" or "value", not {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', filename: Optional[str] = None, scale=None, # LogNorm(), vmin=None, @@ -105,28 +136,45 @@ def plot_regular_mesh_tally_with_geometry( plane_normal: List[float] = [0, 0, 1], rotate_mesh: float = 0, rotate_geometry: float = 0, + required_units=None, + source_strength: float = None ): - slice = dgsp.plot_slice_of_dagmc_geometry( + if required_units is not None: + values = opp.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) + + base_plt = 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, ) - both = plot_regular_mesh_tally( - tally=tally, + plot = plot_regular_mesh_values( + values=value, filename=filename, - scale=scale, # LogNorm(), + scale=scale, vmin=vmin, label=label, - base_plt=slice, + base_plt=base_plt, + extent=extent, x_label=x_label, y_label=y_label, rotate_plot=rotate_mesh, ) - return both + return plot + def plot_regular_mesh_tally( @@ -139,16 +187,69 @@ def plot_regular_mesh_tally( x_label="X [cm]", y_label="Y [cm]", rotate_plot: float = 0, + required_units: str = None, + source_strength: float = None, ): - values = opp.process_dose_tally( - tally=tally, + if required_units is not None: + values = opp.process_tally( + tally=tally, + required_units=required_units, + source_strength=source_strength + ) + else: + values = get_values_from_tally(tally) + + values = reshape_values_to_mesh_shape(tally, values) + + extent = get_tally_extent(tally) + + plot = plot_regular_mesh_values( + values=values[0], + filename=filename, + scale=scale, + vmin=vmin, + label=label, + base_plt=base_plt, + extent=extent, + x_label=x_label, + y_label=y_label, + rotate_plot=rotate_plot, ) + return plot + + +def plot_regular_mesh_dose_tally( + tally, + filename: Optional[str] = None, + scale=None, # LogNorm(), + vmin=None, + label="", + base_plt=None, + x_label="X [cm]", + y_label="Y [cm]", + rotate_plot: float = 0, + required_units='picosievert cm **2 / simulated_particle', + source_strength: float = None, +): + + if required_units is not None: + values = opp.process_dose_tally( + tally=tally, + required_units=required_units, + source_strength=source_strength + ) + else: + values = get_values_from_tally(tally) + + values = reshape_values_to_mesh_shape(tally, values) + extent = get_tally_extent(tally) + if len plot = plot_regular_mesh_values( - values=values, + values=values[0], filename=filename, scale=scale, vmin=vmin, @@ -163,6 +264,56 @@ def plot_regular_mesh_tally( return plot +def 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, +): + + 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, + ) + + both = plot_regular_mesh_dose_tally( + tally=tally, + filename=filename, + scale=scale, # LogNorm(), + vmin=vmin, + label=label, + base_plt=slice, + x_label=x_label, + y_label=y_label, + rotate_plot=rotate_mesh, + required_units=required_units, + source_strength=source_strength + ) + + 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, ):