Merge remote-tracking branch 'upstream/develop' into mg_th5

docs/source/usersguide/input.rst

@@ -121,10 +121,13 @@ multi-group mode.
 ``<cutoff>`` Element
 --------------------
 
-The ``<cutoff>`` element indicates the weight cutoff used below which particles
-undergo Russian roulette. Surviving particles are assigned a user-determined
-weight. Note that weight cutoffs and Russian rouletting are not turned on by
-default. This element has the following attributes/sub-elements:
+The ``<cutoff>`` element indicates two kinds of cutoffs. The first is the weight
+cutoff used below which particles undergo Russian roulette. Surviving particles
+are assigned a user-determined weight. Note that weight cutoffs and Russian
+rouletting are not turned on by default. The second is the energy cutoff which
+is used to kill particles under certain energy. The energy cutoff should not be
+used unless you know particles under the energy are of no importance to results
+you care. This element has the following attributes/sub-elements:
 
   :weight:
     The weight below which particles undergo Russian roulette.
@@ -137,6 +140,11 @@ default. This element has the following attributes/sub-elements:
 
     *Default*: 1.0
 
+  :energy:
+    The energy under which particles will be killed.
+
+    *Default*: 0.0
+
 .. _eigenvalue:
 
 ``<eigenvalue>`` Element

openmc/settings.py

@@ -89,10 +89,14 @@ class Settings(object):
         Seed for the linear congruential pseudorandom number generator
     survival_biasing : bool
         Indicate whether survival biasing is to be used
-    weight : float
-        Weight cutoff below which particle undergo Russian roulette
-    weight_avg : float
-        Weight assigned to particles that are not killed after Russian roulette
+    cutoff : dict
+        Dictionary defining weight cutoff and energy cutoff. The dictionary may
+        have three keys, 'weight', 'weight_avg' and 'energy'. Value for 'weight'
+        should be a float indicating weight cutoff below which particle undergo
+        Russian roulette. Value for 'weight_avg' should be a float indicating
+        weight assigned to particles that are not killed after Russian
+        roulette. Value of energy should be a float indicating energy in MeV
+        below which particle will be killed.
     entropy_dimension : tuple or list
         Number of Shannon entropy mesh cells in the x, y, and z directions,
         respectively
@@ -214,8 +218,7 @@ class Settings(object):
         self._temperature = {}
 
         # Cutoff subelement
-        self._weight = None
-        self._weight_avg = None
+        self._cutoff = None
 
         # Uniform fission source subelement
         self._ufs_dimension = 1
@@ -393,12 +396,8 @@ class Settings(object):
         return self._track
 
     @property
-    def weight(self):
-        return self._weight
-
-    @property
-    def weight_avg(self):
-        return self._weight_avg
+    def cutoff(self):
+        return self._cutoff
 
     @property
     def ufs_dimension(self):
@@ -635,17 +634,29 @@ class Settings(object):
         cv.check_type('survival biasing', survival_biasing, bool)
         self._survival_biasing = survival_biasing
 
-    @weight.setter
-    def weight(self, weight):
-        cv.check_type('weight cutoff', weight, Real)
-        cv.check_greater_than('weight cutoff', weight, 0.0)
-        self._weight = weight
+    @cutoff.setter
+    def cutoff(self, cutoff):
+        if not isinstance(cutoff, Mapping):
+            msg = 'Unable to set cutoff from "{0}" which is not a '\
+                  ' Python dictionary'.format(cutoff)
+            raise ValueError(msg)
+        for key in cutoff:
+            if key == 'weight':
+                cv.check_type('weight cutoff', cutoff['weight'], Real)
+                cv.check_greater_than('weight cutoff', cutoff['weight'], 0.0)
+            elif key == 'weight_avg':
+                cv.check_type('average survival weight', cutoff['weight_avg'],
+                              Real)
+                cv.check_greater_than('average survival weight',
+                                      cutoff['weight_avg'], 0.0)
+            elif key == 'energy':
+                cv.check_type('energy cutoff', cutoff['energy'], Real)
+                cv.check_greater_than('energy cutoff', cutoff['energy'], 0.0)
+            else:
+                msg = 'Unable to set cutoff to "{0}" which is unsupported by '\
+                      'OpenMC'.format(key)
 
-    @weight_avg.setter
-    def weight_avg(self, weight_avg):
-        cv.check_type('average survival weight', weight_avg, Real)
-        cv.check_greater_than('average survival weight', weight_avg, 0.0)
-        self._weight_avg = weight_avg
+        self._cutoff = cutoff
 
     @entropy_dimension.setter
     def entropy_dimension(self, dimension):
@@ -1030,14 +1041,19 @@ class Settings(object):
             element.text = str(self._survival_biasing).lower()
 
     def _create_cutoff_subelement(self):
-        if self._weight is not None:
+        if self._cutoff is not None:
             element = ET.SubElement(self._settings_file, "cutoff")
+            if 'weight' in self._cutoff:
+                subelement = ET.SubElement(element, "weight")
+                subelement.text = str(self._cutoff['weight'])
 
-            subelement = ET.SubElement(element, "weight")
-            subelement.text = str(self._weight)
+            if 'weight_avg' in self._cutoff:
+                subelement = ET.SubElement(element, "weight_avg")
+                subelement.text = str(self._cutoff['weight_avg'])
 
-            subelement = ET.SubElement(element, "weight_avg")
-            subelement.text = str(self._weight_avg)
+            if 'energy' in self._cutoff:
+                subelement = ET.SubElement(element, "energy")
+                subelement.text = str(self._cutoff['energy'])
 
     def _create_entropy_subelement(self):
         if self._entropy_lower_left is not None and \

