diff --git a/.gitignore b/.gitignore index ccdcc64..92d669b 100644 --- a/.gitignore +++ b/.gitignore @@ -136,3 +136,6 @@ _version.py # output xml files from openmc *.xml + +*.html +*.png diff --git a/examples/plot_multiple_materials.py b/examples/plot_multiple_materials.py index b5c17b3..de2248b 100644 --- a/examples/plot_multiple_materials.py +++ b/examples/plot_multiple_materials.py @@ -1,4 +1,4 @@ -from openmc_depletion_plotter import plot_materials +from openmc_depletion_plotter import plot_material import openmc my_mat = openmc.Material() @@ -22,7 +22,7 @@ my_mat_2.set_density('g/cm3', 7.7) my_mat_2.volume = 1 -plot_materials( - [my_mat, my_mat_2], - filenames=['my_mat_1.png', 'my_mat_2.png'] -) +for mat, filename in zip([my_mat, my_mat_2], ['my_mat_1.png', 'my_mat_2.png']): + plotly_figure=plot_material(mat) + plotly_figure.write_image(filename) + diff --git a/examples/plot_single_material.py b/examples/plot_single_material.py index a2b2907..425104e 100644 --- a/examples/plot_single_material.py +++ b/examples/plot_single_material.py @@ -8,7 +8,15 @@ my_mat.add_element('Al', 0.5) my_mat.add_element('B', 0.5) my_mat.add_element('Co', 0.5) +my_mat.add_element('Cs', 0.5) +my_mat.add_nuclide('Fe60', 0.1) +my_mat.add_nuclide('Fe61',0.2) +my_mat.add_nuclide('Fe62',0.3) my_mat.set_density('g/cm3', 7.7) my_mat.volume = 1 -plot_material(my_mat, filename='my_mat_on_isotope_chart.png') +plotly_figure = plot_material(my_mat) + +plotly_figure.write_html('my_mat_on_isotope_chart.html') +plotly_figure.write_image('my_mat_on_isotope_chart.png') +plotly_figure.show() \ No newline at end of file diff --git a/openmc_depletion_plotter/__init__.py b/openmc_depletion_plotter/__init__.py index b44267b..a3ebdd8 100644 --- a/openmc_depletion_plotter/__init__.py +++ b/openmc_depletion_plotter/__init__.py @@ -13,6 +13,11 @@ __all__ = ["__version__"] -from .core import plot_material, plot_materials +from .core import plot_material from .utils import find_most_abundant_nuclides_in_material from .utils import find_most_abundant_nuclides_in_materials +from .utils import get_nuclide_atom_densities_from_materials +from .utils import find_most_active_nuclides_in_material +from .utils import find_most_active_nuclides_in_materials +from .utils import get_nuclide_activities_from_materials +from .utils import get_atoms_from_material diff --git a/openmc_depletion_plotter/core.py b/openmc_depletion_plotter/core.py index 7c30590..acee255 100644 --- a/openmc_depletion_plotter/core.py +++ b/openmc_depletion_plotter/core.py @@ -1,275 +1,111 @@ -import matplotlib.pyplot as plt -import numpy as np + import openmc -import pint -import seaborn as sns -# from matplotlib.colors import Normalize -from openmc.data import ATOMIC_SYMBOL, NATURAL_ABUNDANCE import matplotlib.cm as cm -ureg = pint.UnitRegistry() - -stable_nuclides = list(NATURAL_ABUNDANCE.keys()) - - -def get_atoms_from_material(material): - - if material.volume is None: - msg = "material.volume must be set to find the number of atoms present." - raise ValueError(msg) - - # in units of atom / ( barn cm2 ) - atoms_per_barn_cm2 = material.get_nuclide_atom_densities() - volume = material.volume * ureg.cm**3 - - # print(atoms_per_barn_cm2) - isotopes_and_atoms = [] - for key, value in atoms_per_barn_cm2.items(): - # print(key, value[1]) - atoms_per_b_cm = value[1] * ureg.particle / (ureg.barn * ureg.cm) - atoms = atoms_per_b_cm * volume - # print(key, atoms) - # print(key, atoms.to(ureg.particle)) - - atomic_number, mass_number, _ = openmc.data.zam(key) - - isotopes_and_atoms.append( - ( - atomic_number, - mass_number - atomic_number, - atoms.to(ureg.particle).magnitude, - key, - ) - ) - - # print(atoms_per_barn_cm2) - return isotopes_and_atoms - - -def build_grid_of_nuclides(iterable_of_nuclides): - - grid_width = 200 # protons - grid_height = 200 # neutrons - grid = np.array([[0] * grid_width] * grid_height, dtype=float) - - for atomic_number, neutron_number, atoms, _ in iterable_of_nuclides: - # print(atomic_number,neutron_number,atoms ) - grid[atomic_number][neutron_number] = atoms - - return grid - - -def build_grid_of_stable_nuclides(iterable_of_nuclides): - - grid_width = 200 # protons - grid_height = 200 # neutrons - grid = np.array([[0] * grid_width] * grid_height, dtype=float) - - for atomic_number, neutron_number in iterable_of_nuclides: - # print(atomic_number,neutron_number,atoms ) - grid[atomic_number][ - neutron_number - ] = 0.2 # sets the grey scale from 0 (white) to 1 (black) - - return grid - - -def build_grid_of_annotations(iterable_of_nuclides): - - grid_width = 200 # protons - grid_height = 200 # neutrons - grid = np.array([[0] * grid_width] * grid_height, dtype=str) - - for atomic_number, neutron_number, atoms, name in iterable_of_nuclides: - grid[atomic_number][neutron_number] = name - - return grid - - -def get_neutron_range(material): - nucs = material.get_nuclides() - - neutrons = [] - for nuc in nucs: - proton, proton_plus_neutron, _ = openmc.data.zam(nuc) - neutron = proton_plus_neutron - proton - neutrons.append(neutron) - return [min(neutrons), max(neutrons)] - - -def get_proton_range(material): - nucs = material.get_nuclides() - - protons = [] - for nuc in nucs: - proton, proton_plus_neutron, _ = openmc.data.zam(nuc) - protons.append(proton) - return [min(protons), max(protons)] - - -# Get the colormap and set the under and bad colors -# colMap = cm.gist_rainbow -def make_unstable_cm(): - colMap = cm.get_cmap("viridis", 256) - colMap.set_bad(color="white", alpha=100) - colMap.set_under(color="white", alpha=100) - return colMap - - -def make_stable_cm(): - colMap = cm.get_cmap("gray_r", 256) - colMap.set_bad(color="white", alpha=100) - colMap.set_under(color="white", alpha=100) - return colMap +import plotly.graph_objects as go +import matplotlib.cm as cm +from .utils import get_atoms_from_material +from openmc.data import NATURAL_ABUNDANCE -def plot_materials( - my_mats, - filenames, - neutron_range=None, - proton_range=None, -): +stable_nuclides = list(NATURAL_ABUNDANCE.keys()) - for filename, my_mat in zip(filenames, my_mats): - plot_material( - my_mat, - filename=filename, - neutron_range=neutron_range, - proton_range=proton_range, - ) def plot_material( my_mat, - filename=None, - neutron_range=None, - proton_range=None, - isotopes_label_size=None, + show_all=True + # isotopes_label_size=None, ): - plt.figure().clear() - - stable_nuclides_za = [] - for entry in stable_nuclides: - atomic_number, mass_number, _ = openmc.data.zam(entry) - stable_nuclides_za.append((atomic_number, mass_number - atomic_number)) + xycl = get_atoms_from_material(my_mat) - nuclides = get_atoms_from_material(my_mat) - grid = build_grid_of_nuclides(nuclides) - annots = build_grid_of_annotations(nuclides) - stable_grid = build_grid_of_stable_nuclides(stable_nuclides_za) - - sns.set_theme() - # sns.set_style("whitegrid") - sns.set_style("darkgrid", {"axes.facecolor": ".9"}) - # sns.set_style("white") - - # masked_array = np.ma.array (grid, mask=np.zeros(grid)) - data_masked = np.ma.masked_where(grid != 0, grid) - # plt.imshow(data_masked, interpolation = 'none', vmin = 0) - - ax2 = sns.heatmap( - stable_grid, - square=True, - cmap=make_stable_cm(), - cbar=False, - linewidths=0, - ) - - ax = sns.heatmap( - grid, - linewidths=0.1, - vmin=2.204575507634703e20, - vmax=7.173000749202642e22, - square=True, - cmap=make_unstable_cm(), - cbar_kws={"label": "number of atoms"}, - # annot=True, - # interpolation = 'none', - # mask=data_masked + fig = go.Figure() + y_vals, x_vals, c_vals, l_vals = zip(*xycl) + fig.add_trace( + go.Scatter( + x=x_vals, + y=y_vals, + mode='markers', + name='material', + ) ) - plt.gca().set_facecolor("white") - plt.gca().invert_yaxis() - - if neutron_range is None: - neutron_range = get_neutron_range(my_mat) - neutron_range[0] = max(neutron_range[0] - 1, 0) - neutron_range[1] = ( - neutron_range[1] + 1 + 1 - ) # todo remove this extra +1 which is currently needed - if proton_range is None: - proton_range = get_proton_range(my_mat) - proton_range[0] = max(proton_range[0] - 1, 0) - proton_range[1] = ( - proton_range[1] + 1 + 1 - ) # todo remove this extra +1 which is currently needed - - ax.set_xlim(neutron_range[0], neutron_range[1]) - ax.set_ylim(proton_range[0], proton_range[1]) - ax2.set_xlim(neutron_range[0], neutron_range[1]) - ax2.set_ylim(proton_range[0], proton_range[1]) - - # ax.set_title("Number of atoms in material") - ax.set_ylabel("Number of protons") - ax.set_xlabel("Number of neutrons") - ax.grid(True) - # ax.grid(True, which='both') - # ax.axhline(y=0, color='k') - # ax.axvline(x=0, color='k') - # plt.axvline(x=0, color='k') - - # plt.axis('on') - # ax.axis('on') - # ax2.axis('on') + #print(stable_nuclides) + #print(l_vals) + + + for stable in stable_nuclides: + atomic_number, mass_number, _ = openmc.data.zam(stable) + y = atomic_number + x = mass_number - y + fig.add_shape( + x0=x-0.5, + x1=x+0.5, + y0=y+0.5, + y1=y-0.5, + xref='x', + yref='y', + fillcolor='lightgrey', + line_color="lightgrey", + ) - plt.xticks(rotation=0) + for entry in xycl: + y, x, c, l = entry + + color = cm.viridis(c/max(c_vals))[:3] + scaled_color = (color[0]*255,color[1]*255,color[2]*255) + text_color = f'rgb{scaled_color}' + #print(color, text_color) + + if l in stable_nuclides: + line_color = "lightgrey" + line_width = 1 + else: + line_color = "Black" + line_width = 1 + + + fig.add_shape( + x0=x-0.5, + x1=x+0.5, + y0=y+0.5, + y1=y-0.5, + xref='x', + yref='y', + fillcolor=text_color, + # line_color="LightSeaGreen", + line={ + 'color':line_color, + 'width':line_width, + # 'dash':"dashdot", + } + ) - if isotopes_label_size is not None: + # fig.update_xaxes(rangemode="tozero") - for j in range(grid.shape[1]): - for i in range(grid.shape[0]): - # print(i,j, grid[i, j]) - # text = ax.text(j, i, isotope_chart[i, j], + if show_all: + fig.update(layout_yaxis_range = [0, 93]) + fig.update(layout_xaxis_range = [0, 147]) - if grid[i, j] > 0: + height = 93 + width = 147 + ratio = height / width - text = ax.text( - j + 0.5, - i + 0.66, - f"{ATOMIC_SYMBOL[i]}{i+j}", - ha="center", - va="center", - color="w", - fontdict={"size": isotopes_label_size}, - ) - text = ax.text( - j + 0.5, - i + 0.33, - f"{grid[i, j]:.1e}", - ha="center", - va="center", - color="w", - fontdict={"size": isotopes_label_size}, - ) + else: + fig.update(layout_yaxis_range = [0, max(y_vals)+1]) + fig.update(layout_xaxis_range = [0, max(x_vals)+1]) - if (i, j) in stable_nuclides_za: - text = ax.text( - j + 0.5, - i + 0.66, - f"{ATOMIC_SYMBOL[i]}{i+j}", - ha="center", - va="center", - color="w", - fontdict={"size": isotopes_label_size}, - ) + height = max(y_vals)+1 + width = max(x_vals)+1 + ratio = height / width - plt.axis("on") - plt.tight_layout() - ax.set_axis_on() - ax2.set_axis_on() - if filename: - plt.savefig(filename, dpi=800) - return plt - # plt.show() + fig.update_xaxes(title='protons') + fig.update_yaxes(title='neutrons') + fig.update_layout( + # autosize=True + width=1000, + height=1000*ratio + ) + return fig diff --git a/openmc_depletion_plotter/utils.py b/openmc_depletion_plotter/utils.py index 8ef750d..7b50680 100644 --- a/openmc_depletion_plotter/utils.py +++ b/openmc_depletion_plotter/utils.py @@ -1,6 +1,78 @@ from typing import Iterable import openmc +import numpy as np +import openmc +import pint + + +from openmc.data import ATOMIC_SYMBOL, NATURAL_ABUNDANCE + +ureg = pint.UnitRegistry() + +def find_most_active_nuclides_in_material( + material, + exclude=None, +): + + if exclude is None: + excluded_isotopes = [] + else: + if isinstance(exclude, Iterable): + excluded_isotopes = exclude + elif isinstance(exclude, openmc.Material): + excluded_isotopes = exclude.get_nuclides() + + non_excluded_nucs = {} + for key, value in material.get_nuclide_activity().items(): + if key not in excluded_isotopes: + if key not in non_excluded_nucs.keys(): + if value != 0.: + non_excluded_nucs[key] = value + else: + non_excluded_nucs[key] += value + + sorted_dict = { + k: v + for k, v in sorted( + non_excluded_nucs.items(), key=lambda item: item[1], reverse=True + ) + } + return list(sorted_dict.keys()) + + + +def find_most_active_nuclides_in_materials( + materials, + exclude=None, +): + non_excluded_nucs = {} + if exclude is None: + excluded_isotopes = [] + else: + if isinstance(exclude, Iterable): + excluded_isotopes = exclude + elif isinstance(exclude, openmc.Material): + excluded_isotopes = exclude.get_nuclides() + + for material in materials: + for key, value in material.get_nuclide_activity().items(): + if key not in excluded_isotopes: + if key not in non_excluded_nucs.keys(): + if value != 0.: + non_excluded_nucs[key] = value + else: + non_excluded_nucs[key] += value + + sorted_dict = { + k: v + for k, v in sorted( + non_excluded_nucs.items(), key=lambda item: item[1], reverse=True + ) + } + + return list(sorted_dict.keys()) + def find_most_abundant_nuclides_in_material( material, @@ -19,9 +91,9 @@ def find_most_abundant_nuclides_in_material( for key, value in material.get_nuclide_atom_densities().items(): if key not in excluded_isotopes: if key not in non_excluded_nucs.keys(): - non_excluded_nucs[key] = value[1] + non_excluded_nucs[key] = value else: - non_excluded_nucs[key] += value[1] + non_excluded_nucs[key] += value sorted_dict = { k: v @@ -49,9 +121,9 @@ def find_most_abundant_nuclides_in_materials( for key, value in material.get_nuclide_atom_densities().items(): if key not in excluded_isotopes: if key not in non_excluded_nucs.keys(): - non_excluded_nucs[key] = value[1] + non_excluded_nucs[key] = value else: - non_excluded_nucs[key] += value[1] + non_excluded_nucs[key] += value sorted_dict = { k: v @@ -61,3 +133,144 @@ def find_most_abundant_nuclides_in_materials( } return list(sorted_dict.keys()) + + +def get_nuclide_atom_densities_from_materials(nuclides, materials): + all_nuclides_with_atoms = {} + for isotope in nuclides: + all_quants = [] + for material in materials: + quants = material.get_nuclide_atom_densities() + if isotope in quants.keys(): + quant = quants[isotope] + else: + quant=0. + all_quants.append(quant) + all_nuclides_with_atoms[isotope] = all_quants + return all_nuclides_with_atoms + + +def get_nuclide_activities_from_materials(nuclides, materials): + all_nuclides_with_atoms = {} + for isotope in nuclides: + all_quants = [] + for material in materials: + quants = material.get_nuclide_activity() + if isotope in quants.keys(): + quant = quants[isotope] + else: + quant=0. + all_quants.append(quant) + all_nuclides_with_atoms[isotope] = all_quants + return all_nuclides_with_atoms + + +def get_atoms_from_material(material): + + if material.volume is None: + msg = "material.volume must be set to find the number of atoms present." + raise ValueError(msg) + + # in units of atom / ( barn cm2 ) + atoms_per_barn_cm2 = material.get_nuclide_atom_densities() + volume = material.volume * ureg.cm**3 + + # print(atoms_per_barn_cm2) + isotopes_and_atoms = [] + for key, value in atoms_per_barn_cm2.items(): + # print(key, value[1]) + atoms_per_b_cm = value * ureg.particle / (ureg.barn * ureg.cm) + atoms = atoms_per_b_cm * volume + # print(key, atoms) + # print(key, atoms.to(ureg.particle)) + + atomic_number, mass_number, _ = openmc.data.zam(key) + + isotopes_and_atoms.append( + ( + atomic_number, + mass_number - atomic_number, + atoms.to(ureg.particle).magnitude, + key, + ) + ) + + # print(atoms_per_barn_cm2) + return isotopes_and_atoms + + +def build_grid_of_nuclides(iterable_of_nuclides): + + grid_width = 200 # protons + grid_height = 200 # neutrons + grid = np.array([[0] * grid_width] * grid_height, dtype=float) + + for atomic_number, neutron_number, atoms, _ in iterable_of_nuclides: + # print(atomic_number,neutron_number,atoms ) + grid[atomic_number][neutron_number] = atoms + + return grid + + +def build_grid_of_stable_nuclides(iterable_of_nuclides): + + grid_width = 200 # protons + grid_height = 200 # neutrons + grid = np.array([[0] * grid_width] * grid_height, dtype=float) + + for atomic_number, neutron_number in iterable_of_nuclides: + # print(atomic_number,neutron_number,atoms ) + grid[atomic_number][ + neutron_number + ] = 0.2 # sets the grey scale from 0 (white) to 1 (black) + + return grid + + +def build_grid_of_annotations(iterable_of_nuclides): + + grid_width = 200 # protons + grid_height = 200 # neutrons + grid = np.array([[0] * grid_width] * grid_height, dtype=str) + + for atomic_number, neutron_number, atoms, name in iterable_of_nuclides: + grid[atomic_number][neutron_number] = name + + return grid + + +def get_neutron_range(material): + nucs = material.get_nuclides() + + neutrons = [] + for nuc in nucs: + proton, proton_plus_neutron, _ = openmc.data.zam(nuc) + neutron = proton_plus_neutron - proton + neutrons.append(neutron) + return [min(neutrons), max(neutrons)] + + +def get_proton_range(material): + nucs = material.get_nuclides() + + protons = [] + for nuc in nucs: + proton, proton_plus_neutron, _ = openmc.data.zam(nuc) + protons.append(proton) + return [min(protons), max(protons)] + + +# Get the colormap and set the under and bad colors +# colMap = cm.gist_rainbow +def make_unstable_cm(): + colMap = cm.get_cmap("viridis", 256) + colMap.set_bad(color="white", alpha=100) + colMap.set_under(color="white", alpha=100) + return colMap + + +def make_stable_cm(): + colMap = cm.get_cmap("gray_r", 256) + colMap.set_bad(color="white", alpha=100) + colMap.set_under(color="white", alpha=100) + return colMap \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 4adf21d..9bbf860 100644 --- a/setup.cfg +++ b/setup.cfg @@ -29,7 +29,7 @@ python_requires= >=3.6 install_requires= numpy >= 1.21.1 matplotlib - seaborn + plotly pint diff --git a/tests/test_python_api.py b/tests/test_python_api.py index 1b496f1..ecd6c38 100644 --- a/tests/test_python_api.py +++ b/tests/test_python_api.py @@ -1,4 +1,8 @@ -from openmc_depletion_plotter import find_most_abundant_nuclides_in_material, find_most_abundant_nuclides_in_materials +from openmc_depletion_plotter import find_most_abundant_nuclides_in_material +from openmc_depletion_plotter import find_most_abundant_nuclides_in_materials +from openmc_depletion_plotter import find_most_active_nuclides_in_material +from openmc_depletion_plotter import find_most_active_nuclides_in_materials +from openmc_depletion_plotter import get_nuclide_atom_densities_from_materials import openmc @@ -121,3 +125,57 @@ def test_openmc_material_shared_isotope(): ) assert nucs == ['Fe56'] + + +def test_get_nuclide_atom_densities_from_materials(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 0.5) + my_mat.add_nuclide("Li7", 0.5) + + my_mat_2 = openmc.Material() + my_mat_2.add_nuclide("Fe56", 0.6) + my_mat_2.add_nuclide("Li7", 0.6) + + nucs = get_nuclide_atom_densities_from_materials( + nuclides=['Li6', 'Fe56'], + materials=[my_mat, my_mat_2] + ) + + assert list(nucs.keys()) == ['Li6', 'Fe56'] + + # first material + assert isinstance(nucs['Li6'][0], float) + assert isinstance(nucs['Fe56'][0], float) + + # second material + assert isinstance(nucs['Li6'][1], float) + assert isinstance(nucs['Fe56'][1], float) + + +def test_find_most_active_nuclides_in_material(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 0.5) + my_mat.add_nuclide("Li7", 0.5) + my_mat.add_nuclide("U236", 0.5) # unstable + my_mat.volume = 1 + + assert find_most_active_nuclides_in_material(my_mat) == ['U236'] + + +def test_find_most_active_nuclides_in_material(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 0.5) + my_mat.add_nuclide("Li7", 0.5) + my_mat.add_nuclide("U236", 0.5) # unstable + my_mat.volume = 1 + + my_mat2 = openmc.Material() + my_mat2.add_nuclide("Li6", 0.5) + my_mat2.add_nuclide("Li7", 0.5) + my_mat2.add_nuclide("U238", 0.5) # unstable + my_mat2.volume = 1 + + assert find_most_active_nuclides_in_materials([my_mat, my_mat2]) == ['U236', 'U238'] \ No newline at end of file