Update pincell example adding flux spectrum tally and associated plotting script

examples/python/pincell/build-xml.py

@@ -1,21 +1,11 @@
+from math import log10
+
+import numpy as np
 import openmc
 
 ###############################################################################
-#                      Simulation Input File Parameters
-###############################################################################
+# Create materials for the problem
 
-# OpenMC simulation parameters
-batches = 100
-inactive = 10
-particles = 1000
-
-
-###############################################################################
-#                 Exporting to OpenMC materials.xml file
-###############################################################################
-
-
-# Instantiate some Materials and register the appropriate Nuclides
 uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
 uo2.set_density('g/cm3', 10.29769)
 uo2.add_element('U', 1., enrichment=2.4)
@@ -39,14 +29,12 @@ borated_water.add_element('H', 5.0e-2)
 borated_water.add_element('O', 2.4e-2)
 borated_water.add_s_alpha_beta('c_H_in_H2O')
 
-# Instantiate a Materials collection and export to XML
+# Collect the materials together and export to XML
 materials = openmc.Materials([uo2, helium, zircaloy, borated_water])
 materials.export_to_xml()
 
-
-###############################################################################
-#                 Exporting to OpenMC geometry.xml file
 ###############################################################################
+# Define problem geometry
 
 # Create cylindrical surfaces
 fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
@@ -67,22 +55,23 @@ water = openmc.Cell(fill=borated_water, region=+clad_or & box)
 geometry = openmc.Geometry([fuel, gap, clad, water])
 geometry.export_to_xml()
 
-
-###############################################################################
-#                   Exporting to OpenMC settings.xml file
 ###############################################################################
+# Define problem settings
 
-# Instantiate a Settings object, set all runtime parameters, and export to XML
+# Indicate how many particles to run
 settings = openmc.Settings()
-settings.batches = batches
-settings.inactive = inactive
-settings.particles = particles
+settings.batches = 100
+settings.inactive = 10
+settings.particles = 1000
 
 # Create an initial uniform spatial source distribution over fissionable zones
-bounds = [-pitch/2, -pitch/2, -1, pitch/2, pitch/2, 1]
-uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+lower_left = (-pitch/2, -pitch/2, -1)
+upper_right = (pitch/2, pitch/2, 1)
+uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
 settings.source = openmc.source.Source(space=uniform_dist)
 
+# For source convergence checks, add a mesh that can be used to calculate the
+# Shannon entropy
 entropy_mesh = openmc.RegularMesh()
 entropy_mesh.lower_left = (-fuel_or.r, -fuel_or.r)
 entropy_mesh.upper_right = (fuel_or.r, fuel_or.r)
@@ -90,26 +79,34 @@ entropy_mesh.dimension = (10, 10)
 settings.entropy_mesh = entropy_mesh
 settings.export_to_xml()
 
-
-###############################################################################
-#                   Exporting to OpenMC tallies.xml file
 ###############################################################################
+# Define tallies
 
-# Instantiate a tally mesh
+# Create a mesh that will be used for tallying
 mesh = openmc.RegularMesh()
 mesh.dimension = (100, 100)
 mesh.lower_left = (-pitch/2, -pitch/2)
 mesh.upper_right = (pitch/2, pitch/2)
 
-# Instantiate some tally Filters
-energy_filter = openmc.EnergyFilter((0., 4., 20.e6))
+# Create a mesh filter that can be used in a tally
 mesh_filter = openmc.MeshFilter(mesh)
 
-# Instantiate the Tally
-tally = openmc.Tally()
-tally.filters = [energy_filter, mesh_filter]
-tally.scores = ['flux', 'fission', 'nu-fission']
+# Now use the mesh filter in a tally and indicate what scores are desired
+mesh_tally = openmc.Tally(name="Mesh tally")
+mesh_tally.filters = [mesh_filter]
+mesh_tally.scores = ['flux', 'fission', 'nu-fission']
+
+# Let's also create a tally to get the flux energy spectrum. We start by
+# creating an energy filter
+e_min, e_max = 1e-5, 20.0e6
+groups = 500
+energies = np.logspace(log10(e_min), log10(e_max), groups + 1)
+energy_filter = openmc.EnergyFilter(energies)
+
+spectrum_tally = openmc.Tally(name="Flux spectrum")
+spectrum_tally.filters = [energy_filter]
+spectrum_tally.scores = ['flux']
 
 # Instantiate a Tallies collection and export to XML
-tallies = openmc.Tallies([tally])
+tallies = openmc.Tallies([mesh_tally, spectrum_tally])
 tallies.export_to_xml()

examples/python/pincell/plot_spectrum.py

@@ -0,0 +1,23 @@
+import matplotlib.pyplot as plt
+import openmc
+
+
+# Get results from statepoint
+with openmc.StatePoint('statepoint.100.h5') as sp:
+    t = sp.get_tally(name="Flux spectrum")
+
+    # Get the energies from the energy filter
+    energy_filter = t.filters[0]
+    energies = energy_filter.bins[:, 0]
+
+    # Get the flux values
+    mean = t.get_values(value='mean').ravel()
+    uncertainty = t.get_values(value='std_dev').ravel()
+
+# Plot flux spectrum
+fix, ax = plt.subplots()
+ax.loglog(energies, mean, drawstyle='steps-post')
+ax.set_xlabel('Energy [eV]')
+ax.set_ylabel('Flux')
+ax.grid(True, which='both')
+plt.show()