diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 9e3211819..6327d4106 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -121,10 +121,13 @@ multi-group mode. ```` Element -------------------- -The ```` 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 ```` 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: ```` Element diff --git a/openmc/settings.py b/openmc/settings.py index db6bedca1..d24af185f 100644 --- a/openmc/settings.py +++ b/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 \ diff --git a/src/global.F90 b/src/global.F90 index 1e837b8c0..8049db6ef 100644 --- a/src/global.F90 +++ b/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 ! ============================================================================ diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c78315cc8..2719aa68f 100644 --- a/src/input_xml.F90 +++ b/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 diff --git a/src/physics.F90 b/src/physics.F90 index 232effebf..842540d16 100644 --- a/src/physics.F90 +++ b/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 !=============================================================================== diff --git a/tests/test_energy_cutoff/inputs_true.dat b/tests/test_energy_cutoff/inputs_true.dat new file mode 100644 index 000000000..0630d2acc --- /dev/null +++ b/tests/test_energy_cutoff/inputs_true.dat @@ -0,0 +1 @@ +4dc6a7b131f6757ecc9b06e93d425712f0813c546ed9eaa0aafa4990383b2f3a5b742a5a39ef5f27300d2c861b344e78de2d2c569b9439e7893dc06b03bf3b02 \ No newline at end of file diff --git a/tests/test_energy_cutoff/results_true.dat b/tests/test_energy_cutoff/results_true.dat new file mode 100644 index 000000000..5ad463b77 --- /dev/null +++ b/tests/test_energy_cutoff/results_true.dat @@ -0,0 +1,3 @@ +tally 1: +sum = 0.000000E+00 +sum_sq = 0.000000E+00 diff --git a/tests/test_energy_cutoff/test_energy_cutoff.py b/tests/test_energy_cutoff/test_energy_cutoff.py new file mode 100755 index 000000000..6467a5925 --- /dev/null +++ b/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()