mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
3D lattice capability.
This commit is contained in:
Paul Romano
parent
e34d8383ef
commit
d50703937b
7 changed files with 193 additions and 81 deletions
|
|
@ -93,10 +93,12 @@ module constants
|
|||
|
||||
! Lattice boundary crossings
|
||||
integer, parameter :: &
|
||||
LATTICE_LEFT = 1, & ! Flag for crossing left lattice boundary
|
||||
LATTICE_RIGHT = 2, & ! Flag for crossing right lattice boundary
|
||||
LATTICE_BOTTOM = 3, & ! Flag for crossing bottom lattice boundary
|
||||
LATTICE_TOP = 4 ! Flag for crossing top lattice boundary
|
||||
LATTICE_LEFT = 1, & ! Flag for crossing left (x) lattice boundary
|
||||
LATTICE_RIGHT = 2, & ! Flag for crossing right (x) lattice boundary
|
||||
LATTICE_BACK = 3, & ! Flag for crossing back (y) lattice boundary
|
||||
LATTICE_FRONT = 4, & ! Flag for crossing front (y) lattice boundary
|
||||
LATTICE_BOTTOM = 5, & ! Flag for crossing bottom (z) lattice boundary
|
||||
LATTICE_TOP = 6 ! Flag for crossing top (z) lattice boundary
|
||||
|
||||
! Surface types
|
||||
integer, parameter :: &
|
||||
|
|
|
|||
121
src/geometry.F90
121
src/geometry.F90
|
|
@ -73,8 +73,8 @@ contains
|
|||
integer, optional :: search_cells(:)
|
||||
|
||||
integer :: i ! index over cells
|
||||
integer :: x ! x-index for lattice
|
||||
integer :: y ! y-index for lattice
|
||||
integer :: i_x, i_y, i_z ! indices in lattice
|
||||
integer :: n_x, n_y, n_z ! size of lattice
|
||||
integer :: n ! number of cells to search
|
||||
integer :: index_cell ! index in cells array
|
||||
real(8) :: xyz(3) ! temporary location
|
||||
|
|
@ -163,14 +163,23 @@ contains
|
|||
! Set current lattice
|
||||
lat => lattices(c % fill)
|
||||
|
||||
! determine universe based on lattice position
|
||||
! determine lattice index based on position
|
||||
xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
|
||||
x = ceiling((xyz(1) - lat % x0)/lat % width_x)
|
||||
y = ceiling((xyz(2) - lat % y0)/lat % width_y)
|
||||
i_x = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
|
||||
i_y = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
|
||||
n_x = lat % dimension(1)
|
||||
n_y = lat % dimension(2)
|
||||
if (lat % n_dimension == 3) then
|
||||
i_z = ceiling((xyz(3) - lat % lower_left(3))/lat % width(3))
|
||||
n_z = lat % dimension(3)
|
||||
else
|
||||
i_z = 1
|
||||
n_z = 1
|
||||
end if
|
||||
|
||||
! Check if lattice coordinates are within bounds
|
||||
if (x < 1 .or. x > lat % n_x .or. &
|
||||
y < 1 .or. y > lat % n_y) then
|
||||
if (i_x < 1 .or. i_x > n_x .or. i_y < 1 .or. i_y > n_y .or. &
|
||||
i_z < 1 .or. i_z > n_z) then
|
||||
|
||||
! This condition should only get hit in rare circumstances where a
|
||||
! neutron hits the corner of a lattice. In this case, the neutron
|
||||
|
|
@ -192,20 +201,26 @@ contains
|
|||
|
||||
! adjust local position of particle
|
||||
p % coord % next % xyz(1) = p % coord % xyz(1) - &
|
||||
(lat%x0 + (x-0.5_8)*lat%width_x)
|
||||
(lat % lower_left(1) + (i_x - 0.5_8)*lat % width(1))
|
||||
p % coord % next % xyz(2) = p % coord % xyz(2) - &
|
||||
(lat%y0 + (y-0.5_8)*lat%width_y)
|
||||
p % coord % next % xyz(3) = p % coord % xyz(3)
|
||||
p % coord % next % uvw = p % coord % uvw
|
||||
(lat % lower_left(2) + (i_y - 0.5_8)*lat % width(2))
|
||||
if (lat % n_dimension == 3) then
|
||||
p % coord % next % xyz(3) = p % coord % xyz(3) - &
|
||||
(lat % lower_left(3) + (i_z - 0.5_8)*lat % width(3))
|
||||
else
|
||||
p % coord % next % xyz(3) = p % coord % xyz(3)
|
||||
end if
|
||||
p % coord % next % uvw = p % coord % uvw
|
||||
|
||||
! Move particle to next level
|
||||
p % coord => p % coord % next
|
||||
|
||||
! set particle lattice indices
|
||||
p % coord % lattice = c % fill
|
||||
p % coord % lattice_x = x
|
||||
p % coord % lattice_y = y
|
||||
p % coord % universe = lat % element(x,y)
|
||||
p % coord % lattice_x = i_x
|
||||
p % coord % lattice_y = i_y
|
||||
p % coord % lattice_z = i_z
|
||||
p % coord % universe = lat % element(i_x,i_y,i_z)
|
||||
end if
|
||||
|
||||
call find_cell(found)
|
||||
|
|
@ -517,12 +532,11 @@ contains
|
|||
|
||||
integer, intent(in) :: lattice_crossed
|
||||
|
||||
integer :: i_x ! x index in lattice
|
||||
integer :: i_y ! y index in lattice
|
||||
real(8) :: x0 ! half the width of lattice element
|
||||
real(8) :: y0 ! half the height of lattice element
|
||||
logical :: found ! particle found in cell?
|
||||
type(Lattice), pointer :: lat => null()
|
||||
integer :: i_x, i_y, i_z ! indices in lattice
|
||||
integer :: n_x, n_y, n_z ! size of lattice
|
||||
real(8) :: x0, y0, z0 ! half width of lattice element
|
||||
logical :: found ! particle found in cell?
|
||||
type(Lattice), pointer :: lat => null()
|
||||
|
||||
lat => lattices(p % coord % lattice)
|
||||
|
||||
|
|
@ -534,8 +548,9 @@ contains
|
|||
end if
|
||||
|
||||
if (lat % type == LATTICE_RECT) then
|
||||
x0 = lat % width_x * 0.5_8
|
||||
y0 = lat % width_y * 0.5_8
|
||||
x0 = lat % width(1) * 0.5_8
|
||||
y0 = lat % width(2) * 0.5_8
|
||||
if (lat % n_dimension == 3) z0 = lat % width(3) * 0.5_8
|
||||
|
||||
select case (lattice_crossed)
|
||||
case (LATTICE_LEFT)
|
||||
|
|
@ -548,16 +563,26 @@ contains
|
|||
p % coord % lattice_x = p % coord % lattice_x + 1
|
||||
p % coord % xyz(1) = -x0
|
||||
|
||||
case (LATTICE_BOTTOM)
|
||||
case (LATTICE_BACK)
|
||||
! Move particle to bottom element
|
||||
p % coord % lattice_y = p % coord % lattice_y - 1
|
||||
p % coord % xyz(2) = y0
|
||||
|
||||
case (LATTICE_TOP)
|
||||
case (LATTICE_FRONT)
|
||||
! Move particle to top element
|
||||
p % coord % lattice_y = p % coord % lattice_y + 1
|
||||
p % coord % xyz(2) = -y0
|
||||
|
||||
case (LATTICE_BOTTOM)
|
||||
! Move particle to bottom element
|
||||
p % coord % lattice_z = p % coord % lattice_z - 1
|
||||
p % coord % xyz(3) = z0
|
||||
|
||||
case (LATTICE_TOP)
|
||||
! Move particle to top element
|
||||
p % coord % lattice_z = p % coord % lattice_z + 1
|
||||
p % coord % xyz(3) = -z0
|
||||
|
||||
end select
|
||||
elseif (lat % type == LATTICE_HEX) then
|
||||
! TODO: Add hex lattice support
|
||||
|
|
@ -566,7 +591,16 @@ contains
|
|||
! Check to make sure still in lattice
|
||||
i_x = p % coord % lattice_x
|
||||
i_y = p % coord % lattice_y
|
||||
if (i_x < 1 .or. i_x > lat % n_x .or. i_y < 1 .or. i_y > lat % n_y) then
|
||||
i_z = p % coord % lattice_z
|
||||
n_x = lat % dimension(1)
|
||||
n_y = lat % dimension(2)
|
||||
if (lat % n_dimension == 3) then
|
||||
n_z = lat % dimension(3)
|
||||
else
|
||||
n_z = 1
|
||||
end if
|
||||
if (i_x < 1 .or. i_x > n_x .or. i_y < 1 .or. i_y > n_y .or. &
|
||||
i_z < 1 .or. i_z > n_z) then
|
||||
call deallocate_coord(p % coord0 % next)
|
||||
p % coord => p % coord0
|
||||
|
||||
|
|
@ -579,7 +613,7 @@ contains
|
|||
end if
|
||||
else
|
||||
! Find universe for next lattice element
|
||||
p % coord % universe = lat % element(i_x, i_y)
|
||||
p % coord % universe = lat % element(i_x, i_y, i_z)
|
||||
|
||||
! Find cell in next lattice element
|
||||
call find_cell(found)
|
||||
|
|
@ -1086,8 +1120,8 @@ contains
|
|||
z = coord % xyz(3)
|
||||
|
||||
! determine oncoming edge
|
||||
x0 = sign(lat % width_x * 0.5_8, u)
|
||||
y0 = sign(lat % width_y * 0.5_8, v)
|
||||
x0 = sign(lat % width(1) * 0.5_8, u)
|
||||
y0 = sign(lat % width(2) * 0.5_8, v)
|
||||
|
||||
! left and right sides
|
||||
if (abs(x - x0) < FP_PRECISION) then
|
||||
|
|
@ -1117,7 +1151,7 @@ contains
|
|||
end if
|
||||
end if
|
||||
|
||||
! top and bottom sides
|
||||
! front and back sides
|
||||
if (abs(y - y0) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
elseif (v == ZERO) then
|
||||
|
|
@ -1130,14 +1164,39 @@ contains
|
|||
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
|
||||
dist = d
|
||||
if (v > 0) then
|
||||
lattice_crossed = LATTICE_TOP
|
||||
lattice_crossed = LATTICE_FRONT
|
||||
else
|
||||
lattice_crossed = LATTICE_BOTTOM
|
||||
lattice_crossed = LATTICE_BACK
|
||||
end if
|
||||
final_coord => coord
|
||||
end if
|
||||
end if
|
||||
|
||||
if (lat % n_dimension == 3) then
|
||||
z0 = sign(lat % width(3) * 0.5_8, w)
|
||||
|
||||
! top and bottom sides
|
||||
if (abs(z - z0) < FP_PRECISION) then
|
||||
d = INFINITY
|
||||
elseif (w == ZERO) then
|
||||
d = INFINITY
|
||||
else
|
||||
d = (z0 - z)/w
|
||||
end if
|
||||
|
||||
if (d < dist) then
|
||||
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
|
||||
dist = d
|
||||
if (w > 0) then
|
||||
lattice_crossed = LATTICE_TOP
|
||||
else
|
||||
lattice_crossed = LATTICE_BOTTOM
|
||||
end if
|
||||
final_coord => coord
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
|
||||
elseif (lat % type == LATTICE_HEX) then
|
||||
! TODO: Add hex lattice support
|
||||
end if
|
||||
|
|
|
|||
|
|
@ -22,16 +22,14 @@ module geometry_header
|
|||
!===============================================================================
|
||||
|
||||
type Lattice
|
||||
integer :: id ! Universe number for lattice
|
||||
integer :: type ! Type of lattice (rectangular, hex, etc)
|
||||
integer :: level ! Level of lattice
|
||||
integer :: n_x ! number of lattice cells in x-direction
|
||||
integer :: n_y ! number of lattice cells in y-direction
|
||||
real(8) :: x0 ! x-coordinate of lower-left corner of lattice
|
||||
real(8) :: y0 ! y-coordinate of lower-left corner of lattice
|
||||
real(8) :: width_x ! width of lattice cell
|
||||
real(8) :: width_y ! width of lattice cell
|
||||
integer, allocatable :: element(:,:) ! specified universes
|
||||
integer :: id ! Universe number for lattice
|
||||
integer :: type ! Type of lattice (rectangular, hex, etc)
|
||||
integer :: level ! Level of lattice
|
||||
integer :: n_dimension ! Number of dimensions
|
||||
integer, allocatable :: dimension(:) ! number of cells in each direction
|
||||
real(8), allocatable :: lower_left(:) ! lower-left corner of lattice
|
||||
real(8), allocatable :: width(:) ! width of each lattice cell
|
||||
integer, allocatable :: element(:,:,:) ! specified universes
|
||||
end type Lattice
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -413,11 +413,13 @@ contains
|
|||
|
||||
subroutine adjust_indices()
|
||||
|
||||
integer :: i ! index for various purposes
|
||||
integer :: j ! index for various purposes
|
||||
integer :: k ! loop index for lattices
|
||||
integer :: i_array ! index in surfaces/materials array
|
||||
integer :: id ! user-specified id
|
||||
integer :: i ! index for various purposes
|
||||
integer :: j ! index for various purposes
|
||||
integer :: k ! loop index for lattices
|
||||
integer :: m ! loop index for lattices
|
||||
integer :: n_x, n_y, n_z ! size of lattice
|
||||
integer :: i_array ! index in surfaces/materials array
|
||||
integer :: id ! user-specified id
|
||||
type(Cell), pointer :: c => null()
|
||||
type(Lattice), pointer :: l => null()
|
||||
type(TallyObject), pointer :: t => null()
|
||||
|
|
@ -489,18 +491,29 @@ contains
|
|||
|
||||
do i = 1, n_lattices
|
||||
l => lattices(i)
|
||||
do j = 1, l % n_x
|
||||
do k = 1, l % n_y
|
||||
id = l % element(j,k)
|
||||
if (universe_dict % has_key(id)) then
|
||||
l % element(j,k) = universe_dict % get_key(id)
|
||||
else
|
||||
message = "Invalid universe number " // trim(to_str(id)) &
|
||||
// " specified on lattice " // trim(to_str(l % id))
|
||||
call fatal_error()
|
||||
end if
|
||||
n_x = l % dimension(1)
|
||||
n_y = l % dimension(2)
|
||||
if (l % n_dimension == 3) then
|
||||
n_z = l % dimension(3)
|
||||
else
|
||||
n_z = 1
|
||||
end if
|
||||
|
||||
do m = 1, n_z
|
||||
do k = 1, n_y
|
||||
do j = 1, n_x
|
||||
id = l % element(j,k,m)
|
||||
if (universe_dict % has_key(id)) then
|
||||
l % element(j,k,m) = universe_dict % get_key(id)
|
||||
else
|
||||
message = "Invalid universe number " // trim(to_str(id)) &
|
||||
// " specified on lattice " // trim(to_str(l % id))
|
||||
call fatal_error()
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
|
||||
TALLY_LOOP: do i = 1, n_tallies
|
||||
|
|
|
|||
|
|
@ -567,9 +567,9 @@ contains
|
|||
|
||||
use xml_data_geometry_t
|
||||
|
||||
integer :: i, j, k
|
||||
integer :: i, j, k, m
|
||||
integer :: n
|
||||
integer :: n_x, n_y
|
||||
integer :: n_x, n_y, n_z
|
||||
integer :: universe_num
|
||||
integer :: n_cells_in_univ
|
||||
integer :: coeffs_reqd
|
||||
|
|
@ -895,10 +895,10 @@ contains
|
|||
message = "Lattice must be two or three dimensions."
|
||||
call fatal_error()
|
||||
end if
|
||||
n_x = lattice_(i) % dimension(1)
|
||||
n_y = lattice_(i) % dimension(2)
|
||||
l % n_x = n_x
|
||||
l % n_y = n_y
|
||||
|
||||
l % n_dimension = n
|
||||
allocate(l % dimension(n))
|
||||
l % dimension = lattice_(i) % dimension
|
||||
|
||||
! Read lattice lower-left location
|
||||
if (size(lattice_(i) % dimension) /= size(lattice_(i) % lower_left)) then
|
||||
|
|
@ -906,8 +906,9 @@ contains
|
|||
"the number of entries on <dimension>."
|
||||
call fatal_error()
|
||||
end if
|
||||
l % x0 = lattice_(i) % lower_left(1)
|
||||
l % y0 = lattice_(i) % lower_left(2)
|
||||
|
||||
allocate(l % lower_left(n))
|
||||
l % lower_left = lattice_(i) % lower_left
|
||||
|
||||
! Read lattice widths
|
||||
if (size(lattice_(i) % width) /= size(lattice_(i) % lower_left)) then
|
||||
|
|
@ -915,17 +916,37 @@ contains
|
|||
"the number of entries on <lower_left>."
|
||||
call fatal_error()
|
||||
end if
|
||||
l % width_x = lattice_(i) % width(1)
|
||||
l % width_y = lattice_(i) % width(2)
|
||||
|
||||
allocate(l % width(n))
|
||||
l % width = lattice_(i) % width
|
||||
|
||||
! Copy number of dimensions
|
||||
n_x = l % dimension(1)
|
||||
n_y = l % dimension(2)
|
||||
if (l % n_dimension == 3) then
|
||||
n_z = l % dimension(3)
|
||||
else
|
||||
n_z = 1
|
||||
end if
|
||||
allocate(l % element(n_x, n_y, n_z))
|
||||
|
||||
! Check that number of universes matches size
|
||||
if (size(lattice_(i) % universes) /= n_x*n_y*n_z) then
|
||||
message = "Number of universes on <elements> does not match size of &
|
||||
&lattice " // trim(to_str(l % id)) // "."
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! Read universes
|
||||
allocate(l % element(n_x, n_y))
|
||||
do k = 0, n_y - 1
|
||||
do j = 1, n_x
|
||||
l % element(j, n_y - k) = lattice_(i) % universes(j + k*n_x)
|
||||
do m = 1, n_z
|
||||
do k = 0, n_y - 1
|
||||
do j = 1, n_x
|
||||
l % element(j, n_y - k, m) = lattice_(i) % &
|
||||
universes(j + n_x*k + n_x*n_y*(m-1))
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
! Add lattice to dictionary
|
||||
call lattice_dict % add_key(l % id, i)
|
||||
|
||||
|
|
|
|||
|
|
@ -450,7 +450,10 @@ contains
|
|||
type(Lattice), pointer :: lat
|
||||
integer, optional :: unit
|
||||
|
||||
integer :: i ! loop index
|
||||
integer :: unit_ ! unit to write to
|
||||
character(MAX_LINE_LEN) :: string
|
||||
|
||||
|
||||
! set default unit if not specified
|
||||
if (present(unit)) then
|
||||
|
|
@ -461,12 +464,27 @@ contains
|
|||
|
||||
! Write information about lattice
|
||||
write(unit_,*) 'Lattice ' // to_str(lat % id)
|
||||
write(unit_,*) ' n_x = ' // to_str(lat % n_x)
|
||||
write(unit_,*) ' n_y = ' // to_str(lat % n_y)
|
||||
write(unit_,*) ' x0 = ' // to_str(lat % x0)
|
||||
write(unit_,*) ' y0 = ' // to_str(lat % y0)
|
||||
write(unit_,*) ' width_x = ' // to_str(lat % width_x)
|
||||
write(unit_,*) ' width_y = ' // to_str(lat % width_y)
|
||||
|
||||
! Write dimension of lattice
|
||||
string = ""
|
||||
do i = 1, lat % n_dimension
|
||||
string = trim(string) // ' ' // to_str(lat % dimension(i))
|
||||
end do
|
||||
write(unit_,*) ' Dimension =' // string
|
||||
|
||||
! Write lower-left coordinates of lattice
|
||||
string = ""
|
||||
do i = 1, lat % n_dimension
|
||||
string = trim(string) // ' ' // to_str(lat % lower_left(i))
|
||||
end do
|
||||
write(unit_,*) ' Lower-left =' // string
|
||||
|
||||
! Write width of each lattice cell
|
||||
string = ""
|
||||
do i = 1, lat % n_dimension
|
||||
string = trim(string) // ' ' // to_str(lat % width(i))
|
||||
end do
|
||||
write(unit_,*) ' Width =' // string
|
||||
write(unit_,*)
|
||||
|
||||
end subroutine print_lattice
|
||||
|
|
|
|||
|
|
@ -17,6 +17,7 @@ module particle_header
|
|||
integer :: lattice = NONE
|
||||
integer :: lattice_x = NONE
|
||||
integer :: lattice_y = NONE
|
||||
integer :: lattice_z = NONE
|
||||
|
||||
! Particle position and direction for this level
|
||||
real(8) :: xyz(3)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue