From 355fb6f30d6c0a93d606aa2ad337b31189037709 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Sun, 31 Jul 2016 10:20:50 -0400 Subject: [PATCH 1/2] addressed PR comments --- docs/source/pythonapi/examples/mgxs-part-iv.ipynb | 4 ++-- openmc/mgxs/library.py | 3 +-- openmc/mgxs/mgxs.py | 8 +++++--- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb index 5b3bae495..eaae92f7c 100644 --- a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb +++ b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb @@ -519,9 +519,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now we must specify the type of domain over which we would like the `Library` to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the `Library` supports \"material\" \"cell\", \"universe\", and \"mesh\" domain types. In this simple example, we wish to compute multi-group cross sections only for each material andtherefore will use a \"material\" domain type.\n", + "Now we must specify the type of domain over which we would like the `Library` to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the `Library` supports \"material\" \"cell\", \"universe\", and \"mesh\" domain types. In this simple example, we wish to compute multi-group cross sections only for each material and therefore will use a \"material\" domain type.\n", "\n", - "**Note:** By default, the `Library` class will instantiate `MGXS` objects for each and every domain (material, cell, universe, or mesh cell) in the geometry of interest. However, one may specify a subset of these domains to the `Library.domains` property." + "**Note:** By default, the `Library` class will instantiate `MGXS` objects for each and every domain (material, cell, universe, or mesh) in the geometry of interest. However, one may specify a subset of these domains to the `Library.domains` property." ] }, { diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 50d568a2d..84571ce23 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -56,8 +56,7 @@ class Library(object): The types of cross sections in the library (e.g., ['total', 'scatter']) domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} Domain type for spatial homogenization - domains : Iterable of openmc.Material, openmc.Cell, openmc.Universe or - openmc.Mesh + domains : Iterable of openmc.Material, openmc.Cell, openmc.Universe or openmc.Mesh The spatial domain(s) for which MGXS in the Library are computed correction : {'P0', None} Apply the P0 correction to scattering matrices if set to 'P0' diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index f19b50e18..e1c2e6221 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -439,8 +439,7 @@ class MGXS(object): ---------- mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', 'chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission'} The type of multi-group cross section object to return - domain : openmc.Material or openmc.Cell or openmc.Universe or - openmc.Mesh + domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh The domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} The domain type for spatial homogenization @@ -1491,7 +1490,10 @@ class MGXS(object): distribcell_paths=distribcell_paths) # Remove nuclide column since it is homogeneous and redundant - df.drop('nuclide', axis=1, inplace=True) + if self.domain_type == 'mesh': + df.drop('nuclide', axis=1, level=0, inplace=True) + else: + df.drop('nuclide', axis=1, inplace=True) # If the user requested a specific set of nuclides elif self.by_nuclide and nuclides != 'all': From 5074c9831f440a14c81cd84607a132abd24be3d0 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Tue, 2 Aug 2016 03:05:16 -0400 Subject: [PATCH 2/2] reworked surface current tally indexing and fixed issue in Python API current tallies --- openmc/filter.py | 12 ++++++ openmc/statepoint.py | 2 +- src/cmfd_data.F90 | 91 ++++++++++++++++++++++++++++-------------- src/cmfd_input.F90 | 6 +-- src/constants.F90 | 12 +++--- src/input_xml.F90 | 20 +++------- src/mesh.F90 | 15 ++----- src/output.F90 | 60 +++++----------------------- src/tally.F90 | 90 +++++++++++++++++++++--------------------- src/trigger.F90 | 94 +++----------------------------------------- 10 files changed, 151 insertions(+), 251 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index d6ee70f14..9e05cff63 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -772,6 +772,18 @@ class Filter(object): df.loc[:, self.type + ' low'] = lo_bins df.loc[:, self.type + ' high'] = hi_bins + elif self.type == 'surface': + 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] + df = pd.concat([df, pd.DataFrame({self.type : filter_bins})]) + # universe, material, surface, cell, and cellborn filters else: filter_bins = np.repeat(self.bins, self.stride) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 633774665..10d91d194 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -681,7 +681,7 @@ class StatePoint(object): if tally_filter.type == 'surface': surface_ids = [] for bin in tally_filter.bins: - surface_ids.append(summary.surfaces[bin].id) + surface_ids.append(bin) tally_filter.bins = surface_ids if tally_filter.type in ['cell', 'distribcell']: diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 347351e31..c07569ae5 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -50,9 +50,9 @@ contains subroutine compute_xs() use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, & - FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT, & - OUT_FRONT, IN_TOP, OUT_TOP, CMFD_NOACCEL, ZERO, & - ONE, TINY_BIT + FILTER_SURFACE, OUT_LEFT, OUT_RIGHT, OUT_BACK, & + OUT_FRONT, OUT_BOTTOM, OUT_TOP, CMFD_NOACCEL, & + ZERO, ONE, TINY_BIT use error, only: fatal_error use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes,& matching_bins @@ -227,60 +227,91 @@ contains ! Left surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & - (/ i-1, j, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_RIGHT + (/ i, j, k /)) + matching_bins(i_filter_surf) = OUT_LEFT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum - matching_bins(i_filter_surf) = OUT_RIGHT - score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(2,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) = OUT_RIGHT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming + cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum + end if ! 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) = OUT_LEFT + score_index = sum((matching_bins(1:t%n_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 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_RIGHT - score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum + (/ i, j, k /) ) matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing 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-1, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_FRONT + (/ i, j, k /)) + matching_bins(i_filter_surf) = OUT_BACK score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum - matching_bins(i_filter_surf) = OUT_FRONT - score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(6,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) = OUT_FRONT + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming + cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum + end if ! 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) = OUT_BACK + score_index = sum((matching_bins(1:t%n_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 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_FRONT - score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum + (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing 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-1 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_TOP + (/ i, j, k /)) + matching_bins(i_filter_surf) = OUT_BOTTOM score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum - matching_bins(i_filter_surf) = OUT_TOP - score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(10,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) = OUT_TOP + score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming + cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum + end if ! 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) = OUT_BOTTOM + score_index = sum((matching_bins(1:t%n_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 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_TOP - score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum + (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_TOP score_index = sum((matching_bins(1:t%n_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 f69c09fe1..e63a1fcea 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -526,11 +526,11 @@ contains filters(n_filters) % n_bins = 2 * m % n_dimension allocate(filters(n_filters) % int_bins(2 * m % n_dimension)) if (m % n_dimension == 2) then - filters(n_filters) % int_bins = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, & + filters(n_filters) % int_bins = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, & OUT_FRONT /) elseif (m % n_dimension == 3) then - filters(n_filters) % int_bins = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, & - OUT_FRONT, IN_TOP, OUT_TOP /) + filters(n_filters) % int_bins = (/ OUT_LEFT, OUT_RIGHT, OUT_BACK, & + OUT_FRONT, OUT_BOTTOM, OUT_TOP /) end if t % find_filter(FILTER_SURFACE) = n_filters diff --git a/src/constants.F90 b/src/constants.F90 index a22c9ac05..ebb8322e9 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -356,12 +356,12 @@ module constants ! Tally surface current directions integer, parameter :: & - IN_RIGHT = 1, & - OUT_RIGHT = 2, & - IN_FRONT = 3, & - OUT_FRONT = 4, & - IN_TOP = 5, & - OUT_TOP = 6 + 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 ! Tally trigger types and threshold integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3e5a68712..b94c41a65 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3029,9 +3029,7 @@ contains // " specified on tally " // trim(to_str(t % id))) end if - ! Determine number of bins -- this is assuming that the tally is - ! a volume tally and not a surface current tally. If it is a - ! surface current tally, the number of bins will get reset later + ! Determine number of bins t % filters(j) % n_bins = product(m % dimension) ! Allocate and store index of mesh @@ -3646,10 +3644,6 @@ contains &same tally as surface currents") end if - ! Since the number of bins for the mesh filter was already set - ! assuming it was a volume tally, we need to adjust the number - ! of bins - ! Get index of mesh filter k = t % find_filter(FILTER_MESH) @@ -3663,10 +3657,6 @@ contains i_mesh = t % filters(k) % int_bins(1) m => meshes(i_mesh) - ! We need to increase the dimension by one since we also need - ! currents coming into and out of the boundary mesh cells. - t % filters(k) % n_bins = product(m % dimension + 1) - ! Copy filters to temporary array allocate(filters(t % n_filters + 1)) filters(1:t % n_filters) = t % filters @@ -3682,11 +3672,11 @@ contains allocate(t % filters(t % n_filters) % int_bins(& 2 * m % n_dimension)) if (m % n_dimension == 2) then - t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT /) + t % filters(t % n_filters) % int_bins = (/ OUT_LEFT, & + OUT_RIGHT, OUT_BACK, OUT_FRONT /) elseif (m % n_dimension == 3) then - t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /) + t % filters(t % n_filters) % int_bins = (/ OUT_LEFT, & + OUT_RIGHT, OUT_BACK, OUT_FRONT, OUT_BOTTOM, OUT_TOP /) end if t % find_filter(FILTER_SURFACE) = t % n_filters diff --git a/src/mesh.F90 b/src/mesh.F90 index 4ec84345b..0543eaa9f 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -93,29 +93,20 @@ contains ! use in a TallyObject results array !=============================================================================== - pure function mesh_indices_to_bin(m, ijk, surface_current) result(bin) + pure function mesh_indices_to_bin(m, ijk) result(bin) type(RegularMesh), intent(in) :: m integer, intent(in) :: ijk(:) - logical, intent(in), optional :: surface_current integer :: bin integer :: n_y ! number of mesh cells in y direction integer :: n_z ! number of mesh cells in z direction - if (present(surface_current)) then - n_y = m % dimension(2) + 1 - else - n_y = m % dimension(2) - end if + n_y = m % dimension(2) if (m % n_dimension == 2) then bin = (ijk(1) - 1)*n_y + ijk(2) elseif (m % n_dimension == 3) then - if (present(surface_current)) then - n_z = m % dimension(3) + 1 - else - n_z = m % dimension(3) - end if + n_z = m % dimension(3) bin = (ijk(1) - 1)*n_y*n_z + (ijk(2) - 1)*n_z + ijk(3) end if diff --git a/src/output.F90 b/src/output.F90 index 22733a305..efe75ad8b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1047,31 +1047,17 @@ contains ! Left Surface matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_RIGHT + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = OUT_LEFT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_surf) = OUT_RIGHT - filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from 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 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_RIGHT - filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Right", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) - + mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & @@ -1081,31 +1067,17 @@ contains ! Back Surface matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_FRONT + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = OUT_BACK filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_surf) = OUT_FRONT - filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from 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 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_FRONT - filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Front", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) - + mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & @@ -1115,31 +1087,17 @@ contains ! Bottom Surface matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_TOP + mesh_indices_to_bin(m, (/ i, j, k /)) + matching_bins(i_filter_surf) = OUT_BOTTOM filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) - matching_bins(i_filter_surf) = OUT_TOP - filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from 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 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_TOP - filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Top", & - to_str(t % results(1,filter_index) % sum), & - trim(to_str(t % results(1,filter_index) % sum_sq)) - + mesh_indices_to_bin(m, (/ i, j, k /)) matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') & diff --git a/src/tally.F90 b/src/tally.F90 index 94143184a..5f9268216 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2819,10 +2819,10 @@ contains if (uvw(3) > 0) then do j = ijk0(3), ijk1(3) - 1 ijk0(3) = j - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & @@ -2830,12 +2830,12 @@ contains end if end do else - do j = ijk0(3) - 1, ijk1(3), -1 + do j = ijk0(3), ijk1(3) + 1, -1 ijk0(3) = j - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = IN_TOP + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & @@ -2849,10 +2849,10 @@ contains if (uvw(2) > 0) then do j = ijk0(2), ijk1(2) - 1 ijk0(2) = j - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & @@ -2860,12 +2860,12 @@ contains end if end do else - do j = ijk0(2) - 1, ijk1(2), -1 + do j = ijk0(2), ijk1(2) + 1, -1 ijk0(2) = j - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = IN_FRONT + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & @@ -2879,10 +2879,10 @@ contains if (uvw(1) > 0) then do j = ijk0(1), ijk1(1) - 1 ijk0(1) = j - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & @@ -2890,12 +2890,12 @@ contains end if end do else - do j = ijk0(1) - 1, ijk1(1), -1 + do j = ijk0(1), ijk1(1) + 1, -1 ijk0(1) = j - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = IN_RIGHT + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 !$omp atomic t % results(1, filter_index) % value = & @@ -2946,67 +2946,67 @@ contains if (uvw(1) > 0) then ! Crossing into right mesh cell -- this is treated as outgoing ! current from (i,j,k) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) 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 incoming - ! current in (i-1,j,k) + ! Crossing into left mesh cell -- this is treated as outgoing + ! current in (i,j,k) + 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) + end if ijk0(1) = ijk0(1) - 1 xyz_cross(1) = xyz_cross(1) - m % width(1) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = IN_RIGHT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) - end if 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) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) 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 incoming - ! current in (i,j-1,k) + ! Crossing into back mesh cell -- this is treated as outgoing + ! current in (i,j,k) + 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) + end if ijk0(2) = ijk0(2) - 1 xyz_cross(2) = xyz_cross(2) - m % width(2) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = IN_FRONT - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) - end if 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) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + 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 + 1, .true.) + mesh_indices_to_bin(m, ijk0) 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 incoming - ! current in (i,j,k-1) + ! Crossing into bottom mesh cell -- this is treated as outgoing + ! current in (i,j,k) + 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) + end if ijk0(3) = ijk0(3) - 1 xyz_cross(3) = xyz_cross(3) - m % width(3) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - matching_bins(i_filter_surf) = IN_TOP - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) - end if end if end if diff --git a/src/trigger.F90 b/src/trigger.F90 index ed362154b..7d488acb7 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -324,10 +324,11 @@ contains matching_bins(i_filter_ein) = l end if - ! Left Surface matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_RIGHT + mesh_indices_to_bin(m, (/ i, j, k /) + 1) + + ! Left Surface + matching_bins(i_filter_surf) = OUT_LEFT filter_index = & sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -339,33 +340,7 @@ contains end if trigger % variance = std_dev**2 - matching_bins(i_filter_surf) = OUT_RIGHT - filter_index = & - sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - ! Right Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_RIGHT - filter_index = & - sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - matching_bins(i_filter_surf) = OUT_RIGHT filter_index = & sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 @@ -379,22 +354,7 @@ contains trigger % variance = trigger % std_dev**2 ! Back Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_FRONT - filter_index = & - sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - - - matching_bins(i_filter_surf) = OUT_FRONT + matching_bins(i_filter_surf) = OUT_BACK filter_index = & sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -407,20 +367,6 @@ contains trigger % variance = trigger % std_dev**2 ! Front Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_FRONT - filter_index = & - sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - matching_bins(i_filter_surf) = OUT_FRONT filter_index = & sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 @@ -434,21 +380,7 @@ contains trigger % variance = trigger % std_dev**2 ! Bottom Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) - matching_bins(i_filter_surf) = IN_TOP - filter_index = & - sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - - matching_bins(i_filter_surf) = OUT_TOP + matching_bins(i_filter_surf) = OUT_BOTTOM filter_index = & sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) @@ -461,20 +393,6 @@ contains trigger % variance = trigger % std_dev**2 ! Top Surface - matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) - matching_bins(i_filter_surf) = IN_TOP - filter_index = & - sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1 - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - matching_bins(i_filter_surf) = OUT_TOP filter_index = & sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1