Merge remote-tracking branch 'mitcrpgopenmc/develop' into res_scat

src/tally.F90

@@ -498,7 +498,9 @@ contains
       filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
       if (t % all_nuclides) then
-        call score_all_nuclides(p, i_tally, flux, filter_index)
+        if (p % material /= MATERIAL_VOID) then
+          call score_all_nuclides(p, i_tally, flux, filter_index)
+        end if
       else
 
         NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins
@@ -506,23 +508,27 @@ contains
           i_nuclide = t % nuclide_bins(k)
 
           if (i_nuclide > 0) then
-            ! Get pointer to current material
-            mat => materials(p % material)
+            if (p % material /= MATERIAL_VOID) then
+              ! Get pointer to current material
+              mat => materials(p % material)
 
-            ! Determine if nuclide is actually in material
-            NUCLIDE_MAT_LOOP: do j = 1, mat % n_nuclides
-              ! If index of nuclide matches the j-th nuclide listed in the
-              ! material, break out of the loop
-              if (i_nuclide == mat % nuclide(j)) exit
+              ! Determine if nuclide is actually in material
+              NUCLIDE_MAT_LOOP: do j = 1, mat % n_nuclides
+                ! If index of nuclide matches the j-th nuclide listed in the
+                ! material, break out of the loop
+                if (i_nuclide == mat % nuclide(j)) exit
 
-              ! If we've reached the last nuclide in the material, it means
-              ! the specified nuclide to be tallied is not in this material
-              if (j == mat % n_nuclides) then
-                cycle NUCLIDE_BIN_LOOP
-              end if
-            end do NUCLIDE_MAT_LOOP
+                ! If we've reached the last nuclide in the material, it means
+                ! the specified nuclide to be tallied is not in this material
+                if (j == mat % n_nuclides) then
+                  cycle NUCLIDE_BIN_LOOP
+                end if
+              end do NUCLIDE_MAT_LOOP
 
-            atom_density = mat % atom_density(j)
+              atom_density = mat % atom_density(j)
+            else
+              atom_density = ZERO
+            end if
           end if
 
           ! Determine score for each bin
@@ -1167,33 +1173,38 @@ contains
         filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
 
         if (t % all_nuclides) then
-          ! Score reaction rates for each nuclide in material
-          call score_all_nuclides(p, i_tally, flux, filter_index)
-
+          if (p % material /= MATERIAL_VOID) then
+            ! Score reaction rates for each nuclide in material
+            call score_all_nuclides(p, i_tally, flux, filter_index)
+          end if
         else
           NUCLIDE_BIN_LOOP: do b = 1, t % n_nuclide_bins
             ! Get index of nuclide in nuclides array
             i_nuclide = t % nuclide_bins(b)
 
             if (i_nuclide > 0) then
-              ! Get pointer to current material
-              mat => materials(p % material)
+              if (p % material /= MATERIAL_VOID) then
+                ! Get pointer to current material
+                mat => materials(p % material)
 
-              ! Determine if nuclide is actually in material
-              NUCLIDE_MAT_LOOP: do j = 1, mat % n_nuclides
-                ! If index of nuclide matches the j-th nuclide listed in
-                ! the material, break out of the loop
-                if (i_nuclide == mat % nuclide(j)) exit
+                ! Determine if nuclide is actually in material
+                NUCLIDE_MAT_LOOP: do j = 1, mat % n_nuclides
+                  ! If index of nuclide matches the j-th nuclide listed in
+                  ! the material, break out of the loop
+                  if (i_nuclide == mat % nuclide(j)) exit
 
-                ! If we've reached the last nuclide in the material, it
-                ! means the specified nuclide to be tallied is not in this
-                ! material
-                if (j == mat % n_nuclides) then
-                  cycle NUCLIDE_BIN_LOOP
-                end if
-              end do NUCLIDE_MAT_LOOP
+                  ! If we've reached the last nuclide in the material, it
+                  ! means the specified nuclide to be tallied is not in this
+                  ! material
+                  if (j == mat % n_nuclides) then
+                    cycle NUCLIDE_BIN_LOOP
+                  end if
+                end do NUCLIDE_MAT_LOOP
 
-              atom_density = mat % atom_density(j)
+                atom_density = mat % atom_density(j)
+              else
+                atom_density = ZERO
+              end if
             end if
 
             ! Determine score for each bin

tests/test_tally_nuclides/geometry.xml

@@ -0,0 +1,10 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <surface id="1"  type="sphere" coeffs="0. 0. 0. 1.4" />
+  <surface id="2"  type="sphere" coeffs="0. 0. 0. 5.35" boundary="vacuum" />
+
+  <cell id="1"  material="void" surfaces="-1" />      <!-- cavity -->
+  <cell id="2"  material="1"    surfaces=" 1  -2" />  <!-- Pu Core -->
+
+</geometry>

tests/test_tally_nuclides/materials.xml

@@ -0,0 +1,13 @@
+<?xml version="1.0"?>
+<materials>
+
+  <default_xs>71c</default_xs>
+
+  <!-- Plutonium -->
+  <material id="1">
+    <density value="18.4757" units="g/cm3" />
+    <nuclide name="Pu-239" ao="1.0" />
+  </material>
+
+
+</materials>

tests/test_tally_nuclides/results.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+
+import sys
+import numpy as np
+
+# import statepoint
+sys.path.append('../../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()
+
+# extract tally results and convert to vector
+results0 = sp.tallies[0].results
+results1 = sp.tallies[1].results
+results = np.concatenate((results0.flatten(), results1.flatten()))
+
+# 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 out tally results
+outstr += 'tallies:\n'
+for item in results:
+  outstr += "{0:12.6E}\n".format(item)
+
+# write results to file
+with open('results_test.dat','w') as fh:
+    fh.write(outstr)

tests/test_tally_nuclides/results_true.dat

@@ -0,0 +1,27 @@
+k-combined:
+9.914201E-01 3.310472E-02
+tallies:
+7.251620E+00
+1.056750E+01
+1.649542E+00
+5.462124E-01
+1.598276E+00
+5.127366E-01
+5.602079E+00
+6.310407E+00
+7.251620E+00
+1.056750E+01
+1.649542E+00
+5.462124E-01
+1.598276E+00
+5.127366E-01
+5.602079E+00
+6.310407E+00
+7.251620E+00
+1.056750E+01
+1.649542E+00
+5.462124E-01
+1.598276E+00
+5.127366E-01
+5.602079E+00
+6.310407E+00

tests/test_tally_nuclides/settings.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0"?>
+<settings>
+
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>100</particles>
+  </eigenvalue>
+
+  <source>
+    <space>
+      <type>point</type>
+      <parameters>0.0 0.0 0.0</parameters>
+    </space>
+  </source>
+
+</settings>

tests/test_tally_nuclides/tallies.xml

@@ -0,0 +1,14 @@
+<?xml version="1.0"?>
+<tallies>
+
+  <tally id="1">
+    <nuclides>all</nuclides>
+    <scores>total absorption fission scatter</scores>
+  </tally>
+
+  <tally id="2">
+    <nuclides>Pu-239</nuclides>
+    <scores>total absorption fission scatter</scores>
+  </tally>
+
+</tallies>

tests/test_tally_nuclides/test_tally_nuclides.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()