From 1e639f54d4b2f96299d19d6627d10e0aeedb74a8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Jan 2013 18:37:20 -0500 Subject: [PATCH 01/12] Ability to tally arbitrary MT for analog tallies. --- src/input_xml.F90 | 27 ++++++++++++++++++++++++--- src/output.F90 | 6 +++++- src/tally.F90 | 8 ++++++-- 3 files changed, 35 insertions(+), 6 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 374c1ef97..ddcf4fcc7 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1221,6 +1221,7 @@ contains integer :: n_scores ! number of tot scores after adjusting for Pn tally integer :: n_order ! Scattering order requested integer :: n_order_pos ! Position of Scattering order in score name string + integer :: MT ! user-specified MT for score logical :: file_exists ! does tallies.xml file exist? character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word @@ -1908,9 +1909,29 @@ contains t % score_bins(j) = SCORE_EVENTS case default - message = "Unknown scoring function: " // & - trim(tally_(i) % scores(j)) - call fatal_error() + ! Assume that user has specified an MT number + MT = str_to_int(score_name) + + if (MT /= ERROR_INT) then + ! Specified score was an integer + if (MT > 1) then + t % score_bins(j) = MT + + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + else + message = "Invalid MT on : " // & + trim(tally_(i) % scores(j)) + call fatal_error() + end if + + else + ! Specified score was not an integer + message = "Unknown scoring function: " // & + trim(tally_(i) % scores(j)) + call fatal_error() + end if + end select end do t % n_score_bins = n_scores diff --git a/src/output.F90 b/src/output.F90 index ff3fdf163..b6904e3bf 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1623,7 +1623,11 @@ contains end do k = k + n_order - 1 else - score_name = score_names(abs(t % score_bins(k))) + if (t % score_bins(k) > 0) then + score_name = reaction_name(t % score_bins(k)) + else + score_name = score_names(abs(t % score_bins(k))) + end if write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & diff --git a/src/tally.F90 b/src/tally.F90 index e1c05d8d2..45f886bf9 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -337,8 +337,12 @@ contains score = ONE case default - message = "Invalid score type on tally " // to_str(t % id) // "." - call fatal_error() + ! Any other score is assumed to be a MT number. Thus, we just need + ! to check if it matches the MT number of the event + if (p % event_MT /= score_bin) cycle SCORE_LOOP + + score = last_wgt + end select ! Add score to tally From 5131eb3e0998e87798333a8775d59256087b5efb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Jan 2013 19:08:10 -0500 Subject: [PATCH 02/12] Remove special (n,xn) tallies. Need to fix statepoint.py. --- src/constants.F90 | 15 ++++++--------- src/input_xml.F90 | 6 +++--- src/output.F90 | 9 --------- src/tally.F90 | 15 --------------- src/utils/statepoint.py | 5 ++--- 5 files changed, 11 insertions(+), 39 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 83edf571c..fc51c8bc8 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -252,7 +252,7 @@ module constants EVENT_FISSION = 3 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 17 + integer, parameter :: N_SCORE_TYPES = 14 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -263,14 +263,11 @@ module constants SCORE_TRANSPORT = -7, & ! transport reaction rate SCORE_DIFFUSION = -8, & ! diffusion coefficient SCORE_N_1N = -9, & ! (n,1n) rate - SCORE_N_2N = -10, & ! (n,2n) rate - SCORE_N_3N = -11, & ! (n,3n) rate - SCORE_N_4N = -12, & ! (n,4n) rate - SCORE_ABSORPTION = -13, & ! absorption rate - SCORE_FISSION = -14, & ! fission rate - SCORE_NU_FISSION = -15, & ! neutron production rate - SCORE_CURRENT = -16, & ! partial current - SCORE_EVENTS = -17 ! number of events + SCORE_ABSORPTION = -10, & ! absorption rate + SCORE_FISSION = -11, & ! fission rate + SCORE_NU_FISSION = -12, & ! neutron production rate + SCORE_CURRENT = -13, & ! partial current + SCORE_EVENTS = -14 ! number of events ! Maximum scattering order supported integer, parameter :: SCATT_ORDER_MAX = 10 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ddcf4fcc7..2db157d01 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1824,17 +1824,17 @@ contains ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('n2n') - t % score_bins(j) = SCORE_N_2N + t % score_bins(j) = N_2N ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('n3n') - t % score_bins(j) = SCORE_N_3N + t % score_bins(j) = N_3N ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('n4n') - t % score_bins(j) = SCORE_N_4N + t % score_bins(j) = N_4N ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG diff --git a/src/output.F90 b/src/output.F90 index b6904e3bf..14fed4843 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -809,12 +809,6 @@ contains string = trim(string) // ' diffusion' case (SCORE_N_1N) string = trim(string) // ' n1n' - case (SCORE_N_2N) - string = trim(string) // ' n2n' - case (SCORE_N_3N) - string = trim(string) // ' n3n' - case (SCORE_N_4N) - string = trim(string) // ' n4n' case (SCORE_ABSORPTION) string = trim(string) // ' absorption' case (SCORE_FISSION) @@ -1455,9 +1449,6 @@ contains score_names(abs(SCORE_TRANSPORT)) = "Transport Rate" score_names(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient" score_names(abs(SCORE_N_1N)) = "(n,1n) Rate" - score_names(abs(SCORE_N_2N)) = "(n,2n) Rate" - score_names(abs(SCORE_N_3N)) = "(n,3n) Rate" - score_names(abs(SCORE_N_4N)) = "(n,4n) Rate" score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate" score_names(abs(SCORE_FISSION)) = "Fission Rate" score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" diff --git a/src/tally.F90 b/src/tally.F90 index 45f886bf9..32784149e 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -245,21 +245,6 @@ contains ! All events that reach this point are (n,1n) reactions score = last_wgt - case (SCORE_N_2N) - ! Skip any event that is not (n,2n) - if (p % event_MT /= N_2N) cycle SCORE_LOOP - score = last_wgt - - case (SCORE_N_3N) - ! Skip any event that is not (n,3n) - if (p % event_MT /= N_3N) cycle SCORE_LOOP - score = last_wgt - - case (SCORE_N_4N) - ! Skip any event that is not (n,4n) - if (p % event /= N_4N) cycle SCORE_LOOP - score = last_wgt - case (SCORE_ABSORPTION) if (survival_biasing) then ! No absorption events actually occur if survival biasing is on -- diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 584e83772..1e0acbd74 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -12,9 +12,8 @@ filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', score_types = {-1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter', -5: 'scatter-n', -6: 'scatter-pn', -7: 'transport', - -8: 'diffusion', -9: 'n1n', -10: 'n2n', -11: 'n3n', - -12: 'n4n', -13: 'absorption', -14: 'fission', - -15: 'nu-fission', -16: 'current', -17: 'events'} + -8: 'diffusion', -9: 'n1n', -10: 'absorption', -11: 'fission', + -12: 'nu-fission', -13: 'current', -14: 'events'} class BinaryFile(object): def __init__(self, filename): From b2419bf1503ee3643ed24d70ba063ffc537d2d88 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Jan 2013 19:09:07 -0500 Subject: [PATCH 03/12] Add *.binary files to .gitignore. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index 118a58bf4..8aa7c1f97 100644 --- a/.gitignore +++ b/.gitignore @@ -13,6 +13,9 @@ src/openmc # emacs backups *~ +# OpenMC statepoints +*.binary + # Documentation builds docs/build From a44090d1e7e094a8ea6506170f13f3b35357dc5f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 5 Jan 2013 16:13:11 -0500 Subject: [PATCH 04/12] Ability to tally arbitrary MT for trackelngth tallies. Still need to change score_all_nuclides. --- src/DEPENDENCIES | 1 + src/input_xml.F90 | 9 ---- src/tally.F90 | 134 +++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 128 insertions(+), 16 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 0d2e4e2de..1b9b6258a 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -336,6 +336,7 @@ string.o: constants.o string.o: error.o string.o: global.o +tally.o: ace_header.o tally.o: constants.o tally.o: datatypes_header.o tally.o: error.o diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2db157d01..d9dec9fbc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1826,18 +1826,12 @@ contains case ('n2n') t % score_bins(j) = N_2N - ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG case ('n3n') t % score_bins(j) = N_3N - ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG case ('n4n') t % score_bins(j) = N_4N - ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG case ('absorption') t % score_bins(j) = SCORE_ABSORPTION if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -1916,9 +1910,6 @@ contains ! Specified score was an integer if (MT > 1) then t % score_bins(j) = MT - - ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG else message = "Invalid MT on : " // & trim(tally_(i) % scores(j)) diff --git a/src/tally.F90 b/src/tally.F90 index 32784149e..7e5dd771a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,5 +1,6 @@ module tally + use ace_header, only: Reaction use constants use datatypes_header, only: ListInt use error, only: fatal_error @@ -423,10 +424,15 @@ contains integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins + integer :: l ! loop index for nuclides in material + integer :: m ! loop index for reactions integer :: filter_index ! single index for single bin - integer :: i_nuclide ! index in nuclides array + integer :: i_nuclide ! index in nuclides array (from bins) + integer :: i_nuc ! index in nuclides array (from material) + integer :: i_energy ! index in nuclide energy grid integer :: score_bin ! scoring type, e.g. SCORE_FLUX integer :: score_index ! scoring bin index + real(8) :: f ! interpolation factor real(8) :: flux ! tracklength estimate of flux real(8) :: score ! actual score (e.g., flux*xs) real(8) :: atom_density ! atom density of single nuclide in atom/b-cm @@ -434,6 +440,7 @@ contains type(TallyObject), pointer :: t => null() type(Material), pointer :: mat => null() type(ListInt), pointer :: curr_ptr => null() + type(Reaction), pointer :: rxn => null() ! Determine track-length estimate of flux flux = p % wgt * distance @@ -510,51 +517,164 @@ contains score_bin = t % score_bins(j) if (i_nuclide > 0) then - ! Determine macroscopic nuclide cross section + ! ================================================================ + ! DETERMINE NUCLIDE CROSS SECTION + select case(score_bin) case (SCORE_TOTAL) + ! Total cross section is pre-calculated score = micro_xs(i_nuclide) % total * & atom_density * flux + case (SCORE_SCATTER) + ! Scattering cross section is pre-calculated score = (micro_xs(i_nuclide) % total - & micro_xs(i_nuclide) % absorption) * & atom_density * flux + case (SCORE_ABSORPTION) + ! Absorption cross section is pre-calculated score = micro_xs(i_nuclide) % absorption * & atom_density * flux + case (SCORE_FISSION) + ! Fission cross section is pre-calculated score = micro_xs(i_nuclide) % fission * & atom_density * flux + case (SCORE_NU_FISSION) + ! Nu-fission cross section is pre-calculated score = micro_xs(i_nuclide) % nu_fission * & atom_density * flux + case (SCORE_EVENTS) + ! For number of events, just score unity score = ONE + case default - message = "Invalid score type on tally " // to_str(t % id) // "." - call fatal_error() + ! Any other cross section has to be calculated on-the-fly. For + ! cross sections that are used often (e.g. n2n, ngamma, etc. for + ! depletion), it might make sense to optimize this section or + ! pre-calculate cross sections + + if (score_bin > 1) then + ! Set default score + score = ZERO + + ! TODO: The following search for the matching reaction could + ! be replaced by adding a dictionary on each Nuclide instance + ! of the form {MT: i_reaction, ...} + + REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuclide) % reactions(m) + + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor + + if (i_energy >= rxn % threshold) then + score = ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density * flux + end if + + exit REACTION_LOOP + end if + end do REACTION_LOOP + + else + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end if end select else - ! Determine macroscopic material cross section + ! ================================================================ + ! DETERMINE MATERIAL CROSS SECTION + select case(score_bin) case (SCORE_FLUX) + ! For flux, we need no cross section score = flux + case (SCORE_TOTAL) + ! Total cross section is pre-calculated score = material_xs % total * flux + case (SCORE_SCATTER) + ! Scattering cross section is pre-calculated score = (material_xs % total - material_xs % absorption) * flux + case (SCORE_ABSORPTION) + ! Absorption cross section is pre-calculated score = material_xs % absorption * flux + case (SCORE_FISSION) + ! Fission cross section is pre-calculated score = material_xs % fission * flux + case (SCORE_NU_FISSION) + ! Nu-fission cross section is pre-calculated score = material_xs % nu_fission * flux + case (SCORE_EVENTS) + ! For number of events, just score unity score = ONE + case default - message = "Invalid score type on tally " // to_str(t % id) // "." - call fatal_error() + ! Any other cross section has to be calculated on-the-fly. This + ! is somewhat costly since it requires a loop over each nuclide + ! in a material and each reaction in the nuclide + + if (score_bin > 1) then + ! Set default score + score = ZERO + + ! Get pointer to current material + mat => materials(p % material) + + do l = 1, mat % n_nuclides + ! Get atom density + atom_density = mat % atom_density(l) + + ! Get index in nuclides array + i_nuc = mat % nuclide(l) + + ! TODO: The following search for the matching reaction could + ! be replaced by adding a dictionary on each Nuclide + ! instance of the form {MT: i_reaction, ...} + + do m = 1, nuclides(i_nuc) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuc) % reactions(m) + + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor + + if (i_energy >= rxn % threshold) then + score = score + ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density * flux + end if + + exit + end if + end do + + end do + + else + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end if end select end if From dcca75ce0686518984dedee86914760431e46d91 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 5 Jan 2013 17:21:30 -0500 Subject: [PATCH 05/12] Tally arbitrary MT for tracklength tallies with 'all' nuclides. Closes #130 on github. --- src/tally.F90 | 109 +++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 104 insertions(+), 5 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 7e5dd771a..0ce310a30 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -722,13 +722,17 @@ contains integer :: i ! loop index for nuclides in material integer :: j ! loop index for scoring bin types + integer :: m ! loop index for reactions in nuclide integer :: i_nuclide ! index in nuclides array integer :: score_bin ! type of score, e.g. SCORE_FLUX integer :: score_index ! scoring bin index + integer :: i_energy ! index in nuclide energy grid + real(8) :: f ! interpolation factor real(8) :: score ! actual scoring tally value real(8) :: atom_density ! atom density of single nuclide in atom/b-cm type(TallyObject), pointer :: t => null() type(Material), pointer :: mat => null() + type(Reaction), pointer :: rxn => null() ! Get pointer to tally t => tallies(i_tally) @@ -756,20 +760,61 @@ contains select case(score_bin) case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * atom_density * flux + case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & micro_xs(i_nuclide) % absorption) * atom_density * flux + case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * atom_density * flux + case (SCORE_FISSION) score = micro_xs(i_nuclide) % fission * atom_density * flux + case (SCORE_NU_FISSION) score = micro_xs(i_nuclide) % nu_fission * atom_density * flux + case (SCORE_EVENTS) score = ONE + case default - message = "Invalid score type on tally " // to_str(t % id) // "." - call fatal_error() + ! Any other cross section has to be calculated on-the-fly. For cross + ! sections that are used often (e.g. n2n, ngamma, etc. for depletion), + ! it might make sense to optimize this section or pre-calculate cross + ! sections + + if (score_bin > 1) then + ! Set default score + score = ZERO + + ! TODO: The following search for the matching reaction could be + ! replaced by adding a dictionary on each Nuclide instance of the + ! form {MT: i_reaction, ...} + + REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuclide) % reactions(m) + + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor + + if (i_energy >= rxn % threshold) then + score = ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density * flux + end if + + exit REACTION_LOOP + end if + end do REACTION_LOOP + + else + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end if end select ! Determine scoring bin index based on what the index of the nuclide @@ -785,7 +830,7 @@ contains end do NUCLIDE_LOOP ! ========================================================================== - ! SCORE ALL INDIVIDUAL NUCLIDE REACTION RATES + ! SCORE TOTAL MATERIAL REACTION RATES ! Loop over score types for each bin MATERIAL_SCORE_LOOP: do j = 1, t % n_score_bins @@ -796,21 +841,75 @@ contains select case(score_bin) case (SCORE_FLUX) score = flux + case (SCORE_TOTAL) score = material_xs % total * flux + case (SCORE_SCATTER) score = (material_xs % total - material_xs % absorption) * flux + case (SCORE_ABSORPTION) score = material_xs % absorption * flux + case (SCORE_FISSION) score = material_xs % fission * flux + case (SCORE_NU_FISSION) score = material_xs % nu_fission * flux + case (SCORE_EVENTS) score = ONE + case default - message = "Invalid score type on tally " // to_str(t % id) // "." - call fatal_error() + ! Any other cross section has to be calculated on-the-fly. This is + ! somewhat costly since it requires a loop over each nuclide in a + ! material and each reaction in the nuclide + + if (score_bin > 1) then + ! Set default score + score = ZERO + + ! Get pointer to current material + mat => materials(p % material) + + do i = 1, mat % n_nuclides + ! Get atom density + atom_density = mat % atom_density(i) + + ! Get index in nuclides array + i_nuclide = mat % nuclide(i) + + ! TODO: The following search for the matching reaction could + ! be replaced by adding a dictionary on each Nuclide + ! instance of the form {MT: i_reaction, ...} + + do m = 1, nuclides(i_nuclide) % n_reaction + ! Get pointer to reaction + rxn => nuclides(i_nuclide) % reactions(m) + + ! Check if this is the desired MT + if (score_bin == rxn % MT) then + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor + + if (i_energy >= rxn % threshold) then + score = score + ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density * flux + end if + + exit + end if + end do + + end do + + else + message = "Invalid score type on tally " // to_str(t % id) // "." + call fatal_error() + end if end select ! Determine scoring bin index based on what the index of the nuclide From df35773e89b513275fc82377b608d185742e65a5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 7 Jan 2013 16:19:44 -0500 Subject: [PATCH 06/12] Add score_bins > 0 to print_tally. --- src/output.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/output.F90 b/src/output.F90 index 14fed4843..6b021e03f 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -817,6 +817,8 @@ contains string = trim(string) // ' nu-fission' case (SCORE_CURRENT) string = trim(string) // ' current' + case default + string = trim(string) // ' ' // reaction_name(t % score_bins(j)) end select end do write(unit_,*) ' Scores:' // trim(string) From d2941db6cd5bd4f7ac481de671691b18c993d2b2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 7 Jan 2013 16:48:15 -0500 Subject: [PATCH 07/12] Fix reading of scatt_order in statepoint.py. --- src/utils/statepoint.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 1e0acbd74..4240ae8ab 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -181,10 +181,11 @@ class StatePoint(BinaryFile): t.n_nuclides = n_nuclides t.nuclides = self._get_int(n_nuclides) - # Read score bins + # Read score bins and scattering order t.n_scores = self._get_int()[0] t.scores = [score_types[j] for j in self._get_int(t.n_scores)] - t.scatt_order = [score_types[j] for j in self._get_int(t.n_scores)] + t.scatt_order = self._get_int(t.n_scores) + # Read number of user score bins t.n_user_scores = self._get_int()[0] From 7720f3332741ce4c31796345942599caaa0d2087 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 7 Jan 2013 16:58:29 -0500 Subject: [PATCH 08/12] Add all MTs to score_types in statepoint.py. --- src/utils/statepoint.py | 83 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 79 insertions(+), 4 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 4240ae8ab..36766c550 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -10,10 +10,85 @@ import scipy.stats filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', 5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'} -score_types = {-1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter', - -5: 'scatter-n', -6: 'scatter-pn', -7: 'transport', - -8: 'diffusion', -9: 'n1n', -10: 'absorption', -11: 'fission', - -12: 'nu-fission', -13: 'current', -14: 'events'} +score_types = {-1: 'flux', + -2: 'total', + -3: 'scatter', + -4: 'nu-scatter', + -5: 'scatter-n', + -6: 'scatter-pn', + -7: 'transport', + -8: 'diffusion', + -9: 'n1n', + -10: 'absorption', + -11: 'fission', + -12: 'nu-fission', + -13: 'current', + -14: 'events', + 1: '(n,total)', + 2: '(n,elastic)', + 4: '(n,level)', + 11: '(n,2nd)', + 16: '(n,2n)', + 17: '(n,3n)', + 18: '(n,fission)', + 19: '(n,f)', + 20: '(n,nf)', + 21: '(n,2nf)', + 22: '(n,na)', + 23: '(n,n3a)', + 24: '(n,2na)', + 25: '(n,3na)', + 28: '(n,np)', + 29: '(n,n2a)', + 30: '(n,2n2a)', + 32: '(n,nd)', + 33: '(n,nt)', + 34: '(n,nHe-3)', + 35: '(n,nd2a)', + 36: '(n,nt2a)', + 37: '(n,4n)', + 38: '(n,3nf)', + 41: '(n,2np)', + 42: '(n,3np)', + 44: '(n,n2p)', + 45: '(n,npa)', + 91: '(n,nc)', + 101: '(n,disappear)', + 102: '(n,gamma)', + 103: '(n,p)', + 104: '(n,d)', + 105: '(n,t)', + 106: '(n,3He)', + 107: '(n,a)', + 108: '(n,2a)', + 109: '(n,3a)', + 111: '(n,2p)', + 112: '(n,pa)', + 113: '(n,t2a)', + 114: '(n,d2a)', + 115: '(n,pd)', + 116: '(n,pt)', + 117: '(n,da)', + 201: '(n,Xn)', + 202: '(n,Xgamma)', + 203: '(n,Xp)', + 204: '(n,Xd)', + 205: '(n,Xt)', + 206: '(n,X3He)', + 207: '(n,Xa)', + 444: '(damage)', + 649: '(n,pc)', + 699: '(n,dc)', + 749: '(n,tc)', + 799: '(n,3Hec)', + 849: '(n,tc)'} +score_types.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)}) +score_types.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)}) +score_types.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)}) +score_types.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)}) +score_types.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)}) +score_types.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)}) + class BinaryFile(object): def __init__(self, filename): From 5e4509fa53dcfc0f04a2d82996bfa2548bf48e97 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Jan 2013 16:11:22 -0500 Subject: [PATCH 09/12] Added instructions on using JEFF ACE data from NEA. --- docs/source/usersguide/install.rst | 59 ++++++++++++++++++++++++------ 1 file changed, 48 insertions(+), 11 deletions(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index fe50dd87f..05f08e064 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -142,26 +142,63 @@ This will build an executable named ``openmc``. .. _Cygwin: http://www.cygwin.com/ --------------------------- -Cross-Section Configuration +Cross Section Configuration --------------------------- -In order to run a simulation with OpenMC, you will need cross-section data for -each nuclide in your problem. Since OpenMC uses ACE format cross-sections, you -can use nuclear data distributed with MCNP_ or Serpent_. +In order to run a simulation with OpenMC, you will need cross section data for +each nuclide in your problem. Since OpenMC uses ACE format cross sections, you +can use nuclear data that was processed with NJOY, such as that distributed with +MCNP_ or Serpent_. + +Using JEFF Cross Sections from OECD/NEA +--------------------------------------- + +The NEA_ provides processed ACE data from the JEFF_ nuclear library upon +request. A DVD of the data can be requested here_. To use this data with OpenMC, +the following steps must be taken: + +1. Copy and unzip the data on the DVD to a directory on your computer. +2. In the root directory, a file named ``xsdir``, or some variant thereof, + should be present. This file contains a listing of all the cross sections and + is used by MCNP. This file should be converted to a ``cross_sections.xml`` + file for use with OpenMC. A Python script is provided in the OpenMC + distribution for this purpose: + + .. code-block:: sh + + openmc/src/utils/convert_xsdir.py xsdir31 cross_sections.xml + +3. In the converted ``cross_sections.xml`` file, change the contents of the + element to the absolute path of the directory containing the + actual ACE files. +4. Additionally, you may need to change any occurrences of upper-case "ACE" + within the ``cross_sections.xml`` file to lower-case. +5. Either set the :ref:`cross_sections` in a settings.xml file or the + :envvar:`CROSS_SECTIONS` environment variable to the absolute path of the + ``cross_sections.xml`` file. + +Using Cross Sections from MCNP +------------------------------ To use cross sections distributed with MCNP, change the element in the ``cross_sections.xml`` file in the root directory of the OpenMC distribution -to the location of the MCNP cross-sections. Then, either set the +to the location of the MCNP cross sections. Then, either set the :ref:`cross_sections` in a settings.xml file or the :envvar:`CROSS_SECTIONS` environment variable to the absolute path of the ``cross_sections.xml`` file. -Similarly, to use cross-sections distributed with Serpent, change the - element in the ``cross_sections_serpent.xml`` file in the root -directory of the OpenMC distribution to the location of the Serpent -cross-sections. Then, either set the :ref:`cross_sections` in a settings.xml -file or the :envvar:`CROSS_SECTIONS` environment variable to the absolute path -of the ``cross_sections_serpent.xml`` file. +Using Cross Sections from Serpent +--------------------------------- +To use cross sections distributed with Serpent, change the element +in the ``cross_sections_serpent.xml`` file in the root directory of the OpenMC +distribution to the location of the Serpent cross sections. Then, either set the +:ref:`cross_sections` in a settings.xml file or the :envvar:`CROSS_SECTIONS` +environment variable to the absolute path of the ``cross_sections_serpent.xml`` +file. + +.. _NEA: http://www.oecd-nea.org +.. _JEFF: http://www.oecd-nea.org/dbdata/jeff/ +.. _here: http://www.oecd-nea.org/dbdata/pubs/jeff312-cd.html .. _MCNP: http://mcnp.lanl.gov .. _Serpent: http://montecarlo.vtt.fi From 54ad656d19f9047f6bb8dc1eeb8a2f998bd3fcc7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 9 Jan 2013 18:15:40 -0500 Subject: [PATCH 10/12] Updated Windows install instructions. --- docs/source/usersguide/install.rst | 61 +++++++++++++++++++++++++++--- 1 file changed, 55 insertions(+), 6 deletions(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 05f08e064..971c3bb01 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -126,11 +126,33 @@ This will build an executable named ``openmc``. Compiling on Windows -------------------- -To compile OpenMC on a Windows operating system, you will need to first install -Cygwin_, a Linux-like environment for Windows. When configuring Cygwin, make -sure you install both the gfortran compiler as well as git. Once you have -obtained the source code, run the following commands from within the source code -root directory: +Using Cygwin +------------ + +One option for compiling OpenMC on a Windows operating system is to use Cygwin_, +a Linux-like environment for Windows. You will need to first `install +Cygwin`_. When you are asked to select packages, make sure the following are +selected: + +* Devel: gcc4-core +* Devel: gcc4-fortran +* Devel: make + +If you plan on obtaining the source code directly using git, select the +following packages: + +* Devel: git +* Devel: git-completion (Optional) +* Devel: gitk (Optional) + +In order to use the Python scripts provided with OpenMC, you will also need to +install Python. This can be done within Cygwin or directly in Windows. To +install within Cygwin, select the following packages: + +* Python: python (Version > 2.7 recommended) + +Once you have obtained the source code, run the following commands from within +the source code root directory: .. code-block:: sh @@ -139,7 +161,34 @@ root directory: This will build an executable named ``openmc``. -.. _Cygwin: http://www.cygwin.com/ +.. _Cygwin: http://cygwin.com/ +.. _install Cygwin: http://cygwin.com/setup.exe + +Using MinGW +----------- + +An alternate option for installing OpenMC on Windows is using MinGW_, which +stands for Minimalist GNU for Windows. An executable for installing the MinGW +distribution is available on SourceForge_. When installing MinGW, make sure the +following components are selected: + +* MinGW Compiler Suite: Fortran Compiler +* MSYS Basic System + +Once MinGW is installed, copy the OpenMC source distribution to your MinGW home +directory (usually C:\\MinGW\\msys\\1.0\\home\\YourUsername). Once you have +the source code in place, run the following commands from within the MinGW shell +in the root directory of the OpenMC distribution: + +.. code-block:: sh + + cd src + make + +This will build an executable named ``openmc``. + +.. _MinGW: http://www.mingw.org +.. _SourceForge: http://sourceforge.net/projects/mingw --------------------------- Cross Section Configuration From 221f48af171dec61338c71794f387252eb246a6a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 15 Jan 2013 19:21:13 -0500 Subject: [PATCH 11/12] Added constants.o in OBJECTS. --- src/OBJECTS | 1 + 1 file changed, 1 insertion(+) diff --git a/src/OBJECTS b/src/OBJECTS index 1a0aac382..2405d71a2 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -13,6 +13,7 @@ cmfd_output.o \ cmfd_power_solver.o \ cmfd_prod_operator.o \ cmfd_snes_solver.o \ +constants.o \ cross_section.o \ datatypes.o \ datatypes_header.o \ From 63bfd25d34ecffd250825cf1628d014dbe9d34f0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 16 Jan 2013 15:53:10 -0500 Subject: [PATCH 12/12] Search all coordinate levels for cell filters. Closes #139. --- src/DEPENDENCIES | 1 + src/tally.F90 | 48 ++++++++++++++++++++++++++++++++---------------- 2 files changed, 33 insertions(+), 16 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 1b9b6258a..c9b33be4b 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -345,6 +345,7 @@ tally.o: math.o tally.o: mesh.o tally.o: mesh_header.o tally.o: output.o +tally.o: particle_header.o tally.o: search.o tally.o: string.o tally.o: tally_header.o diff --git a/src/tally.F90 b/src/tally.F90 index 0ce310a30..2305a4594 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,18 +1,20 @@ module tally - use ace_header, only: Reaction + use ace_header, only: Reaction use constants use datatypes_header, only: ListInt - use error, only: fatal_error + use error, only: fatal_error use global - use math, only: t_percentile, calc_pn - use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, & - mesh_indices_to_bin, mesh_intersects - use mesh_header, only: StructuredMesh - use output, only: header - use search, only: binary_search - use string, only: to_str - use tally_header, only: TallyResult, TallyMapItem, TallyMapElement + use math, only: t_percentile, calc_pn + use mesh, only: get_mesh_bin, bin_to_mesh_indices, & + get_mesh_indices, mesh_indices_to_bin, & + mesh_intersects + use mesh_header, only: StructuredMesh + use output, only: header + use particle_header, only: LocalCoord + use search, only: binary_search + use string, only: to_str + use tally_header, only: TallyResult, TallyMapItem, TallyMapElement #ifdef MPI use mpi @@ -963,6 +965,7 @@ contains type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() type(Material), pointer :: mat => null() + type(LocalCoord), pointer :: coord => null() t => tallies(i_tally) t % matching_bins = 1 @@ -1010,9 +1013,15 @@ contains case (FILTER_CELL) ! determine next cell bin - ! TODO: Account for cells in multiple levels when performing this filter - t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + coord => p % coord0 + do while(associated(coord)) + position(FILTER_CELL) = 0 + t % matching_bins(i) = get_next_bin(FILTER_CELL, & + coord % cell, i_tally) + if (t % matching_bins(i) /= NO_BIN_FOUND) exit + coord => coord % next + end do + nullify(coord) case (FILTER_CELLBORN) ! determine next cellborn bin @@ -1237,6 +1246,7 @@ contains real(8) :: E ! particle energy type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() + type(LocalCoord), pointer :: coord => null() found_bin = .true. t => tallies(i_tally) @@ -1264,9 +1274,15 @@ contains case (FILTER_CELL) ! determine next cell bin - ! TODO: Account for cells in multiple levels when performing this filter - t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + coord => p % coord0 + do while(associated(coord)) + position(FILTER_CELL) = 0 + t % matching_bins(i) = get_next_bin(FILTER_CELL, & + coord % cell, i_tally) + if (t % matching_bins(i) /= NO_BIN_FOUND) exit + coord => coord % next + end do + nullify(coord) case (FILTER_CELLBORN) ! determine next cellborn bin