diff --git a/src/tally.F90 b/src/tally.F90 index b731c1ae5e..466f94b1c0 100644 --- a/src/tally.F90 +++ b/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 diff --git a/tests/test_tally_nuclides/geometry.xml b/tests/test_tally_nuclides/geometry.xml new file mode 100644 index 0000000000..65954a7853 --- /dev/null +++ b/tests/test_tally_nuclides/geometry.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/test_tally_nuclides/materials.xml b/tests/test_tally_nuclides/materials.xml new file mode 100644 index 0000000000..2761be30cf --- /dev/null +++ b/tests/test_tally_nuclides/materials.xml @@ -0,0 +1,13 @@ + + + + 71c + + + + + + + + + diff --git a/tests/test_tally_nuclides/results.py b/tests/test_tally_nuclides/results.py new file mode 100644 index 0000000000..3aa0946a5d --- /dev/null +++ b/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) diff --git a/tests/test_tally_nuclides/results_true.dat b/tests/test_tally_nuclides/results_true.dat new file mode 100644 index 0000000000..e429381ddb --- /dev/null +++ b/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 diff --git a/tests/test_tally_nuclides/settings.xml b/tests/test_tally_nuclides/settings.xml new file mode 100644 index 0000000000..b2ddb42483 --- /dev/null +++ b/tests/test_tally_nuclides/settings.xml @@ -0,0 +1,17 @@ + + + + + 10 + 5 + 100 + + + + + point + 0.0 0.0 0.0 + + + + diff --git a/tests/test_tally_nuclides/tallies.xml b/tests/test_tally_nuclides/tallies.xml new file mode 100644 index 0000000000..cf20668c83 --- /dev/null +++ b/tests/test_tally_nuclides/tallies.xml @@ -0,0 +1,14 @@ + + + + + all + total absorption fission scatter + + + + Pu-239 + total absorption fission scatter + + + diff --git a/tests/test_tally_nuclides/test_tally_nuclides.py b/tests/test_tally_nuclides/test_tally_nuclides.py new file mode 100644 index 0000000000..6fdbf87459 --- /dev/null +++ b/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()