Merge pull request #243 from mit-crpg/read_sourcefile

Ability to initialize source based on source point file.

docs/source/usersguide/input.rst

@@ -264,7 +264,9 @@ attributes/sub-elements:
 
   :file:
     If this attribute is given, it indicates that the source is to be read from
-    a binary source file whose path is given by the value of this element
+    a binary source file whose path is given by the value of this element. Note,
+    the number of source sites needs to be the same as the number of particles
+    simulated in a fission source generation.
 
     *Default*: None
 

src/source.F90

@@ -1,16 +1,17 @@
 module source
 
-  use bank_header,     only: Bank
+  use bank_header,      only: Bank
   use constants
-  use error,           only: fatal_error
-  use geometry,        only: find_cell
-  use geometry_header, only: BASE_UNIVERSE
+  use error,            only: fatal_error
+  use geometry,         only: find_cell
+  use geometry_header,  only: BASE_UNIVERSE
   use global
-  use math,            only: maxwell_spectrum, watt_spectrum
-  use output,          only: write_message
-  use particle_header, only: Particle
-  use random_lcg,      only: prn, set_particle_seed
-  use string,          only: to_str
+  use math,             only: maxwell_spectrum, watt_spectrum
+  use output,           only: write_message
+  use output_interface, only: BinaryOutput
+  use particle_header,  only: Particle
+  use random_lcg,       only: prn, set_particle_seed
+  use string,           only: to_str
 
 #ifdef MPI
   use mpi
@@ -28,8 +29,9 @@ contains
 
     integer(8) :: i          ! loop index over bank sites
     integer(8) :: id         ! particle id
-
+    integer(4) :: itmp       ! temporary integer
     type(Bank), pointer :: src => null() ! source bank site
+    type(BinaryOutput) :: sp ! statepoint/source binary file
 
     message = "Initializing source particles..."
     call write_message(6)
@@ -38,8 +40,26 @@ contains
       ! Read the source from a binary file instead of sampling from some
       ! assumed source distribution
 
-      message = 'This feature is currently disabled and will be added back in.'
-      call fatal_error()
+      message = 'Reading source file from ' // trim(path_source) // '...'
+      call write_message(6)
+
+      ! Open the binary file
+      call sp % file_open(path_source, 'r', serial = .false.)
+
+      ! Read the file type
+      call sp % read_data(itmp, "filetype")
+
+      ! Check to make sure this is a source file
+      if (itmp /= FILETYPE_SOURCE) then
+        message = "Specified starting source file not a source file type."
+        call fatal_error()
+      end if
+
+      ! Read in the source bank
+      call sp % read_source_bank()
+
+      ! Close file
+      call sp % file_close()
 
     else
       ! Generation source sites from specified distribution in user input

tests/test_source_file/geometry.xml

@@ -0,0 +1,8 @@
+<?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" surfaces="-1" />
+    
+</geometry>

tests/test_source_file/materials.xml

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

tests/test_source_file/results.py

@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+
+import sys
+
+# 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()
+
+# 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_source_file/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+2.977307E-01 2.848670E-03

tests/test_source_file/settings.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0"?>
+<settings>
+
+  <state_point batches="10" />
+  <source_point separate="true" />
+
+  <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_source_file/test_source_file.py

@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+
+import os
+from subprocess import Popen, STDOUT, PIPE, call
+import filecmp
+from nose_mpi import NoseMPI
+import glob
+
+pwd = os.path.dirname(__file__)
+
+settings1="""<?xml version="1.0"?>
+<settings>
+
+  <state_point batches="10" />
+  <source_point separate="true" />
+
+  <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>
+"""
+
+settings2 = """<?xml version="1.0"?>
+<settings>
+
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+
+  <source>
+    <file> source.10.{0} </file>
+  </source>
+
+</settings>
+"""
+
+def setup(): 
+    os.putenv('PWD', pwd)
+    os.chdir(pwd)
+
+def test_run1():
+    openmc_path = pwd + '/../../src/openmc'
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_statepoint_exists():
+    statepoint = glob.glob(pwd + '/statepoint.10.*')
+    assert len(statepoint) == 1
+    assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5')
+    source = glob.glob(pwd + '/source.10.*')
+    assert len(statepoint) == 1
+    assert source[0].endswith('binary') or source[0].endswith('h5')
+
+def test_run2():
+    openmc_path = pwd + '/../../src/openmc'
+    source = glob.glob(pwd + '/source.10.*')
+    with open('settings.xml','w') as fh:
+        fh.write(settings2.format(source[0].split('.')[2]))
+    if int(NoseMPI.mpi_np) > 0:
+        proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path],
+               stderr=STDOUT, stdout=PIPE)
+    else:
+        proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE)
+    print(proc.communicate()[0])
+    returncode = proc.returncode
+    assert returncode == 0
+
+def test_results():
+    statepoint = glob.glob(pwd + '/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
+
+def teardown():
+    with open('settings.xml','w') as fh:
+        fh.write(settings1)
+    output = glob.glob(pwd + '/statepoint.10.*')
+    output += glob.glob(pwd + '/source.10.*')
+    output.append(pwd + '/results_test.dat')
+    for f in output:
+        if os.path.exists(f):
+            os.remove(f)