diff --git a/src/tally.F90 b/src/tally.F90 index a9acb95d88..b488725415 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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 = &