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

Merge remote-tracking branch 'origin/meshout-filter' into mesh-mgxs

Browse source
This commit is contained in:
Sam Shaner 2016-08-02 08:13:35 -04:00
commit e8819e6a77
13 changed files with 159 additions and 258 deletions
Download patch file Download diff file Expand all files Collapse all files

4
docs/source/pythonapi/examples/mgxs-part-iv.ipynb
View file

@ -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."
]
},
{

12
openmc/filter.py
View file

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

3
openmc/mgxs/library.py
View file

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

8
openmc/mgxs/mgxs.py
View file

@ -438,8 +438,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
@ -1490,7 +1489,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':

2
openmc/statepoint.py
View file

@ -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']:

91
src/cmfd_data.F90
View file

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

6
src/cmfd_input.F90
View file

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

12
src/constants.F90
View file

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

20
src/input_xml.F90
View file

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

15
src/mesh.F90
View file

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

60
src/output.F90
View file

@ -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)') &

90
src/tally.F90
View file

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

94
src/trigger.F90
View file

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