crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Seem to have all the bugs fixed for the harmonic tallying; next up will be documentation.

Browse source
This commit is contained in:
Adam Nelson 2014-04-04 19:37:26 -04:00
parent ba39ff4b63
commit 2703239c16
6 changed files with 253 additions and 257 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/constants.F90
View file

@ -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

103
src/input_xml.F90
View file

@ -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 <scores>: " // &
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

3
src/math.F90
View file

@ -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)

81
src/output.F90
View file

@ -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))

315
src/tally.F90
View file

@ -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)

4
src/utils/statepoint.py
View file

@ -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