Merge pull request #341 from smharper/no_surfaces

Ability to model cells with no surfaces

docs/source/usersguide/input.rst

@@ -757,7 +757,10 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
     is on the negative side of surface 3 and the positive side of surface 5, the
     bounding surfaces would be given as "-3 5".
 
-    *Default*: None
+    .. note:: The surface attribute/element can be omitted to make a cell fill
+              its entire universe.
+
+    *Default*: No surfaces
 
   :rotation:
     If the cell is filled with a universe, this element specifies the angles in

src/geometry.F90

@@ -59,7 +59,8 @@ contains
       end if
     end do SURFACE_LOOP
 
-    ! If we've reached here, then the sense matched on every surface
+    ! If we've reached here, then the sense matched on every surface or there
+    ! are no surfaces.
     in_cell = .true.
 
   end function simple_cell_contains

src/input_xml.F90

@@ -1011,17 +1011,18 @@ contains
         call fatal_error("Cannot specify material and fill simultaneously")
       end if
 
-      ! Check to make sure that surfaces were specified
-      if (.not. check_for_node(node_cell, "surfaces")) then
-        call fatal_error("No surfaces specified for cell " &
-             &// trim(to_str(c % id)))
-      end if
-
       ! Allocate array for surfaces and copy
-      n = get_arraysize_integer(node_cell, "surfaces")
+      if (check_for_node(node_cell, "surfaces")) then
+        n = get_arraysize_integer(node_cell, "surfaces")
+      else
+        n = 0
+      end if
       c % n_surfaces = n
-      allocate(c % surfaces(n))
-      call get_node_array(node_cell, "surfaces", c % surfaces)
+
+      if (n > 0) then
+        allocate(c % surfaces(n))
+        call get_node_array(node_cell, "surfaces", c % surfaces)
+      end if
 
       ! Rotation matrix
       if (check_for_node(node_cell, "rotation")) then

src/relaxng/geometry.rnc

@@ -7,7 +7,7 @@ element geometry {
       (element material { ( xsd:int | "void" ) } | 
         attribute material { ( xsd:int | "void" ) })
     ) &
-    (element surfaces { list { xsd:int+ } } | attribute surfaces { list { xsd:int+ } }) &
+    (element surfaces { list { xsd:int* } } | attribute surfaces { list { xsd:int* } })? &
     (element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? &
     (element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })?
   }*
@@ -23,8 +23,8 @@ element geometry {
 
   & element lattice {
     (element id { xsd:int } | attribute id { xsd:int }) &
-    (element type { ( "rectangular" | "hexagonal" ) } | 
-     attribute type { ( "rectangular" | "hexagonal" ) })? &
+    (element type { ( "rect" | "rectangle" | "rectangular" | "hexagonal" ) } | 
+     attribute type { ( "rect" | "rectangle" | "rectangular" | "hexagonal" ) })? &
     (element dimension { list { xsd:positiveInteger+ } } | 
       attribute dimension { list { xsd:positiveInteger+ } }) &
     (element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) &

tests/test_infinite_cell/geometry.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0"?>
+<geometry>
+  <cell id="11" universe="11" material="1"/>
+  <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">
+    <universes>
+      11 12
+      12 11
+    </universes>
+  </lattice>
+
+
+  <surface id="101" type="z-cylinder" coeffs="0.0 0.0 5.0" boundary="vacuum"/>
+  <cell id="101" universe="0" fill="21" surfaces="-101"/>
+</geometry>

tests/test_infinite_cell/materials.xml

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

tests/test_infinite_cell/results.py

@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+
+import sys
+
+# import statepoint
+sys.path.insert(0, '../../src/utils')
+import statepoint
+
+# read in statepoint file
+if len(sys.argv) > 1:
+    sp = statepoint.StatePoint(sys.argv[1])
+else:
+    sp = statepoint.StatePoint('statepoint.10.binary')
+sp.read_results()
+
+# set up output string
+outstr = ''
+ 
+# write out k-combined
+outstr += 'k-combined:\n'
+outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
+
+# write results to file
+with open('results_test.dat','w') as fh:
+    fh.write(outstr)

tests/test_infinite_cell/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+9.998895E-02 2.846817E-04

tests/test_infinite_cell/settings.xml

@@ -0,0 +1,16 @@
+<?xml version="1.0"?>
+<settings>
+
+  <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_infinite_cell/test_infinite_cell.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+
+import os
+from subprocess import Popen, STDOUT, PIPE, call
+import filecmp
+import glob
+from optparse import OptionParser
+
+parser = OptionParser()
+parser.add_option('--mpi_exec', dest='mpi_exec', default='')
+parser.add_option('--mpi_np', dest='mpi_np', default='3')
+parser.add_option('--exe', dest='exe')
+(opts, args) = parser.parse_args()
+cwd = os.getcwd()
+
+def test_run():
+    if opts.mpi_exec != '':
+        proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0, 'OpenMC did not exit successfully.'
+
+def test_created_statepoint():
+    statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
+    assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
+        'Statepoint file is not a binary or hdf5 file.'
+
+def test_results():
+    statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
+    call(['python', 'results.py', statepoint[0]])
+    compare = filecmp.cmp('results_test.dat', 'results_true.dat')
+    if not compare:
+      os.rename('results_test.dat', 'results_error.dat')
+    assert compare, 'Results do not agree.'
+
+def teardown():
+    output = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
+    output.append(os.path.join(cwd, 'results_test.dat'))
+    for f in output:
+        if os.path.exists(f):
+            os.remove(f)
+
+if __name__ == '__main__':
+
+    # test for openmc executable
+    if opts.exe is None:
+        raise Exception('Must specify OpenMC executable from command line with --exe.')
+
+    # run tests
+    try:
+        test_run()
+        test_created_statepoint()
+        test_results()
+    finally:
+        teardown()