From d50703937bec6d4e57f29f9761ceef20fd86f1e3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 29 Jan 2013 13:38:35 -0500 Subject: [PATCH] 3D lattice capability. --- src/constants.F90 | 10 ++-- src/geometry.F90 | 121 ++++++++++++++++++++++++++++++---------- src/geometry_header.F90 | 18 +++--- src/initialize.F90 | 43 +++++++++----- src/input_xml.F90 | 51 ++++++++++++----- src/output.F90 | 30 ++++++++-- src/particle_header.F90 | 1 + 7 files changed, 193 insertions(+), 81 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index dfe727471e..b7b02d86f3 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 :: & diff --git a/src/geometry.F90 b/src/geometry.F90 index a6c3333a4e..c85c4241e2 100644 --- a/src/geometry.F90 +++ b/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 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 58f7d02878..5dda5e4d85 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -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 !=============================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index ca159f47a7..1df088badd 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a951ade98a..b32a314000 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 ." 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 ." 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 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) diff --git a/src/output.F90 b/src/output.F90 index e708b97bc8..e95e7b8466 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/particle_header.F90 b/src/particle_header.F90 index f17030fd6a..965d637cfd 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -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)