crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Move most of mesh module into mesh_header as type-bound procedures

Browse source
This commit is contained in:
Paul Romano 2017-08-07 13:40:49 -05:00
parent 2864915409
commit d2d70a3d42
11 changed files with 384 additions and 390 deletions
Download patch file Download diff file Expand all files Collapse all files

7
src/cmfd_data.F90
View file

@ -58,7 +58,6 @@ contains
use error, only: fatal_error
use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes, &
filters, filter_matches
use mesh, only: mesh_indices_to_bin
use mesh_header, only: RegularMesh
use string, only: to_str
use tally_header, only: TallyObject
@ -166,7 +165,7 @@ contains
! Get bin number for mesh indices
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m,ijk)
m % get_bin_from_indices(ijk)
! Apply energy in filter
if (energy_filters) then
@ -217,7 +216,7 @@ contains
! Get bin number for mesh indices
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m,ijk)
m % get_bin_from_indices(ijk)
if (energy_filters) then
! Apply energy in filter
@ -266,7 +265,7 @@ contains
! Get the bin for this mesh cell
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, (/ i, j, k /))
m % get_bin_from_indices([ i, j, k ])
score_index = 1
do l = 1, size(t % filter)

4
src/cmfd_execute.F90
View file

@ -218,7 +218,7 @@ contains
use error, only: warning, fatal_error
use global, only: meshes, source_bank, work, n_user_meshes, cmfd
use mesh_header, only: RegularMesh
use mesh, only: count_bank_sites, get_mesh_indices
use mesh, only: count_bank_sites
use message_passing
use string, only: to_str
@ -295,7 +295,7 @@ contains
do i = 1, int(work,4)
! Determine spatial bin
call get_mesh_indices(m, source_bank(i) % xyz, ijk, in_mesh)
call m % get_indices(source_bank(i) % xyz, ijk, in_mesh)
! Determine energy bin
n_groups = size(cmfd % egrid) - 1

334
src/mesh.F90
View file

