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

removed redundant code in tally.F90 and addressed other PR comments

Browse source
This commit is contained in:
Sam Shaner 2016-08-29 10:05:00 -04:00
parent 11dea6f26a
commit 26bdadd79a
3 changed files with 188 additions and 452 deletions
Download patch file Download diff file Expand all files Collapse all files

64
docs/source/pythonapi/examples/tally-arithmetic.ipynb
View file

@ -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": [
"<IPython.core.display.Image object>"
]
@ -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 @@
" <th>0</th>\n",
" <td>total</td>\n",
" <td>(nu-fission / (absorption + current))</td>\n",
" <td>1.02431</td>\n",
" <td>0.00704</td>\n",
" <td>1.007511</td>\n",
" <td>0.007114</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\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 @@
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>(absorption + current)</td>\n",
" <td>0.695303</td>\n",
" <td>0.005091</td>\n",
" <td>0.695851</td>\n",
" <td>0.00511</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\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 @@
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>(absorption + current)</td>\n",
" <td>0.985102</td>\n",
" <td>0.005855</td>\n",
" <td>0.970693</td>\n",
" <td>0.005989</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\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 @@
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>(absorption / (absorption + current))</td>\n",
" <td>0.997407</td>\n",
" <td>0.008492</td>\n",
" <td>0.994827</td>\n",
" <td>0.008481</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\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 @@
" <td>10000</td>\n",
" <td>total</td>\n",
" <td>((((((absorption + current) * nu-fission) * ab...</td>\n",
" <td>1.02431</td>\n",
" <td>0.02062</td>\n",
" <td>1.007511</td>\n",
" <td>0.020364</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\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,

20
openmc/filter.py
View file

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

556
src/tally.F90
View file

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