Merge pull request #359 from smharper/lattice_outside

Replace lattice outside materials with outer universes

docs/source/usersguide/input.rst

@@ -833,11 +833,13 @@ sub-elements:
 
     *Default*: None
 
-  :outside:
-    The unique integer identifier of a material that is to be used to fill all
-    space outside of the lattice. This element is optional.
+  :outer:
+    The unique integer identifier of a universe that will be used to fill all
+    space outside of the lattice.  The universe will be tiled repeatedly as if
+    it were placed in a lattice of infinite size.  This element is optional.
 
-    *Default*: The region outside the defined lattice is treated as void.
+    *Default*: An error will be raised if a particle leaves a lattice with no
+    outer universe.
 
   :universes:
     A list of the universe numbers that fill each cell of the lattice.

src/CMakeLists.txt

@@ -269,6 +269,9 @@ install(PROGRAMS utils/statepoint_histogram.py
 install(PROGRAMS utils/statepoint_meshplot.py
   DESTINATION bin
   RENAME statepoint_meshplot)
+install(PROGRAMS utils/update_inputs.py
+  DESTINATION bin
+  RENAME update_inputs)
 install(FILES ../man/man1/openmc.1 DESTINATION share/man/man1)
 install(FILES ../LICENSE DESTINATION "share/doc/${program}/copyright")
 

src/constants.F90

@@ -127,6 +127,9 @@ module constants
        SURF_CONE_Y = 10, & ! Cone parallel to y-axis
        SURF_CONE_Z = 11    ! Cone parallel to z-axis
 
+  ! Flag to say that the outside of a lattice is not defined
+  integer, parameter :: NO_OUTER_UNIVERSE = -22
+
   ! Maximum number of lost particles
   integer, parameter :: MAX_LOST_PARTICLES = 10
 

src/geometry.F90

@@ -280,17 +280,7 @@ contains
               p % coord % xyz = xyz
 
             else
-
-              ! We're outside the lattice, so treat this as a normal cell with
-              ! the material specified for the outside
-
               outside_lattice = .true.
-              p % last_material = p % material
-              p % material = c % material
-
-              ! We'll still make a new coordinate for the particle, as
-              ! distance_to_boundary will still need to track through lattice
-              ! widths even though there's nothing in them but this material
 
             end if
 
@@ -322,13 +312,13 @@ contains
             if (.not. outside_lattice) then
               p % coord % next % universe = lat % universes(i_x,i_y,i_z)
             else
-
-              ! Set universe as the same for subsequent calls to find_cell
-              p % coord % next % universe = p % coord % universe
-
-              ! Set coord cell for calls to distance_to_boundary
-              p % coord % next % cell = index_cell
-
+              if (lat % outer == NO_OUTER_UNIVERSE) then
+                call fatal_error("A particle is outside latttice " &
+                     &// trim(to_str(lat % id)) // " but the lattice has no &
+                     &defined outer universe.")
+              else
+                p % coord % next % universe = lat % outer
+              end if
             end if
 
             ! Move particle to next level
@@ -336,10 +326,9 @@ contains
 
           end if
 
-          if (.not. outside_lattice) then
-            call find_cell(p, found)
-            if (.not. found) exit
-          end if
+          ! Find in the next lowest coordinate level.
+          call find_cell(p, found)
+          if (.not. found) exit
 
         end if
 

src/geometry_header.F90

@@ -30,7 +30,7 @@ module geometry_header
      real(8), allocatable :: lower_left(:)    ! lower-left corner of lattice
      real(8), allocatable :: width(:)         ! width of each lattice cell
      integer, allocatable :: universes(:,:,:) ! specified universes
-     integer              :: outside          ! material to fill area outside
+     integer              :: outer            ! universe to tile outside the lat
   end type Lattice
 
 !===============================================================================

src/initialize.F90

@@ -549,7 +549,7 @@ contains
     integer :: j             ! index for various purposes
     integer :: k             ! loop index for lattices
     integer :: m             ! loop index for lattices
-    integer :: mid, lid      ! material and lattice IDs
+    integer :: lid           ! lattice IDs
     integer :: n_x, n_y, n_z ! size of lattice
     integer :: i_array       ! index in surfaces/materials array
     integer :: id            ! user-specified id
@@ -607,18 +607,8 @@ contains
           c % fill = universe_dict % get_key(id)
         elseif (lattice_dict % has_key(id)) then
           lid = lattice_dict % get_key(id)
-          mid = lattices(lid) % outside
           c % type = CELL_LATTICE
           c % fill = lid
-          if (mid == MATERIAL_VOID) then
-            c % material = mid
-          else if (material_dict % has_key(mid)) then
-            c % material = material_dict % get_key(mid)
-          else
-            call fatal_error("Could not find material " // trim(to_str(mid)) &
-                 &// " specified on lattice " // trim(to_str(lid)))
-          end if
-
         else
           call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "&
                &// trim(to_str(c % id)) // " is neither a universe nor a &
@@ -654,6 +644,16 @@ contains
         end do
       end do
 
+      if (lat % outer /= NO_OUTER_UNIVERSE) then
+        if (universe_dict % has_key(lat % outer)) then
+          lat % outer = universe_dict % get_key(lat % outer)
+        else
+          call fatal_error("Invalid universe number " &
+               &// trim(to_str(lat % outer)) &
+               &// " specified on lattice " // trim(to_str(lat % id)))
+        end if
+      end if
+
     end do
 
     TALLY_LOOP: do i = 1, n_tallies

src/input_xml.F90

@@ -904,7 +904,6 @@ contains
     integer :: universe_num
     integer :: n_cells_in_univ
     integer :: coeffs_reqd
-    integer :: mid
     integer :: temp_int_array3(3)
     integer, allocatable :: temp_int_array(:)
     real(8) :: phi, theta, psi
@@ -1340,15 +1339,18 @@ contains
       end do
       deallocate(temp_int_array)
 
-      ! Read material for area outside lattice
-      lat % outside = MATERIAL_VOID
+      ! Read outer universe for area outside lattice.
+      lat % outer = NO_OUTER_UNIVERSE
+      if (check_for_node(node_lat, "outer")) then
+        call get_node_value(node_lat, "outer", lat % outer)
+      end if
+
+      ! Check for 'outside' nodes which are no longer supported.
       if (check_for_node(node_lat, "outside")) then
-        call get_node_value(node_lat, "outside", mid)
-        if (mid == 0 .or. mid == MATERIAL_VOID) then
-          lat % outside = MATERIAL_VOID
-        else
-          lat % outside = mid
-        end if
+        call fatal_error("The use of 'outside' in lattices is no longer &
+             &supported.  Instead, use 'outer' which defines a universe rather &
+             &than a material.  The utility openmc/src/utils/update_inputs.py &
+             &can be used automatically replace 'outside' with 'outer'.")
       end if
 
       ! Add lattice to dictionary

src/relaxng/geometry.rnc

@@ -30,6 +30,6 @@ element geometry {
     (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) &
     (element width { list { xsd:double+ } } | attribute width { list { xsd:double+ } }) &
     (element universes { list { xsd:int+ } } | attribute universes { list { xsd:int+ } }) &
-    (element outside { xsd:int } | attribute outside { xsd:int })?
+    (element outer { xsd:int } | attribute outer { xsd:int })?
   }*
 }

src/utils/update_inputs.py

@@ -0,0 +1,232 @@
+#!/usr/bin/env python
+"""Update OpenMC's input XML files to the latest format.
+
+Usage information can be obtained by running 'update_inputs.py --help':
+
+usage: update_lattices.py [-h] IN [IN ...]
+
+Update lattices in geometry.xml files to the latest format. This will remove
+'outside' attributes/elements and replace them with 'outer' attributes. Note
+that this script will not delete the given files; it will append '.original'
+to the given files and write new ones.
+
+positional arguments:
+  IN          Input geometry.xml file(s).
+
+optional arguments:
+  -h, --help  show this help message and exit
+
+"""
+
+from __future__ import print_function
+
+import argparse
+from random import randint
+from shutil import move
+import xml.etree.ElementTree as ET
+
+
+description = "Update OpenMC's input XML files to the latest format."
+epilog = """\
+If any of the given files do not match the most up-to-date formatting, then they
+will be automatically rewritten.  The old out-of-date files will not be deleted;
+they will be moved to a new file with '.original' appended to their name.
+
+Formatting changes that will be made:
+
+geometry.xml:  Lattices containing 'outside' attributes/tags will be replaced
+  with lattices containing 'outer' attributes, and the appropriate
+  cells/universes will be added.
+"""
+
+
+def parse_args():
+    """Read the input files from the commandline."""
+    # Create argument parser.
+    parser = argparse.ArgumentParser(
+         description=description,
+         epilog=epilog,
+         formatter_class=argparse.RawTextHelpFormatter)
+    parser.add_argument('input', metavar='IN', type=str, nargs='+',
+                        help='Input geometry.xml file(s).')
+
+    # Parse and return commandline arguments.
+    return parser.parse_args()
+
+
+def get_universe_ids(geometry_root):
+    """Return a set of universe id numbers."""
+    root = geometry_root
+    out = {0}
+
+    # Get the ids of universes defined by cells.
+    for cell in root.iter('cell'):
+        # Get universe attributes.
+        if 'universe' in cell.attrib:
+            uid = cell.attrib['universe']
+            out.add(int(uid))
+
+        # Get universe elements.
+        elif cell.find('universe') is not None:
+            elem = cell.find('universe')
+            uid = elem.text
+            out.add(int(uid))
+
+    # Get the ids of universes defined by lattices.
+    for lat in root.iter('lattice'):
+        # Get id attributes.
+        if 'id' in lat.attrib:
+            uid = lat.attrib['id']
+            out.add(int(uid))
+
+        # Get id elements.
+        elif lat.find('id') is not None:
+            elem = lat.find('id')
+            uid = elem.text
+            out.add(int(uid))
+
+    return out
+
+
+def get_cell_ids(geometry_root):
+    """Return a set of cell id numbers."""
+    root = geometry_root
+    out = set()
+
+    # Get the ids of universes defined by cells.
+    for cell in root.iter('cell'):
+        # Get id attributes.
+        if 'id' in cell.attrib:
+            cid = cell.attrib['id']
+            out.add(int(cid))
+
+        # Get id elements.
+        elif cell.find('id') is not None:
+            elem = cell.find('id')
+            cid = elem.text
+            out.add(int(cid))
+
+    return out
+
+
+def find_new_id(current_ids, preferred=None):
+    """Return a new id that is not already present in current_ids."""
+    distance_from_preferred = 21
+    max_random_attempts = 10000
+
+    # First, try to find an id near the preferred number.
+    if preferred is not None:
+        assert isinstance(preferred, int)
+        for i in range(1, distance_from_preferred):
+            if (preferred - i not in current_ids) and (preferred - i > 0):
+                return preferred - i
+            if (preferred + i not in current_ids) and (preferred + i > 0):
+                return preferred + i
+
+    # If that was unsuccessful, attempt to randomly guess a new id number.
+    for i in range(max_random_attempts):
+        num = randint(1, 2147483647)
+        if num not in current_inds:
+            return num
+
+    # Raise an error if an id was not found.
+    raise RuntimeError('Could not find a unique id number for a new universe.')
+
+
+def get_lat_id(lattice_element):
+    """Return the id integer of the lattice_element."""
+    assert isinstance(lattice_element, ET.Element)
+    if 'id' in lattice_element.attrib:
+        return int(lattice_element.attrib['id'].strip())
+    elif any([child.tag == 'id' for child in lattice_element]):
+        elem = lattice_element.find('id')
+        return int(elem.text.strip())
+    else:
+        raise RuntimeError('Could not find the id for a lattice.')
+
+
+def pop_lat_outside(lattice_element):
+    """Return lattice's outside material and remove from attributes/elements."""
+    assert isinstance(lattice_element, ET.Element)
+
+    # Check attributes.
+    if 'outside' in lattice_element.attrib:
+        material = lattice_element.attrib['outside'].strip()
+        del lattice_element.attrib['outside']
+
+    # Check subelements.
+    elif any([child.tag == 'outside' for child in lattice_element]):
+        elem = lattice_element.find('outside')
+        material = elem.text.strip()
+        lattice_element.remove(elem)
+
+    # No 'outside' specified.  This means the outside is a void.
+    else:
+        material = 'void'
+
+    return material
+
+
+def update_geometry(geometry_root):
+    """Update the given XML geometry tree. Return True if changes were made."""
+    root = geometry_root
+    was_updated = False
+
+    # Ignore files that do not contain lattices.
+    if all([child.tag != 'lattice' for child in root]): return False
+
+    # Get a set of already-used universe and cell ids.
+    uids = get_universe_ids(root)
+    cids = get_cell_ids(root)
+    taken_ids = uids.union(cids)
+
+    # Update the definitions of each lattice
+    for lat in root.iter('lattice'):
+        # Get the lattice's id.
+        lat_id = get_lat_id(lat)
+
+        # Ignore lattices that have 'outer' specified.
+        if any([child.tag == 'outer' for child in lat]): continue
+
+        # Pop the 'outside' material.
+        material = pop_lat_outside(lat)
+
+        # Get an id number for a new outer universe.  Ideally, the id should
+        # be close to the lattice's id.
+        new_uid = find_new_id(taken_ids, preferred=lat_id)
+        assert new_uid not in taken_ids
+
+        # Add the new universe filled with the old 'outside' material to the
+        # geometry.
+        new_cell = ET.Element('cell')
+        new_cell.attrib['id'] = str(new_uid)
+        new_cell.attrib['universe'] = str(new_uid)
+        new_cell.attrib['material'] = material
+        root.append(new_cell)
+        taken_ids.add(new_uid)
+
+        # Add the new universe to the lattice's 'outer' attribute.
+        lat.attrib['outer'] = str(new_uid)
+
+        was_updated = True
+
+    return was_updated
+
+
+if __name__ == '__main__':
+    args = parse_args()
+    for fname in args.input:
+        # Parse the XML data.
+        tree = ET.parse(fname)
+        root = tree.getroot()
+        was_updated = False
+
+        if root.tag == 'geometry':
+            was_updated = update_geometry(root)
+
+        if was_updated:
+            # Move the original geometry file to preserve it.
+            move(fname, fname + '.original')
+
+            # Write a new geometry file.
+            tree.write(fname)

tests/test_infinite_cell/geometry.xml

@@ -4,12 +4,13 @@
   <cell id="12" universe="12" material="2" surfaces=""/>
 
   <lattice id="21" type="rect" dimension="2 2" lower_left="-2.0 -2.0"
-    width="2.0 2.0" outside="2">
+    width="2.0 2.0" outer="22">
     <universes>
       11 12
       12 11
     </universes>
   </lattice>
+  <cell id="22" universe="22" material="2"/>
 
 
   <surface id="101" type="z-cylinder" coeffs="0.0 0.0 5.0" boundary="vacuum"/>