@ -10,116 +10,6 @@ module mesh
contains
!===============================================================================
! GET_MESH_BIN determines the tally bin for a particle in a structured mesh
!===============================================================================
pure subroutine get_mesh_bin(m, xyz, bin)
type(RegularMesh), intent(in) :: m ! mesh pointer
real(8), intent(in) :: xyz(:) ! coordinates
integer, intent(out) :: bin ! tally bin
integer :: n ! size of mesh
integer :: d ! mesh dimension index
integer :: ijk(3) ! indices in mesh
logical :: in_mesh ! was given coordinate in mesh at all?
! Get number of dimensions
n = m % n_dimension
! Loop over the dimensions of the mesh
do d = 1, n
! Check for cases where particle is outside of mesh
if (xyz(d) < m % lower_left(d)) then
bin = NO_BIN_FOUND
return
elseif (xyz(d) > m % upper_right(d)) then
bin = NO_BIN_FOUND
return
end if
end do
! Determine indices
call get_mesh_indices(m, xyz, ijk, in_mesh)
! Convert indices to bin
if (in_mesh) then
bin = mesh_indices_to_bin(m, ijk)
else
bin = NO_BIN_FOUND
end if
end subroutine get_mesh_bin
!===============================================================================
! GET_MESH_INDICES determines the indices of a particle in a structured mesh
!===============================================================================
pure subroutine get_mesh_indices(m, xyz, ijk, in_mesh)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz(:) ! coordinates to check
integer, intent(out) :: ijk(:) ! indices in mesh
logical, intent(out) :: in_mesh ! were given coords in mesh?
! Find particle in mesh
ijk(:m % n_dimension) = ceiling((xyz(:m % n_dimension) - m % lower_left)/m % width)
! Determine if particle is in mesh
if (any(ijk(:m % n_dimension) < 1) .or. &
any(ijk(:m % n_dimension) > m % dimension)) then
in_mesh = .false.
else
in_mesh = .true.
end if
end subroutine get_mesh_indices
!===============================================================================
! MESH_INDICES_TO_BIN maps (i), (i,j), or (i,j,k) indices to a single bin number
! for use in a TallyObject results array
!===============================================================================
pure function mesh_indices_to_bin(m, ijk) result(bin)
type(RegularMesh), intent(in) :: m
integer, intent(in) :: ijk(:)
integer :: bin
if (m % n_dimension == 1) then
bin = ijk(1)
elseif (m % n_dimension == 2) then
bin = (ijk(2) - 1) * m % dimension(1) + ijk(1)
elseif (m % n_dimension == 3) then
bin = ((ijk(3) - 1) * m % dimension(2) + (ijk(2) - 1)) &
* m % dimension(1) + ijk(1)
end if
end function mesh_indices_to_bin
!===============================================================================
! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to
! (i), (i,j), or (i,j,k) indices
!===============================================================================
pure subroutine bin_to_mesh_indices(m, bin, ijk)
type(RegularMesh), intent(in) :: m
integer, intent(in) :: bin
integer, intent(out) :: ijk(:)
if (m % n_dimension == 1) then
ijk(1) = bin
else if (m % n_dimension == 2) then
ijk(1) = mod(bin - 1, m % dimension(1)) + 1
ijk(2) = (bin - 1)/m % dimension(1) + 1
else if (m % n_dimension == 3) then
ijk(1) = mod(bin - 1, m % dimension(1)) + 1
ijk(2) = mod(bin - 1, m % dimension(1) * m % dimension(2)) &
/ m % dimension(1) + 1
ijk(3) = (bin - 1)/(m % dimension(1) * m % dimension(2)) + 1
end if
end subroutine bin_to_mesh_indices
!===============================================================================
! COUNT_BANK_SITES determines the number of fission bank sites in each cell of a
! given mesh as well as an optional energy group structure. This can be used for
@ -172,7 +62,7 @@ contains
! loop over fission sites and count how many are in each mesh box
FISSION_SITES: do i = 1, n_sites
! determine scoring bin for entropy mesh
call get_mesh_indices(m, bank_array(i) % xyz, ijk, in_mesh)
call m % get_indices(bank_array(i) % xyz, ijk, in_mesh)
! if outside mesh, skip particle
if (.not. in_mesh) then
@ -215,226 +105,4 @@ contains
end subroutine count_bank_sites
!===============================================================================
! MESH_INTERSECTS determines if a line between xyz0 and xyz1 intersects the
! outer boundary of the given mesh. This is important for determining whether a
! track will score to a mesh tally.
!===============================================================================
pure function mesh_intersects_1d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(1)
real(8), intent(in) :: xyz1(1)
logical :: intersects
real(8) :: x0 ! track start point
real(8) :: x1 ! track end point
real(8) :: xm0 ! lower-left coordinates of mesh
real(8) :: xm1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
! Copy coordinates of ending point
x1 = xyz1(1)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
! Set default value for intersects
intersects = .false.
! Check if line intersects left surface
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
intersects = .true.
return
end if
! Check if line intersects right surface
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
intersects = .true.
return
end if
end function mesh_intersects_1d
pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(2)
real(8), intent(in) :: xyz1(2)
logical :: intersects
real(8) :: x0, y0 ! track start point
real(8) :: x1, y1 ! track end point
real(8) :: xi, yi ! track intersection point with mesh
real(8) :: xm0, ym0 ! lower-left coordinates of mesh
real(8) :: xm1, ym1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
y0 = xyz0(2)
! Copy coordinates of ending point
x1 = xyz1(1)
y1 = xyz1(2)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
ym0 = m % lower_left(2)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
ym1 = m % upper_right(2)
! Set default value for intersects
intersects = .false.
! Check if line intersects left surface -- calculate the intersection point
! y
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! x
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection
! point y
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! x
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_2d
pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(3)
real(8), intent(in) :: xyz1(3)
logical :: intersects
real(8) :: x0, y0, z0 ! track start point
real(8) :: x1, y1, z1 ! track end point
real(8) :: xi, yi, zi ! track intersection point with mesh
real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh
real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
y0 = xyz0(2)
z0 = xyz0(3)
! Copy coordinates of ending point
x1 = xyz1(1)
y1 = xyz1(2)
z1 = xyz1(3)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
ym0 = m % lower_left(2)
zm0 = m % lower_left(3)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
ym1 = m % upper_right(2)
zm1 = m % upper_right(3)
! Set default value for intersects
intersects = .false.
! Check if line intersects left surface -- calculate the intersection point
! (y,z)
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! (x,z)
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects bottom surface -- calculate the intersection
! point (x,y)
if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection point
! (y,z)
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! (x,z)
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects top surface -- calculate the intersection point
! (x,y)
if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_3d
end module mesh

