From 0f7cd378dffc0551ace9ba82a1a88f28600784cb Mon Sep 17 00:00:00 2001 From: nhorelik Date: Tue, 9 Apr 2013 13:43:59 -0400 Subject: [PATCH 1/5] Added ability to treat area outside lattices as filled with a material --- src/geometry.F90 | 53 ++++++++++++++++++++++++++++-------- src/geometry_header.F90 | 1 + src/initialize.F90 | 16 ++++++++++- src/input_xml.F90 | 13 +++++++-- src/templates/geometry_t.xml | 1 + 5 files changed, 69 insertions(+), 15 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 64912cb954..9dd771d5d8 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -79,6 +79,8 @@ contains integer :: index_cell ! index in cells array real(8) :: xyz(3) ! temporary location logical :: use_search_cells ! use cells provided as argument + logical :: outside_lattice ! if particle is not inside lattice bounds + logical :: lattice_edge ! if the particle is on a lattice edge type(Cell), pointer :: c ! pointer to cell type(Lattice), pointer :: lat ! pointer to lattice type(Universe), pointer :: univ ! universe to search in @@ -162,6 +164,8 @@ contains ! Set current lattice lat => lattices(c % fill) + + outside_lattice = .false. ! determine lattice index based on position xyz = p % coord % xyz + TINY_BIT * p % coord % uvw @@ -181,19 +185,42 @@ contains 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 - ! may need to be moved diagonally across the lattice. To do so, we - ! remove all lower coordinate levels and then search from universe - ! 0. + ! Check for when particle is on lattice edge + lattice_edge = .false. + if ( abs(xyz(1) - lat % lower_left(1)) < FP_COINCIDENT .or. & + abs(xyz(2) - lat % lower_left(2)) < FP_COINCIDENT) then + lattice_edge = .true. + end if + if (lat % n_dimension == 3) then + if (abs(xyz(3) - lat % lower_left(3)) < FP_COINCIDENT) then + lattice_edge = .true. + end if + end if + + if (lattice_edge) then + + ! In this case the neutron is leaving the lattice, so we move it + ! out, remove all lower coordinate levels and then search from + ! universe 0. + + p % coord => p % coord0 + call deallocate_coord(p % coord % next) - p % coord => p % coord0 - call deallocate_coord(p % coord % next) + ! Reset surface and advance particle a tiny bit + p % surface = NONE + p % coord % xyz = xyz - ! Reset surface and advance particle a tiny bit - p % surface = NONE - p % coord % xyz = xyz + else + + ! We're outside the lattice, so treat this as a normal cell with + ! the material specified for the outside + outside_lattice = .true. + p % last_material = p % material + p % material = c % material + + end if + else ! Create new level of coordinates @@ -223,8 +250,10 @@ contains p % coord % universe = lat % universes(i_x,i_y,i_z) end if - call find_cell(found) - if (.not. found) exit + if (.not. outside_lattice) then + call find_cell(found) + if (.not. found) exit + end if end if ! Found cell so we can return diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index abbad38628..93863865c4 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -30,6 +30,7 @@ module geometry_header 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 !=============================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index 83946165d3..68eecf8a87 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -398,6 +398,7 @@ contains ! increment the index for the cells array within the Universe object and ! then store the index of the Cell object in that array index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1 + univ % cells(index_cell_in_univ(i_univ)) = i end do @@ -417,6 +418,7 @@ contains integer :: j ! index for various purposes integer :: k ! loop index for lattices integer :: m ! loop index for lattices + integer :: mid, lid ! material and lattice ids integer :: n_x, n_y, n_z ! size of lattice integer :: i_array ! index in surfaces/materials array integer :: id ! user-specified id @@ -476,8 +478,20 @@ contains c % type = CELL_FILL c % fill = universe_dict % get_key(id) elseif (lattice_dict % has_key(id)) then + lid = lattice_dict % get_key(id) + mid = lattices(lid) % outside c % type = CELL_LATTICE - c % fill = lattice_dict % get_key(id) + c % fill = lid + if (mid == MATERIAL_VOID) then + c % material = mid + else if (material_dict % has_key(mid)) then + c % material = material_dict % get_key(mid) + else + message = "Could not find material " // trim(to_str(mid)) // & + " specified on lattice " // trim(to_str(lid)) + call fatal_error() + end if + else message = "Specified fill " // trim(to_str(id)) // " on cell " // & trim(to_str(c % id)) // " is neither a universe nor a lattice." diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9d6fe36d23..740b404566 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -578,6 +578,7 @@ contains integer :: universe_num integer :: n_cells_in_univ integer :: coeffs_reqd + integer :: mid real(8) :: phi, theta, psi logical :: file_exists logical :: boundary_exists @@ -607,7 +608,7 @@ contains ! Get number of tags n_cells = size(cell_) - + ! Check for no cells if (n_cells == 0) then message = "No cells found in geometry.xml!" @@ -968,7 +969,15 @@ contains end do end do end do - + + ! Read material for area outside lattice + mid = lattice_(i) % outside + if (mid == 0 .or. mid == MATERIAL_VOID) then + lat % outside = MATERIAL_VOID + else + lat % outside = mid + end if + ! Add lattice to dictionary call lattice_dict % add_key(lat % id, i) diff --git a/src/templates/geometry_t.xml b/src/templates/geometry_t.xml index ae8105a6bd..539a23f6b4 100644 --- a/src/templates/geometry_t.xml +++ b/src/templates/geometry_t.xml @@ -29,6 +29,7 @@ + From 92f88c4406ee2b5c23f8bff318ba33d6089ceb79 Mon Sep 17 00:00:00 2001 From: nhorelik Date: Tue, 9 Apr 2013 13:49:11 -0400 Subject: [PATCH 2/5] Cleaned up whitespace and comments --- src/geometry.F90 | 2 +- src/initialize.F90 | 3 +-- src/input_xml.F90 | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 9dd771d5d8..a8ecf6d01b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -80,7 +80,7 @@ contains real(8) :: xyz(3) ! temporary location logical :: use_search_cells ! use cells provided as argument logical :: outside_lattice ! if particle is not inside lattice bounds - logical :: lattice_edge ! if the particle is on a lattice edge + logical :: lattice_edge ! if particle is on a lattice edge type(Cell), pointer :: c ! pointer to cell type(Lattice), pointer :: lat ! pointer to lattice type(Universe), pointer :: univ ! universe to search in diff --git a/src/initialize.F90 b/src/initialize.F90 index 68eecf8a87..6c25af3c43 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -398,7 +398,6 @@ contains ! increment the index for the cells array within the Universe object and ! then store the index of the Cell object in that array index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1 - univ % cells(index_cell_in_univ(i_univ)) = i end do @@ -418,7 +417,7 @@ contains integer :: j ! index for various purposes integer :: k ! loop index for lattices integer :: m ! loop index for lattices - integer :: mid, lid ! material and lattice ids + integer :: mid, lid ! material and lattice IDs integer :: n_x, n_y, n_z ! size of lattice integer :: i_array ! index in surfaces/materials array integer :: id ! user-specified id diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 740b404566..7c5ad9e0fe 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -608,7 +608,7 @@ contains ! Get number of tags n_cells = size(cell_) - + ! Check for no cells if (n_cells == 0) then message = "No cells found in geometry.xml!" From 9884246221eb076176c07ab160792e319ea4b976 Mon Sep 17 00:00:00 2001 From: nhorelik Date: Wed, 10 Apr 2013 16:30:29 -0400 Subject: [PATCH 3/5] fixed lattice_outside so that distance_to_boundary tracks as expected --- src/geometry.F90 | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index a8ecf6d01b..9dafd75f61 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -164,8 +164,9 @@ contains ! Set current lattice lat => lattices(c % fill) - + outside_lattice = .false. + lattice_edge = .false. ! determine lattice index based on position xyz = p % coord % xyz + TINY_BIT * p % coord % uvw @@ -186,7 +187,6 @@ contains i_z < 1 .or. i_z > n_z) then ! Check for when particle is on lattice edge - lattice_edge = .false. if ( abs(xyz(1) - lat % lower_left(1)) < FP_COINCIDENT .or. & abs(xyz(2) - lat % lower_left(2)) < FP_COINCIDENT) then lattice_edge = .true. @@ -211,17 +211,23 @@ contains p % coord % xyz = xyz else - + ! We're outside the lattice, so treat this as a normal cell with ! the material specified for the outside outside_lattice = .true. p % last_material = p % material p % material = c % material - + + ! We'll still make a new coordinate for the particle, as + ! distance_to_boundary will still need to track through lattice + ! widths even though there's nothing in them but this material + end if - - else + + end if + + if (.not. lattice_edge) then ! Create new level of coordinates allocate(p % coord % next) @@ -247,13 +253,25 @@ contains p % coord % lattice_x = i_x p % coord % lattice_y = i_y p % coord % lattice_z = i_z - p % coord % universe = lat % universes(i_x,i_y,i_z) + if (.not. outside_lattice) then + p % coord % universe = lat % universes(i_x,i_y,i_z) + else + + ! If we define a void universe, we could + p % coord % universe = NONE + + ! Set coord cell for distance_to_boundary + p % coord % cell = index_cell + + end if + end if if (.not. outside_lattice) then call find_cell(found) if (.not. found) exit end if + end if ! Found cell so we can return From 40543168aebff63ce17cd9996f26ebefb9df7155 Mon Sep 17 00:00:00 2001 From: nhorelik Date: Wed, 10 Apr 2013 16:38:35 -0400 Subject: [PATCH 4/5] Changed universe on new coord for lattice outside to be the same as previously --- src/geometry.F90 | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 9dafd75f61..7f5e28ef3c 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -245,26 +245,26 @@ contains 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 = i_x - p % coord % lattice_y = i_y - p % coord % lattice_z = i_z + p % coord % next% lattice = c % fill + p % coord % next% lattice_x = i_x + p % coord % next% lattice_y = i_y + p % coord % next% lattice_z = i_z if (.not. outside_lattice) then - p % coord % universe = lat % universes(i_x,i_y,i_z) + p % coord % next % universe = lat % universes(i_x,i_y,i_z) else - ! If we define a void universe, we could - p % coord % universe = NONE + ! Set universe as the same for subsequent calls to find_cell + p % coord % next % universe = p % coord % universe - ! Set coord cell for distance_to_boundary - p % coord % cell = index_cell + ! Set coord cell for calls to distance_to_boundary + p % coord % next % cell = index_cell end if + ! Move particle to next level + p % coord => p % coord % next + end if if (.not. outside_lattice) then From c9d0bb94f77e15791f926f52bbd9c5de629fc786 Mon Sep 17 00:00:00 2001 From: Nicholas Horelik Date: Thu, 18 Apr 2013 22:01:17 -0300 Subject: [PATCH 5/5] Added check for particle on upper_right of lattice --- src/geometry.F90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 7f5e28ef3c..f0e7f4dd8e 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -78,6 +78,7 @@ contains integer :: n ! number of cells to search integer :: index_cell ! index in cells array real(8) :: xyz(3) ! temporary location + real(8) :: upper_right(3) ! lattice upper_right logical :: use_search_cells ! use cells provided as argument logical :: outside_lattice ! if particle is not inside lattice bounds logical :: lattice_edge ! if particle is on a lattice edge @@ -187,12 +188,21 @@ contains i_z < 1 .or. i_z > n_z) then ! Check for when particle is on lattice edge + upper_right(1) = lat % lower_left(1) + & + lat % width(1) * dble(lat % dimension(1)) + upper_right(2) = lat % lower_left(2) + & + lat % width(2) * dble(lat % dimension(2)) if ( abs(xyz(1) - lat % lower_left(1)) < FP_COINCIDENT .or. & - abs(xyz(2) - lat % lower_left(2)) < FP_COINCIDENT) then + abs(xyz(2) - lat % lower_left(2)) < FP_COINCIDENT .or. & + abs(upper_right(1) - xyz(1)) < FP_COINCIDENT .or. & + abs(upper_right(2) - xyz(2)) < FP_COINCIDENT) then lattice_edge = .true. end if if (lat % n_dimension == 3) then - if (abs(xyz(3) - lat % lower_left(3)) < FP_COINCIDENT) then + upper_right(3) = lat % lower_left(3) + & + lat % width(3) * dble(lat % dimension(3)) + if (abs(xyz(3) - lat % lower_left(3)) < FP_COINCIDENT .or. & + abs(upper_right(3) - xyz(3)) < FP_COINCIDENT) then lattice_edge = .true. end if end if