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

Started transition of lattices to polymorphism

Browse source
This commit is contained in:
Sterling Harper 2014-08-13 16:46:34 +02:00
parent 34bb23d1cc
commit ca4b00e053
7 changed files with 305 additions and 241 deletions
Download patch file Download diff file Expand all files Collapse all files

68
src/geometry.F90
View file

@ -2,7 +2,8 @@ module geometry
use constants
use error, only: fatal_error, warning
use geometry_header, only: Cell, Surface, Universe, Lattice
use geometry_header, only: Cell, Surface, Universe, Lattice, &
&RectLattice, HexLattice
use global
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord, Particle
@ -138,7 +139,7 @@ contains
real(8) :: xyz(3) ! temporary location
logical :: use_search_cells ! use cells provided as argument
type(Cell), pointer, save :: c => null() ! pointer to cell
type(Lattice), pointer, save :: lat => null() ! pointer to lattice
class(Lattice), pointer, save :: lat => null() ! pointer to lattice
type(Universe), pointer, save :: univ => null() ! universe to search in
!$omp threadprivate(c, lat, univ)
@ -222,7 +223,7 @@ contains
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
! Set current lattice
lat => lattices(c % fill)
lat => lattices(c % fill) % obj
! Determine lattice indices
i_xyz = get_lat_indices(lat, &
@ -279,15 +280,17 @@ contains
!===============================================================================
function get_lat_trans(lat, xyz, i_xyz) result(xyz_t)
type(Lattice), pointer, intent(in) :: lat
real(8) , intent(in) :: xyz(3)
integer , intent(in) :: i_xyz(3)
class(Lattice), pointer, intent(in) :: lat
real(8) , intent(in) :: xyz(3)
integer , intent(in) :: i_xyz(3)
real(8) :: xyz_t(3)
integer :: n_rings
if (lat % type == LATTICE_RECT) then
select type(lat)
type is (RectLattice)
xyz_t(1) = xyz(1) - (lat % lower_left(1) + &
&(i_xyz(1) - 0.5_8)*lat % width(1))
xyz_t(2) = xyz(2) - (lat % lower_left(2) + &
@ -299,7 +302,7 @@ contains
xyz_t(3) = xyz(3)
end if
else if (lat % type == LATTICE_HEX) then
type is (HexLattice)
n_rings = lat % dimension(1)
xyz_t(1) = xyz(1) - (lat % lower_left(1) + &
@ -314,7 +317,7 @@ contains
xyz_t(3) = xyz(3)
end if
end if
end select
end function get_lat_trans
@ -324,13 +327,15 @@ contains
!===============================================================================
function is_valid_lat_index(lat, i_xyz) result(is_valid)
type(Lattice), pointer, intent(in) :: lat
integer , intent(in) :: i_xyz(3)
logical :: is_valid
class(Lattice), pointer, intent(in) :: lat
integer , intent(in) :: i_xyz(3)
logical :: is_valid
integer :: n_rings, n_x, n_y, n_z
if (lat % type == LATTICE_RECT) then
select type(lat)
type is (RectLattice)
n_x = lat % dimension(1)
n_y = lat % dimension(2)
if (lat % n_dimension == 3) then
@ -343,7 +348,7 @@ contains
&i_xyz(2) > 0 .and. i_xyz(2) <= n_y .and. &
&i_xyz(3) > 0 .and. i_xyz(3) <= n_z
else if (lat % type == LATTICE_HEX) then
type is (HexLattice)
n_rings = lat % dimension(1)
if (lat % n_dimension == 2) then
@ -358,7 +363,7 @@ contains
&i_xyz(1) + i_xyz(2) < 3*n_rings .and. &
&i_xyz(3) > 0 .and. i_xyz(3) < n_z + 1)
end if
end select
end function is_valid_lat_index
@ -368,9 +373,9 @@ contains
!===============================================================================
function get_lat_indices(lat, xyz, uvw) result(i_xyz)
type(Lattice), pointer, intent(in) :: lat
real(8) , intent(in) :: xyz(3)
real(8) , intent(in) :: uvw(3)
class(Lattice), pointer, intent(in) :: lat
real(8) , intent(in) :: xyz(3)
real(8) , intent(in) :: uvw(3)
integer :: i_xyz(3)
@ -379,8 +384,10 @@ contains
real(8) :: dists(4)
integer :: n_rings
integer :: i, j, k, loc(1)
if (lat % type == LATTICE_RECT) then
select type(lat)
type is (RectLattice)
! Find approximate indices using ceiling division.
i_xyz(1) = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
i_xyz(2) = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
@ -416,7 +423,7 @@ contains
! end if
!end if
else if (lat % type == LATTICE_HEX) then
type is (HexLattice)
! Index z direction.
if (lat % n_dimension == 2) then
i_xyz(3) = ceiling((xyz(3) - lat % lower_left(3)) / lat % width(2))
@ -533,7 +540,7 @@ contains
end if
end if
end select
end function get_lat_indices
@ -849,11 +856,11 @@ contains
real(8) :: x0, y0, z0 ! half width of lattice element
real(8) :: half_pitch ! half pitch of a hexagonal lattice
logical :: found ! particle found in cell?
type(Lattice), pointer, save :: lat => null()
class(Lattice), pointer, save :: lat => null()
type(LocalCoord), pointer :: parent_coord
!$omp threadprivate(lat)
lat => lattices(p % coord % lattice)
lat => lattices(p % coord % lattice) % obj
if (verbosity >= 10 .or. trace) then
message = " Crossing lattice " // trim(to_str(lat % id)) // &
@ -952,7 +959,7 @@ contains
logical :: on_surface ! is particle on surface?
type(Cell), pointer, save :: cl => null()
type(Surface), pointer, save :: surf => null()
type(Lattice), pointer, save :: lat => null()
class(Lattice), pointer, save :: lat => null()
type(LocalCoord), pointer, save :: coord => null()
type(LocalCoord), pointer, save :: final_coord => null()
@ -1413,8 +1420,11 @@ contains
! FIND MINIMUM DISTANCE TO LATTICE SURFACES
LAT_COORD: if (coord % lattice /= NONE) then
lat => lattices(coord % lattice)
LAT_TYPE: if (lat % type == LATTICE_RECT) then
lat => lattices(coord % lattice) % obj
LAT_TYPE: select type(lat)
type is (RectLattice)
! copy local coordinates
x = coord % xyz(1)
y = coord % xyz(2)
@ -1480,7 +1490,7 @@ contains
end if
end if
elseif (lat % type == LATTICE_HEX) then LAT_TYPE
type is (HexLattice) LAT_TYPE
! Copy local coordinates.
x = coord % xyz(1)
y = coord % xyz(2)
@ -1592,7 +1602,7 @@ contains
end if
end if
end if
end if LAT_TYPE
end select LAT_TYPE
! If the lattice boundary is coincident with the parent cell boundary,
! we need to make sure that the lattice is not selected. This is

30
src/geometry_header.F90
View file

@ -21,18 +21,28 @@ module geometry_header
! triangular)
!===============================================================================
type Lattice
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 :: universes(:,:,:) ! specified universes
integer :: outside ! material to fill area outside
type, abstract :: Lattice
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 :: universes(:,:,:) ! specified universes
integer :: outside ! material to fill area outside
end type Lattice
type, extends(Lattice) :: RectLattice
end type RectLattice
type, extends(Lattice) :: HexLattice
end type HexLattice
type LatticeContainer
class(Lattice), allocatable :: obj
end type LatticeContainer
!===============================================================================
! SURFACE type defines a first- or second-order surface that can be used to
! construct closed volumes (cells)

14
src/global.F90
View file

@ -6,7 +6,7 @@ module global
use cmfd_header
use constants
use dict_header, only: DictCharInt, DictIntInt
use geometry_header, only: Cell, Universe, Lattice, Surface
use geometry_header, only: Cell, Universe, Lattice, LatticeContainer, Surface
use material_header, only: Material
use mesh_header, only: StructuredMesh
use plot_header, only: ObjectPlot
@ -26,12 +26,12 @@ module global
! GEOMETRY-RELATED VARIABLES
! Main arrays
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(ObjectPlot),allocatable, target :: plots(:)
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(LatticeContainer), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(ObjectPlot), allocatable, target :: plots(:)
! Size of main arrays
integer :: n_cells ! # of cells

8
src/hdf5_summary.F90
View file

@ -105,7 +105,7 @@ contains
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
type(Lattice), pointer :: lat => null()
class(Lattice), pointer :: lat => null()
! Use H5LT interface to write number of geometry objects
call su % write_data(n_cells, "n_cells", group="geometry")
@ -275,11 +275,11 @@ contains
lat => lattices(i)
! Write lattice type
select case(lat % type)
case (LATTICE_RECT)
select type (lat)
type is (LatticeRect)
call su % write_data("rectangular", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
case (LATTICE_HEX)
type is (LatticeHex)
call su % write_data("hexagonal", "type", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end select

16
src/initialize.F90
View file

@ -7,7 +7,8 @@ module initialize
use energy_grid, only: unionized_grid
use error, only: fatal_error, warning
use geometry, only: neighbor_lists
use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
use global
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict, read_plots_xml
@ -566,7 +567,7 @@ contains
integer :: i_array ! index in surfaces/materials array
integer :: id ! user-specified id
type(Cell), pointer :: c => null()
type(Lattice), pointer :: lat => null()
class(Lattice), pointer :: lat => null()
type(TallyObject), pointer :: t => null()
do i = 1, n_cells
@ -622,7 +623,7 @@ contains
c % fill = universe_dict % get_key(id)
elseif (lattice_dict % has_key(id)) then
lid = lattice_dict % get_key(id)
mid = lattices(lid) % outside
mid = lattices(lid) % obj % outside
c % type = CELL_LATTICE
c % fill = lid
if (mid == MATERIAL_VOID) then
@ -647,9 +648,10 @@ contains
! ADJUST UNIVERSE INDICES FOR EACH LATTICE
do i = 1, n_lattices
lat => lattices(i)
lat => lattices(i) % obj
select type (lat)
if (lat % type == LATTICE_RECT) then
type is (RectLattice)
n_x = lat % dimension(1)
n_y = lat % dimension(2)
if (lat % n_dimension == 3) then
@ -673,7 +675,7 @@ contains
end do
end do
else
type is (HexLattice)
n_rings = lat % dimension(1)
if (lat % n_dimension == 2) then
n_z = lat % dimension(2)
@ -701,7 +703,7 @@ contains
end do
end do
end if
end select
end do

394
src/input_xml.F90
View file

@ -4,7 +4,7 @@ module input_xml
use constants
use dict_header, only: DictIntInt, ElemKeyValueCI
use error, only: fatal_error, warning
use geometry_header, only: Cell, Surface, Lattice
use geometry_header, only: Cell, Surface, Lattice, RectLattice, HexLattice
use global
use list_header, only: ListChar, ListReal
use mesh_header, only: StructuredMesh
@ -827,7 +827,7 @@ contains
integer :: i, j, k, m, i_x, i_a, input_index
integer :: n
integer :: n_x, n_y, n_z, n_rings
integer :: n_x, n_y, n_z, n_rings, n_rlats, n_hlats
integer :: universe_num
integer :: n_cells_in_univ
integer :: coeffs_reqd
@ -839,16 +839,17 @@ contains
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Lattice), pointer :: lat => null()
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
class(Lattice), pointer :: lat => null()
type(Node), pointer :: doc => null()
type(Node), pointer :: node_cell => null()
type(Node), pointer :: node_surf => null()
type(Node), pointer :: node_lat => null()
type(NodeList), pointer :: node_cell_list => null()
type(NodeList), pointer :: node_surf_list => null()
type(NodeList), pointer :: node_lat_list => null()
type(NodeList), pointer :: node_rlat_list => null()
type(NodeList), pointer :: node_hlat_list => null()
! Display output message
message = "Reading geometry XML file..."
@ -1187,17 +1188,21 @@ contains
! READ LATTICES FROM GEOMETRY.XML
! Get pointer to list of XML <lattice>
call get_node_list(doc, "lattice", node_lat_list)
call get_node_list(doc, "lattice", node_rlat_list)
call get_node_list(doc, "hex_lattice", node_hlat_list)
! Allocate lattices array
n_lattices = get_list_size(node_lat_list)
n_rlats = get_list_size(node_rlat_list)
n_hlats = get_list_size(node_hlat_list)
n_lattices = n_rlats + n_hlats
allocate(lattices(n_lattices))
do i = 1, n_lattices
lat => lattices(i)
RECT_LATTICES: do i = 1, n_rlats
allocate(RectLattice::lattices(i) % obj)
lat => lattices(i) % obj
! Get pointer to i-th lattice
call get_list_item(node_lat_list, i, node_lat)
call get_list_item(node_rlat_list, i, node_lat)
! ID of lattice
if (check_for_node(node_lat, "id")) then
@ -1214,28 +1219,11 @@ contains
call fatal_error()
end if
! Read lattice type
word = ''
if (check_for_node(node_lat, "type")) &
call get_node_value(node_lat, "type", word)
call lower_case(word)
select case (trim(word))
case ('rect', 'rectangle', 'rectangular')
lat % type = LATTICE_RECT
case ('hex', 'hexagon', 'hexagonal')
lat % type = LATTICE_HEX
case default
message = "Invalid lattice type: " // trim(word)
call fatal_error()
end select
! Read number of lattice cells in each dimension
n = get_arraysize_integer(node_lat, "dimension")
if (lat % type == LATTICE_RECT .and. n /= 2 .and. n /= 3) then
if (n /= 2 .and. n /= 3) then
message = "Rectangular lattice must be two or three dimensions."
call fatal_error()
else if (lat % type == LATTICE_HEX .and. n /= 1 .and. n/= 2) then
message = "Hexagonal lattice must be one or two dimension(s)."
end if
lat % n_dimension = n
@ -1243,28 +1231,16 @@ contains
call get_node_array(node_lat, "dimension", lat % dimension)
! Read lattice lower-left location
if (lat % type == LATTICE_RECT) then
if (size(lat % dimension) /= &
&get_arraysize_double(node_lat, "lower_left")) then
message = "Number of entries on <lower_left> must be the same as &
&the number of entries on <dimension>."
call fatal_error()
end if
allocate(lat % lower_left(n))
call get_node_array(node_lat, "lower_left", lat % lower_left)
else
if (size(lat % dimension) /= &
&get_arraysize_double(node_lat, "lower_left") - 1) then
message = "Number of entries on <lower_left> must be one greater &
&than the number of entries on <dimension>."
call fatal_error()
end if
allocate(lat % lower_left(n+1))
call get_node_array(node_lat, "lower_left", lat % lower_left)
if (size(lat % dimension) /= &
&get_arraysize_double(node_lat, "lower_left")) then
message = "Number of entries on <lower_left> must be the same as &
&the number of entries on <dimension>."
call fatal_error()
end if
allocate(lat % lower_left(n))
call get_node_array(node_lat, "lower_left", lat % lower_left)
! Read lattice widths
if (size(lat % dimension) /= &
get_arraysize_double(node_lat, "width")) then
@ -1277,145 +1253,35 @@ contains
call get_node_array(node_lat, "width", lat % width)
! Copy number of dimensions
if (lat % type == LATTICE_RECT) then
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
allocate(lat % universes(n_x, n_y, n_z))
n_x = lat % dimension(1)
n_y = lat % dimension(2)
if (lat % n_dimension == 3) then
n_z = lat % dimension(3)
else
n_rings = lat % dimension(1)
if (lat % n_dimension == 2) then
n_z = lat % dimension(2)
else
n_z = 1
end if
allocate(lat % universes(2*n_rings - 1, 2*n_rings - 1, n_z))
n_z = 1
end if
allocate(lat % universes(n_x, n_y, n_z))
! Check that number of universes matches size
n = get_arraysize_integer(node_lat, "universes")
if (lat % type == LATTICE_RECT) then
if (n /= n_x*n_y*n_z) then
message = "Number of universes on <universes> does not match size of &
&lattice " // trim(to_str(lat % id)) // "."
call fatal_error()
end if
else
if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then
message = "Number of universes on <universes> does not match size of &
&lattice " // trim(to_str(lat % id)) // "."
call fatal_error()
end if
if (n /= n_x*n_y*n_z) then
message = "Number of universes on <universes> does not match size of &
&lattice " // trim(to_str(lat % id)) // "."
call fatal_error()
end if
allocate(temp_int_array(n))
call get_node_array(node_lat, "universes", temp_int_array)
! Read universes
if (lat % type == LATTICE_RECT) then
do m = 1, n_z
do k = 0, n_y - 1
do j = 1, n_x
lat % universes(j, n_y - k, m) = &
temp_int_array(j + n_x*k + n_x*n_y*(m-1))
end do
do m = 1, n_z
do k = 0, n_y - 1
do j = 1, n_x
lat % universes(j, n_y - k, m) = &
&temp_int_array(j + n_x*k + n_x*n_y*(m-1))
end do
end do
else
do m = 1, n_z
do k = 1, 2*n_rings - 1
do j = 1, 2*n_rings - 1
lat % universes(j, k, m) = -1
end do
end do
end do
! TODO: The index walk in the k loops can be made a little more
! efficient. They currently sometimes change index values and then
! change them again before use.
! Universes in hexagonal lattices are stored in a manner that represents
! a skewed coordinate system: (x, alpha) rather than (x, y). There is
! no obvious, direct relationship between the order of universes in the
! input and the order that they will be stored in the skewed array so
! the following code walks a set of index values across the skewed array
! in a manner that matches the input order. Note that i_x = 0, i_a = 0
! corresponds to the center of the hexagonal lattice.
input_index = 1
do m = 1, n_z
! Initialize lattice indecies.
i_x = 1
i_a = n_rings - 1
! Map upper triangular region of hexagonal lattice.
do k = 1, n_rings-1
! Walk index to lower-left neighbor of last row start.
i_x = i_x - 1
do j = 1, k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*k
i_a = i_a + k
end do
! Map middle square region of hexagonal lattice.
do k = 1, 2*n_rings - 1
if (mod(k, 2) == 1) then
! Walk index to lower-left neighbor of last row start.
i_x = i_x - 1
else
! Walk index to lower-right neighbor of last row start
i_x = i_x + 1
i_a = i_a - 1
end if
do j = 1, n_rings - mod(k-1, 2)
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*(n_rings - mod(k-1, 2))
i_a = i_a + n_rings - mod(k-1, 2)
end do
! Map lower triangular region of hexagonal lattice.
do k = 1, n_rings-1
! Walk index to lower-right neighbor of last row start.
i_x = i_x + 1
i_a = i_a - 1
do j = 1, n_rings - k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*(n_rings - k)
i_a = i_a + n_rings - k
end do
end do
end if
end do
deallocate(temp_int_array)
! Read material for area outside lattice
@ -1432,7 +1298,183 @@ contains
! Add lattice to dictionary
call lattice_dict % add_key(lat % id, i)
end do
end do RECT_LATTICES
HEX_LATTICES: do i = 1, n_hlats
allocate(HexLattice::lattices(i) % obj)
lat => lattices(i) % obj
! Get pointer to i-th lattice
call get_list_item(node_hlat_list, i, node_lat)
! ID of lattice
if (check_for_node(node_lat, "id")) then
call get_node_value(node_lat, "id", lat % id)
else
message = "Must specify id of lattice in geometry XML file."
call fatal_error()
end if
! Check to make sure 'id' hasn't been used
if (lattice_dict % has_key(lat % id)) then
message = "Two or more lattices use the same unique ID: " // &
to_str(lat % id)
call fatal_error()
end if
! Read number of lattice cells in each dimension
n = get_arraysize_integer(node_lat, "dimension")
if (n /= 1 .and. n/= 2) then
message = "Hexagonal lattice must be one or two dimension(s)."
call fatal_error()
end if
lat % n_dimension = n
allocate(lat % dimension(n))
call get_node_array(node_lat, "dimension", lat % dimension)
! Read lattice lower-left location
if (size(lat % dimension) /= &
&get_arraysize_double(node_lat, "lower_left") - 1) then
message = "Number of entries on <lower_left> must be one greater &
&than the number of entries on <dimension>."
call fatal_error()
end if
allocate(lat % lower_left(n+1))
call get_node_array(node_lat, "lower_left", lat % lower_left)
! Read lattice widths
if (size(lat % dimension) /= &
get_arraysize_double(node_lat, "width")) then
message = "Number of entries on <width> must be the same as &
&the number of entries on <dimension>."
call fatal_error()
end if
allocate(lat % width(n))
call get_node_array(node_lat, "width", lat % width)
! Copy number of dimensions
n_rings = lat % dimension(1)
if (lat % n_dimension == 2) then
n_z = lat % dimension(2)
else
n_z = 1
end if
allocate(lat % universes(2*n_rings - 1, 2*n_rings - 1, n_z))
! Check that number of universes matches size
n = get_arraysize_integer(node_lat, "universes")
if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then
message = "Number of universes on <universes> does not match size of &
&lattice " // trim(to_str(lat % id)) // "."
call fatal_error()
end if
allocate(temp_int_array(n))
call get_node_array(node_lat, "universes", temp_int_array)
! Read universes
! TODO: The index walk in the k loops can be made a little more
! efficient. They currently sometimes change index values and then
! change them again before use.
! Universes in hexagonal lattices are stored in a manner that represents
! a skewed coordinate system: (x, alpha) rather than (x, y). There is
! no obvious, direct relationship between the order of universes in the
! input and the order that they will be stored in the skewed array so
! the following code walks a set of index values across the skewed array
! in a manner that matches the input order. Note that i_x = 0, i_a = 0
! corresponds to the center of the hexagonal lattice.
input_index = 1
do m = 1, n_z
! Initialize lattice indecies.
i_x = 1
i_a = n_rings - 1
! Map upper triangular region of hexagonal lattice.
do k = 1, n_rings-1
! Walk index to lower-left neighbor of last row start.
i_x = i_x - 1
do j = 1, k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*k
i_a = i_a + k
end do
! Map middle square region of hexagonal lattice.
do k = 1, 2*n_rings - 1
if (mod(k, 2) == 1) then
! Walk index to lower-left neighbor of last row start.
i_x = i_x - 1
else
! Walk index to lower-right neighbor of last row start
i_x = i_x + 1
i_a = i_a - 1
end if
do j = 1, n_rings - mod(k-1, 2)
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*(n_rings - mod(k-1, 2))
i_a = i_a + n_rings - mod(k-1, 2)
end do
! Map lower triangular region of hexagonal lattice.
do k = 1, n_rings-1
! Walk index to lower-right neighbor of last row start.
i_x = i_x + 1
i_a = i_a - 1
do j = 1, n_rings - k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
&temp_int_array(input_index)
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*(n_rings - k)
i_a = i_a + n_rings - k
end do
end do
deallocate(temp_int_array)
! Read material for area outside lattice
lat % outside = MATERIAL_VOID
if (check_for_node(node_lat, "outside")) then
call get_node_value(node_lat, "outside", mid)
if (mid == 0 .or. mid == MATERIAL_VOID) then
lat % outside = MATERIAL_VOID
else
lat % outside = mid
end if
end if
! Add lattice to dictionary
call lattice_dict % add_key(lat % id, i)
end do HEX_LATTICES
! Close geometry XML file
call close_xmldoc(doc)

16
src/output.F90
View file

@ -256,7 +256,7 @@ contains
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
class(Lattice), pointer :: l => null()
type(LocalCoord), pointer :: coord => null()
! display type of particle
@ -292,7 +292,7 @@ contains
! Print information on lattice
if (coord % lattice /= NONE) then
l => lattices(coord % lattice)
l => lattices(coord % lattice) % obj
write(ou,*) ' Lattice = ' // trim(to_str(l % id))
write(ou,*) ' Lattice position = (' // trim(to_str(&
p % coord % lattice_x)) // ',' // trim(to_str(&
@ -355,7 +355,7 @@ contains
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
class(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
! Set unit to stdout if not already set
@ -384,7 +384,7 @@ contains
u => universes(c % fill)
write(unit_,*) ' Fill = Universe ' // to_str(u % id)
case (CELL_LATTICE)
l => lattices(c % fill)
l => lattices(c % fill) % obj
write(unit_,*) ' Fill = Lattice ' // to_str(l % id)
end select
@ -471,8 +471,8 @@ contains
subroutine print_lattice(lat, unit)
type(Lattice), pointer :: lat
integer, optional :: unit
class(Lattice), pointer :: lat
integer, optional :: unit
integer :: i ! loop index
integer :: unit_ ! unit to write to
@ -919,7 +919,7 @@ contains
type(Surface), pointer :: s => null()
type(Cell), pointer :: c => null()
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
class(Lattice), pointer :: l => null()
! print summary of surfaces
call header("SURFACE SUMMARY", unit=UNIT_SUMMARY)
@ -946,7 +946,7 @@ contains
if (n_lattices > 0) then
call header("LATTICE SUMMARY", unit=UNIT_SUMMARY)
do i = 1, n_lattices
l => lattices(i)
l => lattices(i) % obj
call print_lattice(l, unit=UNIT_SUMMARY)
end do
end if