Merge pull request #42 from fusion-energy/adding_decay_heat_plot

added plot_decay_heat method

.github/workflows/ci.yml

@@ -44,3 +44,6 @@ jobs:
           # plot_isotope_charts_fisson_irradiation.py
 
           python examples/pulse_schedule/plot_pulse_schedule.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

@@ -0,0 +1,94 @@
+
+import openmc
+import openmc_depletion_plotter
+import openmc.deplete
+import math
+from pathlib import Path
+
+# these should be set to where you have the chain and cross section file saved
+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()
+my_material.add_nuclide('Fe56', 71.17, percent_type='ao')
+my_material.set_density('g/cm3', 7.9)
+
+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()
+
+# 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 = [
+    (5*60, 1.116E+10),
+    (5*60, 0),
+    (5*60, 0),
+    (5*60, 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,
+    timestep_units='s'
+)
+
+integrator.integrate()
+
+results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
+
+plot = results.plot_decay_heat_vs_time(
+    x_scale='log',
+    y_scale='log',
+    excluded_material=my_material,
+    show_top=10
+)
+plot.show()
+
+plot.write_html('decay_heat.html')

pyproject.toml

@@ -34,7 +34,7 @@ write_to = "src/openmc_depletion_plotter/_version.py"
 
 [project.optional-dependencies]
 tests = [
-    "pytest"
+    "pytest",
 ]
 
 [project.scripts]

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/app.py

@@ -89,7 +89,7 @@ def main():
 
         activity_or_atoms = st.sidebar.selectbox(
             label="Plot",
-            options=("activity", "number of atoms"),
+            options=("activity", "number of atoms", "decay heat"),
             index=0,
             key="activity_or_atoms",
             help="",
@@ -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:
@@ -201,6 +209,23 @@ def main():
                 include_total=include_total,
                 path=materials_file.name,
             )
+        elif activity_or_atoms == "decay heat":
+            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,
+                include_total=include_total,
+                # x_axis_title=None,
+                plotting_backend=backend,
+                # units="W/g", # TODO add drop down option
+                # threshold=None,
+                title="Decay heat of nuclides in material",
+                material_index=material_index,
+            )
         elif activity_or_atoms == "number of atoms":
             plot = results.plot_atoms_vs_time(
                 # todo allow no materials to be excluded
@@ -222,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,18 +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 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())
@@ -162,7 +162,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 +174,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(
@@ -213,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'
@@ -265,5 +268,124 @@ 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(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'
+        raise ValueError(msg)
+
+    for key, value in all_nuclides_with_decay_heat.items():
+        if sum(value) != 0:
+            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)
+
+    if include_total:
+        total = find_total_decay_heat_in_materials(
+            materials=all_materials, units=units, exclude=excluded_material
+        )
+        if plotting_backend == "plotly":
+            figure.add_trace(
+                go.Scatter(
+                    mode="lines",
+                    x=time_steps,
+                    y=total,
+                    name="total",
+                    line=dict(dash="longdashdot", color="black"),
+                )
+            )
+        else:
+            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