From b0b40c352d67982ab6130471ff00c75c2ab1cf1c Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 27 Oct 2021 11:54:45 +0100 Subject: [PATCH 1/7] ground work for scaling units --- regular_mesh_plotter/core.py | 87 ++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index 2ff1f32..c52a15a 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -141,10 +141,47 @@ def plot_regular_mesh_tally( rotate_plot: float = 0, ): + values = opp.process_tally( + tally=tally, + ) + + values = reshape_values_to_mesh_shape(tally, values) + + extent = get_tally_extent(tally) + + plot = plot_regular_mesh_values( + values=values, + 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, +): + values = opp.process_dose_tally( tally=tally, ) + values = reshape_values_to_mesh_shape(tally, values) + extent = get_tally_extent(tally) plot = plot_regular_mesh_values( @@ -163,6 +200,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, +): + + 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, + ) + + 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) + reshaped_values = [] + for value in values: + reshaped_value = value.reshape(shape) + reshaped_values.append(reshaped_value) + return reshaped_values + + def get_tally_extent( tally, ): From f80b29dc11e3447450018a1e45ff5e4d78ecb927 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 27 Oct 2021 17:55:06 +0100 Subject: [PATCH 2/7] added unit propergation --- regular_mesh_plotter/__init__.py | 3 ++- regular_mesh_plotter/core.py | 5 +++++ 2 files changed, 7 insertions(+), 1 deletion(-) 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 c52a15a..4a622a3 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -105,6 +105,7 @@ def plot_regular_mesh_tally_with_geometry( plane_normal: List[float] = [0, 0, 1], rotate_mesh: float = 0, rotate_geometry: float = 0, + required_units='1 / simulated_particle' ): slice = dgsp.plot_slice_of_dagmc_geometry( @@ -124,6 +125,7 @@ def plot_regular_mesh_tally_with_geometry( x_label=x_label, y_label=y_label, rotate_plot=rotate_mesh, + required_units=required_units ) return both @@ -139,10 +141,12 @@ def plot_regular_mesh_tally( x_label="X [cm]", y_label="Y [cm]", rotate_plot: float = 0, + required_units='1 / simulated_particle' ): values = opp.process_tally( tally=tally, + required_units=required_units ) values = reshape_values_to_mesh_shape(tally, values) @@ -164,6 +168,7 @@ def plot_regular_mesh_tally( return plot + def plot_regular_mesh_dose_tally( tally, filename: Optional[str] = None, From 45e70aef010c24a86b092c7811c6dfcf1c4d09e4 Mon Sep 17 00:00:00 2001 From: Jonathan Date: Wed, 27 Oct 2021 23:52:41 +0100 Subject: [PATCH 3/7] added required_units and option for no unit conversion --- regular_mesh_plotter/core.py | 31 ++++++++++++++++++++++++------- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index 4a622a3..b5c0155 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -105,7 +105,7 @@ def plot_regular_mesh_tally_with_geometry( plane_normal: List[float] = [0, 0, 1], rotate_mesh: float = 0, rotate_geometry: float = 0, - required_units='1 / simulated_particle' + required_units=None ): slice = dgsp.plot_slice_of_dagmc_geometry( @@ -141,20 +141,29 @@ def plot_regular_mesh_tally( x_label="X [cm]", y_label="Y [cm]", rotate_plot: float = 0, - required_units='1 / simulated_particle' + required_units=None, + #todo add scaling option + #source_strength ): - values = opp.process_tally( - tally=tally, - required_units=required_units - ) + if required_units is not None: + values = opp.process_tally( + tally=tally, + required_units=required_units + ) + else: + data_frame = tally.get_pandas_dataframe() + values = ( + np.array(data_frame["mean"]), + np.array(data_frame["std. dev."]) + ) values = reshape_values_to_mesh_shape(tally, values) extent = get_tally_extent(tally) plot = plot_regular_mesh_values( - values=values, + values=values[0], filename=filename, scale=scale, vmin=vmin, @@ -179,10 +188,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' + #todo add scaling option + #source_strength ): values = opp.process_dose_tally( tally=tally, + required_units=required_units ) values = reshape_values_to_mesh_shape(tally, values) @@ -218,6 +231,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' + #todo add scaling option + #source_strength ): slice = dgsp.plot_slice_of_dagmc_geometry( @@ -237,6 +253,7 @@ def plot_regular_mesh_dose_tally_with_geometry( x_label=x_label, y_label=y_label, rotate_plot=rotate_mesh, + required_units=required_units ) return both From 5a58d7016fcf6391d61d7d8f50b675df11db6930 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 29 Oct 2021 11:43:47 +0100 Subject: [PATCH 4/7] started added second plot --- regular_mesh_plotter/core.py | 127 +++++++++++++++++++++++++---------- 1 file changed, 92 insertions(+), 35 deletions(-) diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index b5c0155..4ba1446 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,6 +91,18 @@ 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 plot_regular_mesh_tally_with_geometry( tally, @@ -105,30 +117,65 @@ 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 + 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) + + 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, + if isinstance(values, np.ndarray): + values = reshape_values_to_mesh_shape(tally, values) + else: + values = reshape_values_to_mesh_shape(tally, values[0]) + + plot = plot_regular_mesh_values( + values=values, 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, - required_units=required_units ) + + elif len(values) == 2: + value_std_dev = reshape_values_to_mesh_shape(tally, values[1]) + plot_std_dev = plot_regular_mesh_values( + values=value_std_dev, + 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_mesh, + ) + return [plot, plot_std_dev] - return both + + + return plot def plot_regular_mesh_tally( @@ -141,22 +188,18 @@ def plot_regular_mesh_tally( x_label="X [cm]", y_label="Y [cm]", rotate_plot: float = 0, - required_units=None, - #todo add scaling option - #source_strength + required_units: str = None, + source_strength: float = None, ): if required_units is not None: values = opp.process_tally( tally=tally, - required_units=required_units + required_units=required_units, + source_strength=source_strength ) else: - data_frame = tally.get_pandas_dataframe() - values = ( - np.array(data_frame["mean"]), - np.array(data_frame["std. dev."]) - ) + values = get_values_from_tally(tally) values = reshape_values_to_mesh_shape(tally, values) @@ -181,6 +224,7 @@ def plot_regular_mesh_tally( def plot_regular_mesh_dose_tally( tally, filename: Optional[str] = None, + filename_std_dev: Optional[str] = None, scale=None, # LogNorm(), vmin=None, label="", @@ -188,22 +232,26 @@ 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' - #todo add scaling option - #source_strength + required_units='picosievert cm **2 / simulated_particle', + source_strength: float = None, ): - values = opp.process_dose_tally( - tally=tally, - required_units=required_units - ) + 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, @@ -231,9 +279,8 @@ 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' - #todo add scaling option - #source_strength + required_units='picosievert cm **2 / simulated_particle', + source_strength: float = None, ): slice = dgsp.plot_slice_of_dagmc_geometry( @@ -253,7 +300,8 @@ def plot_regular_mesh_dose_tally_with_geometry( x_label=x_label, y_label=y_label, rotate_plot=rotate_mesh, - required_units=required_units + required_units=required_units, + source_strength=source_strength ) return both @@ -265,11 +313,20 @@ def reshape_values_to_mesh_shape(tally, values): # 2d mesh has a shape in the form [1, 400, 400] if 1 in shape: shape.remove(1) - reshaped_values = [] - for value in values: - reshaped_value = value.reshape(shape) - reshaped_values.append(reshaped_value) - return reshaped_values + return values.reshape(shape) + + # if isinstance(values, tuple): + # # values is a tupe of tally results and std dev + # reshaped_values = [] + # for value in values: + # reshaped_value = value.reshape(shape) + # reshaped_values.append(reshaped_value) + # return reshaped_values + # else: + # # values contains just the tally results + # # simulations with a sigle batch have just a tally result and no std dev + # return values.reshape(shape) + def get_tally_extent( From e0217609e977f7b0b75fbfc2b4b11b3dbe1691a9 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 29 Oct 2021 16:49:33 +0100 Subject: [PATCH 5/7] adding in option for std dev to one function --- examples/plot_mesh_with_geometry_slice.py | 17 ------- regular_mesh_plotter/core.py | 61 +++++++++-------------- 2 files changed, 23 insertions(+), 55 deletions(-) delete mode 100644 examples/plot_mesh_with_geometry_slice.py 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/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index 4ba1446..287456e 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -104,9 +104,28 @@ def get_values_from_tally(tally): ) 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, @@ -129,6 +148,8 @@ def plot_regular_mesh_tally_with_geometry( ) 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) @@ -139,13 +160,8 @@ def plot_regular_mesh_tally_with_geometry( rotate_plot=rotate_geometry, ) - if isinstance(values, np.ndarray): - values = reshape_values_to_mesh_shape(tally, values) - else: - values = reshape_values_to_mesh_shape(tally, values[0]) - plot = plot_regular_mesh_values( - values=values, + values=value, filename=filename, scale=scale, vmin=vmin, @@ -156,28 +172,11 @@ def plot_regular_mesh_tally_with_geometry( y_label=y_label, rotate_plot=rotate_mesh, ) - - elif len(values) == 2: - value_std_dev = reshape_values_to_mesh_shape(tally, values[1]) - plot_std_dev = plot_regular_mesh_values( - values=value_std_dev, - 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_mesh, - ) - return [plot, plot_std_dev] - - return plot + def plot_regular_mesh_tally( tally, filename: Optional[str] = None, @@ -224,7 +223,6 @@ def plot_regular_mesh_tally( def plot_regular_mesh_dose_tally( tally, filename: Optional[str] = None, - filename_std_dev: Optional[str] = None, scale=None, # LogNorm(), vmin=None, label="", @@ -315,19 +313,6 @@ def reshape_values_to_mesh_shape(tally, values): shape.remove(1) return values.reshape(shape) - # if isinstance(values, tuple): - # # values is a tupe of tally results and std dev - # reshaped_values = [] - # for value in values: - # reshaped_value = value.reshape(shape) - # reshaped_values.append(reshaped_value) - # return reshaped_values - # else: - # # values contains just the tally results - # # simulations with a sigle batch have just a tally result and no std dev - # return values.reshape(shape) - - def get_tally_extent( tally, From 0085f4c517938a67ce624ec04b4f7e2db1ce3aee Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 29 Oct 2021 16:51:33 +0100 Subject: [PATCH 6/7] added draft example for each function --- ...eate_statepoint_file_with_meshes_openmc.py | 94 ++++++++++++++++ ...tatepoint_file_with_meshes_openmc_dagmc.py | 104 ++++++++++++++++++ examples/plot_regular_mesh_dose_tally.py | 25 +++++ ...t_regular_mesh_dose_tally_with_geometry.py | 28 +++++ examples/plot_regular_mesh_tally.py | 23 ++++ .../plot_regular_mesh_tally_with_geometry.py | 30 +++++ examples/plot_regular_mesh_values.py | 23 ++++ .../plot_regular_mesh_values_with_geometry.py | 27 +++++ 8 files changed, 354 insertions(+) create mode 100644 examples/create_statepoint_file_with_meshes_openmc.py create mode 100644 examples/create_statepoint_file_with_meshes_openmc_dagmc.py create mode 100644 examples/plot_regular_mesh_dose_tally.py create mode 100644 examples/plot_regular_mesh_dose_tally_with_geometry.py create mode 100644 examples/plot_regular_mesh_tally.py create mode 100644 examples/plot_regular_mesh_tally_with_geometry.py create mode 100644 examples/plot_regular_mesh_values.py create mode 100644 examples/plot_regular_mesh_values_with_geometry.py 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_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 From e5e0ed1e9c28f26e0cff675d99cbd2750a38a7e2 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Sat, 30 Oct 2021 06:41:59 +0100 Subject: [PATCH 7/7] added v0.1.0 citation --- CITATION.cff | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 CITATION.cff 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"