Added tallies for all types, now to debug.

This commit is contained in:
Adam Nelson 2014-04-04 14:16:26 -04:00
parent 26edd77b45
commit ba39ff4b63

View file

@ -4,7 +4,7 @@ module tally
use constants
use error, only: fatal_error
use global
use math, only: t_percentile, calc_pn
use math, only: t_percentile, calc_pn, calc_rn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, &
get_mesh_indices, mesh_indices_to_bin, &
mesh_intersects_2d, mesh_intersects_3d
@ -41,7 +41,9 @@ contains
integer :: i_tally
integer :: j ! loop index for scoring bins
integer :: k ! loop index for nuclide bins
integer :: n ! loop index for scattering order
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
@ -141,6 +143,36 @@ contains
score = last_wgt / material_xs % total
case (SCORE_FLUX_YN)
! All events score to a flux bin. We actually use a collision
! estimator since there is no way to count 'events' exactly for
! the flux
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
! 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 + &
score * calc_rn(n, p % last_uvw)
!$omp end critical
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_TOTAL)
! All events will score to the total reaction rate. We can just
! use the weight of the particle entering the collision as the
@ -154,6 +186,43 @@ contains
score = last_wgt
end if
case (SCORE_TOTAL_YN)
! All events will score to the total reaction rate. We can just
! use the weight of the particle entering the collision as the
! score
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
! 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 + &
score * calc_rn(n, p % last_uvw)
!$omp end critical
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_SCATTER)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
@ -210,6 +279,37 @@ contains
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_SCATTER_YN)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
j = j + t % moment_order(j)
cycle SCORE_LOOP
end if
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
! 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 of the scattering moment
score = last_wgt * calc_pn(n, mu)
! 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 % last_uvw)
!$omp end critical
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_NU_SCATTER_N)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
@ -246,6 +346,37 @@ contains
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_NU_SCATTER_YN)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
j = j + t % moment_order(j)
cycle SCORE_LOOP
end if
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
! 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 of the scattering moment
score = wgt * calc_pn(n, mu)
! 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 % last_uvw)
!$omp end critical
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_TRANSPORT)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
@ -482,6 +613,10 @@ contains
integer :: k ! loop index for nuclide bins
integer :: l ! loop index for nuclides in material
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)
integer :: i_nuc ! index in nuclides array (from material)
@ -562,10 +697,15 @@ contains
end if
! Determine score for each bin
SCORE_LOOP: do j = 1, t % n_score_bins
j = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
j = j + 1
! determine what type of score bin
score_bin = t % score_bins(j)
! determine scoring bin index
score_index = (k - 1)*t % n_score_bins + j
if (i_nuclide > 0) then
! ================================================================
! DETERMINE NUCLIDE CROSS SECTION
@ -575,11 +715,66 @@ contains
! For flux, we need no cross section
score = flux
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
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, num_n
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_TOTAL)
! Total cross section is pre-calculated
score = micro_xs(i_nuclide) % total * &
atom_density * flux
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
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_SCATTER)
! Scattering cross section is pre-calculated
score = (micro_xs(i_nuclide) % total - &
@ -659,10 +854,64 @@ contains
! For flux, we need no cross section
score = flux
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
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, num_n
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_TOTAL)
! Total cross section is pre-calculated
score = material_xs % total * flux
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
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_SCATTER)
! Scattering cross section is pre-calculated
score = (material_xs % total - material_xs % absorption) * flux
@ -739,9 +988,6 @@ contains
end select
end if
! Determine scoring bin index
score_index = (k - 1)*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
@ -782,6 +1028,10 @@ 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 :: 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
integer :: score_index ! scoring bin index
@ -812,18 +1062,77 @@ contains
atom_density = mat % atom_density(i)
! Loop over score types for each bin
SCORE_LOOP: do j = 1, t % n_score_bins
j = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
j = j + 1
! determine what type of score bin
score_bin = t % score_bins(j)
! Determine scoring bin index based on what the index of the nuclide
! is in the nuclides array
score_index = (i_nuclide - 1)*t % n_score_bins + j
! Determine macroscopic nuclide cross section
select case(score_bin)
case (SCORE_FLUX)
score = flux
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
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, num_n
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_TOTAL)
score = micro_xs(i_nuclide) % total * atom_density * flux
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
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * atom_density * flux
@ -883,10 +1192,6 @@ contains
end if
end select
! Determine scoring bin index based on what the index of the nuclide
! is in the nuclides array
score_index = (i_nuclide - 1)*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
@ -901,18 +1206,78 @@ contains
! SCORE TOTAL MATERIAL REACTION RATES
! Loop over score types for each bin
MATERIAL_SCORE_LOOP: do j = 1, t % n_score_bins
j = 0
MATERIAL_SCORE_LOOP: do q = 1, t % n_user_score_bins
j = j + 1
! determine what type of score bin
score_bin = t % score_bins(j)
! Determine scoring bin index based on what the index of the nuclide
! is in the nuclides array
score_index = n_nuclides_total*t % n_score_bins + j
! Determine macroscopic material cross section
select case(score_bin)
case (SCORE_FLUX)
score = flux
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
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, num_n
! 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
end do
j = j + t % moment_order(j)
cycle MATERIAL_SCORE_LOOP
case (SCORE_TOTAL)
score = material_xs % total * flux
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
! 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
end do
j = j + t % moment_order(j)
cycle MATERIAL_SCORE_LOOP
case (SCORE_SCATTER)
score = (material_xs % total - material_xs % absorption) * flux
@ -983,10 +1348,6 @@ contains
end if
end select
! Determine scoring bin index based on what the index of the nuclide
! is in the nuclides array
score_index = n_nuclides_total*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &
@ -1013,6 +1374,10 @@ contains
integer :: j ! loop index for direction
integer :: k ! loop index for mesh cell crossings
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
integer :: ijk_cross(3) ! indices of mesh cell crossed
@ -1233,18 +1598,79 @@ contains
end if
! Determine score for each bin
SCORE_LOOP: do j = 1, t % n_score_bins
j = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
j = j + 1
! determine what type of score bin
score_bin = t % score_bins(j)
! Determine scoring bin index
score_index = (b - 1)*t % n_score_bins + j
if (i_nuclide > 0) then
! Determine macroscopic nuclide cross section
select case(score_bin)
case (SCORE_FLUX)
score = flux
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
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, num_n
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_TOTAL)
score = micro_xs(i_nuclide) % total * &
atom_density * flux
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
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * &
@ -1273,8 +1699,63 @@ contains
select case(score_bin)
case (SCORE_FLUX)
score = flux
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
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, num_n
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_TOTAL)
score = material_xs % total * flux
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
! 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
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
case (SCORE_SCATTER)
score = (material_xs % total - material_xs % absorption) * flux
case (SCORE_ABSORPTION)
@ -1294,9 +1775,6 @@ contains
end select
end if
! Determine scoring bin index
score_index = (b - 1)*t % n_score_bins + j
! Add score to tally
!$omp critical
t % results(score_index, filter_index) % value = &