Merge pull request #607 from smharper/pyapi_resonance_scattering

Add resonance scattering to the Python API

openmc/settings.py

@@ -9,6 +9,7 @@ import numpy as np
 from openmc.clean_xml import *
 from openmc.checkvalue import (check_type, check_length, check_value,
                                check_greater_than, check_less_than)
+from openmc import Nuclide
 from openmc.source import Source
 
 if sys.version_info[0] >= 3:
@@ -125,6 +126,8 @@ class SettingsFile(object):
         Coordinates of the lower-left point of the UFS mesh
     ufs_upper_right : tuple or list
         Coordinates of the upper-right point of the UFS mesh
+    resonance_scattering : ResonanceScattering or iterable of ResonanceScattering
+        The elastic scattering model to use for resonant isotopes
 
     """
 
@@ -205,6 +208,8 @@ class SettingsFile(object):
         self._run_mode_subelement = None
         self._source_element = None
 
+        self._resonance_scattering = None
+
     @property
     def run_mode(self):
         return self._run_mode
@@ -393,9 +398,13 @@ class SettingsFile(object):
     def dd_count_interactions(self):
         return self._dd_count_interactions
 
+    @property
+    def resonance_scattering(self):
+        return self._resonance_scattering
+
     @run_mode.setter
     def run_mode(self, run_mode):
-        if 'run_mode' not in ['eigenvalue', 'fixed source']:
+        if run_mode not in ['eigenvalue', 'fixed source']:
             msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \
                   'and "fixed source" are supported."'.format(run_mode)
             raise ValueError(msg)
@@ -764,6 +773,16 @@ class SettingsFile(object):
 
         self._dd_count_interactions = interactions
 
+    @resonance_scattering.setter
+    def resonance_scattering(self, res):
+        if isinstance(res, Iterable):
+            check_type('resonance_scattering', res, Iterable,
+                       ResonanceScattering)
+            self._resonance_scattering = res
+        else:
+            check_type('resonance_scattering', res, ResonanceScattering)
+            self._resonance_scattering = [res]
+
     def _create_run_mode_subelement(self):
 
         if self.run_mode == 'eigenvalue':
@@ -1043,6 +1062,17 @@ class SettingsFile(object):
             subelement = ET.SubElement(element, "count_interactions")
             subelement.text = str(self._dd_count_interactions).lower()
 
+    def _create_resonance_scattering_element(self):
+        if self.resonance_scattering is None: return
+
+        element = ET.SubElement(self._settings_file, "resonance_scattering")
+
+        for r in self.resonance_scattering:
+            if r.nuclide.name != r.nuclide_0K.name:
+                raise ValueError("The nuclide and nuclide_0K attributes of "
+                     "a ResonantScattering object must have identical names.")
+            r.create_xml_subelement(element)
+
     def export_to_xml(self):
         """Create a settings.xml file that can be used for a simulation.
 
@@ -1079,6 +1109,7 @@ class SettingsFile(object):
         self._create_track_subelement()
         self._create_ufs_subelement()
         self._create_dd_subelement()
+        self._create_resonance_scattering_element()
 
         # Clean the indentation in the file to be user-readable
         clean_xml_indentation(self._settings_file)
@@ -1087,3 +1118,104 @@ class SettingsFile(object):
         tree = ET.ElementTree(self._settings_file)
         tree.write("settings.xml", xml_declaration=True,
                              encoding='utf-8', method="xml")