src/global.F90

@@ -307,6 +307,7 @@ module global
 
   logical :: survival_biasing = .false.
   real(8) :: weight_cutoff = 0.25_8
+  real(8) :: energy_cutoff = ZERO
   real(8) :: weight_survive = ONE
 
   ! ============================================================================

src/input_xml.F90

@@ -670,8 +670,15 @@ contains
     ! Cutoffs
     if (check_for_node(doc, "cutoff")) then
       call get_node_ptr(doc, "cutoff", node_cutoff)
-      call get_node_value(node_cutoff, "weight", weight_cutoff)
-      call get_node_value(node_cutoff, "weight_avg", weight_survive)
+      if (check_for_node(node_cutoff, "weight")) then
+        call get_node_value(node_cutoff, "weight", weight_cutoff)
+      end if
+      if (check_for_node(node_cutoff, "weight_avg")) then
+        call get_node_value(node_cutoff, "weight_avg", weight_survive)
+      end if
+      if (check_for_node(node_cutoff, "energy")) then
+        call get_node_value(node_cutoff, "energy", energy_cutoff)
+      end if
     end if
 
     ! Particle trace

src/physics.F90

@@ -127,6 +127,13 @@ contains
       if (.not. p % alive) return
     end if
 
+    ! Kill neutron under certain energy
+    if (p % E < energy_cutoff) then
+      p % alive = .false.
+      p % wgt = ZERO
+      p % last_wgt = ZERO
+    end if
+
   end subroutine sample_reaction
 
 !===============================================================================

tests/test_energy_cutoff/inputs_true.dat

@@ -0,0 +1 @@
+4dc6a7b131f6757ecc9b06e93d425712f0813c546ed9eaa0aafa4990383b2f3a5b742a5a39ef5f27300d2c861b344e78de2d2c569b9439e7893dc06b03bf3b02

tests/test_energy_cutoff/results_true.dat

@@ -0,0 +1,3 @@
+tally 1:
+sum = 0.000000E+00
+sum_sq = 0.000000E+00

tests/test_energy_cutoff/test_energy_cutoff.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+
+import os
+import sys
+sys.path.insert(0, os.pardir)
+from testing_harness import PyAPITestHarness
+import openmc
+
+
+class EnergyCutoffTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+        # Set energy cutoff
+        energy_cutoff = 4e-6
+
+        # Material is composed of H-1
+        mat = openmc.Material(material_id=1, name='mat')
+        mat.set_density('atom/b-cm', 0.069335)
+        mat.add_nuclide('H1', 40.0)
+        materials_file = openmc.Materials([mat])
+        materials_file.export_to_xml()
+
+        # Cell is box with reflective boundary
+        x1 = openmc.XPlane(surface_id=1, x0=-1)
+        x2 = openmc.XPlane(surface_id=2, x0=1)
+        y1 = openmc.YPlane(surface_id=3, y0=-1)
+        y2 = openmc.YPlane(surface_id=4, y0=1)
+        z1 = openmc.ZPlane(surface_id=5, z0=-1)
+        z2 = openmc.ZPlane(surface_id=6, z0=1)
+        for surface in [x1, x2, y1, y2, z1, z2]:
+            surface.boundary_type = 'reflective'
+        box = openmc.Cell(cell_id=1, name='box')
+        box.region = +x1 & -x2 & +y1 & -y2 & +z1 & -z2
+        box.fill = mat
+        root = openmc.Universe(universe_id=0, name='root universe')
+        root.add_cell(box)
+        geometry = openmc.Geometry(root)
+        geometry.export_to_xml()
+
+        # Set the running parameters
+        settings_file = openmc.Settings()
+        settings_file.run_mode = 'fixed source'
+        settings_file.batches = 10
+        settings_file.particles = 100
+        settings_file.cutoff = {'energy': energy_cutoff}
+        bounds = [-1, -1, -1, 1, 1, 1]
+        uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
+        watt_dist = openmc.stats.Watt()
+        settings_file.source = openmc.source.Source(space=uniform_dist,
+                                                    energy=watt_dist)
+        settings_file.export_to_xml()
+
+        # Tally flux under energy cutoff
+        tallies = openmc.Tallies()
+        tally = openmc.Tally(1)
+        tally.scores = ['flux']
+        energy_filter = openmc.filter.EnergyFilter((0.0, energy_cutoff))
+        tally.filters = [energy_filter]
+        tallies.append(tally)
+        tallies.export_to_xml()
+
+    def _get_results(self):
+        """Digest info in the statepoint and return as a string."""
+        # Read the statepoint file.
+        sp = openmc.StatePoint(self._sp_name)
+
+        # Write out tally data.
+        outstr = ''
+        t = sp.get_tally()
+        outstr += 'tally {0}:\n'.format(t.id)
+        outstr += 'sum = {0:12.6E}\n'.format(t.sum[0, 0, 0])
+        outstr += 'sum_sq = {0:12.6E}\n'.format(t.sum_sq[0, 0, 0])
+
+        return outstr
+
+    def _cleanup(self):
+        super(EnergyCutoffTestHarness, self)._cleanup()
+        f = os.path.join(os.getcwd(), 'tallies.xml')
+        if os.path.exists(f):
+            os.remove(f)
+
+
+if __name__ == '__main__':
+    harness = EnergyCutoffTestHarness('statepoint.10.h5', True)
+    harness.main()