359
src/mesh_header.F90
View file

@ -1,5 +1,7 @@
module mesh_header
use constants, only: NO_BIN_FOUND
implicit none
!===============================================================================
@ -16,6 +18,363 @@ module mesh_header
real(8), allocatable :: lower_left(:) ! lower-left corner of mesh
real(8), allocatable :: upper_right(:) ! upper-right corner of mesh
real(8), allocatable :: width(:) ! width of each mesh cell
contains
procedure :: get_bin => regular_get_bin
procedure :: get_indices => regular_get_indices
procedure :: get_bin_from_indices => regular_get_bin_from_indices
procedure :: get_indices_from_bin => regular_get_indices_from_bin
procedure :: intersects => regular_intersects
end type RegularMesh
contains
!===============================================================================
! GET_MESH_BIN determines the tally bin for a particle in a structured mesh
!===============================================================================
pure subroutine regular_get_bin(this, xyz, bin)
class(RegularMesh), intent(in) :: this
real(8), intent(in) :: xyz(:) ! coordinates
integer, intent(out) :: bin ! tally bin
integer :: n ! size of mesh
integer :: d ! mesh dimension index
integer :: ijk(3) ! indices in mesh
logical :: in_mesh ! was given coordinate in mesh at all?
! Get number of dimensions
n = this % n_dimension
! Loop over the dimensions of the mesh
do d = 1, n
! Check for cases where particle is outside of mesh
if (xyz(d) < this % lower_left(d)) then
bin = NO_BIN_FOUND
return
elseif (xyz(d) > this % upper_right(d)) then
bin = NO_BIN_FOUND
return
end if
end do
! Determine indices
call this % get_indices(xyz, ijk, in_mesh)
! Convert indices to bin
if (in_mesh) then
bin = this % get_bin_from_indices(ijk)
else
bin = NO_BIN_FOUND
end if
end subroutine regular_get_bin
!===============================================================================
! GET_MESH_INDICES determines the indices of a particle in a structured mesh
!===============================================================================
pure subroutine regular_get_indices(this, xyz, ijk, in_mesh)
class(RegularMesh), intent(in) :: this
real(8), intent(in) :: xyz(:) ! coordinates to check
integer, intent(out) :: ijk(:) ! indices in mesh
logical, intent(out) :: in_mesh ! were given coords in mesh?
! Find particle in mesh
ijk(:this % n_dimension) = ceiling((xyz(:this % n_dimension) - &
this % lower_left)/this % width)
! Determine if particle is in mesh
if (any(ijk(:this % n_dimension) < 1) .or. &
any(ijk(:this % n_dimension) > this % dimension)) then
in_mesh = .false.
else
in_mesh = .true.
end if
end subroutine regular_get_indices
!===============================================================================
! MESH_INDICES_TO_BIN maps (i), (i,j), or (i,j,k) indices to a single bin number
! for use in a TallyObject results array
!===============================================================================
pure function regular_get_bin_from_indices(this, ijk) result(bin)
class(RegularMesh), intent(in) :: this
integer, intent(in) :: ijk(:)
integer :: bin
if (this % n_dimension == 1) then
bin = ijk(1)
elseif (this % n_dimension == 2) then
bin = (ijk(2) - 1) * this % dimension(1) + ijk(1)
elseif (this % n_dimension == 3) then
bin = ((ijk(3) - 1) * this % dimension(2) + (ijk(2) - 1)) &
* this % dimension(1) + ijk(1)
end if
end function regular_get_bin_from_indices
!===============================================================================
! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to
! (i), (i,j), or (i,j,k) indices
!===============================================================================
pure subroutine regular_get_indices_from_bin(this, bin, ijk)
class(RegularMesh), intent(in) :: this
integer, intent(in) :: bin
integer, intent(out) :: ijk(:)
if (this % n_dimension == 1) then
ijk(1) = bin
else if (this % n_dimension == 2) then
ijk(1) = mod(bin - 1, this % dimension(1)) + 1
ijk(2) = (bin - 1)/this % dimension(1) + 1
else if (this % n_dimension == 3) then
ijk(1) = mod(bin - 1, this % dimension(1)) + 1
ijk(2) = mod(bin - 1, this % dimension(1) * this % dimension(2)) &
/ this % dimension(1) + 1
ijk(3) = (bin - 1)/(this % dimension(1) * this % dimension(2)) + 1
end if
end subroutine regular_get_indices_from_bin
!===============================================================================
! MESH_INTERSECTS determines if a line between xyz0 and xyz1 intersects the
! outer boundary of the given mesh. This is important for determining whether a
! track will score to a mesh tally.
!===============================================================================
pure function regular_intersects(this, xyz0, xyz1) result(intersects)
class(RegularMesh), intent(in) :: this
real(8), intent(in) :: xyz0(1)
real(8), intent(in) :: xyz1(1)
logical :: intersects
select case(this % n_dimension)
case (1)
intersects = mesh_intersects_1d(this, xyz0, xyz1)
case (2)
intersects = mesh_intersects_2d(this, xyz0, xyz1)
case (3)
intersects = mesh_intersects_3d(this, xyz0, xyz1)
end select
end function regular_intersects
pure function mesh_intersects_1d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(:)
real(8), intent(in) :: xyz1(:)
logical :: intersects
real(8) :: x0 ! track start point
real(8) :: x1 ! track end point
real(8) :: xm0 ! lower-left coordinates of mesh
real(8) :: xm1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
! Copy coordinates of ending point
x1 = xyz1(1)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
! Set default value for intersects
intersects = .false.
! Check if line intersects left surface
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
intersects = .true.
return
end if
! Check if line intersects right surface
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
intersects = .true.
return
end if
end function mesh_intersects_1d
pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(:)
real(8), intent(in) :: xyz1(:)
logical :: intersects
real(8) :: x0, y0 ! track start point
real(8) :: x1, y1 ! track end point
real(8) :: xi, yi ! track intersection point with mesh
real(8) :: xm0, ym0 ! lower-left coordinates of mesh
real(8) :: xm1, ym1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
y0 = xyz0(2)
! Copy coordinates of ending point
x1 = xyz1(1)
y1 = xyz1(2)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
ym0 = m % lower_left(2)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
ym1 = m % upper_right(2)
! Set default value for intersects
intersects = .false.
! Check if line intersects left surface -- calculate the intersection point
! y
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! x
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection
! point y
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! x
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_2d
pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(:)
real(8), intent(in) :: xyz1(:)
logical :: intersects
real(8) :: x0, y0, z0 ! track start point
real(8) :: x1, y1, z1 ! track end point
real(8) :: xi, yi, zi ! track intersection point with mesh
real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh
real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh
! Copy coordinates of starting point
x0 = xyz0(1)
y0 = xyz0(2)
z0 = xyz0(3)
! Copy coordinates of ending point
x1 = xyz1(1)
y1 = xyz1(2)
z1 = xyz1(3)
! Copy coordinates of mesh lower_left
xm0 = m % lower_left(1)
ym0 = m % lower_left(2)
zm0 = m % lower_left(3)
! Copy coordinates of mesh upper_right
xm1 = m % upper_right(1)
ym1 = m % upper_right(2)
zm1 = m % upper_right(3)
! Set default value for intersects
intersects = .false.
! Check if line intersects left surface -- calculate the intersection point
! (y,z)
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! (x,z)
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects bottom surface -- calculate the intersection
! point (x,y)
if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection point
! (y,z)
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! (x,z)
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects top surface -- calculate the intersection point
! (x,y)
if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_3d
end module mesh_header