+
+
+class ResonanceScattering(object):
+    """Specification of the elastic scattering model for resonant isotopes
+
+    Attributes
+    ----------
+    nuclide : openmc.nuclide.Nuclide
+        The nuclide affected by this resonance scattering treatment.
+    nuclide_0K : openmc.nuclide.Nuclide
+        This should be the same isotope as the nuclide attribute above, but it
+        should have an xs attribute that identifies 0 Kelvin data.
+    method : str
+        The method used to sample outgoing scattering energies.  Valid options
+        are 'ARES', 'CXS' (constant cross section), 'DBRC' (Doppler broadening
+        rejection correction), and 'WCM' (weight correction method).
+    E_min : Real
+        The minimum energy above which the specified method is applied.  By
+        default, CXS will be used below E_min.
+    E_max : Real
+        The maximum energy below which the specified method is applied.  By
+        default, the asymptotic target-at-rest model is applied  above E_max.
+
+    """
+
+    def __init__(self):
+        self._nuclide = None
+        self._nuclide_0K = None
+        self._method = None
+        self._E_min = None
+        self._E_max = None
+
+    @property
+    def nuclide(self):
+        return self._nuclide
+
+    @property
+    def nuclide_0K(self):
+        return self._nuclide_0K
+
+    @property
+    def method(self):
+        return self._method
+
+    @property
+    def E_min(self):
+        return self._E_min
+
+    @property
+    def E_max(self):
+        return self._E_max
+
+    @nuclide.setter
+    def nuclide(self, nuc):
+        check_type('nuclide', nuc, Nuclide)
+        if nuc.zaid == None: raise ValueError("The nuclide must have an "
+             "explicitly defined zaid attribute.")
+        self._nuclide = nuc
+
+    @nuclide_0K.setter
+    def nuclide_0K(self, nuc):
+        check_type('nuclide_0K', nuc, Nuclide)
+        if nuc.zaid == None: raise ValueError("The nuclide_0K must have an "
+             "explicitly defined zaid attribute.")
+        self._nuclide_0K = nuc
+
+    @method.setter
+    def method(self, m):
+        check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM'))
+        self._method = m
+
+    @E_min.setter
+    def E_min(self, E):
+        check_type('E_min', E, Real)
+        check_greater_than('E_min', E, 0, True)
+        self._E_min = E
+
+    @E_max.setter
+    def E_max(self, E):
+        check_type('E_max', E, Real)
+        check_greater_than('E_max', E, 0, True)
+        self._E_max = E
+
+    def create_xml_subelement(self, xml_element):
+        scatterer = ET.SubElement(xml_element, "scatterer")
+        subelement = ET.SubElement(scatterer, 'nuclide')
+        subelement.text = self.nuclide.name
+        if self.method is not None:
+            subelement = ET.SubElement(scatterer, 'method')
+            subelement.text = self.method
+        subelement = ET.SubElement(scatterer, 'xs_label')
+        subelement.text = str(self.nuclide.zaid) + '.' + str(self.nuclide.xs)
+        subelement = ET.SubElement(scatterer, 'xs_label_0K')
+        subelement.text = str(self.nuclide_0K.zaid) + '.' \
+             + str(self.nuclide_0K.xs)
+        if self.E_min is not None:
+            subelement = ET.SubElement(scatterer, 'E_min')
+            subelement.text = str(self.E_min)
+        if self.E_max is not None:
+            subelement = ET.SubElement(scatterer, 'E_max')
+            subelement.text = str(self.E_max)

tests/test_resonance_scattering/geometry.xml

@@ -1,8 +0,0 @@
-<?xml version="1.0"?>
-<geometry>
-
-  <!-- Sphere with radius 10 -->
-  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
-  <cell id="1" material="1" region="-1" />
-    
-</geometry>

tests/test_resonance_scattering/inputs_true.dat

@@ -0,0 +1 @@
+ece83bb075ed8144af89ce7cebf1577dcb2489d2e9ce4afbe61a3e4398837e7a9aaa2ae0cea0a6542f51ca5e0d119b570c675ed1dca0d74237cd5fdce0b606a3

tests/test_resonance_scattering/materials.xml

@@ -1,9 +0,0 @@
-<?xml version="1.0"?>
-<materials>
-
-  <material id="1">
-    <density value="4.5" units="g/cc" />
-    <nuclide name="U-238" xs="71c" ao="1.0" />
-  </material>
-
-</materials>

tests/test_resonance_scattering/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-6.842159E-02 8.481029E-04
+1.440556E+00 6.383274E-02

tests/test_resonance_scattering/settings.xml

