From fc2681d65fcb79849a0f4cb38e49d10ce0dd30c8 Mon Sep 17 00:00:00 2001 From: "Jonathan A. Walsh" Date: Thu, 12 Jun 2014 15:30:38 -0600 Subject: [PATCH 01/18] don't resample ptable xs if diff temperature but same E --- src/ace_header.F90 | 1 + src/cross_section.F90 | 28 ++++++++++++++++++++++------ src/geometry.F90 | 1 + src/material_header.F90 | 1 + 4 files changed, 25 insertions(+), 6 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index e5bf1bb3a4..532a355d3f 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -241,6 +241,7 @@ module ace_header ! Information for URR probability table use logical :: use_ptable ! in URR range with probability tables? + real(8) :: last_prn end type NuclideMicroXS !=============================================================================== diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 7426519cda..80995da861 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -92,9 +92,9 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, mat % first_E) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, mat % first_E) end if ! ======================================================================== @@ -135,11 +135,12 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, first_E) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy + real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid @@ -184,6 +185,8 @@ contains ! Initialize sab treatment to false micro_xs(i_nuclide) % index_sab = NONE micro_xs(i_nuclide) % elastic_sab = ZERO + + ! Initialize URR probability table treatment to false micro_xs(i_nuclide) % use_ptable = .false. ! Initialize nuclide cross-sections to zero @@ -233,7 +236,7 @@ contains if (urr_ptables_on .and. nuc % urr_present) then if (E > nuc % urr_data % energy(1) .and. & E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then - call calculate_urr_xs(i_nuclide, E) + call calculate_urr_xs(i_nuclide, E, first_E) end if end if @@ -334,10 +337,11 @@ contains ! from probability tables !=============================================================================== - subroutine calculate_urr_xs(i_nuclide, E) + subroutine calculate_urr_xs(i_nuclide, E, first_E) integer, intent(in) :: i_nuclide ! index into nuclides array real(8), intent(in) :: E ! energy + real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_energy ! index for energy integer :: i_table ! index for table @@ -370,8 +374,20 @@ contains (urr % energy(i_energy + 1) - urr % energy(i_energy)) ! sample probability table using the cumulative distribution - r = prn() + ! (if we're dealing with an isotope that we've previously encountered at + ! this energy but a different temperature, use the original random number to ! preserve correlation of temperature in probability tables) + if (E /= ZERO .and. E == first_E .and. & + & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & + & + abs(micro_xs(:) % last_E - E)) == ZERO) then + r = micro_xs(minloc(dble(abs(nuclides(:) % zaid & + & - nuclides(i_nuclide) % zaid)) + abs(micro_xs(:) % last_E & + & - E),1)) % last_prn + else + r = prn() + end if + i_table = 1 + do if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit i_table = i_table + 1 diff --git a/src/geometry.F90 b/src/geometry.F90 index 24d3ef8b06..1a2544d28d 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -190,6 +190,7 @@ contains ! set material p % last_material = p % material p % material = c % material + materials(p % material) % first_E = p % E elseif (c % type == CELL_FILL) then ! ==================================================================== diff --git a/src/material_header.F90 b/src/material_header.F90 index cbfd917498..d4c7e81290 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -12,6 +12,7 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + real(8) :: first_E ! particle energy upon entering mat ! S(a,b) data references integer :: n_sab = 0 ! number of S(a,b) tables From ad3519a0efd8a49679e60486f0b2528196130fe9 Mon Sep 17 00:00:00 2001 From: "Jonathan A. Walsh" Date: Thu, 12 Jun 2014 16:48:24 -0600 Subject: [PATCH 02/18] set last prn of a nuclide whenever a prn is needed --- src/cross_section.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 80995da861..2f7bcf9e3a 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -384,6 +384,7 @@ contains & - E),1)) % last_prn else r = prn() + micro_xs(i_nuclide) % last_prn = r end if i_table = 1 From 3bd733bd89f88157870aaba94d4ea3b81ef5d47f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 10 Jul 2014 18:24:53 +0200 Subject: [PATCH 03/18] Added zero-padded statepoint names --- src/state_point.F90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/state_point.F90 b/src/state_point.F90 index 83d768a1ba..e3ce3f9972 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -36,15 +36,26 @@ contains subroutine write_state_point() + character(MAX_FILE_LEN) :: batch_form + character(MAX_FILE_LEN) :: batch_string character(MAX_FILE_LEN) :: filename + integer :: n_digits integer :: i integer :: j integer, allocatable :: temp_array(:) type(TallyObject), pointer :: t => null() + ! Count digits needed to represent largest batch number (for file name). + n_digits = 1 + do while (n_batches / 10**(n_digits) /= 0 & + & .and. n_batches / 10 **(n_digits-1) /= 1) + n_digits = n_digits + 1 + end do + ! Convert the current batch batch to a zero-padded string. + write(batch_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits + write(batch_string, batch_form) current_batch ! Set filename for state point - filename = trim(path_output) // 'statepoint.' // & - trim(to_str(current_batch)) + filename = trim(path_output) // 'statepoint.' // batch_string ! Append appropriate extension #ifdef HDF5 From d1f863b2be48078281e6a56d05b5ae22f1b38093 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 11 Jul 2014 16:04:10 +0200 Subject: [PATCH 04/18] Wrapped batch number zero padding into a function --- src/state_point.F90 | 40 +++++++++++++++++++++++++++------------- 1 file changed, 27 insertions(+), 13 deletions(-) diff --git a/src/state_point.F90 b/src/state_point.F90 index e3ce3f9972..5cb2f5ad7a 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -36,26 +36,14 @@ contains subroutine write_state_point() - character(MAX_FILE_LEN) :: batch_form - character(MAX_FILE_LEN) :: batch_string character(MAX_FILE_LEN) :: filename - integer :: n_digits integer :: i integer :: j integer, allocatable :: temp_array(:) type(TallyObject), pointer :: t => null() - ! Count digits needed to represent largest batch number (for file name). - n_digits = 1 - do while (n_batches / 10**(n_digits) /= 0 & - & .and. n_batches / 10 **(n_digits-1) /= 1) - n_digits = n_digits + 1 - end do - ! Convert the current batch batch to a zero-padded string. - write(batch_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits - write(batch_string, batch_form) current_batch ! Set filename for state point - filename = trim(path_output) // 'statepoint.' // batch_string + filename = trim(path_output) // 'statepoint.' // zero_padded_batch() ! Append appropriate extension #ifdef HDF5 @@ -821,6 +809,32 @@ contains end subroutine load_state_point +!=============================================================================== +! STATE_POINT -- This module handles writing and reading state point +!=============================================================================== + + function zero_padded_batch() + character(MAX_FILE_LEN) :: zero_padded_batch + + character(MAX_FILE_LEN) :: batch_form + character(MAX_FILE_LEN) :: batch_string + character(MAX_FILE_LEN) :: filename + integer :: n_digits + + ! Count digits needed to represent largest batch number. + n_digits = 1 + do while (n_batches / 10**(n_digits) /= 0 & + & .and. n_batches / 10 **(n_digits-1) /= 1) + n_digits = n_digits + 1 + end do + + ! Write a format string of the form In.n where n is n_digits. + write(batch_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits + + ! Format the batch number + write(zero_padded_batch, batch_form) current_batch + end function zero_padded_batch + subroutine read_source ! TODO write this routine ! TODO what if n_particles does not match source bank From 244e9970e216f21820f858c1206c8467e24d50cf Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 11 Jul 2014 16:09:16 +0200 Subject: [PATCH 05/18] Added zero padded filenames for source files Also correcteted statepoint filenames in write_source_point --- src/state_point.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/state_point.F90 b/src/state_point.F90 index 5cb2f5ad7a..cb276e2700 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -304,7 +304,7 @@ contains if (source_separate) then ! Set filename - filename = trim(path_output) // 'source.' // trim(to_str(current_batch)) + filename = trim(path_output) // 'source.' // zero_padded_batch() #ifdef HDF5 filename = trim(filename) // '.h5' #else @@ -324,8 +324,7 @@ contains else ! Set filename for state point - filename = trim(path_output) // 'statepoint.' // & - trim(to_str(current_batch)) + filename = trim(path_output) // 'statepoint.' // zero_padded_batch() #ifdef HDF5 filename = trim(filename) // '.h5' #else From aa2f94ef3485a40968dc78823d6c740cad036552 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 11 Jul 2014 17:29:56 +0200 Subject: [PATCH 06/18] Updated tests for zero padded filenames --- tests/test_sourcepoint_batch/results.py | 2 +- .../test_sourcepoint_batch.py | 2 +- tests/test_sourcepoint_interval/results.py | 2 +- .../test_sourcepoint_interval.py | 2 +- .../test_sourcepoint_restart.py | 14 +++++------ tests/test_statepoint_batch/results.py | 2 +- .../test_statepoint_batch.py | 14 +++++------ .../test_statepoint_interval.py | 16 ++++++------- tests/test_statepoint_restart/results.py | 2 +- .../test_statepoint_restart.py | 24 +++++++++---------- 10 files changed, 40 insertions(+), 40 deletions(-) diff --git a/tests/test_sourcepoint_batch/results.py b/tests/test_sourcepoint_batch/results.py index 9b0f577eb0..96984b2fae 100644 --- a/tests/test_sourcepoint_batch/results.py +++ b/tests/test_sourcepoint_batch/results.py @@ -10,7 +10,7 @@ import statepoint if len(sys.argv) > 1: sp = statepoint.StatePoint(sys.argv[1]) else: - sp = statepoint.StatePoint('statepoint.8.binary') + sp = statepoint.StatePoint('statepoint.08.binary') sp.read_results() sp.read_source() diff --git a/tests/test_sourcepoint_batch/test_sourcepoint_batch.py b/tests/test_sourcepoint_batch/test_sourcepoint_batch.py index a1c89690b8..4258a54de7 100644 --- a/tests/test_sourcepoint_batch/test_sourcepoint_batch.py +++ b/tests/test_sourcepoint_batch/test_sourcepoint_batch.py @@ -30,7 +30,7 @@ def test_statepoint_exists(): 'Statepoint file detected that is not binary or hdf5.' def test_results(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.8.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.08.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: diff --git a/tests/test_sourcepoint_interval/results.py b/tests/test_sourcepoint_interval/results.py index 9b0f577eb0..96984b2fae 100644 --- a/tests/test_sourcepoint_interval/results.py +++ b/tests/test_sourcepoint_interval/results.py @@ -10,7 +10,7 @@ import statepoint if len(sys.argv) > 1: sp = statepoint.StatePoint(sys.argv[1]) else: - sp = statepoint.StatePoint('statepoint.8.binary') + sp = statepoint.StatePoint('statepoint.08.binary') sp.read_results() sp.read_source() diff --git a/tests/test_sourcepoint_interval/test_sourcepoint_interval.py b/tests/test_sourcepoint_interval/test_sourcepoint_interval.py index a583370c35..2a089747f6 100644 --- a/tests/test_sourcepoint_interval/test_sourcepoint_interval.py +++ b/tests/test_sourcepoint_interval/test_sourcepoint_interval.py @@ -30,7 +30,7 @@ def test_statepoint_exists(): 'Statepoint file detected that is not binary or hdf5.' def test_results(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.8.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.08.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: diff --git a/tests/test_sourcepoint_restart/test_sourcepoint_restart.py b/tests/test_sourcepoint_restart/test_sourcepoint_restart.py index d42763be3e..47033b0121 100644 --- a/tests/test_sourcepoint_restart/test_sourcepoint_restart.py +++ b/tests/test_sourcepoint_restart/test_sourcepoint_restart.py @@ -28,7 +28,7 @@ def test_created_statepoint(): assert len(statepoint) == 2, '2 statepoint files must exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint file must either be binary or hdf5.' - sourcepoint = glob.glob(os.path.join(cwd, 'source.7.*')) + sourcepoint = glob.glob(os.path.join(cwd, 'source.07.*')) assert len(sourcepoint) == 1, 'Either multiple or no source files found.' assert sourcepoint[0].endswith('binary') or sourcepoint[0].endswith('h5'),\ 'Source file must either be binary or hdf5.' @@ -43,8 +43,8 @@ def test_results(): os.remove(statepoint[0]) def test_restart_form1(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) - sourcepoint = glob.glob(os.path.join(cwd, 'source.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) + sourcepoint = glob.glob(os.path.join(cwd, 'source.07.*')) if opts.mpi_exec != '': proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, '-r', statepoint[0], sourcepoint[0], cwd], stderr=STDOUT, stdout=PIPE) @@ -70,8 +70,8 @@ def test_results_form1(): os.remove(statepoint[0]) def test_restart_form2(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) - sourcepoint = glob.glob(os.path.join(cwd, 'source.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) + sourcepoint = glob.glob(os.path.join(cwd, 'source.07.*')) if opts.mpi_exec != '': proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, '--restart', statepoint[0], sourcepoint[0], cwd], stderr=STDOUT, stdout=PIPE) @@ -97,8 +97,8 @@ def test_results_form2(): os.remove(statepoint[0]) def test_restart_serial(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) - sourcepoint = glob.glob(os.path.join(cwd, 'source.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) + sourcepoint = glob.glob(os.path.join(cwd, 'source.07.*')) proc = Popen([opts.exe, '--restart', statepoint[0], sourcepoint[0], cwd], stderr=STDOUT, stdout=PIPE) print(proc.communicate()[0]) returncode = proc.returncode diff --git a/tests/test_statepoint_batch/results.py b/tests/test_statepoint_batch/results.py index 31c0ac3caa..112c005167 100644 --- a/tests/test_statepoint_batch/results.py +++ b/tests/test_statepoint_batch/results.py @@ -10,7 +10,7 @@ import statepoint if len(sys.argv) > 1: sp = statepoint.StatePoint(sys.argv[1]) else: - sp = statepoint.StatePoint('statepoint.9.binary') + sp = statepoint.StatePoint('statepoint.09.binary') sp.read_results() # set up output string diff --git a/tests/test_statepoint_batch/test_statepoint_batch.py b/tests/test_statepoint_batch/test_statepoint_batch.py index a594b2c5a9..e1951d23fd 100644 --- a/tests/test_statepoint_batch/test_statepoint_batch.py +++ b/tests/test_statepoint_batch/test_statepoint_batch.py @@ -24,21 +24,21 @@ def test_run(): assert returncode == 0, 'OpenMC did not exit successfully.' def test_statepoints_exist(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.3.*')) - assert len(statepoint) == 1, 'Either multiple or no statepoint.3 files exist.' + statepoint = glob.glob(os.path.join(cwd, 'statepoint.03.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint.03 files exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint.3 file is not a binary or hdf5 file.' - statepoint = glob.glob(os.path.join(cwd, 'statepoint.6.*')) - assert len(statepoint) == 1, 'Either multiple or no statepoint.6 files exist.' + statepoint = glob.glob(os.path.join(cwd, 'statepoint.06.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint.06 files exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint.6 file is not a binary or hdf5 file.' - statepoint = glob.glob(os.path.join(cwd, 'statepoint.9.*')) - assert len(statepoint) == 1, 'Either multiple or no statepoint.9 files exist.' + statepoint = glob.glob(os.path.join(cwd, 'statepoint.09.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint.09 files exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint.9 file is not a binary or hdf5 file.' def test_results(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.9.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.09.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: diff --git a/tests/test_statepoint_interval/test_statepoint_interval.py b/tests/test_statepoint_interval/test_statepoint_interval.py index b7c30d3384..1099566e58 100644 --- a/tests/test_statepoint_interval/test_statepoint_interval.py +++ b/tests/test_statepoint_interval/test_statepoint_interval.py @@ -25,20 +25,20 @@ def test_run(): assert returncode == 0, 'OpenMC did not exit successfully.' def test_statepoints_exist(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.2.*')) - assert len(statepoint) == 1, 'Either multiple or no statepoint.2 files exist.' + statepoint = glob.glob(os.path.join(cwd, 'statepoint.02.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint.02 files exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint.2 file is not a binary or hdf5 file.' - statepoint = glob.glob(os.path.join(cwd, 'statepoint.4.*')) - assert len(statepoint) == 1, 'Either multiple or no statepoint.4 files exist.' + statepoint = glob.glob(os.path.join(cwd, 'statepoint.04.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint.04 files exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint.4 file is not a binary or hdf5 file.' - statepoint = glob.glob(os.path.join(cwd, 'statepoint.6.*')) - assert len(statepoint) == 1, 'Either multiple or no statepoint.6 files exist.' + statepoint = glob.glob(os.path.join(cwd, 'statepoint.06.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint.06 files exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint.6 file is not a binary or hdf5 file.' - statepoint = glob.glob(os.path.join(cwd, 'statepoint.8.*')) - assert len(statepoint) == 1, 'Either multiple or no statepoint.8 files exist.' + statepoint = glob.glob(os.path.join(cwd, 'statepoint.08.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint.08 files exist.' assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ 'Statepoint.8 file is not a binary or hdf5 file.' statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) diff --git a/tests/test_statepoint_restart/results.py b/tests/test_statepoint_restart/results.py index d0b3a55e04..dd9dedd035 100644 --- a/tests/test_statepoint_restart/results.py +++ b/tests/test_statepoint_restart/results.py @@ -11,7 +11,7 @@ import statepoint if len(sys.argv) > 1: sp = statepoint.StatePoint(sys.argv[1]) else: - sp = statepoint.StatePoint('statepoint.7.binary') + sp = statepoint.StatePoint('statepoint.07.binary') sp.read_results() # extract tally results and convert to vector diff --git a/tests/test_statepoint_restart/test_statepoint_restart.py b/tests/test_statepoint_restart/test_statepoint_restart.py index 9473cc6923..a26a94442a 100644 --- a/tests/test_statepoint_restart/test_statepoint_restart.py +++ b/tests/test_statepoint_restart/test_statepoint_restart.py @@ -24,13 +24,13 @@ def test_run(): assert returncode == 0, 'OpenMC did not exit successfully.' def test_created_statepoint(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) 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.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: @@ -38,7 +38,7 @@ def test_results(): assert compare, 'Initial test results do not agree.' def test_restart_form1(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) if opts.mpi_exec != '': proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, '-r', statepoint[0], cwd], stderr=STDOUT, stdout=PIPE) @@ -49,13 +49,13 @@ def test_restart_form1(): assert returncode == 0, 'OpenMC restart 1 did not exit successfully.' def test_created_statepoint_form1(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) 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_form1(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: @@ -63,7 +63,7 @@ def test_results_form1(): assert compare, 'Restart 1 test results do not agree.' def test_restart_form2(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) if opts.mpi_exec != '': proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, '--restart', statepoint[0], cwd], stderr=STDOUT, stdout=PIPE) @@ -74,13 +74,13 @@ def test_restart_form2(): assert returncode == 0, 'OpenMC restart 2 did not exit successfully.' def test_created_statepoint_form2(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) 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_form2(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: @@ -88,20 +88,20 @@ def test_results_form2(): assert compare, 'Restart 2 test results do not agree.' def test_restart_serial(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) proc = Popen([opts.exe, '--restart', statepoint[0], cwd], stderr=STDOUT, stdout=PIPE) print(proc.communicate()[0]) returncode = proc.returncode assert returncode == 0, 'OpenMC restart serial did not exit successfully.' def test_created_statepoint_serial(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) 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_serial(): - statepoint = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + statepoint = glob.glob(os.path.join(cwd, 'statepoint.07.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: @@ -109,7 +109,7 @@ def test_results_serial(): assert compare, 'Restart serial test results do not agree.' def teardown(): - output = glob.glob(os.path.join(cwd, 'statepoint.7.*')) + output = glob.glob(os.path.join(cwd, 'statepoint.07.*')) output.append(os.path.join(cwd, 'results_test.dat')) for f in output: if os.path.exists(f): From 147af675fefc671bb9a4b1fa4fb9090fdae693c3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 11 Jul 2014 17:36:09 +0200 Subject: [PATCH 07/18] Commented zero_padded_batch --- src/state_point.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/state_point.F90 b/src/state_point.F90 index cb276e2700..a4c4375722 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -809,7 +809,9 @@ contains end subroutine load_state_point !=============================================================================== -! STATE_POINT -- This module handles writing and reading state point +! ZERO_PADDED_BATCH -- This function returns the current batch number as a +! string formatted into a zero padded integer. For example, batch 3 out of 100 +! will be returned as '003'. !=============================================================================== function zero_padded_batch() From df261e9a92a6125d7e88dc52aa128ab29c610630 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 18 Jul 2014 15:27:36 +0200 Subject: [PATCH 08/18] Moved zero-padding features to string.F90 --- src/state_point.F90 | 40 +++++++---------------------------- src/string.F90 | 51 ++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 58 insertions(+), 33 deletions(-) diff --git a/src/state_point.F90 b/src/state_point.F90 index a4c4375722..8aff2ec364 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -16,7 +16,7 @@ module state_point use error, only: fatal_error, warning use global use output, only: write_message, time_stamp - use string, only: to_str + use string, only: to_str, zero_padded, count_digits use output_interface use tally_header, only: TallyObject @@ -43,7 +43,8 @@ contains type(TallyObject), pointer :: t => null() ! Set filename for state point - filename = trim(path_output) // 'statepoint.' // zero_padded_batch() + filename = trim(path_output) // 'statepoint.' // & + & zero_padded(current_batch, count_digits(n_batches)) ! Append appropriate extension #ifdef HDF5 @@ -304,7 +305,9 @@ contains if (source_separate) then ! Set filename - filename = trim(path_output) // 'source.' // zero_padded_batch() + filename = trim(path_output) // 'source.' // & + & zero_padded(current_batch, count_digits(n_batches)) + #ifdef HDF5 filename = trim(filename) // '.h5' #else @@ -324,7 +327,8 @@ contains else ! Set filename for state point - filename = trim(path_output) // 'statepoint.' // zero_padded_batch() + filename = trim(path_output) // 'statepoint.' // & + & zero_padded(current_batch, count_digits(n_batches)) #ifdef HDF5 filename = trim(filename) // '.h5' #else @@ -808,34 +812,6 @@ contains end subroutine load_state_point -!=============================================================================== -! ZERO_PADDED_BATCH -- This function returns the current batch number as a -! string formatted into a zero padded integer. For example, batch 3 out of 100 -! will be returned as '003'. -!=============================================================================== - - function zero_padded_batch() - character(MAX_FILE_LEN) :: zero_padded_batch - - character(MAX_FILE_LEN) :: batch_form - character(MAX_FILE_LEN) :: batch_string - character(MAX_FILE_LEN) :: filename - integer :: n_digits - - ! Count digits needed to represent largest batch number. - n_digits = 1 - do while (n_batches / 10**(n_digits) /= 0 & - & .and. n_batches / 10 **(n_digits-1) /= 1) - n_digits = n_digits + 1 - end do - - ! Write a format string of the form In.n where n is n_digits. - write(batch_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits - - ! Format the batch number - write(zero_padded_batch, batch_form) current_batch - end function zero_padded_batch - subroutine read_source ! TODO write this routine ! TODO what if n_particles does not match source bank diff --git a/src/string.F90 b/src/string.F90 index e322894b55..d04dbd0532 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -1,7 +1,7 @@ module string use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL - use error, only: warning + use error, only: fatal_error, warning use global, only: message implicit none @@ -184,6 +184,39 @@ contains end subroutine upper_case +!=============================================================================== +! ZERO_PADDED returns a string of the input integer padded with zeros to the +! desired number of digits. Do not include the sign in n_digits for negative +! integers. +!=============================================================================== + +function zero_padded(num, n_digits) result(str) + integer, intent(in) :: num + integer, intent(in) :: n_digits + character(11) :: str + + character(8) :: zp_form + + ! Make sure n_digits is reasonable. 10 digits is the maximum needed for the + ! largest integer(4). + if (n_digits > 10) then + message = 'int4_to_padstr called with unreasonably large n_digits (>10)' + call fatal_error() + end if + + ! Write a format string of the form '(In.m)' where n is the max width and + ! m is the min width. If a sign is present, then n must be one greater + ! than m. + if (num < 0) then + write(zp_form, '("(I", I0, ".", I0, ")")') n_digits+1, n_digits + else + write(zp_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits + end if + + ! Format the number. + write(str, zp_form) num +end function zero_padded + !=============================================================================== ! IS_NUMBER determines whether a string of characters is all 0-9 characters !=============================================================================== @@ -268,6 +301,22 @@ contains end function ends_with +!=============================================================================== +! COUNT_DIGITS returns the number of digits needed to represent the input +! integer. +!=============================================================================== + + function count_digits(num) result(n_digits) + integer, intent(in) :: num + integer :: n_digits + + n_digits = 1 + do while (num / 10**(n_digits) /= 0 .and. abs(num / 10 **(n_digits-1)) /= 1) + n_digits = n_digits + 1 + end do + + end function count_digits + !=============================================================================== ! INT4_TO_STR converts an integer(4) to a string. !=============================================================================== From e99f964d475d2023511a937f2ac3f9d26a97bf1b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 29 Jul 2014 10:25:46 +0200 Subject: [PATCH 09/18] Fixed zero_padded error message typo --- src/string.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/string.F90 b/src/string.F90 index d04dbd0532..c53645c15e 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -200,7 +200,7 @@ function zero_padded(num, n_digits) result(str) ! Make sure n_digits is reasonable. 10 digits is the maximum needed for the ! largest integer(4). if (n_digits > 10) then - message = 'int4_to_padstr called with unreasonably large n_digits (>10)' + message = 'zero_padded called with an unreasonably large n_digits (>10)' call fatal_error() end if From 0c998286ad92b5e7e890ec8c28275841188e70e5 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 30 Jul 2014 15:07:03 +0200 Subject: [PATCH 10/18] Fixed overflow error in count_digits --- src/string.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/string.F90 b/src/string.F90 index c53645c15e..f78b94d07f 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -311,7 +311,10 @@ end function zero_padded integer :: n_digits n_digits = 1 - do while (num / 10**(n_digits) /= 0 .and. abs(num / 10 **(n_digits-1)) /= 1) + do while (num / 10**(n_digits) /= 0 .and. abs(num / 10 **(n_digits-1)) /= 1& + &.and. n_digits /= 10) + ! Note that 10 digits is the maximum needed to represent an integer(4) so + ! the loop automatically exits when n_digits = 10. n_digits = n_digits + 1 end do From 3d38155179383ab2ab284dc7aa3a393e7e07874a Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Wed, 30 Jul 2014 19:06:31 -0600 Subject: [PATCH 11/18] fixed comment that was over 80 characters --- src/cross_section.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 2f7bcf9e3a..815e6e90d7 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -374,8 +374,10 @@ contains (urr % energy(i_energy + 1) - urr % energy(i_energy)) ! sample probability table using the cumulative distribution + ! (if we're dealing with an isotope that we've previously encountered at - ! this energy but a different temperature, use the original random number to ! preserve correlation of temperature in probability tables) + ! this energy but a different temperature, use the original random number to + ! preserve correlation of temperature in probability tables) if (E /= ZERO .and. E == first_E .and. & & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & & + abs(micro_xs(:) % last_E - E)) == ZERO) then From 560e979836afb8d0cb135af565da09f235c16d0f Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Wed, 30 Jul 2014 22:10:49 -0600 Subject: [PATCH 12/18] check for void material to avoid out of bounds --- src/geometry.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 1a2544d28d..03eddb5ec4 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -190,7 +190,8 @@ contains ! set material p % last_material = p % material p % material = c % material - materials(p % material) % first_E = p % E + if(p % material /= MATERIAL_VOID) & + & materials(p % material) % first_E = p % E elseif (c % type == CELL_FILL) then ! ==================================================================== From 60127faddca151c328f397012c118ae2afd13760 Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Thu, 31 Jul 2014 21:22:14 -0600 Subject: [PATCH 13/18] removed vestigial check, variable --- src/cross_section.F90 | 14 ++++++-------- src/geometry.F90 | 2 -- src/material_header.F90 | 1 - 3 files changed, 6 insertions(+), 11 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 815e6e90d7..cb97e5512a 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -92,9 +92,9 @@ contains ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, mat % first_E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, mat % first_E) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E) end if ! ======================================================================== @@ -135,12 +135,11 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, first_E) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy - real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid @@ -236,7 +235,7 @@ contains if (urr_ptables_on .and. nuc % urr_present) then if (E > nuc % urr_data % energy(1) .and. & E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then - call calculate_urr_xs(i_nuclide, E, first_E) + call calculate_urr_xs(i_nuclide, E) end if end if @@ -337,11 +336,10 @@ contains ! from probability tables !=============================================================================== - subroutine calculate_urr_xs(i_nuclide, E, first_E) + subroutine calculate_urr_xs(i_nuclide, E) integer, intent(in) :: i_nuclide ! index into nuclides array real(8), intent(in) :: E ! energy - real(8), intent(in) :: first_E ! particle energy when it entered this mat integer :: i_energy ! index for energy integer :: i_table ! index for table @@ -378,7 +376,7 @@ contains ! (if we're dealing with an isotope that we've previously encountered at ! this energy but a different temperature, use the original random number to ! preserve correlation of temperature in probability tables) - if (E /= ZERO .and. E == first_E .and. & + if (E /= ZERO .and. & & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & & + abs(micro_xs(:) % last_E - E)) == ZERO) then r = micro_xs(minloc(dble(abs(nuclides(:) % zaid & diff --git a/src/geometry.F90 b/src/geometry.F90 index 03eddb5ec4..24d3ef8b06 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -190,8 +190,6 @@ contains ! set material p % last_material = p % material p % material = c % material - if(p % material /= MATERIAL_VOID) & - & materials(p % material) % first_E = p % E elseif (c % type == CELL_FILL) then ! ==================================================================== diff --git a/src/material_header.F90 b/src/material_header.F90 index d4c7e81290..cbfd917498 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -12,7 +12,6 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm - real(8) :: first_E ! particle energy upon entering mat ! S(a,b) data references integer :: n_sab = 0 ! number of S(a,b) tables From 3e0870ac96988e5408741f05a356f1648de62f9d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 5 Aug 2014 22:49:44 -0400 Subject: [PATCH 14/18] Don't use PWD environment variable when setting path to input files. Closes #302 on github. --- src/initialize.F90 | 28 ++++++++++------------------ 1 file changed, 10 insertions(+), 18 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index c14a8a02a7..8d465ef7d7 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -211,13 +211,13 @@ contains call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) - ! Adjust displacements + ! Adjust displacements bank_disp = bank_disp - bank_disp(1) ! Define MPI_BANK for fission sites bank_blocks = (/ 1, 3, 3, 1 /) bank_types = (/ MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) - call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & + call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) @@ -308,13 +308,9 @@ contains integer :: argc ! number of command line arguments integer :: last_flag ! index of last flag integer :: filetype - character(MAX_FILE_LEN) :: pwd ! present working directory character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments type(BinaryOutput) :: sp - ! Get working directory - call GET_ENVIRONMENT_VARIABLE("PWD", pwd) - ! Check number of command line arguments and allocate argv argc = COMMAND_ARGUMENT_COUNT() @@ -414,7 +410,7 @@ contains ! Read number of threads i = i + 1 -#ifdef _OPENMP +#ifdef _OPENMP ! Read and set number of OpenMP threads n_threads = int(str_to_int(argv(i)), 4) if (n_threads < 1) then @@ -454,16 +450,12 @@ contains ! Determine directory where XML input files are if (argc > 0 .and. last_flag < argc) then path_input = argv(last_flag + 1) - ! Need to add working directory if the given path is a relative path - if (.not. starts_with(path_input, "/")) then - path_input = trim(pwd) // "/" // trim(path_input) - end if else - path_input = pwd + path_input = '' end if ! Add slash at end of directory if it isn't there - if (.not. ends_with(path_input, "/")) then + if (.not. ends_with(path_input, "/") .and. len_trim(path_input) > 0) then path_input = trim(path_input) // "/" end if @@ -563,7 +555,7 @@ contains integer :: m ! loop index for lattices integer :: mid, lid ! material and lattice IDs integer :: n_x, n_y, n_z ! size of lattice - integer :: i_array ! index in surfaces/materials array + integer :: i_array ! index in surfaces/materials array integer :: id ! user-specified id type(Cell), pointer :: c => null() type(Lattice), pointer :: lat => null() @@ -634,7 +626,7 @@ contains " specified on lattice " // trim(to_str(lid)) call fatal_error() end if - + else message = "Specified fill " // trim(to_str(id)) // " on cell " // & trim(to_str(c % id)) // " is neither a universe nor a lattice." @@ -805,7 +797,7 @@ contains sum_percent = sum_percent + x*awr end do sum_percent = ONE / sum_percent - mat % density = -mat % density * N_AVOGADRO & + mat % density = -mat % density * N_AVOGADRO & / MASS_NEUTRON * sum_percent end if @@ -866,7 +858,7 @@ contains ! Allocate source bank allocate(source_bank(work), STAT=alloc_err) - ! Check for allocation errors + ! Check for allocation errors if (alloc_err /= 0) then message = "Failed to allocate source bank." call fatal_error() @@ -894,7 +886,7 @@ contains allocate(fission_bank(3*work), STAT=alloc_err) #endif - ! Check for allocation errors + ! Check for allocation errors if (alloc_err /= 0) then message = "Failed to allocate fission bank." call fatal_error() From 590fd55b1787d076aa9b0c249f2662d301d04c1f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 6 Aug 2014 15:38:17 +0200 Subject: [PATCH 15/18] Fixed tests for zero_padded_statepoint --- src/CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ee5f882b7e..0025ba79e7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -342,9 +342,9 @@ foreach(test ${TESTS}) # Handle restart tests separately if(${test} MATCHES "test_statepoint_restart") - set(RESTART_FILE statepoint.7.h5) + set(RESTART_FILE statepoint.07.h5) elseif(${test} MATCHES "test_sourcepoint_restart") - set(RESTART_FILE statepoint.7.h5 source.7.h5) + set(RESTART_FILE statepoint.07.h5 source.07.h5) elseif(${test} MATCHES "test_particle_restart") set(RESTART_FILE particle_12_192.h5) else(${test} MATCHES "test_statepoint_restart") @@ -355,9 +355,9 @@ foreach(test ${TESTS}) # Handle restart tests separately if(${test} MATCHES "test_statepoint_restart") - set(RESTART_FILE statepoint.7.binary) + set(RESTART_FILE statepoint.07.binary) elseif(${test} MATCHES "test_sourcepoint_restart") - set(RESTART_FILE statepoint.7.binary source.7.binary) + set(RESTART_FILE statepoint.07.binary source.07.binary) elseif(${test} MATCHES "test_particle_restart") set(RESTART_FILE particle_12_192.binary) else(${test} MATCHES "test_statepoint_restart") From e6e4898045adc982d4a7e4b89d0818bb8eb7107b Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Sat, 9 Aug 2014 20:18:38 -0600 Subject: [PATCH 16/18] use linked list to find multiple instances of a single nuclide --- src/ace.F90 | 28 ++++++++++++++++++++++++++++ src/ace_header.F90 | 4 ++++ src/cross_section.F90 | 32 +++++++++++++++++++++----------- src/initialize.F90 | 5 ++++- 4 files changed, 57 insertions(+), 12 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index adff30dbeb..58f7fd5707 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -7,6 +7,7 @@ module ace use error, only: fatal_error, warning use fission, only: nu_total use global + use list_header, only: ListElemInt use material_header, only: Material use output, only: write_message use set_header, only: SetChar @@ -1428,4 +1429,31 @@ contains end function get_real +!=============================================================================== +! SAME_NUCLIDE_LIST creates a linked list for each nuclide containing the +! indices in the nuclides array of all other instances of that nuclide. For +! example, the same nuclide may exist at multiple temperatures resulting +! in multiple entries in the nuclides array for a single zaid number. +!=============================================================================== + + subroutine same_nuclide_list() + + integer :: i ! index in nuclides array + integer :: j ! index in nuclides array + type(ListElemInt), pointer :: nuc_list => null() ! pointer to nuclide list + + do i = 1, n_nuclides_total + allocate(nuclides(i) % nuc_list) + nuc_list => nuclides(i) % nuc_list + do j = 1, n_nuclides_total + if (nuclides(i) % zaid == nuclides(j) % zaid) then + nuc_list % data = j + allocate(nuc_list % next) + nuc_list => nuc_list % next + end if + end do + end do + + end subroutine same_nuclide_list + end module ace diff --git a/src/ace_header.F90 b/src/ace_header.F90 index fe645683be..2067030c1e 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -2,6 +2,7 @@ module ace_header use constants, only: MAX_FILE_LEN use endf_header, only: Tab1 + use list_header, only: ListElemInt implicit none @@ -95,6 +96,9 @@ module ace_header real(8) :: awr ! weight of nucleus in neutron masses real(8) :: kT ! temperature in MeV (k*T) + ! Linked list of indices in nuclides array of instances of this same nuclide + type(ListElemInt), pointer :: nuc_list => null() + ! Energy grid information integer :: n_grid ! # of nuclide grid points integer, allocatable :: grid_index(:) ! pointers to union grid diff --git a/src/cross_section.F90 b/src/cross_section.F90 index cb97e5512a..cc810dc5e6 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -5,6 +5,7 @@ module cross_section use error, only: fatal_error use fission, only: nu_total use global + use list_header, only: ListElemInt use material_header, only: Material use particle_header, only: Particle use random_lcg, only: prn @@ -349,9 +350,11 @@ contains real(8) :: capture ! (n,gamma) cross section real(8) :: fission ! fission cross section real(8) :: inelastic ! inelastic cross section - type(UrrData), pointer, save :: urr => null() - type(Nuclide), pointer, save :: nuc => null() - type(Reaction), pointer, save :: rxn => null() + logical :: same_nuc ! do we know the xs for this nuclide at this energy? + type(UrrData), pointer, save :: urr => null() + type(Nuclide), pointer, save :: nuc => null() + type(Reaction), pointer, save :: rxn => null() + type(ListElemInt), pointer :: nuc_list => null() !$omp threadprivate(urr, nuc, rxn) micro_xs(i_nuclide) % use_ptable = .true. @@ -373,15 +376,22 @@ contains ! sample probability table using the cumulative distribution - ! (if we're dealing with an isotope that we've previously encountered at + ! if we're dealing with a nuclide that we've previously encountered at ! this energy but a different temperature, use the original random number to - ! preserve correlation of temperature in probability tables) - if (E /= ZERO .and. & - & minval(dble(abs(nuclides(:) % zaid - nuclides(i_nuclide) % zaid)) & - & + abs(micro_xs(:) % last_E - E)) == ZERO) then - r = micro_xs(minloc(dble(abs(nuclides(:) % zaid & - & - nuclides(i_nuclide) % zaid)) + abs(micro_xs(:) % last_E & - & - E),1)) % last_prn + ! preserve correlation of temperature in probability tables + same_nuc = .false. + nuc_list => nuc % nuc_list + do + if (E /= ZERO .and. E == micro_xs(nuc_list % data) % last_E) then + same_nuc = .true. + exit + end if + nuc_list => nuc_list % next + if (.not. associated(nuc_list % next)) exit + end do + + if (same_nuc) then + r = micro_xs(nuc_list % data) % last_prn else r = prn() micro_xs(i_nuclide) % last_prn = r diff --git a/src/initialize.F90 b/src/initialize.F90 index c14a8a02a7..a5c64101f9 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1,6 +1,6 @@ module initialize - use ace, only: read_xs + use ace, only: read_xs, same_nuclide_list use bank_header, only: Bank use constants use dict_header, only: DictIntInt, ElemKeyValueII @@ -104,6 +104,9 @@ contains call read_xs() call time_read_xs % stop() + ! Create linked lists for multiple instances of the same nuclide + call same_nuclide_list() + ! Construct unionized energy grid from cross-sections if (grid_method == GRID_UNION) then call time_unionize % start() From 251f71deeb7b7947802345a5b8c53e646cad96f8 Mon Sep 17 00:00:00 2001 From: Jon Walsh Date: Sat, 9 Aug 2014 21:26:22 -0600 Subject: [PATCH 17/18] make the linked list threadprivate --- src/cross_section.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 23df55acd5..18a23a6b0f 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -356,7 +356,7 @@ contains type(Nuclide), pointer, save :: nuc => null() type(Reaction), pointer, save :: rxn => null() type(ListElemInt), pointer :: nuc_list => null() -!$omp threadprivate(urr, nuc, rxn) +!$omp threadprivate(urr, nuc, rxn, nuc_list) micro_xs(i_nuclide) % use_ptable = .true. From 7ac0a72200f6d093322fd82a106a1e8acd76560c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Aug 2014 19:38:44 -0400 Subject: [PATCH 18/18] Fix restart for test_particle_restart (missed this when merging walshjon:band-samp). --- src/CMakeLists.txt | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0025ba79e7..6db242a090 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -167,28 +167,28 @@ if(petsc) list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_FLAPACK_LIB}) endif() - # If libdl wasn't found, search /usr/lib64 + # If libdl wasn't found, search /usr/lib64 if(PETSC_DL_LIB STREQUAL "PETSC_DL_LIB-NOTFOUND") find_library(PETSC_DL_LIB libdl.so /usr/lib64) list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_DL_LIB-NOTFOUND) list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_DL_LIB}) endif() - # If libm wasn't found, search /usr/lib64 + # If libm wasn't found, search /usr/lib64 if(PETSC_M_LIB STREQUAL "PETSC_M_LIB-NOTFOUND") find_library(PETSC_M_LIB libm.so /usr/lib64) list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_M_LIB-NOTFOUND) list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_M_LIB}) endif() - # If libpthread wasn't found, search /usr/lib64 + # If libpthread wasn't found, search /usr/lib64 if(PETSC_PTHREAD_LIB STREQUAL "PETSC_PTHREAD_LIB-NOTFOUND") find_library(PETSC_PTHREAD_LIB libpthread.so /usr/lib64) list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_PTHREAD_LIB-NOTFOUND) list(INSERT PETSC_PACKAGE_LIBS 0 ${PETSC_PTHREAD_LIB}) endif() - # If librt wasn't found, search /usr/lib64 + # If librt wasn't found, search /usr/lib64 if(PETSC_RT_LIB STREQUAL "PETSC_RT_LIB-NOTFOUND") find_library(PETSC_RT_LIB librt.so /usr/lib64) list(REMOVE_ITEM PETSC_PACKAGE_LIBS PETSC_RT_LIB-NOTFOUND) @@ -213,7 +213,7 @@ execute_process(COMMAND git rev-parse HEAD if(GIT_SHA1_SUCCESS EQUAL 0) add_definitions(-DGIT_SHA1="${GIT_SHA1}") endif() - + #=============================================================================== # FoX Fortran XML Library #=============================================================================== @@ -300,7 +300,7 @@ foreach(test ${TESTS}) # Get test information get_filename_component(TEST_NAME ${test} NAME) - get_filename_component(TEST_PATH ${test} PATH) + get_filename_component(TEST_PATH ${test} PATH) # Check for running standard tests (no valgrind, no gcov) if(NOT ${MEM_CHECK} AND NOT ${COVERAGE}) @@ -309,15 +309,15 @@ foreach(test ${TESTS}) if (${MPI_ENABLED}) # Preform a parallel test - add_test(NAME ${TEST_NAME} - WORKING_DIRECTORY ${TEST_PATH} + add_test(NAME ${TEST_NAME} + WORKING_DIRECTORY ${TEST_PATH} COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $ --mpi_exec $ENV{MPI_DIR}/bin/mpiexec) else(${MPI_ENABLED}) # Perform a serial test - add_test(NAME ${TEST_NAME} + add_test(NAME ${TEST_NAME} WORKING_DIRECTORY ${TEST_PATH} COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $) @@ -346,7 +346,7 @@ foreach(test ${TESTS}) elseif(${test} MATCHES "test_sourcepoint_restart") set(RESTART_FILE statepoint.07.h5 source.07.h5) elseif(${test} MATCHES "test_particle_restart") - set(RESTART_FILE particle_12_192.h5) + set(RESTART_FILE particle_12_885.h5) else(${test} MATCHES "test_statepoint_restart") message(FATAL_ERROR "Restart test ${test} not recognized") endif(${test} MATCHES "test_statepoint_restart") @@ -359,7 +359,7 @@ foreach(test ${TESTS}) elseif(${test} MATCHES "test_sourcepoint_restart") set(RESTART_FILE statepoint.07.binary source.07.binary) elseif(${test} MATCHES "test_particle_restart") - set(RESTART_FILE particle_12_192.binary) + set(RESTART_FILE particle_12_885.binary) else(${test} MATCHES "test_statepoint_restart") message(FATAL_ERROR "Restart test ${test} not recognized") endif(${test} MATCHES "test_statepoint_restart")