diff --git a/src/openmc_depletion_plotter/__init__.py b/src/openmc_depletion_plotter/__init__.py index fe32711..f77bc81 100644 --- a/src/openmc_depletion_plotter/__init__.py +++ b/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 diff --git a/src/openmc_depletion_plotter/results.py b/src/openmc_depletion_plotter/results.py index 8d42f9b..101165b 100644 --- a/src/openmc_depletion_plotter/results.py +++ b/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 diff --git a/src/openmc_depletion_plotter/utils.py b/src/openmc_depletion_plotter/utils.py index 51d2156..1ebd552 100644 --- a/src/openmc_depletion_plotter/utils.py +++ b/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: diff --git a/tests/test_plots.py b/tests/test_plots.py new file mode 100644 index 0000000..8fb85a7 --- /dev/null +++ b/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) diff --git a/tests/test_python_api.py b/tests/test_python_api.py index ed16805..076c58e 100644 --- a/tests/test_python_api.py +++ b/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 diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 0000000..076c58e --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,270 @@ +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_atoms_from_materials +from openmc_depletion_plotter import get_decay_heat_from_materials +import openmc + + +def test_find_nuclides_iron(): + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + + nucs = find_most_abundant_nuclides_in_material(material=my_mat) + + assert nucs == ["Fe56", "Fe54", "Fe57", "Fe58"] + + +def test_find_nuclides_lithium(): + + my_mat = openmc.Material() + my_mat.add_element("Li", 1) + + nucs = find_most_abundant_nuclides_in_material(material=my_mat) + + assert nucs == ["Li7", "Li6"] + + +def test_find_nuclides_lithium_enriched(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 1) + my_mat.add_nuclide("Li7", 0.5) + + nucs = find_most_abundant_nuclides_in_material(material=my_mat) + + assert nucs == ["Li6", "Li7"] + + +def test_find_nuclides_iron_single_exclusion(): + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + + nucs = find_most_abundant_nuclides_in_material(material=my_mat, exclude=["Fe56"]) + + assert nucs == ["Fe54", "Fe57", "Fe58"] + + +def test_find_nuclides_iron_multiple_exclusion(): + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + + nucs = find_most_abundant_nuclides_in_material( + material=my_mat, exclude=["Fe56", "Fe57"] + ) + + assert nucs == ["Fe54", "Fe58"] + + +def test_two_identical_materials(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 1) + my_mat.add_nuclide("Li7", 0.5) + my_mat.volume = 1 + + my_mat_2 = openmc.Material() + my_mat_2.add_nuclide("Li6", 1) + my_mat_2.add_nuclide("Li7", 0.5) + my_mat_2.volume = 1 + + nucs = find_most_abundant_nuclides_in_materials( + materials=[my_mat, my_mat_2], + ) + + assert nucs == ["Li6", "Li7"] + + +def test_two_similar_materials(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 0.5) + my_mat.add_nuclide("Li7", 0.5) + my_mat.volume = 1 + + my_mat_2 = openmc.Material() + my_mat_2.add_nuclide("Li6", 0.6) + my_mat_2.add_nuclide("Li7", 0.7) + my_mat_2.volume = 1 + + nucs = find_most_abundant_nuclides_in_materials( + materials=[my_mat, my_mat_2], + ) + + assert nucs == ["Li7", "Li6"] + + +def test_openmc_material(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 0.5) + my_mat.add_nuclide("Li7", 0.5) + my_mat.volume = 1 + + my_mat_2 = openmc.Material() + my_mat_2.add_nuclide("Fe56", 0.6) + my_mat_2.add_nuclide("Be9", 0.6) + my_mat_2.volume = 1 + + nucs = find_most_abundant_nuclides_in_materials( + materials=[my_mat, my_mat_2], exclude=my_mat + ) + + assert nucs == ["Fe56", "Be9"] + + +def test_openmc_material_shared_isotope(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 0.5) + my_mat.add_nuclide("Li7", 0.5) + my_mat.volume = 1 + + my_mat_2 = openmc.Material() + my_mat_2.add_nuclide("Fe56", 0.6) + my_mat_2.add_nuclide("Li7", 0.6) + my_mat_2.volume = 1 + + nucs = find_most_abundant_nuclides_in_materials( + materials=[my_mat, my_mat_2], exclude=my_mat + ) + + assert nucs == ["Fe56"] + + +def test_get_nuclide_atoms_from_materials(): + + my_mat = openmc.Material() + my_mat.add_nuclide("Li6", 0.5) + my_mat.add_nuclide("Li7", 0.5) + my_mat.volume = 1 + + my_mat_2 = openmc.Material() + my_mat_2.add_nuclide("Fe56", 0.6) + my_mat_2.add_nuclide("Li7", 0.6) + my_mat_2.volume = 1 + + nucs = get_nuclide_atoms_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(material=my_mat, units="Bq") == [ + "U236" + ] + assert find_most_active_nuclides_in_material(material=my_mat, units="Bq/g") == [ + "U236" + ] + + +def test_find_most_active_nuclides_in_materials(): + + 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( + materials=[my_mat, my_mat2], units="Bq" + ) == ["U236", "U238"] + 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