3
src/output.F90
View file

@ -12,7 +12,6 @@ module output
use global
use math, only: t_percentile
use mesh_header, only: RegularMesh
use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices
use message_passing, only: master, n_procs
use nuclide_header
use particle_header, only: LocalCoord, Particle
@ -1055,7 +1054,7 @@ contains
do i = 1, n_cells
! Get the indices for this cell
call bin_to_mesh_indices(m, i, ijk)
call m % get_indices_from_bin(i, ijk)
matches(i_filter_mesh) % bins % data(1) = i
! Write the header for this cell

3
src/physics.F90
View file

@ -8,7 +8,6 @@ module physics
use global
use material_header, only: Material
use math
use mesh, only: get_mesh_indices
use message_passing
use nuclide_header
use output, only: write_message
@ -1090,7 +1089,7 @@ contains
if (ufs) then
! Determine indices on ufs mesh for current location
call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh)
call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")

3
src/physics_mg.F90
View file

@ -8,7 +8,6 @@ module physics_mg
use material_header, only: Material
use math, only: rotate_angle
use mgxs_header, only: Mgxs, MgxsContainer
use mesh, only: get_mesh_indices
use message_passing
use output, only: write_message
use particle_header, only: Particle
@ -196,7 +195,7 @@ contains
if (ufs) then
! Determine indices on ufs mesh for current location
call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh)
call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)

