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