From bc387186ba6b08ce87685dbd7116e29310ce0f74 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 23 Aug 2016 09:41:53 -0400 Subject: [PATCH 01/12] changed partial current tallies to net current tallies --- examples/xml/pincell/geometry.xml | 14 +- examples/xml/pincell/settings.xml | 12 +- examples/xml/pincell/tallies.xml | 11 +- src/cmfd_data.F90 | 28 +-- src/cmfd_input.F90 | 5 +- src/constants.F90 | 12 +- src/input_xml.F90 | 5 +- src/output.F90 | 24 +- src/tally.F90 | 351 +++++++++++++++++++++++++----- src/trigger.F90 | 12 +- 10 files changed, 356 insertions(+), 118 deletions(-) diff --git a/examples/xml/pincell/geometry.xml b/examples/xml/pincell/geometry.xml index f67f9e74c..f4e1bf1dd 100644 --- a/examples/xml/pincell/geometry.xml +++ b/examples/xml/pincell/geometry.xml @@ -8,16 +8,16 @@ --> - - - + + + - - - - + + + + diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml index 443af9cde..9402b979a 100644 --- a/examples/xml/pincell/settings.xml +++ b/examples/xml/pincell/settings.xml @@ -3,7 +3,7 @@ - 100 + 20 10 1000 @@ -14,8 +14,8 @@ - -0.62992 -0.62992 -1. - 0.62992 0.62992 1. + -8.62992 -8.62992 -1. + 8.62992 8.62992 1. @@ -24,9 +24,9 @@ bounds for a mesh over which the Shannon entropy should be calculated. The extent in the z direction is made arbitrarily large. --> - -0.39218 -0.39218 -1.e50 - 0.39218 0.39218 1.e50 + -3.39218 -3.39218 -1.e50 + 3.39218 3.39218 1.e50 10 10 1 - \ No newline at end of file + diff --git a/examples/xml/pincell/tallies.xml b/examples/xml/pincell/tallies.xml index 73242b913..2e39c083d 100644 --- a/examples/xml/pincell/tallies.xml +++ b/examples/xml/pincell/tallies.xml @@ -2,9 +2,9 @@ - 100 100 1 - -0.62992 -0.62992 -1.e50 - 0.62992 0.62992 1.e50 + 2 2 1 + -8.62992 -8.62992 -1.e50 + 8.62992 8.62992 1.e50 @@ -13,4 +13,9 @@ flux fission nu-fission + + + current + + diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index d74c1b17e..0a1490bb9 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -51,8 +51,8 @@ contains subroutine compute_xs() use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, & - FILTER_SURFACE, OUT_LEFT, OUT_RIGHT, OUT_BACK, & - OUT_FRONT, OUT_BOTTOM, OUT_TOP, CMFD_NOACCEL, & + FILTER_SURFACE, LEFT, RIGHT, BACK, & + FRONT, BOTTOM, TOP, CMFD_NOACCEL, & ZERO, ONE, TINY_BIT use error, only: fatal_error use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes,& @@ -237,7 +237,7 @@ contains ! Left surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_LEFT + matching_bins(i_filter_surf) = LEFT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t%stride) + 1 ! outgoing cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum @@ -245,7 +245,7 @@ contains if (i > 1) then matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i-1, j, k /)) - matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_surf) = RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! incoming cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum @@ -255,7 +255,7 @@ contains if (i < nx) then matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i+1, j, k /) ) - matching_bins(i_filter_surf) = OUT_LEFT + matching_bins(i_filter_surf) = LEFT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! incoming cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum @@ -263,7 +263,7 @@ contains matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) ) - matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_surf) = RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! outgoing cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum @@ -271,7 +271,7 @@ contains ! Back surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_BACK + matching_bins(i_filter_surf) = BACK score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! outgoing cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum @@ -279,7 +279,7 @@ contains if (j > 1) then matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j-1, k /)) - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! incoming cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum @@ -289,7 +289,7 @@ contains if (j < ny) then matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j+1, k /)) - matching_bins(i_filter_surf) = OUT_BACK + matching_bins(i_filter_surf) = BACK score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! incoming cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum @@ -297,7 +297,7 @@ contains matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! outgoing cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum @@ -305,7 +305,7 @@ contains ! Bottom surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_BOTTOM + matching_bins(i_filter_surf) = BOTTOM score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! outgoing cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum @@ -313,7 +313,7 @@ contains if (k > 1) then matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k-1 /)) - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! incoming cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum @@ -323,7 +323,7 @@ contains if (k < nz) then matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k+1 /)) - matching_bins(i_filter_surf) = OUT_BOTTOM + matching_bins(i_filter_surf) = BOTTOM score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! incoming cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum @@ -331,7 +331,7 @@ contains matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 ! outgoing cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 5fe44ad43..77e8f4348 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -534,10 +534,9 @@ contains filt % n_bins = 2 * m % n_dimension allocate(filt % surfaces(2 * m % n_dimension)) if (m % n_dimension == 2) then - filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT /) + filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT /) elseif (m % n_dimension == 3) then - filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT, & - OUT_BOTTOM, OUT_TOP /) + filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT, BOTTOM, TOP /) end if end select t % find_filter(FILTER_SURFACE) = n_filters diff --git a/src/constants.F90 b/src/constants.F90 index 6d5f2fc10..ceb68b9e2 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -359,12 +359,12 @@ module constants ! Tally surface current directions integer, parameter :: & - OUT_LEFT = 1, & ! x min - OUT_RIGHT = 2, & ! x max - OUT_BACK = 3, & ! y min - OUT_FRONT = 4, & ! y max - OUT_BOTTOM = 5, & ! z min - OUT_TOP = 6 ! z max + LEFT = 1, & ! x min + RIGHT = 2, & ! x max + BACK = 3, & ! y min + FRONT = 4, & ! y max + BOTTOM = 5, & ! z min + TOP = 6 ! z max ! Tally trigger types and threshold integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a60d8196e..8de198386 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3701,10 +3701,9 @@ contains filt % n_bins = 2 * m % n_dimension allocate(filt % surfaces(2 * m % n_dimension)) if (m % n_dimension == 2) then - filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT /) + filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT /) elseif (m % n_dimension == 3) then - filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT,& - OUT_BOTTOM, OUT_TOP /) + filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT, BOTTOM, TOP /) end if end select t % find_filter(FILTER_SURFACE) = size(t % filters) diff --git a/src/output.F90 b/src/output.F90 index 9e23f11d8..7aacd6713 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1037,66 +1037,66 @@ contains ! Left Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_LEFT + matching_bins(i_filter_surf) = LEFT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Left", & + "Net Current on Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Right Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_surf) = RIGHT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Right", & + "Net Current on Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Back Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_BACK + matching_bins(i_filter_surf) = BACK filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Back", & + "Net Current on Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Front Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = FRONT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Front", & + "Net Current on Front", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Bottom Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_BOTTOM + matching_bins(i_filter_surf) = BOTTOM filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Bottom", & + "Net Current on Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Top Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = TOP filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Top", & + "Net Current on Top", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) end do diff --git a/src/tally.F90 b/src/tally.F90 index 3c3dc6a69..6bb36d9ce 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2367,7 +2367,7 @@ contains call get_mesh_indices(m, xyz0, ijk0(:m % n_dimension), start_in_mesh) call get_mesh_indices(m, xyz1, ijk1(:m % n_dimension), end_in_mesh) - ! Check to if start or end is in mesh -- if not, check if track still + ! Check to see if start or end is in mesh -- if not, check if track still ! intersects with mesh if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then if (m % n_dimension == 2) then @@ -2407,8 +2407,10 @@ contains if (uvw(3) > 0) then do j = ijk0(3), ijk1(3) - 1 ijk0(3) = j + + ! OUT_TOP if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = TOP matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2417,12 +2419,14 @@ contains t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if - end do - else - do j = ijk0(3), ijk1(3) + 1, -1 - ijk0(3) = j - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BOTTOM + + ! IN_BOTTOM + if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 0 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) < m % dimension(3)) then + ijk0(3) = ijk0(3) + 1 + matching_bins(i_filter_surf) = BOTTOM matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2430,6 +2434,40 @@ contains !$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt + ijk0(3) = ijk0(3) - 1 + end if + end do + else + do j = ijk0(3), ijk1(3) + 1, -1 + ijk0(3) = j + + ! OUT_BOTTOM + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then + matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + end if + + ! IN_TOP + if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) > 1 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) <= m % dimension(3) + 1) then + ijk0(3) = ijk0(3) - 1 + matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + ijk0(3) = ijk0(3) + 1 end if end do end if @@ -2439,8 +2477,10 @@ contains if (uvw(2) > 0) then do j = ijk0(2), ijk1(2) - 1 ijk0(2) = j + + ! OUT_FRONT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = FRONT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2449,12 +2489,14 @@ contains t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if - end do - else - do j = ijk0(2), ijk1(2) + 1, -1 - ijk0(2) = j - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BACK + + ! IN_BACK + if (ijk0(1) >= 1 .and. ijk0(2) >= 0 .and. ijk0(3) >= 1 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) < m % dimension(2) .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(2) = ijk0(2) + 1 + matching_bins(i_filter_surf) = BACK matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2462,6 +2504,40 @@ contains !$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt + ijk0(2) = ijk0(2) - 1 + end if + end do + else + do j = ijk0(2), ijk1(2) + 1, -1 + ijk0(2) = j + + ! OUT_BACK + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then + matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + end if + + ! IN_FRONT + if (ijk0(1) >= 1 .and. ijk0(2) > 1 .and. ijk0(3) >= 1 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) + 1 .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(2) = ijk0(2) - 1 + matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + ijk0(2) = ijk0(2) + 1 end if end do end if @@ -2471,8 +2547,10 @@ contains if (uvw(1) > 0) then do j = ijk0(1), ijk1(1) - 1 ijk0(1) = j + + ! OUT_RIGHT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_surf) = RIGHT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2481,12 +2559,14 @@ contains t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt end if - end do - else - do j = ijk0(1), ijk1(1) + 1, -1 - ijk0(1) = j - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_LEFT + + ! IN_LEFT + if (ijk0(1) >= 0 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & + ijk0(1) < m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(1) = ijk0(1) + 1 + matching_bins(i_filter_surf) = LEFT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2494,6 +2574,40 @@ contains !$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt + ijk0(1) = ijk0(1) - 1 + end if + end do + else + do j = ijk0(1), ijk1(1) + 1, -1 + ijk0(1) = j + + ! OUT_LEFT + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then + matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + end if + + ! IN_RIGHT + if (ijk0(1) > 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & + ijk0(1) <= m % dimension(1) + 1 .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(1) = ijk0(1) - 1 + matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + ijk0(1) = ijk0(1) + 1 end if end do end if @@ -2533,94 +2647,215 @@ contains distance = minval(d) - ! Now use the minimum distance and diretion of the particle to + ! Now use the minimum distance and direction of the particle to ! determine which surface was crossed if (distance == d(1)) then if (uvw(1) > 0) then - ! Crossing into right mesh cell -- this is treated as outgoing - ! current from (i,j,k) + + ! OUT_RIGHT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_surf) = RIGHT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if + + ! IN_LEFT + if (ijk0(1) >= 0 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & + ijk0(1) < m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(1) = ijk0(1) + 1 + matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + ijk0(1) = ijk0(1) - 1 + end if + ijk0(1) = ijk0(1) + 1 xyz_cross(1) = xyz_cross(1) + m % width(1) else - ! Crossing into left mesh cell -- this is treated as outgoing - ! current in (i,j,k) + + ! OUT_LEFT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_LEFT + matching_bins(i_filter_surf) = LEFT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt end if + + ! IN_RIGHT + if (ijk0(1) > 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & + ijk0(1) <= m % dimension(1) + 1 .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(1) = ijk0(1) - 1 + matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + ijk0(1) = ijk0(1) + 1 + end if + ijk0(1) = ijk0(1) - 1 xyz_cross(1) = xyz_cross(1) - m % width(1) end if elseif (distance == d(2)) then if (uvw(2) > 0) then - ! Crossing into front mesh cell -- this is treated as outgoing - ! current in (i,j,k) + + ! OUT_FRONT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = FRONT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if + + ! IN_BACK + if (ijk0(1) >= 1 .and. ijk0(2) >= 0 .and. ijk0(3) >= 1 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) < m % dimension(2) .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(2) = ijk0(2) + 1 + matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + ijk0(2) = ijk0(2) - 1 + end if + ijk0(2) = ijk0(2) + 1 xyz_cross(2) = xyz_cross(2) + m % width(2) else - ! Crossing into back mesh cell -- this is treated as outgoing - ! current in (i,j,k) + + ! OUT_BACK if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BACK + matching_bins(i_filter_surf) = BACK matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt end if + + ! IN_FRONT + if (ijk0(1) >= 1 .and. ijk0(2) > 1 .and. ijk0(3) >= 1 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) + 1 .and. & + ijk0(3) <= m % dimension(3)) then + ijk0(2) = ijk0(2) - 1 + matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + ijk0(2) = ijk0(2) + 1 + end if + ijk0(2) = ijk0(2) - 1 xyz_cross(2) = xyz_cross(2) - m % width(2) end if else if (distance == d(3)) then if (uvw(3) > 0) then - ! Crossing into top mesh cell -- this is treated as outgoing - ! current in (i,j,k) + + ! OUT_TOP if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = TOP matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if + + ! IN_BOTTOM + if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 0 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) < m % dimension(3)) then + ijk0(3) = ijk0(3) + 1 + matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + ijk0(3) = ijk0(3) - 1 + end if + ijk0(3) = ijk0(3) + 1 xyz_cross(3) = xyz_cross(3) + m % width(3) else - ! Crossing into bottom mesh cell -- this is treated as outgoing - ! current in (i,j,k) + + ! OUT_BOTTOM if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BOTTOM + matching_bins(i_filter_surf) = BOTTOM matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt end if + + ! IN_TOP + if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) > 1 .and. & + ijk0(1) <= m % dimension(1) .and. & + ijk0(2) <= m % dimension(2) .and. & + ijk0(3) <= m % dimension(3) + 1) then + ijk0(3) = ijk0(3) - 1 + matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value - p % wgt + ijk0(3) = ijk0(3) + 1 + end if + ijk0(3) = ijk0(3) - 1 xyz_cross(3) = xyz_cross(3) - m % width(3) end if end if - ! Determine scoring index - if (matching_bins(i_filter_surf) > 0) then - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 - - ! Check for errors - if (filter_index <= 0 .or. filter_index > & - t % total_filter_bins) then - call fatal_error("Score index outside range.") - end if - - ! Add to surface current tally -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - ! Calculate new coordinates xyz0 = xyz0 + distance * uvw end do diff --git a/src/trigger.F90 b/src/trigger.F90 index 4af03a2de..6f886cd8e 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -332,7 +332,7 @@ contains mesh_indices_to_bin(m, (/ i, j, k /)) ! Left Surface - matching_bins(i_filter_surf) = OUT_LEFT + matching_bins(i_filter_surf) = LEFT filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -345,7 +345,7 @@ contains trigger % variance = std_dev**2 ! Right Surface - matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_surf) = RIGHT filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -358,7 +358,7 @@ contains trigger % variance = trigger % std_dev**2 ! Back Surface - matching_bins(i_filter_surf) = OUT_BACK + matching_bins(i_filter_surf) = BACK filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -371,7 +371,7 @@ contains trigger % variance = trigger % std_dev**2 ! Front Surface - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = FRONT filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -384,7 +384,7 @@ contains trigger % variance = trigger % std_dev**2 ! Bottom Surface - matching_bins(i_filter_surf) = OUT_BOTTOM + matching_bins(i_filter_surf) = BOTTOM filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -397,7 +397,7 @@ contains trigger % variance = trigger % std_dev**2 ! Top Surface - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = TOP filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) From 93862cea249e4441beb714e719120da3bd2b7d0a Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 23 Aug 2016 11:33:23 -0400 Subject: [PATCH 02/12] changed current tallying so in/out partial currents are tallies separately --- src/cmfd_data.F90 | 122 ++++++++++++++++---------------------------- src/cmfd_header.F90 | 2 +- src/cmfd_input.F90 | 7 ++- src/constants.F90 | 18 ++++--- src/input_xml.F90 | 11 ++-- src/output.F90 | 85 +++++++++++++++++++++++++----- src/tally.F90 | 72 +++++++++++++------------- src/trigger.F90 | 12 ++--- 8 files changed, 185 insertions(+), 144 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 0a1490bb9..befdc9c46 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -51,8 +51,9 @@ contains subroutine compute_xs() use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, & - FILTER_SURFACE, LEFT, RIGHT, BACK, & - FRONT, BOTTOM, TOP, CMFD_NOACCEL, & + FILTER_SURFACE, OUT_LEFT, OUT_RIGHT, OUT_BACK, & + OUT_FRONT, OUT_BOTTOM, OUT_TOP, IN_LEFT, IN_RIGHT, & + IN_BACK, IN_FRONT, IN_BOTTOM, IN_TOP, CMFD_NOACCEL, & ZERO, ONE, TINY_BIT use error, only: fatal_error use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes,& @@ -234,106 +235,74 @@ contains matching_bins(i_filter_ein) = ng - h + 1 end if - ! Left surface + ! Get the bin for this mesh cell matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /)) - matching_bins(i_filter_surf) = LEFT + + ! Left surface + matching_bins(i_filter_surf) = OUT_LEFT score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t%stride) + 1 ! outgoing + * t % stride) + 1 cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum - if (i > 1) then - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i-1, j, k /)) - matching_bins(i_filter_surf) = RIGHT - score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! incoming - cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum - end if + matching_bins(i_filter_surf) = IN_LEFT + score_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum ! Right surface - if (i < nx) then - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i+1, j, k /) ) - matching_bins(i_filter_surf) = LEFT - score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! incoming - cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum - end if - - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j, k /) ) - matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! outgoing + * t % stride) + 1 + cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum + + matching_bins(i_filter_surf) = IN_RIGHT + score_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum ! Back surface - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j, k /)) - matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_surf) = OUT_BACK score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! outgoing + * t % stride) + 1 cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum - if (j > 1) then - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j-1, k /)) - matching_bins(i_filter_surf) = FRONT - score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! incoming - cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum - end if + matching_bins(i_filter_surf) = IN_BACK + score_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum ! Front surface - if (j < ny) then - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j+1, k /)) - matching_bins(i_filter_surf) = BACK - score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! incoming - cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum - end if - - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j, k /)) - matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! outgoing + * t % stride) + 1 + cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum + + matching_bins(i_filter_surf) = IN_FRONT + score_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum ! Bottom surface - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j, k /)) - matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_surf) = OUT_BOTTOM score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! outgoing + * t % stride) + 1 cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum - if (k > 1) then - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j, k-1 /)) - matching_bins(i_filter_surf) = TOP - score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! incoming - cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum - end if + matching_bins(i_filter_surf) = IN_BOTTOM + score_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum ! Top surface - if (k < nz) then - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j, k+1 /)) - matching_bins(i_filter_surf) = BOTTOM - score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! incoming - cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum - end if - - matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i, j, k /)) - matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_surf) = OUT_TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 ! outgoing + * t % stride) + 1 + cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum + + matching_bins(i_filter_surf) = IN_TOP + score_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum end if TALLY @@ -478,7 +447,6 @@ contains ! Get leakage leakage = ZERO LEAK: do l = 1, 3 - leakage = leakage + ((cmfd % current(4*l,g,i,j,k) - & cmfd % current(4*l-1,g,i,j,k))) - & ((cmfd % current(4*l-2,g,i,j,k) - & diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index e743cc928..7d253c308 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -126,7 +126,7 @@ contains if (.not. allocated(this % hxyz)) allocate(this % hxyz(3,nx,ny,nz)) ! Allocate surface currents - if (.not. allocated(this % current)) allocate(this % current(12,ng,nx,ny,nz)) + if (.not. allocated(this % current)) allocate(this % current(6,ng,nx,ny,nz)) ! Allocate source distributions if (.not. allocated(this % cmfd_src)) allocate(this % cmfd_src(ng,nx,ny,nz)) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 77e8f4348..50f078a18 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -534,9 +534,12 @@ contains filt % n_bins = 2 * m % n_dimension allocate(filt % surfaces(2 * m % n_dimension)) if (m % n_dimension == 2) then - filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT /) + filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT, & + IN_LEFT, IN_RIGHT, IN_BACK, IN_FRONT /) elseif (m % n_dimension == 3) then - filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT, BOTTOM, TOP /) + filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT, & + OUT_BOTTOM, OUT_TOP, IN_LEFT, IN_RIGHT, IN_BACK, IN_FRONT, & + IN_BOTTOM, IN_TOP /) end if end select t % find_filter(FILTER_SURFACE) = n_filters diff --git a/src/constants.F90 b/src/constants.F90 index ceb68b9e2..481ca8e03 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -359,12 +359,18 @@ module constants ! Tally surface current directions integer, parameter :: & - LEFT = 1, & ! x min - RIGHT = 2, & ! x max - BACK = 3, & ! y min - FRONT = 4, & ! y max - BOTTOM = 5, & ! z min - TOP = 6 ! z max + OUT_LEFT = 1, & ! x min + OUT_RIGHT = 2, & ! x max + OUT_BACK = 3, & ! y min + OUT_FRONT = 4, & ! y max + OUT_BOTTOM = 5, & ! z min + OUT_TOP = 6, & ! z max + IN_LEFT = 7, & ! x min + IN_RIGHT = 8, & ! x max + IN_BACK = 9, & ! y min + IN_FRONT = 10, & ! y max + IN_BOTTOM = 11, & ! z min + IN_TOP = 12 ! z max ! Tally trigger types and threshold integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8de198386..bcdf310e8 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3698,12 +3698,15 @@ contains allocate(SurfaceFilter :: t % filters(n_filters) % obj) select type (filt => t % filters(size(t % filters)) % obj) type is (SurfaceFilter) - filt % n_bins = 2 * m % n_dimension - allocate(filt % surfaces(2 * m % n_dimension)) + filt % n_bins = 4 * m % n_dimension + allocate(filt % surfaces(4 * m % n_dimension)) if (m % n_dimension == 2) then - filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT /) + filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT, & + IN_LEFT, IN_RIGHT, IN_BACK, IN_FRONT /) elseif (m % n_dimension == 3) then - filt % surfaces = (/ LEFT, RIGHT, BACK, FRONT, BOTTOM, TOP /) + filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT, & + OUT_BOTTOM, OUT_TOP, IN_LEFT, IN_RIGHT, IN_BACK, & + IN_FRONT, IN_BOTTOM, IN_TOP /) end if end select t % find_filter(FILTER_SURFACE) = size(t % filters) diff --git a/src/output.F90 b/src/output.F90 index 7aacd6713..eae654787 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1034,43 +1034,85 @@ contains matching_bins(i_filter_ein))) end if + ! Get the bin for this mesh cell + ! Left Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_surf) = OUT_LEFT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Net Current on Left", & + "Outgoing Current on Left", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) + + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = IN_LEFT + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & + "Incoming Current on Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Right Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Net Current on Right", & + "Outgoing Current on Right", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) + + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = IN_RIGHT + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & + "Incoming Current on Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Back Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_surf) = OUT_BACK filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Net Current on Back", & + "Outgoing Current on Back", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) + + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = IN_BACK + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & + "Incoming Current on Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Front Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_surf) = OUT_FRONT + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & + "Net Current on Front", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) + + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & @@ -1081,26 +1123,45 @@ contains ! Bottom Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_surf) = OUT_BOTTOM filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Net Current on Bottom", & + "Outgoing Current on Bottom", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) + + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = IN_BOTTOM + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & + "Incoming Current on Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) ! Top Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) - matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Net Current on Top", & + "Outgoing Current on Top", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) + + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = IN_TOP + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & + "Incoming Current on Top", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) end do - end do end do end do diff --git a/src/tally.F90 b/src/tally.F90 index 6bb36d9ce..7f743d134 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2410,7 +2410,7 @@ contains ! OUT_TOP if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_surf) = OUT_TOP matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2426,7 +2426,7 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) < m % dimension(3)) then ijk0(3) = ijk0(3) + 1 - matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_surf) = IN_BOTTOM matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2443,14 +2443,14 @@ contains ! OUT_BOTTOM if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_surf) = OUT_BOTTOM matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt end if ! IN_TOP @@ -2459,14 +2459,14 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) <= m % dimension(3) + 1) then ijk0(3) = ijk0(3) - 1 - matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_surf) = IN_TOP matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt ijk0(3) = ijk0(3) + 1 end if end do @@ -2480,7 +2480,7 @@ contains ! OUT_FRONT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_surf) = OUT_FRONT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2496,7 +2496,7 @@ contains ijk0(2) < m % dimension(2) .and. & ijk0(3) <= m % dimension(3)) then ijk0(2) = ijk0(2) + 1 - matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_surf) = IN_BACK matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2513,14 +2513,14 @@ contains ! OUT_BACK if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_surf) = OUT_BACK matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt end if ! IN_FRONT @@ -2529,14 +2529,14 @@ contains ijk0(2) <= m % dimension(2) + 1 .and. & ijk0(3) <= m % dimension(3)) then ijk0(2) = ijk0(2) - 1 - matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_surf) = IN_FRONT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt ijk0(2) = ijk0(2) + 1 end if end do @@ -2550,7 +2550,7 @@ contains ! OUT_RIGHT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_surf) = OUT_RIGHT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2566,7 +2566,7 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) <= m % dimension(3)) then ijk0(1) = ijk0(1) + 1 - matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_surf) = IN_LEFT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2583,14 +2583,14 @@ contains ! OUT_LEFT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_surf) = OUT_LEFT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt end if ! IN_RIGHT @@ -2599,14 +2599,14 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) <= m % dimension(3)) then ijk0(1) = ijk0(1) - 1 - matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_surf) = IN_RIGHT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt ijk0(1) = ijk0(1) + 1 end if end do @@ -2655,7 +2655,7 @@ contains ! OUT_RIGHT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_surf) = OUT_RIGHT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2671,7 +2671,7 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) <= m % dimension(3)) then ijk0(1) = ijk0(1) + 1 - matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_surf) = IN_LEFT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2688,14 +2688,14 @@ contains ! OUT_LEFT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_surf) = OUT_LEFT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt end if ! IN_RIGHT @@ -2704,14 +2704,14 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) <= m % dimension(3)) then ijk0(1) = ijk0(1) - 1 - matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_surf) = IN_RIGHT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt ijk0(1) = ijk0(1) + 1 end if @@ -2723,7 +2723,7 @@ contains ! OUT_FRONT if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_surf) = OUT_FRONT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2739,7 +2739,7 @@ contains ijk0(2) < m % dimension(2) .and. & ijk0(3) <= m % dimension(3)) then ijk0(2) = ijk0(2) + 1 - matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_surf) = IN_BACK matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2756,14 +2756,14 @@ contains ! OUT_BACK if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_surf) = OUT_BACK matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt end if ! IN_FRONT @@ -2772,14 +2772,14 @@ contains ijk0(2) <= m % dimension(2) + 1 .and. & ijk0(3) <= m % dimension(3)) then ijk0(2) = ijk0(2) - 1 - matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_surf) = IN_FRONT matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt ijk0(2) = ijk0(2) + 1 end if @@ -2791,7 +2791,7 @@ contains ! OUT_TOP if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_surf) = OUT_TOP matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2807,7 +2807,7 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) < m % dimension(3)) then ijk0(3) = ijk0(3) + 1 - matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_surf) = IN_BOTTOM matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & @@ -2824,14 +2824,14 @@ contains ! OUT_BOTTOM if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_surf) = OUT_BOTTOM matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt end if ! IN_TOP @@ -2840,14 +2840,14 @@ contains ijk0(2) <= m % dimension(2) .and. & ijk0(3) <= m % dimension(3) + 1) then ijk0(3) = ijk0(3) - 1 - matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_surf) = IN_TOP matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & - t % results(1, filter_index) % value - p % wgt + t % results(1, filter_index) % value + p % wgt ijk0(3) = ijk0(3) + 1 end if diff --git a/src/trigger.F90 b/src/trigger.F90 index 6f886cd8e..4af03a2de 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -332,7 +332,7 @@ contains mesh_indices_to_bin(m, (/ i, j, k /)) ! Left Surface - matching_bins(i_filter_surf) = LEFT + matching_bins(i_filter_surf) = OUT_LEFT filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -345,7 +345,7 @@ contains trigger % variance = std_dev**2 ! Right Surface - matching_bins(i_filter_surf) = RIGHT + matching_bins(i_filter_surf) = OUT_RIGHT filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -358,7 +358,7 @@ contains trigger % variance = trigger % std_dev**2 ! Back Surface - matching_bins(i_filter_surf) = BACK + matching_bins(i_filter_surf) = OUT_BACK filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -371,7 +371,7 @@ contains trigger % variance = trigger % std_dev**2 ! Front Surface - matching_bins(i_filter_surf) = FRONT + matching_bins(i_filter_surf) = OUT_FRONT filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -384,7 +384,7 @@ contains trigger % variance = trigger % std_dev**2 ! Bottom Surface - matching_bins(i_filter_surf) = BOTTOM + matching_bins(i_filter_surf) = OUT_BOTTOM filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -397,7 +397,7 @@ contains trigger % variance = trigger % std_dev**2 ! Top Surface - matching_bins(i_filter_surf) = TOP + matching_bins(i_filter_surf) = OUT_TOP filter_index = & sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) From ba42547eb3a0bb06c83094f5b21f26530eaf3098 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 23 Aug 2016 15:26:24 -0400 Subject: [PATCH 03/12] added tally-arithmetic.ipynb with six-factor formula --- .../pythonapi/examples/tally-arithmetic.ipynb | 533 ++++++++++++------ openmc/filter.py | 18 +- 2 files changed, 360 insertions(+), 191 deletions(-) diff --git a/docs/source/pythonapi/examples/tally-arithmetic.ipynb b/docs/source/pythonapi/examples/tally-arithmetic.ipynb index dfd493f16..088af07b3 100644 --- a/docs/source/pythonapi/examples/tally-arithmetic.ipynb +++ b/docs/source/pythonapi/examples/tally-arithmetic.ipynb @@ -136,8 +136,8 @@ "max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective')\n", "min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective')\n", "max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective')\n", - "min_z = openmc.ZPlane(z0=-0.63, boundary_type='reflective')\n", - "max_z = openmc.ZPlane(z0=+0.63, boundary_type='reflective')" + "min_z = openmc.ZPlane(z0=-100., boundary_type='vacuum')\n", + "max_z = openmc.ZPlane(z0=+100., boundary_type='vacuum')" ] }, { @@ -264,7 +264,7 @@ "settings_file.output = {'tallies': True}\n", "\n", "# Create an initial uniform spatial source distribution over fissionable zones\n", - "bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n", + "bounds = [-0.63, -0.63, -100., 0.63, 0.63, 100.]\n", "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", "settings_file.source = openmc.source.Source(space=uniform_dist)\n", "\n", @@ -339,7 +339,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AHFwInLqDpadAAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDctMjJUMjE6Mzk6\nNDYtMDU6MDBOOEOsAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA3LTIyVDIxOjM5OjQ2LTA1OjAw\nP2X7EAAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AIFw8YNfjajoIAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDgtMjNUMTU6MjQ6\nNTMtMDQ6MDANSReSAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA4LTIzVDE1OjI0OjUzLTA0OjAw\nfBSvLgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] @@ -410,7 +410,32 @@ "tally.filters.append(energy_filter)\n", "tally.scores = ['absorption', 'total']\n", "tally.nuclides = [o16, h1]\n", - "tallies_file.append(tally)" + "tallies_file.append(tally)\n", + "\n", + "# Instantiate a tally mesh \n", + "mesh = openmc.Mesh(mesh_id=1)\n", + "mesh.type = 'regular'\n", + "mesh.dimension = [1, 1, 1]\n", + "mesh.lower_left = [-0.63, -0.63, -100.]\n", + "mesh.width = [1.26, 1.26, 200.]\n", + "mesh_filter = openmc.Filter(type='mesh', bins=[mesh.id])\n", + "mesh_filter.mesh = mesh\n", + "\n", + "# Instantiate thermal, fast, and total leakage tallies\n", + "leak = openmc.Tally(name='leakage')\n", + "leak.filters = [mesh_filter]\n", + "leak.scores = ['current']\n", + "tallies_file.append(leak)\n", + "\n", + "thermal_leak = openmc.Tally(name='thermal leakage')\n", + "thermal_leak.filters = [mesh_filter, openmc.Filter(type='energy', bins=[0., 0.625e-6])]\n", + "thermal_leak.scores = ['current']\n", + "tallies_file.append(thermal_leak)\n", + "\n", + "fast_leak = openmc.Tally(name='fast leakage')\n", + "fast_leak.filters = [mesh_filter, openmc.Filter(type='energy', bins=[0.625e-6, 20.])]\n", + "fast_leak.scores = ['current']\n", + "tallies_file.append(fast_leak)" ] }, { @@ -484,12 +509,12 @@ "outputs": [], "source": [ "# Instantiate energy filter to illustrate Tally slicing\n", - "energy_filter = openmc.Filter(type='energy', bins=np.logspace(np.log10(1e-8), np.log10(20), 10))\n", + "fine_energy_filter = openmc.Filter(type='energy', bins=np.logspace(np.log10(1e-8), np.log10(20), 10))\n", "\n", "# Instantiate flux Tally in moderator and fuel\n", "tally = openmc.Tally(name='need-to-slice')\n", "tally.filters = [openmc.Filter(type='cell', bins=[fuel_cell.id, moderator_cell.id])]\n", - "tally.filters.append(energy_filter)\n", + "tally.filters.append(fine_energy_filter)\n", "tally.scores = ['nu-fission', 'scatter']\n", "tally.nuclides = [h1, u238]\n", "tallies_file.append(tally)" @@ -541,9 +566,10 @@ "\n", " Copyright: 2011-2016 Massachusetts Institute of Technology\n", " License: http://openmc.readthedocs.io/en/latest/license.html\n", - " Version: 0.7.1\n", - " Git SHA1: 3d68c07625e33cd64188df03ee03e9c31b3d4b74\n", - " Date/Time: 2016-07-22 21:39:46\n", + " Version: 0.8.0\n", + " Git SHA1: 93862cea249e4441beb714e719120da3bd2b7d0a\n", + " Date/Time: 2016-08-23 15:24:53\n", + " MPI Processes: 1\n", "\n", " ===========================================================================\n", " ========================> INITIALIZATION <=========================\n", @@ -553,12 +579,12 @@ " Reading geometry XML file...\n", " Reading cross sections XML file...\n", " Reading materials XML file...\n", - " Reading U235.71c from /home/romano/openmc/data/nndc_hdf5/U235_71c.h5\n", - " Reading U238.71c from /home/romano/openmc/data/nndc_hdf5/U238_71c.h5\n", - " Reading O16.71c from /home/romano/openmc/data/nndc_hdf5/O16_71c.h5\n", - " Reading H1.71c from /home/romano/openmc/data/nndc_hdf5/H1_71c.h5\n", - " Reading B10.71c from /home/romano/openmc/data/nndc_hdf5/B10_71c.h5\n", - " Reading Zr90.71c from /home/romano/openmc/data/nndc_hdf5/Zr90_71c.h5\n", + " Reading U235.71c from /Users/sam/git/openmc-sam/data/nndc_hdf5/U235_71c.h5\n", + " Reading U238.71c from /Users/sam/git/openmc-sam/data/nndc_hdf5/U238_71c.h5\n", + " Reading O16.71c from /Users/sam/git/openmc-sam/data/nndc_hdf5/O16_71c.h5\n", + " Reading H1.71c from /Users/sam/git/openmc-sam/data/nndc_hdf5/H1_71c.h5\n", + " Reading B10.71c from /Users/sam/git/openmc-sam/data/nndc_hdf5/B10_71c.h5\n", + " Reading Zr90.71c from /Users/sam/git/openmc-sam/data/nndc_hdf5/Zr90_71c.h5\n", " Maximum neutron transport energy: 20.0000 MeV for U235.71c\n", " Reading tallies XML file...\n", " Building neighboring cells lists for each surface...\n", @@ -570,26 +596,26 @@ "\n", " Bat./Gen. k Average k \n", " ========= ======== ==================== \n", - " 1/1 1.03471 \n", - " 2/1 1.03257 \n", - " 3/1 1.00600 \n", - " 4/1 1.04547 \n", - " 5/1 1.02287 \n", - " 6/1 1.05752 \n", - " 7/1 1.04283 1.05017 +/- 0.00734\n", - " 8/1 1.05189 1.05074 +/- 0.00428\n", - " 9/1 1.01645 1.04217 +/- 0.00909\n", - " 10/1 1.04978 1.04369 +/- 0.00721\n", - " 11/1 1.03459 1.04218 +/- 0.00608\n", - " 12/1 1.04019 1.04189 +/- 0.00514\n", - " 13/1 1.05985 1.04414 +/- 0.00499\n", - " 14/1 1.02111 1.04158 +/- 0.00509\n", - " 15/1 1.04774 1.04219 +/- 0.00459\n", - " 16/1 1.00733 1.03902 +/- 0.00523\n", - " 17/1 1.02224 1.03763 +/- 0.00497\n", - " 18/1 1.03263 1.03724 +/- 0.00459\n", - " 19/1 1.01611 1.03573 +/- 0.00451\n", - " 20/1 1.04692 1.03648 +/- 0.00426\n", + " 1/1 0.96168 \n", + " 2/1 0.96651 \n", + " 3/1 1.00678 \n", + " 4/1 0.98773 \n", + " 5/1 1.01883 \n", + " 6/1 1.02959 \n", + " 7/1 0.99859 1.01409 +/- 0.01550\n", + " 8/1 1.03441 1.02086 +/- 0.01123\n", + " 9/1 1.06097 1.03089 +/- 0.01279\n", + " 10/1 1.06094 1.03690 +/- 0.01159\n", + " 11/1 1.04687 1.03856 +/- 0.00961\n", + " 12/1 1.02982 1.03731 +/- 0.00821\n", + " 13/1 1.03520 1.03705 +/- 0.00712\n", + " 14/1 0.99508 1.03239 +/- 0.00782\n", + " 15/1 1.03973 1.03312 +/- 0.00703\n", + " 16/1 1.03807 1.03357 +/- 0.00638\n", + " 17/1 1.03091 1.03335 +/- 0.00583\n", + " 18/1 1.01421 1.03188 +/- 0.00556\n", + " 19/1 0.99339 1.02913 +/- 0.00583\n", + " 20/1 1.04827 1.03040 +/- 0.00558\n", " Creating state point statepoint.20.h5...\n", "\n", " ===========================================================================\n", @@ -599,28 +625,28 @@ "\n", " =======================> TIMING STATISTICS <=======================\n", "\n", - " Total time for initialization = 3.5600E-01 seconds\n", - " Reading cross sections = 2.3400E-01 seconds\n", - " Total time in simulation = 1.8333E+01 seconds\n", - " Time in transport only = 1.8325E+01 seconds\n", - " Time in inactive batches = 2.6950E+00 seconds\n", - " Time in active batches = 1.5638E+01 seconds\n", - " Time synchronizing fission bank = 1.0000E-03 seconds\n", - " Sampling source sites = 0.0000E+00 seconds\n", + " Total time for initialization = 5.4300E-01 seconds\n", + " Reading cross sections = 3.5400E-01 seconds\n", + " Total time in simulation = 1.8555E+01 seconds\n", + " Time in transport only = 1.8464E+01 seconds\n", + " Time in inactive batches = 2.4600E+00 seconds\n", + " Time in active batches = 1.6095E+01 seconds\n", + " Time synchronizing fission bank = 1.1000E-02 seconds\n", + " Sampling source sites = 4.0000E-03 seconds\n", " SEND/RECV source sites = 1.0000E-03 seconds\n", " Time accumulating tallies = 0.0000E+00 seconds\n", - " Total time for finalization = 1.0000E-03 seconds\n", - " Total time elapsed = 1.8711E+01 seconds\n", - " Calculation Rate (inactive) = 4638.22 neutrons/second\n", - " Calculation Rate (active) = 2398.00 neutrons/second\n", + " Total time for finalization = 3.0000E-03 seconds\n", + " Total time elapsed = 1.9120E+01 seconds\n", + " Calculation Rate (inactive) = 5081.30 neutrons/second\n", + " Calculation Rate (active) = 2329.92 neutrons/second\n", "\n", " ============================> RESULTS <============================\n", "\n", - " k-effective (Collision) = 1.03296 +/- 0.00669\n", - " k-effective (Track-length) = 1.03648 +/- 0.00426\n", - " k-effective (Absorption) = 1.03431 +/- 0.00702\n", - " Combined k-effective = 1.03621 +/- 0.00456\n", - " Leakage Fraction = 0.00000 +/- 0.00000\n", + " k-effective (Collision) = 1.02791 +/- 0.00553\n", + " k-effective (Track-length) = 1.03040 +/- 0.00558\n", + " k-effective (Absorption) = 1.02011 +/- 0.00491\n", + " Combined k-effective = 1.02461 +/- 0.00398\n", + " Leakage Fraction = 0.01677 +/- 0.00109\n", "\n" ] }, @@ -674,8 +700,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We have a tally of the total fission rate and the total absorption rate, so we can calculate k-infinity as:\n", - "$$k_\\infty = \\frac{\\langle \\nu \\Sigma_f \\phi \\rangle}{\\langle \\Sigma_a \\phi \\rangle}$$\n", + "We have a tally of the total fission rate and the total absorption rate, so we can calculate k-eff as:\n", + "$$k_{eff} = \\frac{\\langle \\nu \\Sigma_f \\phi \\rangle}{\\langle \\Sigma_a \\phi \\rangle + \\langle L \\rangle}$$\n", "In this notation, $\\langle \\cdot \\rangle^a_b$ represents an OpenMC that is integrated over region $a$ and energy range $b$. If $a$ or $b$ is not reported, it means the value represents an integral over all space or all energy, respectively." ] }, @@ -704,17 +730,17 @@ " \n", " 0\n", " total\n", - " (nu-fission / absorption)\n", - " 1.038387\n", - " 0.006141\n", + " (nu-fission / (absorption + current))\n", + " 1.02431\n", + " 0.00704\n", " \n", " \n", "\n", "" ], "text/plain": [ - " nuclide score mean std. dev.\n", - "0 total (nu-fission / absorption) 1.04e+00 6.14e-03" + " nuclide score mean std. dev.\n", + "0 total (nu-fission / (absorption + current)) 1.02e+00 7.04e-03" ] }, "execution_count": 24, @@ -723,10 +749,17 @@ } ], "source": [ - "# Compute k-infinity using tally arithmetic\n", + "# Get the fission and absorption rate tallies\n", "fiss_rate = sp.get_tally(name='fiss. rate')\n", "abs_rate = sp.get_tally(name='abs. rate')\n", - "keff = fiss_rate / abs_rate\n", + "\n", + "# Get the leakage tally\n", + "leak = sp.get_tally(name='leakage')\n", + "leak = leak.summation(filter_type='surface', remove_filter=True)\n", + "leak = leak.summation(filter_type='mesh', remove_filter=True)\n", + "\n", + "# Compute k-infinity using tally arithmetic\n", + "keff = fiss_rate / (abs_rate + leak)\n", "keff.get_pandas_dataframe()" ] }, @@ -734,9 +767,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Notice that even though the neutron production rate and absorption rate are separate tallies, we still get a first-order estimate of the uncertainty on the quotient of them automatically!\n", + "Notice that even though the neutron production rate, absorption rate, and current are separate tallies, we still get a first-order estimate of the uncertainty on the quotient of them automatically!\n", "\n", - "Often in textbooks you'll see k-infinity represented using the four-factor formula $$k_\\infty = p \\epsilon f \\eta.$$ Let's analyze each of these factors, starting with the resonance escape probability which is defined as $$p=\\frac{\\langle\\Sigma_a\\phi\\rangle_T}{\\langle\\Sigma_a\\phi\\rangle}$$ where the subscript $T$ means thermal energies." + "Often in textbooks you'll see k-eff represented using the six-factor formula $$k_{eff} = p \\epsilon f \\eta P_{FNL} P_{TNL}.$$ Let's analyze each of these factors, starting with the resonance escape probability which is defined as $$p=\\frac{\\langle\\Sigma_a\\phi\\rangle_T + \\langle L \\rangle_T}{\\langle\\Sigma_a\\phi\\rangle + \\langle L \\rangle_T}$$ where the subscript $T$ means thermal energies." ] }, { @@ -768,17 +801,20 @@ " 0.0\n", " 6.250000e-07\n", " total\n", - " absorption\n", - " 0.693337\n", - " 0.004109\n", + " (absorption + current)\n", + " 0.695303\n", + " 0.005091\n", " \n", " \n", "\n", "" ], "text/plain": [ - " energy low [MeV] energy high [MeV] nuclide score mean std. dev.\n", - "0 0.00e+00 6.25e-07 total absorption 6.93e-01 4.11e-03" + " energy low [MeV] energy high [MeV] nuclide score \\\n", + "0 0.00e+00 6.25e-07 total (absorption + current) \n", + "\n", + " mean std. dev. \n", + "0 6.95e-01 5.09e-03 " ] }, "execution_count": 25, @@ -789,7 +825,10 @@ "source": [ "# Compute resonance escape probability using tally arithmetic\n", "therm_abs_rate = sp.get_tally(name='therm. abs. rate')\n", - "res_esc = therm_abs_rate / abs_rate\n", + "thermal_leak = sp.get_tally(name='thermal leakage')\n", + "thermal_leak = thermal_leak.summation(filter_type='surface', remove_filter=True)\n", + "thermal_leak = thermal_leak.summation(filter_type='mesh', remove_filter=True)\n", + "res_esc = (therm_abs_rate + thermal_leak) / (abs_rate + thermal_leak)\n", "res_esc.get_pandas_dataframe()" ] }, @@ -831,8 +870,8 @@ " 6.250000e-07\n", " total\n", " nu-fission\n", - " 1.203042\n", - " 0.0076\n", + " 1.202639\n", + " 0.010348\n", " \n", " \n", "\n", @@ -840,7 +879,7 @@ ], "text/plain": [ " energy low [MeV] energy high [MeV] nuclide score mean std. dev.\n", - "0 0.00e+00 6.25e-07 total nu-fission 1.20e+00 7.60e-03" + "0 0.00e+00 6.25e-07 total nu-fission 1.20e+00 1.03e-02" ] }, "execution_count": 26, @@ -896,8 +935,8 @@ " 10000\n", " total\n", " absorption\n", - " 0.748413\n", - " 0.004723\n", + " 0.749349\n", + " 0.006731\n", " \n", " \n", "\n", @@ -905,10 +944,10 @@ ], "text/plain": [ " energy low [MeV] energy high [MeV] cell nuclide score mean \\\n", - "0 0.00e+00 6.25e-07 10000 total absorption 7.48e-01 \n", + "0 0.00e+00 6.25e-07 10000 total absorption 7.49e-01 \n", "\n", " std. dev. \n", - "0 4.72e-03 " + "0 6.73e-03 " ] }, "execution_count": 27, @@ -927,7 +966,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The final factor is the number of fission neutrons produced per absorption in fuel, calculated as $$\\eta = \\frac{\\langle \\nu\\Sigma_f\\phi \\rangle_T}{\\langle \\Sigma_a \\phi \\rangle^F_T}$$" + "The next factor is the number of fission neutrons produced per absorption in fuel, calculated as $$\\eta = \\frac{\\langle \\nu\\Sigma_f\\phi \\rangle_T}{\\langle \\Sigma_a \\phi \\rangle^F_T}$$" ] }, { @@ -962,8 +1001,8 @@ " 10000\n", " total\n", " (nu-fission / absorption)\n", - " 1.663385\n", - " 0.011253\n", + " 1.663736\n", + " 0.015707\n", " \n", " \n", "\n", @@ -974,7 +1013,7 @@ "0 0.00e+00 6.25e-07 10000 total \n", "\n", " score mean std. dev. \n", - "0 (nu-fission / absorption) 1.66e+00 1.13e-02 " + "0 (nu-fission / absorption) 1.66e+00 1.57e-02 " ] }, "execution_count": 28, @@ -992,7 +1031,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now we can calculate $k_\\infty$ using the product of the factors form the four-factor formula." + "There are two leakage factors to account for fast and thermal leakage. The fast non-leakage probability is computed as $$P_{FNL} = \\frac{\\langle \\Sigma_a\\phi \\rangle + \\langle L \\rangle_T}{\\langle \\Sigma_a \\phi \\rangle + \\langle L \\rangle}$$" ] }, { @@ -1001,6 +1040,130 @@ "metadata": { "collapsed": false }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
energy low [MeV]energy high [MeV]nuclidescoremeanstd. dev.
00.06.250000e-07total(absorption + current)0.9851020.005855
\n", + "
" + ], + "text/plain": [ + " energy low [MeV] energy high [MeV] nuclide score \\\n", + "0 0.00e+00 6.25e-07 total (absorption + current) \n", + "\n", + " mean std. dev. \n", + "0 9.85e-01 5.86e-03 " + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p_fnl = (abs_rate + thermal_leak) / (abs_rate + leak)\n", + "p_fnl.get_pandas_dataframe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The final factor is the thermalnon-leakage probability and is computed as $$P_{TNL} = \\frac{\\langle \\Sigma_a\\phi \\rangle_T}{\\langle \\Sigma_a \\phi \\rangle_T + \\langle L \\rangle_T}$$" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
energy low [MeV]energy high [MeV]nuclidescoremeanstd. dev.
00.06.250000e-07total(absorption / (absorption + current))0.9974070.008492
\n", + "
" + ], + "text/plain": [ + " energy low [MeV] energy high [MeV] nuclide \\\n", + "0 0.00e+00 6.25e-07 total \n", + "\n", + " score mean std. dev. \n", + "0 (absorption / (absorption + current)) 9.97e-01 8.49e-03 " + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p_tnl = therm_abs_rate / (therm_abs_rate + thermal_leak)\n", + "p_tnl.get_pandas_dataframe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can calculate $k_{eff}$ using the product of the factors form the four-factor formula." + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -1026,9 +1189,9 @@ " 6.250000e-07\n", " 10000\n", " total\n", - " (((absorption * nu-fission) * absorption) * (n...\n", - " 1.038387\n", - " 0.01316\n", + " ((((((absorption + current) * nu-fission) * ab...\n", + " 1.02431\n", + " 0.02062\n", " \n", " \n", "\n", @@ -1039,16 +1202,16 @@ "0 0.00e+00 6.25e-07 10000 total \n", "\n", " score mean std. dev. \n", - "0 (((absorption * nu-fission) * absorption) * (n... 1.04e+00 1.32e-02 " + "0 ((((((absorption + current) * nu-fission) * ab... 1.02e+00 2.06e-02 " ] }, - "execution_count": 29, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "keff = res_esc * fast_fiss * therm_util * eta\n", + "keff = res_esc * fast_fiss * therm_util * eta * p_fnl * p_tnl\n", "keff.get_pandas_dataframe()" ] }, @@ -1063,7 +1226,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 32, "metadata": { "collapsed": false, "scrolled": true @@ -1079,7 +1242,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 33, "metadata": { "collapsed": false }, @@ -1109,8 +1272,8 @@ " 6.250000e-07\n", " (U238 / total)\n", " (nu-fission / flux)\n", - " 6.636968e-07\n", - " 4.132875e-09\n", + " 6.662479e-07\n", + " 6.039323e-09\n", " \n", " \n", " 1\n", @@ -1119,8 +1282,8 @@ " 6.250000e-07\n", " (U238 / total)\n", " (scatter / flux)\n", - " 2.099856e-01\n", - " 1.232455e-03\n", + " 2.099897e-01\n", + " 1.843251e-03\n", " \n", " \n", " 2\n", @@ -1129,8 +1292,8 @@ " 6.250000e-07\n", " (U235 / total)\n", " (nu-fission / flux)\n", - " 3.552458e-01\n", - " 2.252681e-03\n", + " 3.568130e-01\n", + " 3.255144e-03\n", " \n", " \n", " 3\n", @@ -1139,8 +1302,8 @@ " 6.250000e-07\n", " (U235 / total)\n", " (scatter / flux)\n", - " 5.554345e-03\n", - " 3.265385e-05\n", + " 5.555326e-03\n", + " 4.893022e-05\n", " \n", " \n", " 4\n", @@ -1149,8 +1312,8 @@ " 2.000000e+01\n", " (U238 / total)\n", " (nu-fission / flux)\n", - " 7.126668e-03\n", - " 5.296883e-05\n", + " 7.215044e-03\n", + " 4.968448e-05\n", " \n", " \n", " 5\n", @@ -1159,8 +1322,8 @@ " 2.000000e+01\n", " (U238 / total)\n", " (scatter / flux)\n", - " 2.277460e-01\n", - " 1.003558e-03\n", + " 2.273966e-01\n", + " 8.969811e-04\n", " \n", " \n", " 6\n", @@ -1169,8 +1332,8 @@ " 2.000000e+01\n", " (U235 / total)\n", " (nu-fission / flux)\n", - " 8.010911e-03\n", - " 6.802256e-05\n", + " 7.969615e-03\n", + " 5.374119e-05\n", " \n", " \n", " 7\n", @@ -1179,8 +1342,8 @@ " 2.000000e+01\n", " (U235 / total)\n", " (scatter / flux)\n", - " 3.367794e-03\n", - " 1.443644e-05\n", + " 3.362798e-03\n", + " 1.286767e-05\n", " \n", " \n", "\n", @@ -1198,17 +1361,17 @@ "7 10000 6.25e-07 2.00e+01 (U235 / total) \n", "\n", " score mean std. dev. \n", - "0 (nu-fission / flux) 6.64e-07 4.13e-09 \n", - "1 (scatter / flux) 2.10e-01 1.23e-03 \n", - "2 (nu-fission / flux) 3.55e-01 2.25e-03 \n", - "3 (scatter / flux) 5.55e-03 3.27e-05 \n", - "4 (nu-fission / flux) 7.13e-03 5.30e-05 \n", - "5 (scatter / flux) 2.28e-01 1.00e-03 \n", - "6 (nu-fission / flux) 8.01e-03 6.80e-05 \n", - "7 (scatter / flux) 3.37e-03 1.44e-05 " + "0 (nu-fission / flux) 6.66e-07 6.04e-09 \n", + "1 (scatter / flux) 2.10e-01 1.84e-03 \n", + "2 (nu-fission / flux) 3.57e-01 3.26e-03 \n", + "3 (scatter / flux) 5.56e-03 4.89e-05 \n", + "4 (nu-fission / flux) 7.22e-03 4.97e-05 \n", + "5 (scatter / flux) 2.27e-01 8.97e-04 \n", + "6 (nu-fission / flux) 7.97e-03 5.37e-05 \n", + "7 (scatter / flux) 3.36e-03 1.29e-05 " ] }, - "execution_count": 31, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" } @@ -1227,7 +1390,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 34, "metadata": { "collapsed": false }, @@ -1236,11 +1399,11 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[[ 6.63696783e-07]\n", - " [ 3.55245846e-01]]\n", + "[[[ 6.66247898e-07]\n", + " [ 3.56812954e-01]]\n", "\n", - " [[ 7.12666800e-03]\n", - " [ 8.01091088e-03]]]\n" + " [[ 7.21504433e-03]\n", + " [ 7.96961502e-03]]]\n" ] } ], @@ -1259,7 +1422,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 35, "metadata": { "collapsed": false }, @@ -1268,9 +1431,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[[ 0.00555435]]\n", + "[[[ 0.00555533]]\n", "\n", - " [[ 0.00336779]]]\n" + " [[ 0.0033628 ]]]\n" ] } ], @@ -1283,7 +1446,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 36, "metadata": { "collapsed": false }, @@ -1292,8 +1455,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[[ 0.22774598]\n", - " [ 0.00336779]]]\n" + "[[[ 0.22739657]\n", + " [ 0.0033628 ]]]\n" ] } ], @@ -1314,7 +1477,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 37, "metadata": { "collapsed": false }, @@ -1345,7 +1508,7 @@ " U238\n", " nu-fission\n", " 0.000002\n", - " 7.473789e-09\n", + " 1.057199e-08\n", " \n", " \n", " 1\n", @@ -1354,8 +1517,8 @@ " 6.250000e-07\n", " U235\n", " nu-fission\n", - " 0.861547\n", - " 4.131310e-03\n", + " 0.856784\n", + " 5.730044e-03\n", " \n", " \n", " 2\n", @@ -1364,8 +1527,8 @@ " 2.000000e+01\n", " U238\n", " nu-fission\n", - " 0.082356\n", - " 5.560461e-04\n", + " 0.082495\n", + " 5.176027e-04\n", " \n", " \n", " 3\n", @@ -1374,8 +1537,8 @@ " 2.000000e+01\n", " U235\n", " nu-fission\n", - " 0.092574\n", - " 7.315442e-04\n", + " 0.091123\n", + " 5.574052e-04\n", " \n", " \n", "\n", @@ -1383,19 +1546,19 @@ ], "text/plain": [ " cell energy low [MeV] energy high [MeV] nuclide score mean \\\n", - "0 10000 0.00e+00 6.25e-07 U238 nu-fission 1.61e-06 \n", - "1 10000 0.00e+00 6.25e-07 U235 nu-fission 8.62e-01 \n", - "2 10000 6.25e-07 2.00e+01 U238 nu-fission 8.24e-02 \n", - "3 10000 6.25e-07 2.00e+01 U235 nu-fission 9.26e-02 \n", + "0 10000 0.00e+00 6.25e-07 U238 nu-fission 1.60e-06 \n", + "1 10000 0.00e+00 6.25e-07 U235 nu-fission 8.57e-01 \n", + "2 10000 6.25e-07 2.00e+01 U238 nu-fission 8.25e-02 \n", + "3 10000 6.25e-07 2.00e+01 U235 nu-fission 9.11e-02 \n", "\n", " std. dev. \n", - "0 7.47e-09 \n", - "1 4.13e-03 \n", - "2 5.56e-04 \n", - "3 7.32e-04 " + "0 1.06e-08 \n", + "1 5.73e-03 \n", + "2 5.18e-04 \n", + "3 5.57e-04 " ] }, - "execution_count": 35, + "execution_count": 37, "metadata": {}, "output_type": "execute_result" } @@ -1408,7 +1571,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 38, "metadata": { "collapsed": false }, @@ -1438,8 +1601,8 @@ " 1.080060e-07\n", " H1\n", " scatter\n", - " 4.599225\n", - " 0.015973\n", + " 4.547947\n", + " 0.028000\n", " \n", " \n", " 1\n", @@ -1448,8 +1611,8 @@ " 1.166529e-06\n", " H1\n", " scatter\n", - " 2.037260\n", - " 0.011236\n", + " 2.003068\n", + " 0.008587\n", " \n", " \n", " 2\n", @@ -1458,8 +1621,8 @@ " 1.259921e-05\n", " H1\n", " scatter\n", - " 1.662552\n", - " 0.010280\n", + " 1.647225\n", + " 0.011136\n", " \n", " \n", " 3\n", @@ -1468,8 +1631,8 @@ " 1.360790e-04\n", " H1\n", " scatter\n", - " 1.872201\n", - " 0.012136\n", + " 1.831367\n", + " 0.010196\n", " \n", " \n", " 4\n", @@ -1478,8 +1641,8 @@ " 1.469734e-03\n", " H1\n", " scatter\n", - " 2.080459\n", - " 0.013155\n", + " 2.039613\n", + " 0.008059\n", " \n", " \n", " 5\n", @@ -1488,8 +1651,8 @@ " 1.587401e-02\n", " H1\n", " scatter\n", - " 2.154996\n", - " 0.011975\n", + " 2.137523\n", + " 0.012885\n", " \n", " \n", " 6\n", @@ -1498,8 +1661,8 @@ " 1.714488e-01\n", " H1\n", " scatter\n", - " 2.218740\n", - " 0.008528\n", + " 2.170725\n", + " 0.012669\n", " \n", " \n", " 7\n", @@ -1508,8 +1671,8 @@ " 1.851749e+00\n", " H1\n", " scatter\n", - " 2.010517\n", - " 0.009187\n", + " 2.002724\n", + " 0.010768\n", " \n", " \n", " 8\n", @@ -1518,8 +1681,8 @@ " 2.000000e+01\n", " H1\n", " scatter\n", - " 0.372022\n", - " 0.003196\n", + " 0.371624\n", + " 0.002959\n", " \n", " \n", "\n", @@ -1527,29 +1690,29 @@ ], "text/plain": [ " cell energy low [MeV] energy high [MeV] nuclide score mean \\\n", - "0 10002 1.00e-08 1.08e-07 H1 scatter 4.60e+00 \n", - "1 10002 1.08e-07 1.17e-06 H1 scatter 2.04e+00 \n", - "2 10002 1.17e-06 1.26e-05 H1 scatter 1.66e+00 \n", - "3 10002 1.26e-05 1.36e-04 H1 scatter 1.87e+00 \n", - "4 10002 1.36e-04 1.47e-03 H1 scatter 2.08e+00 \n", - "5 10002 1.47e-03 1.59e-02 H1 scatter 2.15e+00 \n", - "6 10002 1.59e-02 1.71e-01 H1 scatter 2.22e+00 \n", - "7 10002 1.71e-01 1.85e+00 H1 scatter 2.01e+00 \n", + "0 10002 1.00e-08 1.08e-07 H1 scatter 4.55e+00 \n", + "1 10002 1.08e-07 1.17e-06 H1 scatter 2.00e+00 \n", + "2 10002 1.17e-06 1.26e-05 H1 scatter 1.65e+00 \n", + "3 10002 1.26e-05 1.36e-04 H1 scatter 1.83e+00 \n", + "4 10002 1.36e-04 1.47e-03 H1 scatter 2.04e+00 \n", + "5 10002 1.47e-03 1.59e-02 H1 scatter 2.14e+00 \n", + "6 10002 1.59e-02 1.71e-01 H1 scatter 2.17e+00 \n", + "7 10002 1.71e-01 1.85e+00 H1 scatter 2.00e+00 \n", "8 10002 1.85e+00 2.00e+01 H1 scatter 3.72e-01 \n", "\n", " std. dev. \n", - "0 1.60e-02 \n", - "1 1.12e-02 \n", - "2 1.03e-02 \n", - "3 1.21e-02 \n", - "4 1.32e-02 \n", - "5 1.20e-02 \n", - "6 8.53e-03 \n", - "7 9.19e-03 \n", - "8 3.20e-03 " + "0 2.80e-02 \n", + "1 8.59e-03 \n", + "2 1.11e-02 \n", + "3 1.02e-02 \n", + "4 8.06e-03 \n", + "5 1.29e-02 \n", + "6 1.27e-02 \n", + "7 1.08e-02 \n", + "8 2.96e-03 " ] }, - "execution_count": 36, + "execution_count": 38, "metadata": {}, "output_type": "execute_result" } @@ -1565,21 +1728,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 2", "language": "python", - "name": "python3" + "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 3 + "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" + "pygments_lexer": "ipython2", + "version": "2.7.12" } }, "nbformat": 4, diff --git a/openmc/filter.py b/openmc/filter.py index 9b4eb67dd..771d40401 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -781,12 +781,18 @@ class Filter(object): filter_bins = np.repeat(self.bins, self.stride) tile_factor = data_size / len(filter_bins) filter_bins = np.tile(filter_bins, tile_factor) - filter_bins = [x if x != 1 else 'x-min' for x in filter_bins] - filter_bins = [x if x != 2 else 'x-max' for x in filter_bins] - filter_bins = [x if x != 3 else 'y-min' for x in filter_bins] - filter_bins = [x if x != 4 else 'y-max' for x in filter_bins] - filter_bins = [x if x != 5 else 'z-min' for x in filter_bins] - filter_bins = [x if x != 6 else 'z-max' for x in filter_bins] + filter_bins = [x if x != 1 else 'x-min out' for x in filter_bins] + filter_bins = [x if x != 2 else 'x-max out' for x in filter_bins] + filter_bins = [x if x != 3 else 'y-min out' for x in filter_bins] + filter_bins = [x if x != 4 else 'y-max out' for x in filter_bins] + filter_bins = [x if x != 5 else 'z-min out' for x in filter_bins] + filter_bins = [x if x != 6 else 'z-max out' for x in filter_bins] + filter_bins = [x if x != 7 else 'x-min in' for x in filter_bins] + filter_bins = [x if x != 8 else 'x-max in' for x in filter_bins] + filter_bins = [x if x != 9 else 'y-min in' for x in filter_bins] + filter_bins = [x if x != 10 else 'y-max in' for x in filter_bins] + filter_bins = [x if x != 11 else 'z-min in' for x in filter_bins] + filter_bins = [x if x != 12 else 'z-max in' for x in filter_bins] df = pd.concat([df, pd.DataFrame({self.type : filter_bins})]) # universe, material, surface, cell, and cellborn filters From ac17ca4e612babd0e6509b2d94a3d4319dce9896 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 23 Aug 2016 15:27:34 -0400 Subject: [PATCH 04/12] reverted pincell xml example --- examples/xml/pincell/geometry.xml | 14 +++++++------- examples/xml/pincell/settings.xml | 12 ++++++------ examples/xml/pincell/tallies.xml | 11 +++-------- 3 files changed, 16 insertions(+), 21 deletions(-) diff --git a/examples/xml/pincell/geometry.xml b/examples/xml/pincell/geometry.xml index f4e1bf1dd..f67f9e74c 100644 --- a/examples/xml/pincell/geometry.xml +++ b/examples/xml/pincell/geometry.xml @@ -8,16 +8,16 @@ --> - - - + + + - - - - + + + + diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml index 9402b979a..443af9cde 100644 --- a/examples/xml/pincell/settings.xml +++ b/examples/xml/pincell/settings.xml @@ -3,7 +3,7 @@ - 20 + 100 10 1000 @@ -14,8 +14,8 @@ - -8.62992 -8.62992 -1. - 8.62992 8.62992 1. + -0.62992 -0.62992 -1. + 0.62992 0.62992 1. @@ -24,9 +24,9 @@ bounds for a mesh over which the Shannon entropy should be calculated. The extent in the z direction is made arbitrarily large. --> - -3.39218 -3.39218 -1.e50 - 3.39218 3.39218 1.e50 + -0.39218 -0.39218 -1.e50 + 0.39218 0.39218 1.e50 10 10 1 - + \ No newline at end of file diff --git a/examples/xml/pincell/tallies.xml b/examples/xml/pincell/tallies.xml index 2e39c083d..73242b913 100644 --- a/examples/xml/pincell/tallies.xml +++ b/examples/xml/pincell/tallies.xml @@ -2,9 +2,9 @@ - 2 2 1 - -8.62992 -8.62992 -1.e50 - 8.62992 8.62992 1.e50 + 100 100 1 + -0.62992 -0.62992 -1.e50 + 0.62992 0.62992 1.e50 @@ -13,9 +13,4 @@ flux fission nu-fission - - - current - - From 28f50a635a6365bbf365858e07046680b82798f4 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 23 Aug 2016 15:31:44 -0400 Subject: [PATCH 05/12] removed unnecessary repeated filter mesh bin look up in output.F90 --- src/cmfd_header.F90 | 2 +- src/output.F90 | 26 ++------------------------ 2 files changed, 3 insertions(+), 25 deletions(-) diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 7d253c308..e743cc928 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -126,7 +126,7 @@ contains if (.not. allocated(this % hxyz)) allocate(this % hxyz(3,nx,ny,nz)) ! Allocate surface currents - if (.not. allocated(this % current)) allocate(this % current(6,ng,nx,ny,nz)) + if (.not. allocated(this % current)) allocate(this % current(12,ng,nx,ny,nz)) ! Allocate source distributions if (.not. allocated(this % cmfd_src)) allocate(this % cmfd_src(ng,nx,ny,nz)) diff --git a/src/output.F90 b/src/output.F90 index eae654787..2ee213e8e 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1035,10 +1035,10 @@ contains end if ! Get the bin for this mesh cell - - ! Left Surface matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, (/ i, j, k /)) + + ! Left Surface matching_bins(i_filter_surf) = OUT_LEFT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1047,8 +1047,6 @@ contains to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = IN_LEFT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1058,8 +1056,6 @@ contains trim(to_str(t % results(1,filter_index) % sum_sq)) ! Right Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1068,8 +1064,6 @@ contains to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1079,8 +1073,6 @@ contains trim(to_str(t % results(1,filter_index) % sum_sq)) ! Back Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_BACK filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1089,8 +1081,6 @@ contains to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = IN_BACK filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1100,8 +1090,6 @@ contains trim(to_str(t % results(1,filter_index) % sum_sq)) ! Front Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1110,8 +1098,6 @@ contains to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1121,8 +1107,6 @@ contains trim(to_str(t % results(1,filter_index) % sum_sq)) ! Bottom Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_BOTTOM filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1131,8 +1115,6 @@ contains to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = IN_BOTTOM filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1142,8 +1124,6 @@ contains trim(to_str(t % results(1,filter_index) % sum_sq)) ! Top Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 @@ -1152,8 +1132,6 @@ contains to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = IN_TOP filter_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 From 11dea6f26aaa07794a484b8d009a75cce36a8784 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 23 Aug 2016 16:01:33 -0400 Subject: [PATCH 06/12] fixed error in cmfd surface filter allocation and updated cmfd test results --- src/cmfd_data.F90 | 12 +- src/cmfd_input.F90 | 14 +- tests/test_cmfd_feed/results_true.dat | 528 ++++++++++++++++++++++ tests/test_cmfd_nofeed/results_true.dat | 528 ++++++++++++++++++++++ tests/test_score_current/results_true.dat | 2 +- 5 files changed, 1070 insertions(+), 14 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index befdc9c46..e1bf1f9c3 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -251,12 +251,12 @@ contains cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum ! Right surface - matching_bins(i_filter_surf) = OUT_RIGHT + matching_bins(i_filter_surf) = IN_RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum - matching_bins(i_filter_surf) = IN_RIGHT + matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum @@ -273,12 +273,12 @@ contains cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum ! Front surface - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = IN_FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum - matching_bins(i_filter_surf) = IN_FRONT + matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum @@ -295,12 +295,12 @@ contains cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum ! Top surface - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = IN_TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum - matching_bins(i_filter_surf) = IN_TOP + matching_bins(i_filter_surf) = OUT_TOP score_index = sum((matching_bins(1:size(t % filters)) - 1) & * t % stride) + 1 cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 50f078a18..327b55b42 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -531,15 +531,15 @@ contains allocate(SurfaceFilter :: filters(n_filters) % obj) select type(filt => filters(n_filters) % obj) type is(SurfaceFilter) - filt % n_bins = 2 * m % n_dimension - allocate(filt % surfaces(2 * m % n_dimension)) + filt % n_bins = 4 * m % n_dimension + allocate(filt % surfaces(4 * m % n_dimension)) if (m % n_dimension == 2) then - filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT, & - IN_LEFT, IN_RIGHT, IN_BACK, IN_FRONT /) + filt % surfaces = (/ OUT_LEFT, IN_LEFT, IN_RIGHT, OUT_RIGHT, & + OUT_BACK, IN_BACK, IN_FRONT, OUT_FRONT /) elseif (m % n_dimension == 3) then - filt % surfaces = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, OUT_FRONT, & - OUT_BOTTOM, OUT_TOP, IN_LEFT, IN_RIGHT, IN_BACK, IN_FRONT, & - IN_BOTTOM, IN_TOP /) + filt % surfaces = (/ OUT_LEFT, IN_LEFT, IN_RIGHT, OUT_RIGHT, & + OUT_BACK, IN_BACK, IN_FRONT, OUT_FRONT, & + OUT_BOTTOM, IN_BOTTOM, IN_TOP, OUT_TOP /) end if end select t % find_filter(FILTER_SURFACE) = n_filters diff --git a/tests/test_cmfd_feed/results_true.dat b/tests/test_cmfd_feed/results_true.dat index 04103129d..5fb87e5a3 100644 --- a/tests/test_cmfd_feed/results_true.dat +++ b/tests/test_cmfd_feed/results_true.dat @@ -136,6 +136,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +0.000000E+00 +0.000000E+00 +5.514939E+00 +1.528899E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 5.514939E+00 1.528899E+00 5.032131E+00 @@ -148,6 +160,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +2.770358E+00 +3.879191E-01 +7.294002E+00 +2.675589E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 7.294002E+00 2.675589E+00 7.036008E+00 @@ -160,6 +184,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +5.032131E+00 +1.275040E+00 +8.668860E+00 +3.776102E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 8.668860E+00 3.776102E+00 8.352414E+00 @@ -172,6 +208,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +7.036008E+00 +2.490719E+00 +9.345868E+00 +4.380719E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 9.345868E+00 4.380719E+00 9.093766E+00 @@ -184,6 +232,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +8.352414E+00 +3.501945E+00 +9.223771E+00 +4.270119E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 9.223771E+00 4.270119E+00 9.219150E+00 @@ -196,6 +256,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +9.093766E+00 +4.158282E+00 +8.530966E+00 +3.651778E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 8.530966E+00 3.651778E+00 8.690373E+00 @@ -208,6 +280,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +9.219150E+00 +4.264346E+00 +7.204424E+00 +2.604203E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 7.204424E+00 2.604203E+00 7.513640E+00 @@ -220,6 +304,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +8.690373E+00 +3.785262E+00 +5.326721E+00 +1.426975E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 5.326721E+00 1.426975E+00 5.661144E+00 @@ -232,6 +328,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +7.513640E+00 +2.833028E+00 +2.847310E+00 +4.090440E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 2.847310E+00 4.090440E-01 3.025812E+00 @@ -244,6 +352,426 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +5.661144E+00 +1.607138E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 diff --git a/tests/test_cmfd_nofeed/results_true.dat b/tests/test_cmfd_nofeed/results_true.dat index 54825f5db..91da28751 100644 --- a/tests/test_cmfd_nofeed/results_true.dat +++ b/tests/test_cmfd_nofeed/results_true.dat @@ -136,6 +136,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +0.000000E+00 +0.000000E+00 +5.555000E+00 +1.551579E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 5.555000E+00 1.551579E+00 5.095000E+00 @@ -148,6 +160,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +2.833000E+00 +4.078910E-01 +7.271000E+00 +2.659755E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 7.271000E+00 2.659755E+00 7.026000E+00 @@ -160,6 +184,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +5.095000E+00 +1.310819E+00 +8.577000E+00 +3.703215E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 8.577000E+00 3.703215E+00 8.572000E+00 @@ -172,6 +208,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +7.026000E+00 +2.486552E+00 +9.393000E+00 +4.422429E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 9.393000E+00 4.422429E+00 9.261000E+00 @@ -184,6 +232,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +8.572000E+00 +3.680852E+00 +9.265000E+00 +4.305625E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 9.265000E+00 4.305625E+00 9.303000E+00 @@ -196,6 +256,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +9.261000E+00 +4.304411E+00 +8.535000E+00 +3.659395E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 8.535000E+00 3.659395E+00 8.693000E+00 @@ -208,6 +280,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +9.303000E+00 +4.350791E+00 +7.104000E+00 +2.544182E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 7.104000E+00 2.544182E+00 7.334000E+00 @@ -220,6 +304,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +8.693000E+00 +3.799545E+00 +5.168000E+00 +1.344390E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 5.168000E+00 1.344390E+00 5.416000E+00 @@ -232,6 +328,18 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +7.334000E+00 +2.700052E+00 +2.724000E+00 +3.745680E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 2.724000E+00 3.745680E-01 2.960000E+00 @@ -244,6 +352,426 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 +5.416000E+00 +1.471086E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 diff --git a/tests/test_score_current/results_true.dat b/tests/test_score_current/results_true.dat index 68dbe069c..cf139cce2 100644 --- a/tests/test_score_current/results_true.dat +++ b/tests/test_score_current/results_true.dat @@ -1 +1 @@ -c2921f159dac64099862c1cd9c6d421c977991f621f954f893ec1351cfcea6794ca2c98c9c2dcc3411f1b8dac91ec83bd895788ba33179222c450a7df1d64f1e \ No newline at end of file +6e432700c1fb8641d106471f1fd19a0bf0f00f8b03a97131f2d73732c5a8eea48edf910c7b27b83058a7914c47c29dc7e7899b7dbb83485309b4f7d9f20471d2 \ No newline at end of file From 26bdadd79aac3712450d8c0612ac3edcb68e720f Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Mon, 29 Aug 2016 10:05:00 -0400 Subject: [PATCH 07/12] removed redundant code in tally.F90 and addressed other PR comments --- .../pythonapi/examples/tally-arithmetic.ipynb | 64 +- openmc/filter.py | 20 +- src/tally.F90 | 556 +++++------------- 3 files changed, 188 insertions(+), 452 deletions(-) diff --git a/docs/source/pythonapi/examples/tally-arithmetic.ipynb b/docs/source/pythonapi/examples/tally-arithmetic.ipynb index 088af07b3..b4b41f8d2 100644 --- a/docs/source/pythonapi/examples/tally-arithmetic.ipynb +++ b/docs/source/pythonapi/examples/tally-arithmetic.ipynb @@ -339,7 +339,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AIFw8YNfjajoIAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDgtMjNUMTU6MjQ6\nNTMtMDQ6MDANSReSAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA4LTIzVDE1OjI0OjUzLTA0OjAw\nfBSvLgAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AIHQoDEQSZg0YAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDgtMjlUMTA6MDM6\nMTctMDQ6MDAuYJb8AAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA4LTI5VDEwOjAzOjE3LTA0OjAw\nXz0uQAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] @@ -567,8 +567,8 @@ " Copyright: 2011-2016 Massachusetts Institute of Technology\n", " License: http://openmc.readthedocs.io/en/latest/license.html\n", " Version: 0.8.0\n", - " Git SHA1: 93862cea249e4441beb714e719120da3bd2b7d0a\n", - " Date/Time: 2016-08-23 15:24:53\n", + " Git SHA1: 11dea6f26aaa07794a484b8d009a75cce36a8784\n", + " Date/Time: 2016-08-29 10:03:18\n", " MPI Processes: 1\n", "\n", " ===========================================================================\n", @@ -625,20 +625,20 @@ "\n", " =======================> TIMING STATISTICS <=======================\n", "\n", - " Total time for initialization = 5.4300E-01 seconds\n", - " Reading cross sections = 3.5400E-01 seconds\n", - " Total time in simulation = 1.8555E+01 seconds\n", - " Time in transport only = 1.8464E+01 seconds\n", - " Time in inactive batches = 2.4600E+00 seconds\n", - " Time in active batches = 1.6095E+01 seconds\n", - " Time synchronizing fission bank = 1.1000E-02 seconds\n", - " Sampling source sites = 4.0000E-03 seconds\n", - " SEND/RECV source sites = 1.0000E-03 seconds\n", + " Total time for initialization = 4.6700E-01 seconds\n", + " Reading cross sections = 2.8600E-01 seconds\n", + " Total time in simulation = 1.9770E+01 seconds\n", + " Time in transport only = 1.9753E+01 seconds\n", + " Time in inactive batches = 2.4920E+00 seconds\n", + " Time in active batches = 1.7278E+01 seconds\n", + " Time synchronizing fission bank = 0.0000E+00 seconds\n", + " Sampling source sites = 0.0000E+00 seconds\n", + " SEND/RECV source sites = 0.0000E+00 seconds\n", " Time accumulating tallies = 0.0000E+00 seconds\n", - " Total time for finalization = 3.0000E-03 seconds\n", - " Total time elapsed = 1.9120E+01 seconds\n", - " Calculation Rate (inactive) = 5081.30 neutrons/second\n", - " Calculation Rate (active) = 2329.92 neutrons/second\n", + " Total time for finalization = 2.0000E-03 seconds\n", + " Total time elapsed = 2.0257E+01 seconds\n", + " Calculation Rate (inactive) = 5016.05 neutrons/second\n", + " Calculation Rate (active) = 2170.39 neutrons/second\n", "\n", " ============================> RESULTS <============================\n", "\n", @@ -731,8 +731,8 @@ " 0\n", " total\n", " (nu-fission / (absorption + current))\n", - " 1.02431\n", - " 0.00704\n", + " 1.007511\n", + " 0.007114\n", " \n", " \n", "\n", @@ -740,7 +740,7 @@ ], "text/plain": [ " nuclide score mean std. dev.\n", - "0 total (nu-fission / (absorption + current)) 1.02e+00 7.04e-03" + "0 total (nu-fission / (absorption + current)) 1.01e+00 7.11e-03" ] }, "execution_count": 24, @@ -802,8 +802,8 @@ " 6.250000e-07\n", " total\n", " (absorption + current)\n", - " 0.695303\n", - " 0.005091\n", + " 0.695851\n", + " 0.00511\n", " \n", " \n", "\n", @@ -814,7 +814,7 @@ "0 0.00e+00 6.25e-07 total (absorption + current) \n", "\n", " mean std. dev. \n", - "0 6.95e-01 5.09e-03 " + "0 6.96e-01 5.11e-03 " ] }, "execution_count": 25, @@ -1064,8 +1064,8 @@ " 6.250000e-07\n", " total\n", " (absorption + current)\n", - " 0.985102\n", - " 0.005855\n", + " 0.970693\n", + " 0.005989\n", " \n", " \n", "\n", @@ -1076,7 +1076,7 @@ "0 0.00e+00 6.25e-07 total (absorption + current) \n", "\n", " mean std. dev. \n", - "0 9.85e-01 5.86e-03 " + "0 9.71e-01 5.99e-03 " ] }, "execution_count": 29, @@ -1093,7 +1093,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The final factor is the thermalnon-leakage probability and is computed as $$P_{TNL} = \\frac{\\langle \\Sigma_a\\phi \\rangle_T}{\\langle \\Sigma_a \\phi \\rangle_T + \\langle L \\rangle_T}$$" + "The final factor is the thermal non-leakage probability and is computed as $$P_{TNL} = \\frac{\\langle \\Sigma_a\\phi \\rangle_T}{\\langle \\Sigma_a \\phi \\rangle_T + \\langle L \\rangle_T}$$" ] }, { @@ -1126,8 +1126,8 @@ " 6.250000e-07\n", " total\n", " (absorption / (absorption + current))\n", - " 0.997407\n", - " 0.008492\n", + " 0.994827\n", + " 0.008481\n", " \n", " \n", "\n", @@ -1138,7 +1138,7 @@ "0 0.00e+00 6.25e-07 total \n", "\n", " score mean std. dev. \n", - "0 (absorption / (absorption + current)) 9.97e-01 8.49e-03 " + "0 (absorption / (absorption + current)) 9.95e-01 8.48e-03 " ] }, "execution_count": 30, @@ -1190,8 +1190,8 @@ " 10000\n", " total\n", " ((((((absorption + current) * nu-fission) * ab...\n", - " 1.02431\n", - " 0.02062\n", + " 1.007511\n", + " 0.020364\n", " \n", " \n", "\n", @@ -1202,7 +1202,7 @@ "0 0.00e+00 6.25e-07 10000 total \n", "\n", " score mean std. dev. \n", - "0 ((((((absorption + current) * nu-fission) * ab... 1.02e+00 2.06e-02 " + "0 ((((((absorption + current) * nu-fission) * ab... 1.01e+00 2.04e-02 " ] }, "execution_count": 31, diff --git a/openmc/filter.py b/openmc/filter.py index 771d40401..4961ce725 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -17,6 +17,13 @@ _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup'] +_CURRENT_NAMES = {1: 'x-min out', 2: 'x-max out', + 3: 'y-min out', 4: 'y-max out', + 5: 'z-min out', 6: 'z-max out', + 7: 'x-min in', 8: 'x-max in', + 9: 'y-min in', 10: 'y-max in', + 11: 'z-min in', 12: 'z-max in'} + class Filter(object): """A filter used to constrain a tally to a specific criterion, e.g. only tally events when the particle is in a certain cell and energy range. @@ -781,18 +788,7 @@ class Filter(object): filter_bins = np.repeat(self.bins, self.stride) tile_factor = data_size / len(filter_bins) filter_bins = np.tile(filter_bins, tile_factor) - filter_bins = [x if x != 1 else 'x-min out' for x in filter_bins] - filter_bins = [x if x != 2 else 'x-max out' for x in filter_bins] - filter_bins = [x if x != 3 else 'y-min out' for x in filter_bins] - filter_bins = [x if x != 4 else 'y-max out' for x in filter_bins] - filter_bins = [x if x != 5 else 'z-min out' for x in filter_bins] - filter_bins = [x if x != 6 else 'z-max out' for x in filter_bins] - filter_bins = [x if x != 7 else 'x-min in' for x in filter_bins] - filter_bins = [x if x != 8 else 'x-max in' for x in filter_bins] - filter_bins = [x if x != 9 else 'y-min in' for x in filter_bins] - filter_bins = [x if x != 10 else 'y-max in' for x in filter_bins] - filter_bins = [x if x != 11 else 'z-min in' for x in filter_bins] - filter_bins = [x if x != 12 else 'z-max in' for x in filter_bins] + filter_bins = [_CURRENT_NAMES[x] for x in filter_bins] df = pd.concat([df, pd.DataFrame({self.type : filter_bins})]) # universe, material, surface, cell, and cellborn filters diff --git a/src/tally.F90 b/src/tally.F90 index 7f743d134..e34fecd14 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2321,6 +2321,9 @@ contains integer :: i_tally integer :: j ! loop indices integer :: k ! loop indices + integer :: d1 ! dimension index + integer :: d2 ! dimension index + integer :: d3 ! dimension index integer :: ijk0(3) ! indices of starting coordinates integer :: ijk1(3) ! indices of ending coordinates integer :: n_cross ! number of surface crossings @@ -2398,231 +2401,96 @@ contains ! ======================================================================= ! SPECIAL CASES WHERE TWO INDICES ARE THE SAME - x_same = (ijk0(1) == ijk1(1)) - y_same = (ijk0(2) == ijk1(2)) - z_same = (ijk0(3) == ijk1(3)) + ! Loop over the dimensions + do d1 = 1, 3 - if (x_same .and. y_same) then - ! Only z crossings - if (uvw(3) > 0) then - do j = ijk0(3), ijk1(3) - 1 - ijk0(3) = j + ! Get the other dimensions + d2 = mod(d1, 3) + 1 + d3 = mod(d1 + 1, 3) + 1 - ! OUT_TOP - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_TOP - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + ! If cell index in dimension d1 and d2 are the same and the cells is + ! within the mesh bounds + if (ijk0(d1) == ijk1(d1) .and. ijk0(d2) == ijk1(d2) .and. & + ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1) .and. & + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then + + ! Only d3 crossings + if (uvw(d3) > 0) then + + ! Loop over d3 cells + do j = ijk0(d3), ijk1(d3) - 1 + + ! Outward current on d3 max surface + if (j >= 1 .and. j <= m % dimension(d3)) then + ijk0(d3) = j + matching_bins(i_filter_surf) = d3 * 2 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + end if - ! IN_BOTTOM - if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 0 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) < m % dimension(3)) then - ijk0(3) = ijk0(3) + 1 - matching_bins(i_filter_surf) = IN_BOTTOM - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + ! Inward current on d3 min surface + if (j >= 0 .and. j < m % dimension(d3)) then + ijk0(d3) = j + 1 + matching_bins(i_filter_surf) = d3 * 2 + 5 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(3) = ijk0(3) - 1 - end if - end do - else - do j = ijk0(3), ijk1(3) + 1, -1 - ijk0(3) = j + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + end if + end do + else - ! OUT_BOTTOM - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BOTTOM - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if + do j = ijk0(d3), ijk1(d3) + 1, -1 - ! IN_TOP - if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) > 1 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) <= m % dimension(3) + 1) then - ijk0(3) = ijk0(3) - 1 - matching_bins(i_filter_surf) = IN_TOP - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + + ! Outward current on d3 min surface + if (j >= 1 .and. j <= m % dimension(d3)) then + ijk0(d3) = j + matching_bins(i_filter_surf) = d3 * 2 - 1 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(3) = ijk0(3) + 1 - end if - end do + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + end if + + ! Inward current on d3 max surface + if (j > 1 .and. j <= m % dimension(d3) + 1) then + ijk0(d3) = j - 1 + matching_bins(i_filter_surf) = d3 * 2 + 6 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + end if + end do + end if + cycle end if - cycle - elseif (x_same .and. z_same) then - ! Only y crossings - if (uvw(2) > 0) then - do j = ijk0(2), ijk1(2) - 1 - ijk0(2) = j - - ! OUT_FRONT - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_FRONT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_BACK - if (ijk0(1) >= 1 .and. ijk0(2) >= 0 .and. ijk0(3) >= 1 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) < m % dimension(2) .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(2) = ijk0(2) + 1 - matching_bins(i_filter_surf) = IN_BACK - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(2) = ijk0(2) - 1 - end if - end do - else - do j = ijk0(2), ijk1(2) + 1, -1 - ijk0(2) = j - - ! OUT_BACK - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BACK - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_FRONT - if (ijk0(1) >= 1 .and. ijk0(2) > 1 .and. ijk0(3) >= 1 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) + 1 .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(2) = ijk0(2) - 1 - matching_bins(i_filter_surf) = IN_FRONT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(2) = ijk0(2) + 1 - end if - end do - end if - cycle - elseif (y_same .and. z_same) then - ! Only x crossings - if (uvw(1) > 0) then - do j = ijk0(1), ijk1(1) - 1 - ijk0(1) = j - - ! OUT_RIGHT - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_RIGHT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_LEFT - if (ijk0(1) >= 0 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & - ijk0(1) < m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(1) = ijk0(1) + 1 - matching_bins(i_filter_surf) = IN_LEFT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(1) = ijk0(1) - 1 - end if - end do - else - do j = ijk0(1), ijk1(1) + 1, -1 - ijk0(1) = j - - ! OUT_LEFT - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_LEFT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_RIGHT - if (ijk0(1) > 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & - ijk0(1) <= m % dimension(1) + 1 .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(1) = ijk0(1) - 1 - matching_bins(i_filter_surf) = IN_RIGHT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(1) = ijk0(1) + 1 - end if - end do - end if - cycle - end if + end do ! ======================================================================= ! GENERIC CASE ! Bounding coordinates - do j = 1, 3 - if (uvw(j) > 0) then - xyz_cross(j) = m % lower_left(j) + ijk0(j) * m % width(j) + do d1 = 1, 3 + if (uvw(d1) > 0) then + xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1) else - xyz_cross(j) = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) + xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1) end if end do @@ -2634,11 +2502,11 @@ contains ! special case where the cosine of the angle is zero since this would ! result in a divide-by-zero. - do j = 1, 3 - if (uvw(j) == 0) then - d(j) = INFINITY + do d1 = 1, 3 + if (uvw(d1) == 0) then + d(d1) = INFINITY else - d(j) = (xyz_cross(j) - xyz0(j))/uvw(j) + d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1) end if end do @@ -2650,211 +2518,83 @@ contains ! Now use the minimum distance and direction of the particle to ! determine which surface was crossed - if (distance == d(1)) then - if (uvw(1) > 0) then + ! Loop over the dimensions + do d1 = 1, 3 - ! OUT_RIGHT - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_RIGHT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + ! Get the other dimensions + d2 = mod(d1, 3) + 1 + d3 = mod(d1 + 1, 3) + 1 + + if (distance == d(d1)) then + + ! Check dimension d2 and d3 indices + if (ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & + ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then + + if (uvw(d1) > 0) then + + ! Outward current on d1 max surface + if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1)) then + matching_bins(i_filter_surf) = d1 * 2 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + end if - ! IN_LEFT - if (ijk0(1) >= 0 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & - ijk0(1) < m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(1) = ijk0(1) + 1 - matching_bins(i_filter_surf) = IN_LEFT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + ! Inward current on d1 min surface + if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1)) then + ijk0(d1) = ijk0(d1) + 1 + matching_bins(i_filter_surf) = d1 * 2 + 5 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(1) = ijk0(1) - 1 - end if + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + ijk0(d1) = ijk0(d1) - 1 + end if - ijk0(1) = ijk0(1) + 1 - xyz_cross(1) = xyz_cross(1) + m % width(1) - else + ijk0(d1) = ijk0(d1) + 1 + xyz_cross(d1) = xyz_cross(d1) + m % width(d1) + else - ! OUT_LEFT - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_LEFT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + ! Outward current on d1 min surface + if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1)) then + matching_bins(i_filter_surf) = d1 * 2 - 1 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + end if - ! IN_RIGHT - if (ijk0(1) > 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 1 .and. & - ijk0(1) <= m % dimension(1) + 1 .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(1) = ijk0(1) - 1 - matching_bins(i_filter_surf) = IN_RIGHT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + ! Inward current on d1 max surface + if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1) then + ijk0(d1) = ijk0(d1) - 1 + matching_bins(i_filter_surf) = d1 * 2 + 6 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(1) = ijk0(1) + 1 - end if + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + ijk0(d1) = ijk0(d1) + 1 + end if - ijk0(1) = ijk0(1) - 1 - xyz_cross(1) = xyz_cross(1) - m % width(1) + ijk0(d1) = ijk0(d1) - 1 + xyz_cross(d1) = xyz_cross(d1) - m % width(d1) + end if + end if end if - elseif (distance == d(2)) then - if (uvw(2) > 0) then - - ! OUT_FRONT - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_FRONT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_BACK - if (ijk0(1) >= 1 .and. ijk0(2) >= 0 .and. ijk0(3) >= 1 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) < m % dimension(2) .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(2) = ijk0(2) + 1 - matching_bins(i_filter_surf) = IN_BACK - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(2) = ijk0(2) - 1 - end if - - ijk0(2) = ijk0(2) + 1 - xyz_cross(2) = xyz_cross(2) + m % width(2) - else - - ! OUT_BACK - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BACK - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_FRONT - if (ijk0(1) >= 1 .and. ijk0(2) > 1 .and. ijk0(3) >= 1 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) + 1 .and. & - ijk0(3) <= m % dimension(3)) then - ijk0(2) = ijk0(2) - 1 - matching_bins(i_filter_surf) = IN_FRONT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(2) = ijk0(2) + 1 - end if - - ijk0(2) = ijk0(2) - 1 - xyz_cross(2) = xyz_cross(2) - m % width(2) - end if - else if (distance == d(3)) then - if (uvw(3) > 0) then - - ! OUT_TOP - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_TOP - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_BOTTOM - if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) >= 0 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) < m % dimension(3)) then - ijk0(3) = ijk0(3) + 1 - matching_bins(i_filter_surf) = IN_BOTTOM - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(3) = ijk0(3) - 1 - end if - - ijk0(3) = ijk0(3) + 1 - xyz_cross(3) = xyz_cross(3) + m % width(3) - else - - ! OUT_BOTTOM - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = OUT_BOTTOM - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! IN_TOP - if (ijk0(1) >= 1 .and. ijk0(2) >= 1 .and. ijk0(3) > 1 .and. & - ijk0(1) <= m % dimension(1) .and. & - ijk0(2) <= m % dimension(2) .and. & - ijk0(3) <= m % dimension(3) + 1) then - ijk0(3) = ijk0(3) - 1 - matching_bins(i_filter_surf) = IN_TOP - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(3) = ijk0(3) + 1 - end if - - ijk0(3) = ijk0(3) - 1 - xyz_cross(3) = xyz_cross(3) - m % width(3) - end if - end if + end do ! Calculate new coordinates xyz0 = xyz0 + distance * uvw From 329f7d452a91c4276d598490e6877f1249ffedfe Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Mon, 29 Aug 2016 12:17:51 -0400 Subject: [PATCH 08/12] fixed bug in using cycle in nested do loop in tally.F90 --- src/tally.F90 | 149 +++++++++++++++++++++++++------------------------- 1 file changed, 76 insertions(+), 73 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index e34fecd14..b9d945778 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2410,19 +2410,17 @@ contains ! If cell index in dimension d1 and d2 are the same and the cells is ! within the mesh bounds - if (ijk0(d1) == ijk1(d1) .and. ijk0(d2) == ijk1(d2) .and. & - ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1) .and. & - ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then + if (ijk0(d1) == ijk1(d1) .and. ijk0(d2) == ijk1(d2)) then ! Only d3 crossings if (uvw(d3) > 0) then ! Loop over d3 cells do j = ijk0(d3), ijk1(d3) - 1 + ijk0(d3) = j ! Outward current on d3 max surface - if (j >= 1 .and. j <= m % dimension(d3)) then - ijk0(d3) = j + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then matching_bins(i_filter_surf) = d3 * 2 matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) @@ -2434,7 +2432,9 @@ contains end if ! Inward current on d3 min surface - if (j >= 0 .and. j < m % dimension(d3)) then + if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1) .and. & + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & + j >= 0 .and. j < m % dimension(d3)) then ijk0(d3) = j + 1 matching_bins(i_filter_surf) = d3 * 2 + 5 matching_bins(i_filter_mesh) = & @@ -2444,16 +2444,16 @@ contains !$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt + ijk0(d3) = j end if end do else do j = ijk0(d3), ijk1(d3) + 1, -1 - + ijk0(d3) = j ! Outward current on d3 min surface - if (j >= 1 .and. j <= m % dimension(d3)) then - ijk0(d3) = j + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then matching_bins(i_filter_surf) = d3 * 2 - 1 matching_bins(i_filter_mesh) = & mesh_indices_to_bin(m, ijk0) @@ -2465,7 +2465,9 @@ contains end if ! Inward current on d3 max surface - if (j > 1 .and. j <= m % dimension(d3) + 1) then + if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1) .and. & + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & + j > 1 .and. j <= m % dimension(d3) + 1) then ijk0(d3) = j - 1 matching_bins(i_filter_surf) = d3 * 2 + 6 matching_bins(i_filter_mesh) = & @@ -2475,10 +2477,11 @@ contains !$omp atomic t % results(1, filter_index) % value = & t % results(1, filter_index) % value + p % wgt + ijk0(d3) = j end if end do end if - cycle + cycle TALLY_LOOP end if end do @@ -2525,73 +2528,73 @@ contains d2 = mod(d1, 3) + 1 d3 = mod(d1 + 1, 3) + 1 + ! Check distance and dimension d2 and d3 indices if (distance == d(d1)) then - ! Check dimension d2 and d3 indices - if (ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & - ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then + if (uvw(d1) > 0) then - if (uvw(d1) > 0) then - - ! Outward current on d1 max surface - if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1)) then - matching_bins(i_filter_surf) = d1 * 2 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 + ! Outward current on d1 max surface + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then + matching_bins(i_filter_surf) = d1 * 2 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 !$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! Inward current on d1 min surface - if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1)) then - ijk0(d1) = ijk0(d1) + 1 - matching_bins(i_filter_surf) = d1 * 2 + 5 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(d1) = ijk0(d1) - 1 - end if - - ijk0(d1) = ijk0(d1) + 1 - xyz_cross(d1) = xyz_cross(d1) + m % width(d1) - else - - ! Outward current on d1 min surface - if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1)) then - matching_bins(i_filter_surf) = d1 * 2 - 1 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! Inward current on d1 max surface - if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1) then - ijk0(d1) = ijk0(d1) - 1 - matching_bins(i_filter_surf) = d1 * 2 + 6 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(d1) = ijk0(d1) + 1 - end if - - ijk0(d1) = ijk0(d1) - 1 - xyz_cross(d1) = xyz_cross(d1) - m % width(d1) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if + + ! Inward current on d1 min surface + if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1) .and. & + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & + ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then + ijk0(d1) = ijk0(d1) + 1 + matching_bins(i_filter_surf) = d1 * 2 + 5 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + ijk0(d1) = ijk0(d1) - 1 + end if + + ijk0(d1) = ijk0(d1) + 1 + xyz_cross(d1) = xyz_cross(d1) + m % width(d1) + else + + ! Outward current on d1 min surface + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then + matching_bins(i_filter_surf) = d1 * 2 - 1 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + end if + + ! Inward current on d1 max surface + if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and. & + ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & + ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then + ijk0(d1) = ijk0(d1) - 1 + matching_bins(i_filter_surf) = d1 * 2 + 6 + matching_bins(i_filter_mesh) = & + mesh_indices_to_bin(m, ijk0) + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 +!$omp atomic + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt + ijk0(d1) = ijk0(d1) + 1 + end if + + ijk0(d1) = ijk0(d1) - 1 + xyz_cross(d1) = xyz_cross(d1) - m % width(d1) end if end if end do From bc49c74d88aa5794a38a77afdd3fe98e725bfea7 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Mon, 29 Aug 2016 12:20:41 -0400 Subject: [PATCH 09/12] updated tally-arithmetic.ipynb --- .../pythonapi/examples/tally-arithmetic.ipynb | 62 +++++++++---------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/docs/source/pythonapi/examples/tally-arithmetic.ipynb b/docs/source/pythonapi/examples/tally-arithmetic.ipynb index b4b41f8d2..2fecc3c29 100644 --- a/docs/source/pythonapi/examples/tally-arithmetic.ipynb +++ b/docs/source/pythonapi/examples/tally-arithmetic.ipynb @@ -339,7 +339,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AIHQoDEQSZg0YAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDgtMjlUMTA6MDM6\nMTctMDQ6MDAuYJb8AAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA4LTI5VDEwOjAzOjE3LTA0OjAw\nXz0uQAAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AIHQwUARggawYAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDgtMjlUMTI6MjA6\nMDEtMDQ6MDBSHRTRAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA4LTI5VDEyOjIwOjAxLTA0OjAw\nI0CsbQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] @@ -567,8 +567,8 @@ " Copyright: 2011-2016 Massachusetts Institute of Technology\n", " License: http://openmc.readthedocs.io/en/latest/license.html\n", " Version: 0.8.0\n", - " Git SHA1: 11dea6f26aaa07794a484b8d009a75cce36a8784\n", - " Date/Time: 2016-08-29 10:03:18\n", + " Git SHA1: 26bdadd79aac3712450d8c0612ac3edcb68e720f\n", + " Date/Time: 2016-08-29 12:20:02\n", " MPI Processes: 1\n", "\n", " ===========================================================================\n", @@ -625,20 +625,20 @@ "\n", " =======================> TIMING STATISTICS <=======================\n", "\n", - " Total time for initialization = 4.6700E-01 seconds\n", - " Reading cross sections = 2.8600E-01 seconds\n", - " Total time in simulation = 1.9770E+01 seconds\n", - " Time in transport only = 1.9753E+01 seconds\n", - " Time in inactive batches = 2.4920E+00 seconds\n", - " Time in active batches = 1.7278E+01 seconds\n", - " Time synchronizing fission bank = 0.0000E+00 seconds\n", + " Total time for initialization = 4.6800E-01 seconds\n", + " Reading cross sections = 2.8200E-01 seconds\n", + " Total time in simulation = 1.8453E+01 seconds\n", + " Time in transport only = 1.8432E+01 seconds\n", + " Time in inactive batches = 2.5060E+00 seconds\n", + " Time in active batches = 1.5947E+01 seconds\n", + " Time synchronizing fission bank = 2.0000E-03 seconds\n", " Sampling source sites = 0.0000E+00 seconds\n", - " SEND/RECV source sites = 0.0000E+00 seconds\n", - " Time accumulating tallies = 0.0000E+00 seconds\n", - " Total time for finalization = 2.0000E-03 seconds\n", - " Total time elapsed = 2.0257E+01 seconds\n", - " Calculation Rate (inactive) = 5016.05 neutrons/second\n", - " Calculation Rate (active) = 2170.39 neutrons/second\n", + " SEND/RECV source sites = 1.0000E-03 seconds\n", + " Time accumulating tallies = 1.0000E-03 seconds\n", + " Total time for finalization = 3.0000E-03 seconds\n", + " Total time elapsed = 1.8947E+01 seconds\n", + " Calculation Rate (inactive) = 4988.03 neutrons/second\n", + " Calculation Rate (active) = 2351.54 neutrons/second\n", "\n", " ============================> RESULTS <============================\n", "\n", @@ -731,8 +731,8 @@ " 0\n", " total\n", " (nu-fission / (absorption + current))\n", - " 1.007511\n", - " 0.007114\n", + " 1.02431\n", + " 0.00704\n", " \n", " \n", "\n", @@ -740,7 +740,7 @@ ], "text/plain": [ " nuclide score mean std. dev.\n", - "0 total (nu-fission / (absorption + current)) 1.01e+00 7.11e-03" + "0 total (nu-fission / (absorption + current)) 1.02e+00 7.04e-03" ] }, "execution_count": 24, @@ -802,8 +802,8 @@ " 6.250000e-07\n", " total\n", " (absorption + current)\n", - " 0.695851\n", - " 0.00511\n", + " 0.695303\n", + " 0.005091\n", " \n", " \n", "\n", @@ -814,7 +814,7 @@ "0 0.00e+00 6.25e-07 total (absorption + current) \n", "\n", " mean std. dev. \n", - "0 6.96e-01 5.11e-03 " + "0 6.95e-01 5.09e-03 " ] }, "execution_count": 25, @@ -1064,8 +1064,8 @@ " 6.250000e-07\n", " total\n", " (absorption + current)\n", - " 0.970693\n", - " 0.005989\n", + " 0.985102\n", + " 0.005855\n", " \n", " \n", "\n", @@ -1076,7 +1076,7 @@ "0 0.00e+00 6.25e-07 total (absorption + current) \n", "\n", " mean std. dev. \n", - "0 9.71e-01 5.99e-03 " + "0 9.85e-01 5.86e-03 " ] }, "execution_count": 29, @@ -1126,8 +1126,8 @@ " 6.250000e-07\n", " total\n", " (absorption / (absorption + current))\n", - " 0.994827\n", - " 0.008481\n", + " 0.997407\n", + " 0.008492\n", " \n", " \n", "\n", @@ -1138,7 +1138,7 @@ "0 0.00e+00 6.25e-07 total \n", "\n", " score mean std. dev. \n", - "0 (absorption / (absorption + current)) 9.95e-01 8.48e-03 " + "0 (absorption / (absorption + current)) 9.97e-01 8.49e-03 " ] }, "execution_count": 30, @@ -1190,8 +1190,8 @@ " 10000\n", " total\n", " ((((((absorption + current) * nu-fission) * ab...\n", - " 1.007511\n", - " 0.020364\n", + " 1.02431\n", + " 0.02062\n", " \n", " \n", "\n", @@ -1202,7 +1202,7 @@ "0 0.00e+00 6.25e-07 10000 total \n", "\n", " score mean std. dev. \n", - "0 ((((((absorption + current) * nu-fission) * ab... 1.01e+00 2.04e-02 " + "0 ((((((absorption + current) * nu-fission) * ab... 1.02e+00 2.06e-02 " ] }, "execution_count": 31, From f63cad37bf8ff9c41aa30f8b17b92e301639a1b7 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 30 Aug 2016 08:09:03 -0400 Subject: [PATCH 10/12] removed special case for tally surface current --- src/tally.F90 | 109 ++++---------------------------------------------- 1 file changed, 7 insertions(+), 102 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index b9d945778..d5af4ddc1 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2320,7 +2320,6 @@ contains integer :: i integer :: i_tally integer :: j ! loop indices - integer :: k ! loop indices integer :: d1 ! dimension index integer :: d2 ! dimension index integer :: d3 ! dimension index @@ -2340,9 +2339,6 @@ contains real(8) :: filt_score ! score applied by filters logical :: start_in_mesh ! particle's starting xyz in mesh? logical :: end_in_mesh ! particle's ending xyz in mesh? - logical :: x_same ! same starting/ending x index (i) - logical :: y_same ! same starting/ending y index (j) - logical :: z_same ! same starting/ending z index (k) type(TallyObject), pointer :: t type(RegularMesh), pointer :: m @@ -2398,96 +2394,6 @@ contains if (matching_bins(i_filter_energy) == NO_BIN_FOUND) cycle end if - ! ======================================================================= - ! SPECIAL CASES WHERE TWO INDICES ARE THE SAME - - ! Loop over the dimensions - do d1 = 1, 3 - - ! Get the other dimensions - d2 = mod(d1, 3) + 1 - d3 = mod(d1 + 1, 3) + 1 - - ! If cell index in dimension d1 and d2 are the same and the cells is - ! within the mesh bounds - if (ijk0(d1) == ijk1(d1) .and. ijk0(d2) == ijk1(d2)) then - - ! Only d3 crossings - if (uvw(d3) > 0) then - - ! Loop over d3 cells - do j = ijk0(d3), ijk1(d3) - 1 - ijk0(d3) = j - - ! Outward current on d3 max surface - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = d3 * 2 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! Inward current on d3 min surface - if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1) .and. & - ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & - j >= 0 .and. j < m % dimension(d3)) then - ijk0(d3) = j + 1 - matching_bins(i_filter_surf) = d3 * 2 + 5 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(d3) = j - end if - end do - else - - do j = ijk0(d3), ijk1(d3) + 1, -1 - ijk0(d3) = j - - ! Outward current on d3 min surface - if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = d3 * 2 - 1 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - end if - - ! Inward current on d3 max surface - if (ijk0(d1) >= 1 .and. ijk0(d1) <= m % dimension(d1) .and. & - ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. & - j > 1 .and. j <= m % dimension(d3) + 1) then - ijk0(d3) = j - 1 - matching_bins(i_filter_surf) = d3 * 2 + 6 - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0) - filter_index = sum((matching_bins(1:size(t % filters)) - 1) & - * t % stride) + 1 -!$omp atomic - t % results(1, filter_index) % value = & - t % results(1, filter_index) % value + p % wgt - ijk0(d3) = j - end if - end do - end if - cycle TALLY_LOOP - end if - end do - - ! ======================================================================= - ! GENERIC CASE - ! Bounding coordinates do d1 = 1, 3 if (uvw(d1) > 0) then @@ -2497,14 +2403,13 @@ contains end if end do - do k = 1, n_cross + do j = 1, n_cross ! Reset scoring bin index matching_bins(i_filter_surf) = 0 ! Calculate distance to each bounding surface. We need to treat ! special case where the cosine of the angle is zero since this would ! result in a divide-by-zero. - do d1 = 1, 3 if (uvw(d1) == 0) then d(d1) = INFINITY @@ -2514,13 +2419,10 @@ contains end do ! Determine the closest bounding surface of the mesh cell by - ! calculating the minimum distance - + ! calculating the minimum distance. Then use the minimum distance and + ! direction of the particle to determine which surface was crossed. distance = minval(d) - ! Now use the minimum distance and direction of the particle to - ! determine which surface was crossed - ! Loop over the dimensions do d1 = 1, 3 @@ -2528,9 +2430,10 @@ contains d2 = mod(d1, 3) + 1 d3 = mod(d1 + 1, 3) + 1 - ! Check distance and dimension d2 and d3 indices + ! Check whether distance is the shortest distance if (distance == d(d1)) then + ! Check whether particle is moving in positive d1 direction if (uvw(d1) > 0) then ! Outward current on d1 max surface @@ -2563,6 +2466,8 @@ contains ijk0(d1) = ijk0(d1) + 1 xyz_cross(d1) = xyz_cross(d1) + m % width(d1) + + ! The particle is moving in the negative d1 direction else ! Outward current on d1 min surface From 368e5cc87a5ad705c992f9b907f262163c960483 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 31 Aug 2016 07:24:31 -0400 Subject: [PATCH 11/12] fixed issue in mdgxs.get_pandas_dataframe(...) --- openmc/mgxs/mdgxs.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py index 5e5d840b6..7a76a6c9d 100644 --- a/openmc/mgxs/mdgxs.py +++ b/openmc/mgxs/mdgxs.py @@ -343,7 +343,7 @@ class MDGXS(MGXS): # Construct a collection of the nuclides to retrieve from the xs tally if self.by_nuclide: if nuclides == 'all' or nuclides == 'sum' or nuclides == ['sum']: - query_nuclides = self.get_all_nuclides() + query_nuclides = self.get_nuclides() else: query_nuclides = nuclides else: @@ -553,7 +553,7 @@ class MDGXS(MGXS): # Construct a collection of the nuclides to report if self.by_nuclide: if nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() elif nuclides == 'sum': nuclides = ['sum'] else: @@ -751,25 +751,27 @@ class MDGXS(MGXS): if self.by_nuclide and nuclides == 'sum': # Use tally summation to sum across all nuclides - query_nuclides = self.get_all_nuclides() - xs_tally = self.xs_tally.summation(nuclides=query_nuclides) + query_nuclides = [nuclides] + xs_tally = self.xs_tally.summation(nuclides=self.get_nuclides()) df = xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) # Remove nuclide column since it is homogeneous and redundant if self.domain_type == 'mesh': - df.drop('nuclide', axis=1, level=0, inplace=True) + df.drop('sum(nuclide)', axis=1, level=0, inplace=True) else: - df.drop('nuclide', axis=1, inplace=True) + df.drop('sum(nuclide)', axis=1, inplace=True) # If the user requested a specific set of nuclides elif self.by_nuclide and nuclides != 'all': + query_nuclides = nuclides xs_tally = self.xs_tally.get_slice(nuclides=nuclides) df = xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) # If the user requested all nuclides, keep nuclide column in dataframe else: + query_nuclides = nuclides df = self.xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths) @@ -781,7 +783,7 @@ class MDGXS(MGXS): # Override energy groups bounds with indices all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int) - all_groups = np.repeat(all_groups, self.num_nuclides) + all_groups = np.repeat(all_groups, len(query_nuclides)) if 'energy low [MeV]' in df and 'energyout low [MeV]' in df: df.rename(columns={'energy low [MeV]': 'group in'}, inplace=True) @@ -1239,7 +1241,7 @@ class ChiDelayed(MDGXS): delayed_nu_fission_out = self.tallies['delayed-nu-fission-out'] # Sum out all nuclides - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() delayed_nu_fission_in = delayed_nu_fission_in.summation\ (nuclides=nuclides) delayed_nu_fission_out = delayed_nu_fission_out.summation\ @@ -1261,7 +1263,7 @@ class ChiDelayed(MDGXS): # Get chi delayed for all nuclides in the domain elif nuclides == 'all': - nuclides = self.get_all_nuclides() + nuclides = self.get_nuclides() xs = self.xs_tally.get_values(filters=filters, filter_bins=filter_bins, nuclides=nuclides, value=value) From fc1d856cf85cbd9918f0f5c2c9be887afa9abd3a Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 1 Sep 2016 10:29:48 -0400 Subject: [PATCH 12/12] fixed bug in MDGXS.get_pandas_dataframe(...) --- openmc/mgxs/mdgxs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py index 7a76a6c9d..8941d4c6f 100644 --- a/openmc/mgxs/mdgxs.py +++ b/openmc/mgxs/mdgxs.py @@ -771,7 +771,7 @@ class MDGXS(MGXS): # If the user requested all nuclides, keep nuclide column in dataframe else: - query_nuclides = nuclides + query_nuclides = self.nuclides df = self.xs_tally.get_pandas_dataframe( distribcell_paths=distribcell_paths)