From 2703239c168f74e8fa264ad6c39e0755d28444e8 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 4 Apr 2014 19:37:26 -0400 Subject: [PATCH] Seem to have all the bugs fixed for the harmonic tallying; next up will be documentation. --- src/constants.F90 | 4 + src/input_xml.F90 | 103 +++++++------ src/math.F90 | 3 +- src/output.F90 | 81 ++++++----- src/tally.F90 | 315 ++++++++++++++++++---------------------- src/utils/statepoint.py | 4 +- 6 files changed, 253 insertions(+), 257 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index e39a5344ff..a80b26d278 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -302,6 +302,10 @@ module constants MOMENT_N_STRS(2) = (/ "scatter- ", & "nu-scatter- "/) + ! Location in MOMENT_STRS where the YN data begins + integer, parameter :: YN_LOC = 3 + + ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 724d528719..7f5afb0f75 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1708,9 +1708,10 @@ contains integer :: n ! size of arrays in mesh specification integer :: n_words ! number of words read integer :: n_filters ! number of filters - integer :: n_new ! number of new scores to add based on Pn tally - integer :: n_scores ! number of tot scores after adjusting for Pn tally - integer :: n_order ! Scattering order requested + integer :: n_new ! number of new scores to add based on Yn/Pn tally + integer :: n_scores ! number of tot scores after adjusting for Yn/Pn tally + integer :: n_bins ! Total new bins for this score + integer :: n_order ! Moment order requested integer :: n_order_pos ! Position of Scattering order in score name string integer :: MT ! user-specified MT for score integer :: iarray3(3) ! temporary integer array @@ -2265,7 +2266,7 @@ contains n_new = 0 do j = 1, n_words call lower_case(sarray(j)) - ! Find if scores(j) is of the form 'moment-p' present in + ! Find if scores(j) is of the form 'moment-p' or 'moment-y' present in ! MOMENT_STRS(:) ! If so, check the order, store if OK, then reset the number to 'n' score_name = trim(sarray(j)) @@ -2286,7 +2287,14 @@ contains sarray(j) = trim(MOMENT_STRS(imomstr)) // & trim(to_str(MAX_ANG_ORDER)) end if - n_new = n_new + n_order + ! Find total number of bins for this case + if (imomstr >= YN_LOC) then + n_bins = (n_order + 1)**2 + else + n_bins = n_order + 1 + end if + ! We subtract one since n_words already included + n_new = n_new + n_bins - 1 exit end if end do @@ -2320,28 +2328,37 @@ contains n_order = MAX_ANG_ORDER end if score_name = trim(MOMENT_STRS(imomstr)) // "n" - exit - end if - end do - do imomstr = 1, size(MOMENT_N_STRS) - if (starts_with(score_name,trim(MOMENT_N_STRS(imomstr)))) then - n_order_pos = scan(score_name,'0123456789') - n_order = int(str_to_int( & - score_name(n_order_pos:(len_trim(score_name)))),4) - if (n_order > MAX_ANG_ORDER) then - ! User requested too many orders; throw a warning and set to the - ! maximum order. - ! The above scheme will essentially take the absolute value - message = "Invalid scattering order of " // trim(to_str(n_order)) // & - " requested. Setting to the maximum permissible value, " // & - trim(to_str(MAX_ANG_ORDER)) - call warning() - n_order = MAX_ANG_ORDER + ! Find total number of bins for this case + if (imomstr >= YN_LOC) then + n_bins = (n_order + 1)**2 + else + n_bins = n_order + 1 end if - score_name = trim(MOMENT_N_STRS(imomstr)) // "n" exit end if end do + ! Now check the Moment_N_Strs, but only if we werent successful above + if (imomstr > size(MOMENT_STRS)) then + do imomstr = 1, size(MOMENT_N_STRS) + if (starts_with(score_name,trim(MOMENT_N_STRS(imomstr)))) then + n_order_pos = scan(score_name,'0123456789') + n_order = int(str_to_int( & + score_name(n_order_pos:(len_trim(score_name)))),4) + if (n_order > MAX_ANG_ORDER) then + ! User requested too many orders; throw a warning and set to the + ! maximum order. + ! The above scheme will essentially take the absolute value + message = "Invalid scattering order of " // trim(to_str(n_order)) // & + " requested. Setting to the maximum permissible value, " // & + trim(to_str(MAX_ANG_ORDER)) + call warning() + n_order = MAX_ANG_ORDER + end if + score_name = trim(MOMENT_N_STRS(imomstr)) // "n" + exit + end if + end do + end if select case (trim(score_name)) case ('flux') @@ -2370,9 +2387,9 @@ contains call fatal_error() end if - t % score_bins(j : j + n_order) = SCORE_FLUX_YN - t % moment_order(j : j + n_order) = (n_order + 1)**2 - j = j + (n_order + 1)**2 + t % score_bins(j : j + n_bins - 1) = SCORE_FLUX_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 case ('total') t % score_bins(j) = SCORE_TOTAL @@ -2389,9 +2406,9 @@ contains call fatal_error() end if - t % score_bins(j : j + n_order) = SCORE_TOTAL_YN - t % moment_order(j : j + n_order) = (n_order + 1)**2 - j = j + (n_order + 1)**2 + t % score_bins(j : j + n_bins - 1) = SCORE_TOTAL_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 case ('scatter') t % score_bins(j) = SCORE_SCATTER @@ -2424,30 +2441,30 @@ contains case ('scatter-pn') t % estimator = ESTIMATOR_ANALOG ! Setup P0:Pn - t % score_bins(j : j + n_order) = SCORE_SCATTER_PN - t % moment_order(j : j + n_order) = n_order - j = j + n_order + t % score_bins(j : j + n_bins - 1) = SCORE_SCATTER_PN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 case ('nu-scatter-pn') t % estimator = ESTIMATOR_ANALOG ! Setup P0:Pn - t % score_bins(j : j + n_order) = SCORE_NU_SCATTER_PN - t % moment_order(j : j + n_order) = n_order - j = j + n_order + t % score_bins(j : j + n_bins - 1) = SCORE_NU_SCATTER_PN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 case ('scatter-yn') t % estimator = ESTIMATOR_ANALOG ! Setup P0:Pn - t % score_bins(j : j + n_order) = SCORE_SCATTER_YN - t % moment_order(j : j + n_order) = (n_order + 1)**2 - j = j + (n_order + 1)**2 + t % score_bins(j : j + n_bins - 1) = SCORE_SCATTER_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 case ('nu-scatter-yn') t % estimator = ESTIMATOR_ANALOG ! Setup P0:Pn - t % score_bins(j : j + n_order) = SCORE_NU_SCATTER_YN - t % moment_order(j : j + n_order) = (n_order + 1)**2 - j = j + (n_order + 1)**2 + t % score_bins(j : j + n_bins - 1) = SCORE_NU_SCATTER_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 case('transport') t % score_bins(j) = SCORE_TRANSPORT @@ -2558,14 +2575,14 @@ contains t % score_bins(j) = MT else message = "Invalid MT on : " // & - trim(sarray(j)) + trim(sarray(l)) call fatal_error() end if else ! Specified score was not an integer message = "Unknown scoring function: " // & - trim(sarray(j)) + trim(sarray(l)) call fatal_error() end if diff --git a/src/math.F90 b/src/math.F90 index aaaa6ac6e9..b3046ccc4a 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -181,7 +181,8 @@ contains if (uvw(1) == ZERO) then phi = ZERO else - phi = atan(uvw(2) / uvw(1)) +! phi = atan(uvw(2) / uvw(1)) + phi = atan2(uvw(2), uvw(1)) end if w2m1 = (ONE - w**2) diff --git a/src/output.F90 b/src/output.F90 index 6c431b2b95..f445b17eaf 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1657,7 +1657,8 @@ contains integer :: score_index ! scoring bin index integer :: i_nuclide ! index in nuclides array integer :: i_listing ! index in xs_listings array - integer :: n_order ! loop index for scattering orders + integer :: n_order ! loop index for moment orders + integer :: nm_order ! loop index for Ynm moment orders real(8) :: t_value ! t-values for confidence intervals real(8) :: alpha ! significance level for CI character(MAX_FILE_LEN) :: filename ! name of output file @@ -1852,7 +1853,7 @@ contains to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) end do - k = k + n_order - 1 + k = k + t % moment_order(k) else if (t % score_bins(k) == SCORE_SCATTER_YN) then score_name = "Scattering Rate" write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & @@ -1860,15 +1861,17 @@ contains to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) do n_order = 1, t % moment_order(k) - score_index = score_index + 1 - score_name = 'Y' // trim(to_str(n_order)) // & - ' Scattering Moment' - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + do nm_order = -n_order, n_order + score_index = score_index + 1 + score_name = 'Y' // trim(to_str(n_order)) // ',' // & + trim(to_str(nm_order)) // ' Scattering Moment' + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) % sum_sq)) + end do end do - k = k + n_order - 1 + k = k + (t % moment_order(k) + 1)**2 - 1 else if (t % score_bins(k) == SCORE_NU_SCATTER_N) then if (t % moment_order(k) == 0) then score_name = "Nu-Scattering Rate" @@ -1895,7 +1898,7 @@ contains to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) end do - k = k + n_order - 1 + k = k + t % moment_order(k) else if (t % score_bins(k) == SCORE_NU_SCATTER_YN) then score_name = "Nu-Scattering Rate" write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & @@ -1903,15 +1906,17 @@ contains to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) do n_order = 1, t % moment_order(k) - score_index = score_index + 1 - score_name = 'Y' // trim(to_str(n_order)) // & - ' Nu-Scattering Moment' - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + do nm_order = -n_order, n_order + score_index = score_index + 1 + score_name = 'Y' // trim(to_str(n_order)) // ',' // & + trim(to_str(nm_order)) // ' Nu-Scatter Moment' + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) % sum_sq)) + end do end do - k = k + n_order - 1 + k = k + (t % moment_order(k) + 1)**2 - 1 else if (t % score_bins(k) == SCORE_FLUX_YN) then score_name = "Flux Moment" write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & @@ -1919,31 +1924,35 @@ contains to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) do n_order = 1, t % moment_order(k) - score_index = score_index + 1 - score_name = 'Y' // trim(to_str(n_order)) // & - ' Flux Moment' - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + do nm_order = -n_order, n_order + score_index = score_index + 1 + score_name = 'Y' // trim(to_str(n_order)) // ',' // & + trim(to_str(nm_order)) // ' Flux Moment' + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) % sum_sq)) + end do end do - k = k + n_order - 1 - else if (t % score_bins(k) == SCORE_NU_SCATTER_YN) then + k = k + (t % moment_order(k) + 1)**2 - 1 + else if (t % score_bins(k) == SCORE_TOTAL_YN) then score_name = "Total Reaction Moment" write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) do n_order = 1, t % moment_order(k) - score_index = score_index + 1 - score_name = 'Y' // trim(to_str(n_order)) // & - ' Total Reaction Moment' - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + do nm_order = -n_order, n_order + score_index = score_index + 1 + score_name = 'Y' // trim(to_str(n_order)) // ',' // & + trim(to_str(nm_order)) // ' Total Reaction Moment' + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) % sum_sq)) + end do end do - k = k + n_order - 1 + k = k + (t % moment_order(k) + 1)**2 - 1 else if (t % score_bins(k) > 0) then score_name = reaction_name(t % score_bins(k)) diff --git a/src/tally.F90 b/src/tally.F90 index b488725415..edb4d72aba 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -43,7 +43,6 @@ contains integer :: k ! loop index for nuclide bins integer :: n ! loop index for legendre order integer :: num_nm ! Number of N,M orders in harmonic - integer :: num_n ! Number of N orders in harmonic integer :: l ! scoring bin loop index, allowing for changing ! position during the loop integer :: filter_index ! single index for single bin @@ -150,27 +149,26 @@ contains score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! get the score + score = last_wgt / material_xs % total + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! get the score - score = last_wgt / material_xs % total ! multiply score by the angular flux moments and store !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & score * calc_rn(n, p % last_uvw) !$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TOTAL) @@ -193,34 +191,32 @@ contains score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! get the score + if (survival_biasing) then + ! We need to account for the fact that some weight was already + ! absorbed + score = last_wgt + p % absorb_wgt + else + score = last_wgt + end if + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! get the score - if (survival_biasing) then - ! We need to account for the fact that some weight was already - ! absorbed - score = last_wgt + p % absorb_wgt - else - score = last_wgt - end if - ! multiply score by the angular flux moments and store !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & score * calc_rn(n, p % last_uvw) !$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_SCATTER) @@ -288,11 +284,10 @@ contains score_index = score_index - 1 ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) @@ -302,12 +297,12 @@ contains ! multiply score by the angular flux moments and store !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & score * calc_rn(n, p % last_uvw) !$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_NU_SCATTER_N) @@ -355,11 +350,10 @@ contains score_index = score_index - 1 ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) @@ -369,12 +363,12 @@ contains ! multiply score by the angular flux moments and store !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & score * calc_rn(n, p % last_uvw) !$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TRANSPORT) @@ -615,7 +609,6 @@ contains integer :: m ! loop index for reactions integer :: n ! loop index for legendre order integer :: num_nm ! Number of N,M orders in harmonic - integer :: num_n ! Number of N orders in harmonic integer :: q ! loop index for scoring bins integer :: filter_index ! single index for single bin integer :: i_nuclide ! index in nuclides array (from bins) @@ -718,28 +711,26 @@ contains case (SCORE_FLUX_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! For flux, we need no cross section + score = flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! For flux, we need no cross section - score = flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TOTAL) @@ -750,29 +741,23 @@ contains case (SCORE_TOTAL_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! Total cross section is pre-calculated - score = micro_xs(i_nuclide) % total * & - atom_density * flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_SCATTER) @@ -857,28 +842,26 @@ contains case (SCORE_FLUX_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! For flux, we need no cross section + score = flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! For flux, we need no cross section - score = flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TOTAL) @@ -888,28 +871,26 @@ contains case (SCORE_TOTAL_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! Total cross section is pre-calculated + score = material_xs % total * flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! Total cross section is pre-calculated - score = material_xs % total * flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_SCATTER) @@ -1030,7 +1011,6 @@ contains integer :: m ! loop index for reactions in nuclide integer :: n ! loop index for legendre order integer :: num_nm ! Number of N,M orders in harmonic - integer :: num_n ! Number of N orders in harmonic integer :: q ! loop index for scoring bins integer :: i_nuclide ! index in nuclides array integer :: score_bin ! type of score, e.g. SCORE_FLUX @@ -1080,28 +1060,26 @@ contains case (SCORE_FLUX_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! For flux, we need no cross section + score = flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! For flux, we need no cross section - score = flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TOTAL) @@ -1110,27 +1088,25 @@ contains case (SCORE_TOTAL_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + score = micro_xs(i_nuclide) % total * atom_density * flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - score = micro_xs(i_nuclide) % total * atom_density * flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_SCATTER) @@ -1224,28 +1200,26 @@ contains case (SCORE_FLUX_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! For flux, we need no cross section + score = flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! For flux, we need no cross section - score = flux - ! multiply score by the angular flux moments and store !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) !$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle MATERIAL_SCORE_LOOP case (SCORE_TOTAL) @@ -1254,28 +1228,26 @@ contains case (SCORE_TOTAL_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! Total cross section is pre-calculated + score = material_xs % total * flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! Total cross section is pre-calculated - score = material_xs % total * flux - ! multiply score by the angular flux moments and store !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) !$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle MATERIAL_SCORE_LOOP case (SCORE_SCATTER) @@ -1376,7 +1348,6 @@ contains integer :: b ! loop index for nuclide bins integer :: n ! loop index for legendre order integer :: num_nm ! Number of N,M orders in harmonic - integer :: num_n ! Number of N orders in harmonic integer :: q ! loop index for scoring bins integer :: ijk0(3) ! indices of starting coordinates integer :: ijk1(3) ! indices of ending coordinates @@ -1616,27 +1587,25 @@ contains case (SCORE_FLUX_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + score = flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - score = flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TOTAL) @@ -1646,29 +1615,27 @@ contains case (SCORE_TOTAL_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! Total cross section is pre-calculated + score = micro_xs(i_nuclide) % total * & + atom_density * flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! Total cross section is pre-calculated - score = micro_xs(i_nuclide) % total * & - atom_density * flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_SCATTER) @@ -1703,27 +1670,25 @@ contains case (SCORE_FLUX_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + score = flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - score = flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TOTAL) @@ -1732,28 +1697,26 @@ contains case (SCORE_TOTAL_YN) score_index = score_index - 1 - ! Calculate the number of moments from t % moment_order - num_n = int(sqrt(real(t % moment_order(j),8))) - 1 + ! Total cross section is pre-calculated + score = material_xs % total * flux + num_nm = 1 ! Find the order for a collection of requested moments ! and store the moment contribution of each - do n = 0, num_n + do n = 0, t % moment_order(j) ! determine scoring bin index score_index = score_index + num_nm ! Update number of total n,m bins for this n (m = [-n: n]) num_nm = 2 * n + 1 - ! Total cross section is pre-calculated - score = material_xs % total * flux - ! multiply score by the angular flux moments and store - !$omp critical - t % results(score_index: score_index + num_nm, filter_index) % value = & - t % results(score_index: score_index + num_nm, filter_index) % value + & - score * calc_rn(n, p % coord0 % uvw) - !$omp end critical +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical end do - j = j + t % moment_order(j) + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_SCATTER) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 888e8bf2eb..3799aa7253 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -6,6 +6,8 @@ from collections import OrderedDict import numpy as np import scipy.stats +REVISION_STATEPOINT = 12 + filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', 5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'} @@ -157,7 +159,7 @@ class StatePoint(object): # Read statepoint revision self.revision = self._get_int(path='revision')[0] - if self.revision != 11: + if self.revision != REVISION_STATEPOINT: raise Exception('Statepoint Revision is not consistent.') # Read OpenMC version