added decay heat to app

.github/workflows/ci.yml

@@ -44,4 +44,6 @@ jobs:
           # plot_isotope_charts_fisson_irradiation.py
 
           python examples/pulse_schedule/plot_pulse_schedule.py
-          python examples/decay_heat_vs_time/decay_heat_vs_time.py
+
+          # simulation requires more nuclear data that available on the CI
+          # python examples/decay_heat_vs_time/decay_heat_vs_time.py

examples/decay_heat_vs_time/decay_heat_vs_time.py

@@ -8,8 +8,8 @@
 from pathlib import Path
 
 # these should be set to where you have the chain and cross section file saved
-openmc.config['cross_sections'] = '/home/jshimwell/openmc_depletion_plotter/tests/cross_sections.xml'
-openmc.config['chain_file'] = '/home/jshimwell/openmc_depletion_plotter/tests/chain-nndc-b7.1.xml'
+openmc.config['cross_sections'] = Path(__file__).parents[2]/'tests'/'cross_sections.xml'
+openmc.config['chain_file'] = Path(__file__).parents[2]/'tests'/'chain-nndc-b7.1.xml'
 
 # makes a simple material
 my_material = openmc.Material()

src/openmc_depletion_plotter/app.py

@@ -142,6 +142,14 @@ def main():
                 key="activity_units",
                 help="",
             )
+        if activity_or_atoms == "decay heat":
+            decay_heat_units = st.sidebar.selectbox(
+                label="Decay heat units",
+                options=('W', 'W/g', 'W/cm3'),
+                index=0,
+                key="decay_heat_units",
+                help="",
+            )
             # todo horizontal_lines
 
         if number_of_depleted_materials == 1:
@@ -202,9 +210,11 @@ def main():
                 path=materials_file.name,
             )
         elif activity_or_atoms == "decay heat":
-            plot = plot_decay_heat_vs_time(
+            print('x_scale line 213 ',x_scale)
+            plot = results.plot_decay_heat_vs_time(
                 excluded_material=material_to_exclude,
                 time_units=time_units,
+                units=decay_heat_units,
                 show_top=show_top,
                 x_scale=x_scale,
                 y_scale=y_scale,
@@ -237,7 +247,6 @@ def main():
                 'activity_or_atoms must be either "activity" or "number of atoms" to plot'
             )
 
-
         if backend == "matplotlib":
             st.pyplot(plot)
         else:

src/openmc_depletion_plotter/results.py

@@ -1,19 +1,18 @@
+import matplotlib.pyplot as plt
+import numpy as np
 import openmc
 import openmc.deplete
-from .utils import add_scale_buttons
-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
-from .utils import get_nuclide_atoms_from_materials
-from .utils import find_most_abundant_nuclides_in_materials
-from .utils import find_total_nuclides_in_materials
 import plotly.graph_objects as go
-import numpy as np
-import matplotlib.pyplot as plt
+
+from .utils import (add_scale_buttons, create_base_plot,
+                    find_most_abundant_nuclides_in_materials,
+                    find_most_active_nuclides_in_materials,
+                    find_total_activity_in_materials,
+                    find_total_decay_heat_in_materials,
+                    find_total_nuclides_in_materials,
+                    get_decay_heat_from_materials,
+                    get_nuclide_activities_from_materials,
+                    get_nuclide_atoms_from_materials, stable_nuclides)
 
 lots_of_nuclides = []
 elements = list(openmc.data.ATOMIC_SYMBOL.values())
@@ -215,6 +214,8 @@ def plot_atoms_vs_time(
         fig = plt.figure()
         plt.xlabel(x_axis_title)
         plt.ylabel("Number of atoms")
+        plt.xscale(x_scale)
+        plt.yscale(y_scale)
         plt.title(title)
     else:
         msg = 'only "plotly" and "matplotlib" plotting_backend are supported. {plotting_backend} is not an option'
@@ -330,7 +331,9 @@ def plot_decay_heat_vs_time(
         plt.cla
         fig = plt.figure()
         plt.xlabel(x_axis_title)
-        plt.ylabel("Number of atoms")
+        plt.ylabel(f"Decay heat [{units}]")
+        plt.xscale(x_scale)
+        plt.yscale(y_scale)
         plt.title(title)
     else:
         msg = 'only "plotly" and "matplotlib" plotting_backend are supported. {plotting_backend} is not an option'
@@ -358,27 +361,22 @@ def plot_decay_heat_vs_time(
             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
-        # )
+        total = find_total_decay_heat_in_materials(
+            materials=all_materials, units=units, 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"),
-            #     )
-            # )
+            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')
+            plt.plot(time_steps, total, 'k--', label='total')
 
     if plotting_backend == "plotly":
         return figure

src/openmc_depletion_plotter/utils.py

@@ -201,7 +201,7 @@ def find_total_activity_in_materials(
     materials,
     exclude=None,
 ):
-    non_excluded_nucs = []
+
     if exclude is None:
         excluded_isotopes = []
     else:
@@ -224,6 +224,34 @@ def find_total_activity_in_materials(
     return materials_activities
 
 
+def find_total_decay_heat_in_materials(
+    units,
+    materials,
+    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()
+
+    materials_decay_heat = []
+    for material in materials:
+        material_activity = 0
+        heat = material.get_decay_heat(by_nuclide=True, units=units)
+
+        for key, value in heat.items():
+            if key not in excluded_isotopes:
+                material_activity += value
+
+        materials_decay_heat.append(material_activity)
+
+    return materials_decay_heat
+
+
 def find_most_abundant_nuclides_in_material(
     material,
     exclude=None,