6
src/plot.F90
View file

@ -10,8 +10,6 @@ module plot
use geometry_header, only: Cell, root_universe
use global
use hdf5_interface
use mesh, only: get_mesh_indices
use mesh_header, only: RegularMesh
use output, only: write_message, time_stamp
use particle_header, only: LocalCoord, Particle
use plot_header
@ -241,8 +239,8 @@ contains
width = xyz_ur_plot - xyz_ll_plot
associate (m => pl % meshlines_mesh)
call get_mesh_indices(m, xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh)
call get_mesh_indices(m, xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh)
call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh)
call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh)
! sweep through all meshbins on this plane and draw borders
do i = ijk_ll(outer), ijk_ur(outer)

24
src/tally.F90
View file

@ -9,10 +9,6 @@ module tally
use geometry_header
use global
use math, only: t_percentile, calc_pn, calc_rn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, &
get_mesh_indices, mesh_indices_to_bin, &
mesh_intersects_1d, mesh_intersects_2d, &
mesh_intersects_3d
use mesh_header, only: RegularMesh
use message_passing
use output, only: header
@ -3259,19 +3255,13 @@ contains
n_dim = m % n_dimension
! Determine indices for starting and ending location
call get_mesh_indices(m, xyz0, ijk0, start_in_mesh)
call get_mesh_indices(m, xyz1, ijk1, end_in_mesh)
call m % get_indices(xyz0, ijk0, start_in_mesh)
call m % get_indices(xyz1, ijk1, end_in_mesh)
! Check to see if start or end is in mesh -- if not, check if track still
! intersects with mesh
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
if (n_dim == 1) then
if (.not. mesh_intersects_1d(m, xyz0, xyz1)) cycle
else if (n_dim == 2) then
if (.not. mesh_intersects_2d(m, xyz0, xyz1)) cycle
else
if (.not. mesh_intersects_3d(m, xyz0, xyz1)) cycle
end if
if (.not. m % intersects(xyz0, xyz1)) cycle
end if
! Calculate number of surface crossings
@ -3350,7 +3340,7 @@ contains
all(ijk0(:n_dim) <= m % dimension)) then
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 1
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
m % get_bin_from_indices(ijk0)
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &
@ -3389,7 +3379,7 @@ contains
ijk0(d1) = ijk0(d1) + 1
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 2
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
m % get_bin_from_indices(ijk0)
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &
@ -3412,7 +3402,7 @@ contains
all(ijk0(:n_dim) <= m % dimension)) then
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 3
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
m % get_bin_from_indices(ijk0)
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &
@ -3451,7 +3441,7 @@ contains
ijk0(d1) = ijk0(d1) - 1
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4
filter_matches(i_filter_mesh) % bins % data(1) = &
mesh_indices_to_bin(m, ijk0)
m % get_bin_from_indices(ijk0)
filter_index = 1
do k = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % &

