diff --git a/examples/xml/pincell/geometry.xml b/examples/xml/pincell/geometry.xml
index f67f9e74c2..f4e1bf1dd4 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 443af9cde2..9402b979a9 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 73242b9136..2e39c083d4 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 d74c1b17e0..0a1490bb97 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 5fe44ad436..77e8f4348c 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 6d5f2fc107..ceb68b9e2d 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 a60d8196ea..8de1983865 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 9e23f11d88..7aacd67136 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 3c3dc6a693..6bb36d9ce1 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 4af03a2dea..6f886cd8e0 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)