added plot_decay_heat method

src/openmc_depletion_plotter/__init__.py

@@ -25,6 +25,7 @@
 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_decay_heat_from_materials
 from .utils import get_atoms_from_material
 from .utils import create_base_plot
 from .utils import add_stables

src/openmc_depletion_plotter/results.py

@@ -4,6 +4,7 @@
 from .utils import find_most_active_nuclides_in_materials
 from .utils import get_nuclide_activities_from_materials
 from .utils import get_nuclide_activities_from_materials
+from .utils import get_decay_heat_from_materials
 from .utils import find_total_activity_in_materials
 from .utils import stable_nuclides
 from .utils import create_base_plot
@@ -162,7 +163,7 @@ def plot_atoms_vs_time(
     x_axis_title=None,
     plotting_backend="plotly",
     threshold=None,
-    title="Number of of nuclides in material",
+    title="Number of nuclides in material",
     material_index=0,  # zero index as one material in problem
     path="materials.xml",
 ):
@@ -174,6 +175,7 @@ def plot_atoms_vs_time(
         materials = self.export_to_materials(nuc_with_data=lots_of_nuclides,burnup_index=counter, path=path)[
             material_index
         ]
+
         all_materials.append(materials)
 
     most_abundant = find_most_abundant_nuclides_in_materials(
@@ -265,5 +267,126 @@ def plot_atoms_vs_time(
         return plt
 
 
+def plot_decay_heat_vs_time(
+    self,
+    excluded_material=None,
+    time_units="d",
+    show_top=None,
+    x_scale="linear",
+    y_scale="linear",
+    include_total=False,
+    x_axis_title=None,
+    plotting_backend="plotly",
+    units="W/g",
+    threshold=None,
+    title="Decay heat of nuclides in material",
+    material_index=0,  # zero index as one material in problem
+    path="materials.xml",
+):
+
+    time_steps = self.get_times(time_units=time_units)
+
+    all_materials = []
+    for counter, step in enumerate(time_steps):
+        materials = self.export_to_materials(nuc_with_data=lots_of_nuclides,burnup_index=counter, path=path)[
+            material_index
+        ]
+
+        all_materials.append(materials)
+
+    most_abundant = find_most_abundant_nuclides_in_materials(
+        materials=all_materials, exclude=excluded_material
+    )
+    if show_top is not None:
+        nuclides = most_abundant[:show_top]
+    else:
+        nuclides = most_abundant
+
+    all_nuclides_with_decay_heat = get_decay_heat_from_materials(
+        nuclides=nuclides, materials=all_materials, units=units
+    )
+
+    if x_axis_title is None:
+        x_axis_title = f"Time [{time_units}]"
+    if plotting_backend == "plotly":
+
+        figure = create_base_plot(
+            x_title=x_axis_title,
+            y_title=f"Decay heat [{units}]",
+            title=title,
+            x_scale=x_scale,
+            y_scale=y_scale,
+        )
+
+        if threshold:
+            total = (
+                find_total_nuclides_in_materials(
+                    all_materials, exclude=excluded_material
+                ),
+            )
+            figure.update_layout(yaxis_range=[threshold, max(total)])
+        add_scale_buttons(figure, x_scale, y_scale)
+    elif plotting_backend == "matplotlib":
+        plt.cla
+        fig = plt.figure()
+        plt.xlabel(x_axis_title)
+        plt.ylabel("Number of atoms")
+        plt.title(title)
+    else:
+        msg = 'only "plotly" and "matplotlib" plotting_backend are supported. {plotting_backend} is not an option'
+        raise ValueError(msg)
+
+    for key, value in all_nuclides_with_decay_heat.items():
+        if threshold:
+            value = np.array(value)
+            value[value < threshold] = np.nan
+        if key in stable_nuclides:
+            name = key + " stable"
+        else:
+            name = key
+        if plotting_backend == "plotly":
+            figure.add_trace(
+                go.Scatter(
+                    mode="lines",
+                    x=time_steps,
+                    y=value,
+                    name=name,
+                    # line=dict(shape="hv", width=0),
+                )
+            )
+        else:
+            plt.plot(time_steps, value, label=key)
+
+    # todo
+    if include_total:
+        print('Total decay heat not implemented yet')
+        # TODO make find_total_decay_heat_in_materials
+        # total = find_total_nuclides_in_materials(
+        #     all_materials, exclude=excluded_material
+        # )
+        if plotting_backend == "plotly":
+            pass
+            # figure.add_trace(
+            #     go.Scatter(
+            #         mode="lines",
+            #         x=time_steps,
+            #         y=total,
+            #         name="total",
+            #         line=dict(dash="longdashdot", color="black"),
+            #     )
+            # )
+        else:
+            pass
+            # plt.plot(time_steps, total, 'k--', label='total')
+
+    if plotting_backend == "plotly":
+        return figure
+    else:
+        plt.legend(bbox_to_anchor=(1.25, 1.0))
+        plt.tight_layout()
+        return plt 
+
+
 openmc.deplete.Results.plot_activity_vs_time = plot_activity_vs_time
 openmc.deplete.Results.plot_atoms_vs_time = plot_atoms_vs_time
+openmc.deplete.Results.plot_decay_heat_vs_time = plot_decay_heat_vs_time

src/openmc_depletion_plotter/utils.py

@@ -312,6 +312,24 @@ def find_most_abundant_nuclides_in_materials(
     return list(sorted_dict.keys())
 
 
+def get_decay_heat_from_materials(nuclides, materials, units):
+    all_nuclides_with_decay_heat = {}
+    for isotope in nuclides:
+        all_quants = []
+        for material in materials:
+            quants = material.get_decay_heat(
+                units=units,
+                by_nuclide=True
+            )
+            if isotope in quants.keys():
+                quant = quants[isotope]
+            else:
+                quant = 0.0
+            all_quants.append(quant)
+        all_nuclides_with_decay_heat[isotope] = all_quants
+    return all_nuclides_with_decay_heat
+
+
 def get_nuclide_atoms_from_materials(nuclides, materials):
     all_nuclides_with_atoms = {}
     for isotope in nuclides:
@@ -332,7 +350,10 @@ def get_nuclide_activities_from_materials(nuclides, materials, units):
     for isotope in nuclides:
         all_quants = []
         for material in materials:
-            quants = material.get_activity(by_nuclide=True, units=units)  # units in Bq
+            quants = material.get_activity(
+                by_nuclide=True,  # units in Bq
+                units=units
+            )
             if isotope in quants.keys():
                 quant = quants[isotope]
             else:

tests/test_plots.py

@@ -0,0 +1,111 @@
+import openmc
+import openmc_depletion_plotter
+import plotly
+import openmc.deplete
+import math
+
+
+#TODO plot
+# openmc.deplete.PredictorIntegrator plot_pulse_schedule
+
+# openmc.deplete.Results.plot_atoms_vs_time()
+# openmc.deplete.Results.plot_activity_vs_time()
+
+def test_default_isotope_charts():
+    my_mat = openmc.Material()
+    my_mat.add_element("Fe", 1)
+    my_mat.add_nuclide("Co60", 1)
+    my_mat.set_density('g/cm3', 1)
+    my_mat.volume = 1
+
+    plot = my_mat.plot_isotope_chart_of_atoms()
+    assert isinstance(plot, plotly.graph_objs._figure.Figure)
+
+    plot = my_mat.plot_isotope_chart_of_activity()
+    assert isinstance(plot, plotly.graph_objs._figure.Figure)
+
+
+def test_default_time_plots():
+
+    # makes a simple material from Silver
+    my_material = openmc.Material() 
+    my_material.add_element('Ag', 1, percent_type='ao')
+    my_material.set_density('g/cm3', 10.49)
+
+    sphere_radius = 100
+    volume_of_sphere = (4/3) * math.pi * math.pow(sphere_radius, 3)
+    my_material.volume = volume_of_sphere  # a volume is needed so openmc can find the number of atoms in the cell/material
+    my_material.depletable = True  # depletable = True is needed to tell openmc to update the material with each time step
+
+    materials = openmc.Materials([my_material])
+    materials.export_to_xml()
+
+    # GEOMETRY
+
+    # surfaces
+    sph1 = openmc.Sphere(r=sphere_radius, boundary_type='vacuum')
+
+    # cells, makes a simple sphere cell
+    shield_cell = openmc.Cell(region=-sph1)
+    shield_cell.fill = my_material
+
+    # sets the geometry to the universe that contains just the one cell
+    universe = openmc.Universe(cells=[shield_cell])
+    geometry = openmc.Geometry(universe)
+
+    # creates a 14MeV neutron point source
+    source = openmc.Source()
+    source.space = openmc.stats.Point((0, 0, 0))
+    source.angle = openmc.stats.Isotropic()
+    source.energy = openmc.stats.Discrete([14e6], [1])
+    source.particles = 'neutron'
+
+    # SETTINGS
+
+    # Instantiate a Settings object
+    settings = openmc.Settings()
+    settings.batches = 2
+    settings.inactive = 0
+    settings.particles = 10000
+    settings.source = source
+    settings.run_mode = 'fixed source'
+
+    model = openmc.model.Model(geometry, materials, settings)
+
+    operator = openmc.deplete.CoupledOperator(
+        model=model,
+        normalization_mode="source-rate",  # set for fixed source simulation, otherwise defaults to fission simulation
+        dilute_initial=0,  # set to zero to avoid adding small amounts of isotopes, defaults to adding small amounts of fissionable isotopes
+        reduce_chain=True,  # reduced to only the isotopes present in depletable materials and their possible progeny
+        reduce_chain_level=5,
+    )
+
+    # We define timesteps together with the source rate to make it clearer
+    timesteps_and_source_rates = [
+        (24, 1e20),
+        (24, 0),
+        (24, 0),
+    ]
+
+    # Uses list Python comprehension to get the timesteps and source_rates separately
+    timesteps = [item[0] for item in timesteps_and_source_rates]
+    source_rates = [item[1] for item in timesteps_and_source_rates]
+
+    integrator = openmc.deplete.PredictorIntegrator(
+        operator=operator,
+        timesteps=timesteps,
+        source_rates=source_rates
+    )
+
+    integrator.integrate()
+
+    results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
+
+    plot = results.plot_atoms_vs_time(excluded_material=my_material)
+    assert isinstance(plot, plotly.graph_objs._figure.Figure)
+
+    plot = results.plot_activity_vs_time()
+    assert isinstance(plot, plotly.graph_objs._figure.Figure)
+
+    plot = results.plot_decay_heat_vs_time()
+    assert isinstance(plot, plotly.graph_objs._figure.Figure)

tests/test_python_api.py

@@ -3,6 +3,7 @@
 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_atoms_from_materials
+from openmc_depletion_plotter import get_decay_heat_from_materials
 import openmc
 
 
@@ -198,3 +199,72 @@ def test_find_most_active_nuclides_in_materials():
     assert find_most_active_nuclides_in_materials(
         materials=[my_mat, my_mat2], units="Bq/g"
     ) == ["U236", "U238"]
+
+
+def test_get_nuclide_atoms_from_materials():
+
+    my_mat = openmc.Material()
+    my_mat.add_nuclide("Mn57", 0.5)
+    my_mat.add_nuclide("Mn56", 0.5)
+    my_mat.add_nuclide("Fe53", 0.5)
+    my_mat.volume = 1
+
+    my_mat_2 = openmc.Material()
+    my_mat_2.add_nuclide("Mn57", 0.6)
+    my_mat_2.add_nuclide("Li7", 0.6)
+    my_mat_2.volume = 1
+
+    heat = get_decay_heat_from_materials(
+        nuclides=["Mn57"],
+        materials=[my_mat, my_mat_2],
+        units='W/g'
+    )
+
+    assert list(heat.keys()) == ['Mn57']
+    assert len(heat['Mn57']) == 2
+
+    # one isotope, one material
+    heat = get_decay_heat_from_materials(
+        nuclides=["Mn57"],
+        materials=[my_mat],
+        units='W/g'
+    )
+    assert list(heat.keys()) == ['Mn57']
+    assert len(heat['Mn57']) == 1
+    assert heat['Mn57'][0] != 0
+
+    # one isotope, two materials
+    heat = get_decay_heat_from_materials(
+        nuclides=["Mn57"],
+        materials=[my_mat, my_mat_2],
+        units='W/g'
+    )
+    assert list(heat.keys()) == ['Mn57']
+    assert len(heat['Mn57']) == 2
+    assert heat['Mn57'][0] != 0
+    assert heat['Mn57'][1] != 0
+
+    # two isotopes, one materials
+    heat = get_decay_heat_from_materials(
+        nuclides=["Mn57", "Mn56"],
+        materials=[my_mat],
+        units='W/g'
+    )
+    assert sorted(list(heat.keys())) == sorted(['Mn57', 'Mn56'])
+    assert len(heat['Mn57']) == 1
+    assert heat['Mn57'][0] != 0
+    assert heat['Mn56'][0] != 0
+
+    # two isotopes, two materials
+    heat = get_decay_heat_from_materials(
+        nuclides=["Mn57", "Mn56"],
+        materials=[my_mat, my_mat_2],
+        units='W/g'
+    )
+    assert sorted(list(heat.keys())) == sorted(['Mn57', 'Mn56'])
+    assert len(heat['Mn57']) == 2
+    assert len(heat['Mn56']) == 2
+    assert heat['Mn57'][0] != 0.
+    assert heat['Mn56'][0] != 0.
+    assert heat['Mn57'][1] != 0.
+    assert heat['Mn56'][1] == 0.  # Mn56 is not in 2nd material