28
src/tally_filter.F90
View file

@ -7,10 +7,6 @@ module tally_filter
use global
use hdf5_interface
use mesh_header, only: RegularMesh
use mesh, only: get_mesh_bin, bin_to_mesh_indices, &
get_mesh_indices, mesh_indices_to_bin, &
mesh_intersects_1d, mesh_intersects_2d, &
mesh_intersects_3d
use particle_header, only: Particle
use string, only: to_str
use tally_filter_header, only: TallyFilter, TallyFilterContainer, &
@ -266,7 +262,7 @@ contains
if (estimator /= ESTIMATOR_TRACKLENGTH) then
! If this is an analog or collision tally, then there can only be one
! valid mesh bin.
call get_mesh_bin(m, p % coord(1) % xyz, bin)
call m % get_bin(p % coord(1) % xyz, bin)
if (bin /= NO_BIN_FOUND) then
call match % bins % push_back(bin)
call match % weights % push_back(ONE)
@ -288,25 +284,13 @@ contains
xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
! Determine indices for starting and ending location.
call get_mesh_indices(m, xyz0, ijk0(:n), start_in_mesh)
call get_mesh_indices(m, xyz1, ijk1(:n), end_in_mesh)
call m % get_indices(xyz0, ijk0(:n), start_in_mesh)
call m % get_indices(xyz1, ijk1(:n), end_in_mesh)
! If this is the first iteration of the filter loop, check if the track
! intersects any part of the mesh.
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
if (n == 1) then
if (.not. mesh_intersects_1d(m, xyz0, xyz1)) then
return
end if
else if (n == 2) then
if (.not. mesh_intersects_2d(m, xyz0, xyz1)) then
return
end if
else
if (.not. mesh_intersects_3d(m, xyz0, xyz1)) then
return
end if
end if
if (.not. m % intersects(xyz0, xyz1)) return
end if
! ========================================================================
@ -393,7 +377,7 @@ contains
end if
! Assign the next tally bin and the score.
bin = mesh_indices_to_bin(m, ijk0(:n))
bin = m % get_bin_from_indices(ijk0(:n))
call match % bins % push_back(bin)
call match % weights % push_back(distance / total_distance)
@ -440,7 +424,7 @@ contains
m => meshes(this % mesh)
allocate(ijk(m % n_dimension))
call bin_to_mesh_indices(m, bin, ijk)
call m % get_indices_from_bin(bin, ijk)
if (m % n_dimension == 1) then
label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
elseif (m % n_dimension == 2) then

3
src/trigger.F90
View file

@ -11,7 +11,6 @@ module trigger
use global
use string, only: to_str
use output, only: warning, write_message
use mesh, only: bin_to_mesh_indices
use mesh_header, only: RegularMesh
use message_passing, only: master
use trigger_header, only: TriggerObject
@ -344,7 +343,7 @@ contains
do i = 1, n_cells
! Get the indices for this cell
call bin_to_mesh_indices(m, i, ijk)
call m % get_indices_from_bin(i, ijk)
filter_matches(i_filter_mesh) % bins % data(1) = i
do l = 1, n