@@ -1,27 +0,0 @@
-<?xml version="1.0"?>
-<settings>
-
-  <resonance_scattering>
-    <scatterer>
-      <nuclide>U-238</nuclide>
-      <method>cxs</method>
-      <xs_label>92238.71c</xs_label>
-      <xs_label_0K>92238.71c</xs_label_0K>
-      <E_min>5.0e-6</E_min>
-      <E_max>40.0e-6</E_max>
-    </scatterer>
-  </resonance_scattering>
-
-  <eigenvalue>
-    <batches>10</batches>
-    <inactive>5</inactive>
-    <particles>1000</particles>
-  </eigenvalue>
-
-  <source>
-    <space type="box">
-      <parameters>-4 -4 -4  4  4  4</parameters>
-    </space>
-  </source>
-
-</settings>

tests/test_resonance_scattering/test_resonance_scattering.py

@@ -3,9 +3,81 @@
 import os
 import sys
 sys.path.insert(0, os.pardir)
-from testing_harness import TestHarness
+from testing_harness import PyAPITestHarness
+import openmc
+
+
+class ResonanceScatteringTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+        # Materials
+        mat = openmc.Material(material_id=1)
+        mat.set_density('g/cc', 1.0)
+        mat.add_nuclide('U-238', 1.0)
+        mat.add_nuclide('U-235', 0.02)
+        mat.add_nuclide('Pu-239', 0.02)
+        mat.add_nuclide('H-1', 20.0)
+
+        mats_file = openmc.MaterialsFile()
+        mats_file.default_xs = '71c'
+        mats_file.add_material(mat)
+        mats_file.export_to_xml()
+
+        # Geometry
+        dumb_surface = openmc.XPlane(x0=100)
+        dumb_surface.boundary_type = 'reflective'
+
+        c1 = openmc.Cell(cell_id=1)
+        c1.fill = mat
+        c1.region = -dumb_surface
+
+        root_univ = openmc.Universe(universe_id=0)
+        root_univ.add_cell(c1)
+
+        geometry = openmc.Geometry()
+        geometry.root_universe = root_univ
+        geo_file = openmc.GeometryFile()
+        geo_file.geometry = geometry
+        geo_file.export_to_xml()
+
+        # Settings
+        nuclide = openmc.Nuclide('U-238', '71c')
+        nuclide.zaid = 92238
+        res_scatt_dbrc = openmc.ResonanceScattering()
+        res_scatt_dbrc.nuclide = nuclide
+        res_scatt_dbrc.nuclide_0K = nuclide # This is a bad idea! Just for tests
+        res_scatt_dbrc.method = 'DBRC'
+        res_scatt_dbrc.E_min = 1e-6
+        res_scatt_dbrc.E_max = 210e-6
+
+        nuclide = openmc.Nuclide('U-235', '71c')
+        nuclide.zaid = 92235
+        res_scatt_wcm = openmc.ResonanceScattering()
+        res_scatt_wcm.nuclide = nuclide
+        res_scatt_wcm.nuclide_0K = nuclide
+        res_scatt_wcm.method = 'WCM'
+        res_scatt_wcm.E_min = 1e-6
+        res_scatt_wcm.E_max = 210e-6
+
+        nuclide = openmc.Nuclide('Pu-239', '71c')
+        nuclide.zaid = 94239
+        res_scatt_ares = openmc.ResonanceScattering()
+        res_scatt_ares.nuclide = nuclide
+        res_scatt_ares.nuclide_0K = nuclide
+        res_scatt_ares.method = 'ARES'
+        res_scatt_ares.E_min = 1e-6
+        res_scatt_ares.E_max = 210e-6
+
+        sets_file = openmc.SettingsFile()
+        sets_file.batches = 10
+        sets_file.inactive = 5
+        sets_file.particles = 1000
+        sets_file.source = openmc.source.Source(
+             space=openmc.stats.Box([-4, -4, -4], [4, 4, 4]))
+        sets_file.resonance_scattering = [res_scatt_dbrc, res_scatt_wcm,
+             res_scatt_ares]
+        sets_file.export_to_xml()
 
 
 if __name__ == '__main__':
-    harness = TestHarness('statepoint.10.*')
+    harness = ResonanceScatteringTestHarness('statepoint.10.*')
     harness.main()