From a843d245b6ed2892c8322c27eabb457deb09e278 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 1 Dec 2011 17:55:25 -0500 Subject: [PATCH] Rewrote all geometry routines to allow for nested universes. --- src/DEPENDENCIES | 1 + src/constants.f90 | 3 + src/geometry.f90 | 1092 ++++++++++++++++++--------------------- src/mpi_routines.f90 | 13 +- src/output.f90 | 42 +- src/particle_header.f90 | 91 +++- src/physics.f90 | 94 ++-- src/plot.f90 | 189 +++---- src/source.f90 | 17 +- src/tally.f90 | 10 +- 10 files changed, 772 insertions(+), 780 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 178f95a961..85c0318589 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -140,6 +140,7 @@ output.o: string.o output.o: tally_header.o particle_header.o: constants.o +particle_header.o: geometry_header.o physics.o: constants.o physics.o: cross_section_header.o diff --git a/src/constants.f90 b/src/constants.f90 index 3e272e7c70..332df302ee 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -261,6 +261,9 @@ module constants ! ============================================================================ ! MISCELLANEOUS CONSTANTS + ! indicates that an array index hasn't been set + integer, parameter :: NONE = 0 + ! Codes for read errors -- better hope these numbers are never used in an ! input file! integer, parameter :: ERROR_INT = -huge(0) diff --git a/src/geometry.f90 b/src/geometry.f90 index dbb238fb48..ee81c9d9b1 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -6,7 +6,7 @@ module geometry use geometry_header, only: Cell, Surface, Universe, Lattice use global use output, only: write_message - use particle_header, only: Particle + use particle_header, only: Particle, LocalCoord, deallocate_coord use string, only: int_to_str use tally, only: score_surface_current @@ -62,7 +62,7 @@ contains ! Compare sense of point to specified sense specified_sense = sign(1,expression(i)) - actual_sense = sense(surf, p % xyz_local) + actual_sense = sense(surf, p % coord % xyz) if (actual_sense == specified_sense) then expression(i) = 1 else @@ -90,77 +90,116 @@ contains ! as it's within the geometry !=============================================================================== - recursive subroutine find_cell(univ, p, found) + recursive subroutine find_cell(p, found, search_cells) - type(Universe), pointer :: univ ! universe to search in - type(Particle), pointer :: p ! pointer to particle - logical, intent(inout) :: found ! particle found? + type(Particle), pointer :: p + logical, intent(inout) :: found + integer, optional :: search_cells(:) - integer :: i ! index over cells - integer :: x, y - type(Cell), pointer :: c ! pointer to cell - type(Lattice), pointer :: lat ! pointer to lattice - type(Universe), pointer :: lower_univ ! if particle is in lower universe, - ! use this pointer to call recursively + integer :: i ! index over cells + integer :: x ! x-index for lattice + integer :: y ! y-index for lattice + integer :: n ! number of cells to search + integer :: index_cell ! index in cells array + logical :: use_search_cells ! use cells provided as argument + type(Cell), pointer :: c ! pointer to cell + type(Lattice), pointer :: lat ! pointer to lattice + type(Universe), pointer :: univ ! universe to search in - found = .false. + ! set size of list to search + if (present(search_cells)) then + use_search_cells = .true. + n = size(search_cells) + else + use_search_cells = .false. + univ => universes(p % coord % universe) + n = univ % n_cells + end if - ! determine what region in - do i = 1, univ % n_cells - c => cells(univ % cells(i)) - - if (cell_contains(c, p)) then - ! If this cell contains a universe or lattice, search for the particle - ! in that universe/lattice - if (c % type == CELL_NORMAL) then - found = .true. - - ! set particle attributes - p % cell = univ % cells(i) - p % universe = dict_get_key(universe_dict, univ % id) - p % material = c % material - exit - elseif (c % type == CELL_FILL) then - lower_univ => universes(c % fill) - call find_cell(lower_univ, p, found) - if (found) then - exit - else - message = "Could not locate particle in universe: " - call fatal_error() - end if - elseif (c % type == CELL_LATTICE) then - ! Set current lattice - lat => lattices(c % fill) - p % lattice = c % fill - - ! determine universe based on lattice position - x = ceiling((p%xyz(1) - lat%x0)/lat%width_x) - y = ceiling((p%xyz(2) - lat%y0)/lat%width_y) - lower_univ => universes(lat % element(x,y)) - - ! adjust local position of particle - p%xyz_local(1) = p%xyz(1) - (lat%x0 + (x-0.5)*lat%width_x) - p%xyz_local(2) = p%xyz(2) - (lat%y0 + (y-0.5)*lat%width_y) - p%xyz_local(3) = p%xyz(3) - - ! set particle lattice indices - p % index_x = x - p % index_y = y - - call find_cell(lower_univ, p, found) - if (found) then - exit - else - message = "Could not locate particle in lattice: " & - & // int_to_str(lat % id) - call fatal_error() - end if - end if + do i = 1, n + ! select cells based on whether we are searching a universe or a provided + ! list of cells (this would be for lists of neighbor cells) + if (use_search_cells) then + index_cell = search_cells(i) + else + index_cell = univ % cells(i) end if + ! get pointer to cell + c => cells(index_cell) + + ! print *, 'searching cell ', c % id, cell_contains(c,p) + + if (cell_contains(c, p)) then + ! Set cell on this level + p % coord % cell = index_cell + + if (c % type == CELL_NORMAL) then + ! ================================================================= + ! AT LOWEST UNIVERSE, TERMINATE SEARCH + + ! set material + p % material = c % material + + elseif (c % type == CELL_FILL) then + ! ================================================================= + ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL + + ! Create new level of coordinates + p % in_lower_universe = .true. + allocate(p % coord % next) + + ! Move particle to next level and set universe + p % coord => p % coord % next + p % coord % universe = c % fill + + call find_cell(p, found) + if (.not. found) exit + + elseif (c % type == CELL_LATTICE) then + ! ================================================================= + ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL + + ! Set current lattice + lat => lattices(c % fill) + + ! determine universe based on lattice position + x = ceiling((p % coord % xyz(1) - lat % x0)/lat % width_x) + y = ceiling((p % coord % xyz(2) - lat % y0)/lat % width_y) + + ! Create new level of coordinates + p % in_lower_universe = .true. + allocate(p % coord % next) + + ! adjust local position of particle + p % coord % next % xyz(1) = p % coord % xyz(1) - & + (lat%x0 + (x-0.5_8)*lat%width_x) + 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 + + ! 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) + + call find_cell(p, found) + if (.not. found) exit + end if + + ! Found cell so we can return + found = .true. + return + end if end do + found = .false. + end subroutine find_cell !=============================================================================== @@ -173,10 +212,6 @@ contains type(Particle), pointer :: p integer, intent(in) :: last_cell ! last cell particle was in - integer :: i ! index of neighbors - integer :: index_cell ! index in cells array - integer :: i_x ! x index in lattice - integer :: i_y ! y index in lattice real(8) :: x ! x-x0 for sphere real(8) :: y ! y-y0 for sphere real(8) :: z ! z-z0 for sphere @@ -191,9 +226,6 @@ contains real(8) :: norm ! "norm" of surface normal logical :: found ! particle found in universe? type(Surface), pointer :: surf => null() - type(Cell), pointer :: c => null() - type(Lattice), pointer :: lat => null() - type(Universe), pointer :: lower_univ => null() surf => surfaces(abs(p % surface)) if (verbosity >= 10) then @@ -215,7 +247,7 @@ contains ! TODO: Find a better solution to score surface currents than physically ! moving the particle forward slightly - p % xyz = p % xyz + TINY_BIT * p % uvw + p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw call score_surface_current(p) ! Display message @@ -229,25 +261,31 @@ contains ! ======================================================================= ! PARTICLE REFLECTS FROM SURFACE + ! Do not handle reflective boundary conditions on lower universes + if (p % in_lower_universe) then + message = "Cannot reflect particle off surface in a lower universe." + call fatal_error() + end if + ! Score surface currents since reflection causes the direction of the ! particle to change -- artificially move the particle slightly back in ! case the surface crossing in coincident with a mesh boundary - p % xyz = p % xyz - TINY_BIT * p % uvw + p % coord0 % xyz = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw call score_surface_current(p) - p % xyz = p % xyz + TINY_BIT * p % uvw + p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw ! Copy particle's direction cosines - u = p % uvw(1) - v = p % uvw(2) - w = p % uvw(3) + u = p % coord0 % uvw(1) + v = p % coord0 % uvw(2) + w = p % coord0 % uvw(3) select case (surf%type) case (SURF_PX) - p % uvw = (/ -u, v, w /) + p % coord0 % uvw = (/ -u, v, w /) case (SURF_PY) - p % uvw = (/ u, -v, w /) + p % coord0 % uvw = (/ u, -v, w /) case (SURF_PZ) - p % uvw = (/ u, v, -w /) + p % coord0 % uvw = (/ u, v, -w /) case (SURF_PLANE) ! Find surface coefficients and norm of vector normal to surface n1 = surf % coeffs(1) @@ -262,11 +300,11 @@ contains w = w - 2*dot_prod*n3/norm ! Set vector - p % uvw = (/ u, v, w /) + p % coord0 % uvw = (/ u, v, w /) case (SURF_CYL_X) ! Find y-y0, z-z0 and dot product of direction and surface normal - y = p % xyz(2) - surf % coeffs(1) - z = p % xyz(3) - surf % coeffs(2) + y = p % coord0 % xyz(2) - surf % coeffs(1) + z = p % coord0 % xyz(3) - surf % coeffs(2) R = surf % coeffs(3) dot_prod = v*y + w*z @@ -275,11 +313,11 @@ contains w = w - 2*dot_prod*z/(R*R) ! Set vector - p % uvw = (/ u, v, w /) + p % coord0 % uvw = (/ u, v, w /) case (SURF_CYL_Y) ! Find x-x0, z-z0 and dot product of direction and surface normal - x = p % xyz(1) - surf % coeffs(1) - z = p % xyz(3) - surf % coeffs(2) + x = p % coord0 % xyz(1) - surf % coeffs(1) + z = p % coord0 % xyz(3) - surf % coeffs(2) R = surf % coeffs(3) dot_prod = u*x + w*z @@ -288,11 +326,11 @@ contains w = w - 2*dot_prod*z/(R*R) ! Set vector - p % uvw = (/ u, v, w /) + p % coord0 % uvw = (/ u, v, w /) case (SURF_CYL_Z) ! Find x-x0, y-y0 and dot product of direction and surface normal - x = p % xyz(1) - surf % coeffs(1) - y = p % xyz(2) - surf % coeffs(2) + x = p % coord0 % xyz(1) - surf % coeffs(1) + y = p % coord0 % xyz(2) - surf % coeffs(2) R = surf % coeffs(3) dot_prod = u*x + v*y @@ -301,13 +339,13 @@ contains v = v - 2*dot_prod*y/(R*R) ! Set vector - p % uvw = (/ u, v, w /) + p % coord0 % uvw = (/ u, v, w /) case (SURF_SPHERE) ! Find x-x0, y-y0, z-z0 and dot product of direction and surface ! normal - x = p % xyz(1) - surf % coeffs(1) - y = p % xyz(2) - surf % coeffs(2) - z = p % xyz(3) - surf % coeffs(3) + x = p % coord0 % xyz(1) - surf % coeffs(1) + y = p % coord0 % xyz(2) - surf % coeffs(2) + z = p % coord0 % xyz(3) - surf % coeffs(3) R = surf % coeffs(4) dot_prod = u*x + v*y + w*z @@ -317,7 +355,7 @@ contains w = w - 2*dot_prod*z/(R*R) ! Set vector - p % uvw = (/ u, v, w /) + p % coord0 % uvw = (/ u, v, w /) case default message = "Reflection not supported for surface " // & trim(int_to_str(surf % id)) @@ -325,11 +363,11 @@ contains end select ! Reassign particle's cell and surface - p % cell = last_cell + p % coord0 % cell = last_cell p % surface = -p % surface ! Set previous coordinate going slightly past surface crossing - p % last_xyz = p % xyz + TINY_BIT * p % uvw + p % last_xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw ! Diagnostic message if (verbosity >= 10) then @@ -345,115 +383,30 @@ contains if (p % surface > 0 .and. allocated(surf % neighbor_pos)) then ! If coming from negative side of surface, search all the neighboring ! cells on the positive side - do i = 1, size(surf % neighbor_pos) - index_cell = surf % neighbor_pos(i) - c => cells(index_cell) - if (cell_contains(c, p)) then - if (c % type == CELL_FILL) then - lower_univ => universes(c % fill) - call find_cell(lower_univ, p, found) - if (.not. found) then - message = "Could not locate particle in universe: " - call fatal_error() - end if - elseif (c % type == CELL_LATTICE) then - ! Set current lattice - lat => lattices(c % fill) - p % lattice = c % fill + + call find_cell(p, found, surf % neighbor_pos) + if (found) return - ! determine universe based on lattice position - i_x = ceiling((p%xyz(1) - lat%x0)/lat%width_x) - i_y = ceiling((p%xyz(2) - lat%y0)/lat%width_y) - lower_univ => universes(lat % element(i_x,i_y)) - - ! adjust local position of particle - p%xyz_local(1) = p%xyz(1) - (lat%x0 + (i_x-0.5)*lat%width_x) - p%xyz_local(2) = p%xyz(2) - (lat%y0 + (i_y-0.5)*lat%width_y) - p%xyz_local(3) = p%xyz(3) - - ! set particle lattice indices - p % index_x = i_x - p % index_y = i_y - - call find_cell(lower_univ, p, found) - if (.not. found) then - message = "Could not locate particle in lattice: " // & - trim(int_to_str(lat % id)) - call fatal_error() - end if - else - ! set current pointers - p % cell = index_cell - p % material = c % material - end if - return - end if - end do elseif (p % surface < 0 .and. allocated(surf % neighbor_neg)) then ! If coming from positive side of surface, search all the neighboring ! cells on the negative side - do i = 1, size(surf % neighbor_neg) - index_cell = surf % neighbor_neg(i) - c => cells(index_cell) - if (cell_contains(c, p)) then - if (c % type == CELL_FILL) then - lower_univ => universes(c % fill) - call find_cell(lower_univ, p, found) - if (.not. found) then - message = "Could not locate particle in universe: " - call fatal_error() - end if - elseif (c % type == CELL_LATTICE) then - ! Set current lattice - lat => lattices(c % fill) - p % lattice = c % fill + + call find_cell(p, found, surf % neighbor_neg) + if (found) return - ! determine universe based on lattice position - i_x = ceiling((p%xyz(1) - lat%x0)/lat%width_x) - i_y = ceiling((p%xyz(2) - lat%y0)/lat%width_y) - lower_univ => universes(lat % element(i_x,i_y)) - - ! adjust local position of particle - p%xyz_local(1) = p%xyz(1) - (lat%x0 + (i_x-0.5)*lat%width_x) - p%xyz_local(2) = p%xyz(2) - (lat%y0 + (i_y-0.5)*lat%width_y) - p%xyz_local(3) = p%xyz(3) - - ! set particle lattice indices - p % index_x = i_x - p % index_y = i_y - - call find_cell(lower_univ, p, found) - if (.not. found) then - message = "Could not locate particle in lattice: " // & - trim(int_to_str(lat % id)) - call fatal_error() - end if - else - ! set current pointers - p % cell = index_cell - p % material = c % material - end if - return - end if - end do end if ! ========================================================================== ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS - do i = 1, size(cells) - c => cells(i) - if (cell_contains(c, p)) then - p % cell = i - p % material = c % material - return - end if - end do + call find_cell(p, found) ! Couldn't find next cell anywhere! - message = "After particle crossed surface " // trim(int_to_str(p%surface)) & - // ", it could not be located in any cell and it did not leak." - call fatal_error() + if (.not. found) then + message = "After particle crossed surface " // trim(int_to_str(p%surface)) & + // ", it could not be located in any cell and it did not leak." + call fatal_error() + end if end subroutine cross_surface @@ -480,25 +433,24 @@ contains 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 - type(Universe), pointer :: univ + type(Lattice), pointer :: lat => null() - lat => lattices(p % lattice) + lat => lattices(p % coord % lattice) if (verbosity >= 10) then message = " Crossing lattice " // int_to_str(lat % id) call write_message() end if - u = p % uvw(1) - v = p % uvw(2) + u = p % coord % uvw(1) + v = p % coord % uvw(2) if (lat % type == LATTICE_RECT) then - x = p % xyz_local(1) - y = p % xyz_local(2) - z = p % xyz_local(3) - x0 = lat % width_x * 0.5 - y0 = lat % width_y * 0.5 + x = p % coord % xyz(1) + y = p % coord % xyz(2) + z = p % coord % xyz(3) + x0 = lat % width_x * 0.5_8 + y0 = lat % width_y * 0.5_8 dist = INFINITY @@ -523,23 +475,23 @@ contains dist = min(d_left, d_right, d_top, d_bottom) if (dist == d_left) then ! Move particle to left element - p % index_x = p % index_x - 1 - p % xyz_local(1) = x0 + p % coord % lattice_x = p % coord % lattice_x - 1 + p % coord % xyz(1) = x0 elseif (dist == d_right) then ! Move particle to right element - p % index_x = p % index_x + 1 - p % xyz_local(1) = -x0 + p % coord % lattice_x = p % coord % lattice_x + 1 + p % coord % xyz(1) = -x0 elseif (dist == d_bottom) then ! Move particle to bottom element - p % index_y = p % index_y - 1 - p % xyz_local(2) = y0 + p % coord % lattice_y = p % coord % lattice_y - 1 + p % coord % xyz(2) = y0 elseif (dist == d_top) then ! Move particle to top element - p % index_y = p % index_y + 1 - p % xyz_local(2) = -y0 + p % coord % lattice_y = p % coord % lattice_y + 1 + p % coord % xyz(2) = -y0 end if elseif (lat % type == LATTICE_HEX) then @@ -547,8 +499,8 @@ contains end if ! Check to make sure still in lattice - i_x = p % index_x - i_y = p % index_y + i_x = p % coord % lattice_x + i_y = p % coord % lattice_y if (i_x < 1 .or. i_x > lat % n_x) then message = "Reached edge of lattice " // trim(int_to_str(lat % id)) // & " at position (" // trim(int_to_str(i_x)) // "," // & @@ -562,12 +514,13 @@ contains end if ! Find universe for next lattice element - univ => universes(lat % element(i_x,i_y)) + p % coord % universe = lat % element(i_x, i_y) ! Find cell in next lattice element - call find_cell(univ, p, found) + call find_cell(p, found) if (.not. found) then - message = "Could not locate particle in universe: " + message = "Could not locate particle in universe: " // & + int_to_str(universes(p % coord % universe) % id) call fatal_error() end if @@ -579,20 +532,15 @@ contains ! that has a parent cell, also include the surfaces of the edge of the universe. !=============================================================================== - subroutine distance_to_boundary(p, dist, surf, in_lattice) + subroutine distance_to_boundary(p, dist, surface_crossed, lattice_crossed) type(Particle), pointer :: p real(8), intent(out) :: dist - integer, intent(out) :: surf - logical, intent(out) :: in_lattice + integer, intent(out) :: surface_crossed + logical, intent(out) :: lattice_crossed - integer, allocatable :: expression(:) ! copy of surface list integer :: i ! index for surface in cell - integer :: n_surf ! total number of surfaces to check - integer :: n1 ! number of surfaces in current cell - integer :: n2 ! number of surfaces in parent cell integer :: index_surf ! index in surfaces array (with sign) - integer :: current_surf ! current surface real(8) :: x,y,z ! particle coordinates real(8) :: u,v,w ! particle directions real(8) :: d ! evaluated distance @@ -602,378 +550,372 @@ contains real(8) :: a,b,c,k ! quadratic equation coefficients real(8) :: quad ! discriminant of quadratic equation logical :: on_surface ! is particle on surface? - type(Cell), pointer :: cell_p => null() - type(Cell), pointer :: parent_p => null() - type(Surface), pointer :: surf_p => null() - type(Lattice), pointer :: lat => null() - - if (p % cell == CELL_VOID) then - n_surf = n_surfaces - allocate(expression(n_surfaces)) - expression = (/ (i, i=1, n_surfaces) /) - current_surf = 0 - else - cell_p => cells(p%cell) - - current_surf = p%surface - - ! determine number of surfaces to check - n1 = cell_p % n_surfaces - n2 = 0 - if (cell_p % parent > 0) then - parent_p => cells(cell_p % parent) - n2 = parent_p % n_surfaces - end if - n_surf = n1 + n2 - - ! allocate space and assign expression - allocate(expression(n_surf)) - expression(1:n1) = cell_p % surfaces - if (cell_p % parent > 0) then - expression(n1+1:n1+n2) = parent_p % surfaces - end if - end if - - u = p % uvw(1) - v = p % uvw(2) - w = p % uvw(3) - - ! ========================================================================== - ! FIND MINIMUM DISTANCE TO SURFACES IN CELL AND PARENT + type(Cell), pointer :: cl => null() + type(Surface), pointer :: surf => null() + type(Lattice), pointer :: lat => null() + type(LocalCoord), pointer :: coord => null() + type(LocalCoord), pointer :: final_coord => null() + ! inialize distance to infinity (huge) dist = INFINITY - do i = 1, n_surf - if (i <= n1) then - ! in local cell, so use xyz_local - x = p % xyz_local(1) - y = p % xyz_local(2) - z = p % xyz_local(3) - else - ! in parent cell, so use xyz - x = p % xyz(1) - y = p % xyz(2) - z = p % xyz(3) - end if - ! check for coincident surface -- note that we can't skip the calculation - ! in general because a particle could be on one side of a cylinder and - ! still have a positive distance to the other + ! Get pointer to top-level coordinates + coord => p % coord0 - index_surf = expression(i) - if (index_surf == current_surf) then - on_surface = .true. - else - on_surface = .false. - end if + ! Loop over each universe level + LEVEL_LOOP: do while(associated(coord)) - ! check for operators - index_surf = abs(index_surf) - if (index_surf >= OP_DIFFERENCE) cycle + ! get pointer to cell on this level + cl => cells(coord % cell) - surf_p => surfaces(index_surf) - select case (surf_p%type) - case (SURF_PX) - if (on_surface .or. u == ZERO) then - d = INFINITY + ! copy directional cosines + u = coord % uvw(1) + v = coord % uvw(2) + w = coord % uvw(3) + + ! ======================================================================= + ! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL + + SURFACE_LOOP: do i = 1, cl % n_surfaces + + ! copy local coordinates of particle + x = coord % xyz(1) + y = coord % xyz(2) + z = coord % xyz(3) + + ! check for coincident surface -- note that we can't skip the + ! calculation in general because a particle could be on one side of a + ! cylinder and still have a positive distance to the other + + index_surf = cl % surfaces(i) + if (index_surf == p % surface) then + on_surface = .true. else - x0 = surf_p % coeffs(1) - d = (x0 - x)/u - if (d < ZERO) d = INFINITY - end if - - case (SURF_PY) - if (on_surface .or. v == ZERO) then - d = INFINITY - else - y0 = surf_p % coeffs(1) - d = (y0 - y)/v - if (d < ZERO) d = INFINITY - end if - - case (SURF_PZ) - if (on_surface .or. w == ZERO) then - d = INFINITY - else - z0 = surf_p % coeffs(1) - d = (z0 - z)/w - if (d < ZERO) d = INFINITY - end if - - case (SURF_PLANE) - A = surf_p % coeffs(1) - B = surf_p % coeffs(2) - C = surf_p % coeffs(3) - D = surf_p % coeffs(4) - - tmp = A*u + B*v + C*w - if (on_surface .or. tmp == ZERO) then - d = INFINITY - else - d = -(A*x + B*y + C*w - D)/tmp - if (d < ZERO) d = INFINITY + on_surface = .false. end if - case (SURF_CYL_X) - a = ONE - u*u ! v^2 + w^2 - if (a == ZERO) then - d = INFINITY - else - y0 = surf_p % coeffs(1) - z0 = surf_p % coeffs(2) - r = surf_p % coeffs(3) + ! check for operators + index_surf = abs(index_surf) + if (index_surf >= OP_DIFFERENCE) cycle - y = y - y0 - z = z - z0 - k = y*v + z*w - c = y*y + z*z - r*r - quad = k*k - a*c + ! get pointer to surface + surf => surfaces(index_surf) - if (quad < ZERO) then - ! no intersection with cylinder - - d = INFINITY - - elseif (on_surface) then - ! particle is on the cylinder, thus one distance is - ! positive/negative and the other is zero. The sign of k - ! determines if we are facing in or out - - if (k >= ZERO) then - d = INFINITY - else - d = (-k + sqrt(quad))/a - end if - - elseif (c < ZERO) then - ! particle is inside the cylinder, thus one distance must be - ! negative and one must be positive. The positive distance will - ! be the one with negative sign on sqrt(quad) - - d = (-k + sqrt(quad))/a - - else - ! particle is outside the cylinder, thus both distances are - ! either positive or negative. If positive, the smaller distance - ! is the one with positive sign on sqrt(quad) - - d = (-k - sqrt(quad))/a - if (d < ZERO) d = INFINITY - - end if - end if - - case (SURF_CYL_Y) - a = ONE - v*v ! u^2 + w^2 - if (a == ZERO) then - d = INFINITY - else - x0 = surf_p % coeffs(1) - z0 = surf_p % coeffs(2) - r = surf_p % coeffs(3) - - x = x - x0 - z = z - z0 - k = x*u + z*w - c = x*x + z*z - r*r - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cylinder - - d = INFINITY - - elseif (on_surface) then - ! particle is on the cylinder, thus one distance is - ! positive/negative and the other is zero. The sign of k - ! determines if we are facing in or out - - if (k >= ZERO) then - d = INFINITY - else - d = (-k + sqrt(quad))/a - end if - - elseif (c < ZERO) then - ! particle is inside the cylinder, thus one distance must be - ! negative and one must be positive. The positive distance will - ! be the one with negative sign on sqrt(quad) - - d = (-k + sqrt(quad))/a - - else - ! particle is outside the cylinder, thus both distances are - ! either positive or negative. If positive, the smaller distance - ! is the one with positive sign on sqrt(quad) - - d = (-k - sqrt(quad))/a - if (d < ZERO) d = INFINITY - - end if - end if - - case (SURF_CYL_Z) - a = ONE - w*w ! u^2 + v^2 - if (a == ZERO) then - d = INFINITY - else - x0 = surf_p % coeffs(1) - y0 = surf_p % coeffs(2) - r = surf_p % coeffs(3) - - x = x - x0 - y = y - y0 - k = x*u + y*v - c = x*x + y*y - r*r - quad = k*k - a*c - - if (quad < ZERO) then - ! no intersection with cylinder - - d = INFINITY - - elseif (on_surface) then - ! particle is on the cylinder, thus one distance is - ! positive/negative and the other is zero. The sign of k - ! determines if we are facing in or out - - if (k >= ZERO) then - d = INFINITY - else - d = (-k + sqrt(quad))/a - end if - - elseif (c < ZERO) then - ! particle is inside the cylinder, thus one distance must be - ! negative and one must be positive. The positive distance will - ! be the one with negative sign on sqrt(quad) - - d = (-k + sqrt(quad))/a - - else - ! particle is outside the cylinder, thus both distances are - ! either positive or negative. If positive, the smaller distance - ! is the one with positive sign on sqrt(quad) - - d = (-k - sqrt(quad))/a - if (d <= ZERO) d = INFINITY - - end if - end if - - case (SURF_SPHERE) - x0 = surf_p % coeffs(1) - y0 = surf_p % coeffs(2) - z0 = surf_p % coeffs(3) - r = surf_p % coeffs(4) - - x = x - x0 - y = y - y0 - z = z - z0 - k = x*u + y*v + z*w - c = x*x + y*y + z*z - r*r - quad = k*k - c - - if (quad < ZERO) then - ! no intersection with sphere - - d = INFINITY - - elseif (on_surface) then - ! particle is on the sphere, thus one distance is positive/negative - ! and the other is zero. The sign of k determines if we are facing - ! in or out - - if (k >= ZERO) then + select case (surf % type) + case (SURF_PX) + if (on_surface .or. u == ZERO) then d = INFINITY else + x0 = surf % coeffs(1) + d = (x0 - x)/u + if (d < ZERO) d = INFINITY + end if + + case (SURF_PY) + if (on_surface .or. v == ZERO) then + d = INFINITY + else + y0 = surf % coeffs(1) + d = (y0 - y)/v + if (d < ZERO) d = INFINITY + end if + + case (SURF_PZ) + if (on_surface .or. w == ZERO) then + d = INFINITY + else + z0 = surf % coeffs(1) + d = (z0 - z)/w + if (d < ZERO) d = INFINITY + end if + + case (SURF_PLANE) + A = surf % coeffs(1) + B = surf % coeffs(2) + C = surf % coeffs(3) + D = surf % coeffs(4) + + tmp = A*u + B*v + C*w + if (on_surface .or. tmp == ZERO) then + d = INFINITY + else + d = -(A*x + B*y + C*w - D)/tmp + if (d < ZERO) d = INFINITY + end if + + case (SURF_CYL_X) + a = ONE - u*u ! v^2 + w^2 + if (a == ZERO) then + d = INFINITY + else + y0 = surf % coeffs(1) + z0 = surf % coeffs(2) + r = surf % coeffs(3) + + y = y - y0 + z = z - z0 + k = y*v + z*w + c = y*y + z*z - r*r + quad = k*k - a*c + + if (quad < ZERO) then + ! no intersection with cylinder + + d = INFINITY + + elseif (on_surface) then + ! particle is on the cylinder, thus one distance is + ! positive/negative and the other is zero. The sign of k + ! determines if we are facing in or out + + if (k >= ZERO) then + d = INFINITY + else + d = (-k + sqrt(quad))/a + end if + + elseif (c < ZERO) then + ! particle is inside the cylinder, thus one distance must be + ! negative and one must be positive. The positive distance + ! will be the one with negative sign on sqrt(quad) + + d = (-k + sqrt(quad))/a + + else + ! particle is outside the cylinder, thus both distances are + ! either positive or negative. If positive, the smaller + ! distance is the one with positive sign on sqrt(quad) + + d = (-k - sqrt(quad))/a + if (d < ZERO) d = INFINITY + + end if + end if + + case (SURF_CYL_Y) + a = ONE - v*v ! u^2 + w^2 + if (a == ZERO) then + d = INFINITY + else + x0 = surf % coeffs(1) + z0 = surf % coeffs(2) + r = surf % coeffs(3) + + x = x - x0 + z = z - z0 + k = x*u + z*w + c = x*x + z*z - r*r + quad = k*k - a*c + + if (quad < ZERO) then + ! no intersection with cylinder + + d = INFINITY + + elseif (on_surface) then + ! particle is on the cylinder, thus one distance is + ! positive/negative and the other is zero. The sign of k + ! determines if we are facing in or out + + if (k >= ZERO) then + d = INFINITY + else + d = (-k + sqrt(quad))/a + end if + + elseif (c < ZERO) then + ! particle is inside the cylinder, thus one distance must be + ! negative and one must be positive. The positive distance + ! will be the one with negative sign on sqrt(quad) + + d = (-k + sqrt(quad))/a + + else + ! particle is outside the cylinder, thus both distances are + ! either positive or negative. If positive, the smaller + ! distance is the one with positive sign on sqrt(quad) + + d = (-k - sqrt(quad))/a + if (d < ZERO) d = INFINITY + + end if + end if + + case (SURF_CYL_Z) + a = ONE - w*w ! u^2 + v^2 + if (a == ZERO) then + d = INFINITY + else + x0 = surf % coeffs(1) + y0 = surf % coeffs(2) + r = surf % coeffs(3) + + x = x - x0 + y = y - y0 + k = x*u + y*v + c = x*x + y*y - r*r + quad = k*k - a*c + + if (quad < ZERO) then + ! no intersection with cylinder + + d = INFINITY + + elseif (on_surface) then + ! particle is on the cylinder, thus one distance is + ! positive/negative and the other is zero. The sign of k + ! determines if we are facing in or out + + if (k >= ZERO) then + d = INFINITY + else + d = (-k + sqrt(quad))/a + end if + + elseif (c < ZERO) then + ! particle is inside the cylinder, thus one distance must be + ! negative and one must be positive. The positive distance + ! will be the one with negative sign on sqrt(quad) + + d = (-k + sqrt(quad))/a + + else + ! particle is outside the cylinder, thus both distances are + ! either positive or negative. If positive, the smaller + ! distance is the one with positive sign on sqrt(quad) + + d = (-k - sqrt(quad))/a + if (d <= ZERO) d = INFINITY + + end if + end if + + case (SURF_SPHERE) + x0 = surf % coeffs(1) + y0 = surf % coeffs(2) + z0 = surf % coeffs(3) + r = surf % coeffs(4) + + x = x - x0 + y = y - y0 + z = z - z0 + k = x*u + y*v + z*w + c = x*x + y*y + z*z - r*r + quad = k*k - c + + if (quad < ZERO) then + ! no intersection with sphere + + d = INFINITY + + elseif (on_surface) then + ! particle is on the sphere, thus one distance is + ! positive/negative and the other is zero. The sign of k + ! determines if we are facing in or out + + if (k >= ZERO) then + d = INFINITY + else + d = -k + sqrt(quad) + end if + + elseif (c < ZERO) then + ! particle is inside the sphere, thus one distance must be + ! negative and one must be positive. The positive distance will + ! be the one with negative sign on sqrt(quad) + d = -k + sqrt(quad) + + else + ! particle is outside the sphere, thus both distances are either + ! positive or negative. If positive, the smaller distance is the + ! one with positive sign on sqrt(quad) + + d = -k - sqrt(quad) + if (d < ZERO) d = INFINITY + end if - elseif (c < ZERO) then - ! particle is inside the sphere, thus one distance must be negative - ! and one must be positive. The positive distance will be the one - ! with negative sign on sqrt(quad) + case (SURF_GQ) + message = "Surface distance not yet implement for general quadratic." + call fatal_error() - d = -k + sqrt(quad) - - else - ! particle is outside the sphere, thus both distances are either - ! positive or negative. If positive, the smaller distance is the - ! one with positive sign on sqrt(quad) - - d = -k - sqrt(quad) - if (d < ZERO) d = INFINITY - - end if - - case (SURF_GQ) - message = "Surface distance not yet implement for general quadratic." - call fatal_error() - - end select - - ! Check is calculated distance is new minimum - if (d < dist) then - dist = d - surf = -expression(i) - end if - end do - - ! ========================================================================== - ! FIND MINIMUM DISTANCE TO LATTICE SURFACES - - in_lattice = .false. - if (p % lattice > 0) then - lat => lattices(p % lattice) - if (lat % type == LATTICE_RECT) then - x = p % xyz_local(1) - y = p % xyz_local(2) - z = p % xyz_local(3) - x0 = lat % width_x * 0.5_8 - y0 = lat % width_y * 0.5_8 - - ! left and right sides - if (u > 0) then - d = (x0 - x)/u - else - d = -(x + x0)/u - end if - - ! If the lattice boundary is coincident with the parent cell boundary, - ! we need to make sure that the lattice is not selected. This is - ! complicated by the fact that floating point may determine that one - ! is closer than the other (can't check direct equality). Thus, the - ! logic here checks whether the relative difference is within floating - ! point precision. - - if (d < dist) then - if (abs(d - dist)/dist >= FP_PRECISION) then - dist = d - in_lattice = .true. - end if - end if - - ! top and bottom sides - if (v > 0) then - d = (y0 - y)/v - else - d = -(y + y0)/v - end if + end select + ! Check is calculated distance is new minimum if (d < dist) then - if (abs(d - dist)/dist >= FP_PRECISION) then - dist = d - in_lattice = .true. - end if + dist = d + surface_crossed = -cl % surfaces(i) + final_coord => coord end if - elseif (lat % type == LATTICE_HEX) then - ! TODO: Add hex lattice support - end if - end if + end do SURFACE_LOOP - ! deallocate expression - deallocate(expression) + ! ======================================================================= + ! FIND MINIMUM DISTANCE TO LATTICE SURFACES + + lattice_crossed = .false. + + if (p % coord % lattice /= NONE) then + lat => lattices(p % coord % lattice) + if (lat % type == LATTICE_RECT) then + x = p % coord % xyz(1) + y = p % coord % xyz(2) + z = p % coord % xyz(3) + x0 = lat % width_x * 0.5_8 + y0 = lat % width_y * 0.5_8 + + ! left and right sides + if (u > 0) then + d = (x0 - x)/u + else + d = -(x + x0)/u + end if + + ! If the lattice boundary is coincident with the parent cell boundary, + ! we need to make sure that the lattice is not selected. This is + ! complicated by the fact that floating point may determine that one + ! is closer than the other (can't check direct equality). Thus, the + ! logic here checks whether the relative difference is within floating + ! point precision. + + if (d < dist) then + if (abs(d - dist)/dist >= FP_PRECISION) then + dist = d + lattice_crossed = .true. + final_coord => coord + end if + end if + + ! top and bottom sides + if (v > 0) then + d = (y0 - y)/v + else + d = -(y + y0)/v + end if + + if (d < dist) then + if (abs(d - dist)/dist >= FP_PRECISION) then + dist = d + lattice_crossed = .true. + final_coord => coord + end if + end if + + elseif (lat % type == LATTICE_HEX) then + ! TODO: Add hex lattice support + end if + end if + + coord => coord % next + + end do LEVEL_LOOP + + ! Remove coordinates for any lower levels + call deallocate_coord(final_coord % next) + + ! Move particle to appropriate coordinate level + p % coord => final_coord end subroutine distance_to_boundary diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90 index 2ea05e23f5..a523694772 100644 --- a/src/mpi_routines.f90 +++ b/src/mpi_routines.f90 @@ -347,16 +347,15 @@ contains i_source = i_start + i - 1 p => source_bank(i_source) - p % xyz = temp_bank(i) % xyz - p % xyz_local = temp_bank(i) % xyz - p % last_xyz = temp_bank(i) % xyz - p % uvw = temp_bank(i) % uvw - p % E = temp_bank(i) % E - p % last_E = p % E - ! set defaults call initialize_particle(p) + p % coord % xyz = temp_bank(i) % xyz + p % coord % uvw = temp_bank(i) % uvw + p % last_xyz = temp_bank(i) % xyz + p % E = temp_bank(i) % E + p % last_E = temp_bank(i) % E + end do end subroutine copy_from_bank diff --git a/src/output.f90 b/src/output.f90 index 3f52da8a05..720683ef6e 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -210,44 +210,40 @@ contains case default write(ou,*) 'Unknown Particle ' // int_to_str(p % id) end select - write(ou,*) ' x = ' // real_to_str(p % xyz(1)) - write(ou,*) ' y = ' // real_to_str(p % xyz(2)) - write(ou,*) ' z = ' // real_to_str(p % xyz(3)) - write(ou,*) ' x local = ' // real_to_str(p % xyz_local(1)) - write(ou,*) ' y local = ' // real_to_str(p % xyz_local(2)) - write(ou,*) ' z local = ' // real_to_str(p % xyz_local(3)) - write(ou,*) ' u = ' // real_to_str(p % uvw(1)) - write(ou,*) ' v = ' // real_to_str(p % uvw(2)) - write(ou,*) ' w = ' // real_to_str(p % uvw(3)) + write(ou,*) ' x = ' // real_to_str(p % coord0 % xyz(1)) + write(ou,*) ' y = ' // real_to_str(p % coord0 % xyz(2)) + write(ou,*) ' z = ' // real_to_str(p % coord0 % xyz(3)) + write(ou,*) ' x local = ' // real_to_str(p % coord % xyz(1)) + write(ou,*) ' y local = ' // real_to_str(p % coord % xyz(2)) + write(ou,*) ' z local = ' // real_to_str(p % coord % xyz(3)) + write(ou,*) ' u = ' // real_to_str(p % coord0 % uvw(1)) + write(ou,*) ' v = ' // real_to_str(p % coord0 % uvw(2)) + write(ou,*) ' w = ' // real_to_str(p % coord0 % uvw(3)) write(ou,*) ' Weight = ' // real_to_str(p % wgt) write(ou,*) ' Energy = ' // real_to_str(p % E) - write(ou,*) ' x index = ' // int_to_str(p % index_x) - write(ou,*) ' y index = ' // int_to_str(p % index_y) + write(ou,*) ' x index = ' // int_to_str(p % coord % lattice_x) + write(ou,*) ' y index = ' // int_to_str(p % coord % lattice_y) write(ou,*) ' IE = ' // int_to_str(p % IE) write(ou,*) ' Interpolation factor = ' // real_to_str(p % interp) - if (p % cell > 0) then - c => cells(p % cell) + if (p % coord % cell /= NONE) then + c => cells(p % coord % cell) write(ou,*) ' Cell = ' // int_to_str(c % id) else write(ou,*) ' Cell not determined' end if - if (p % surface > 0) then + if (p % surface /= NONE) then s => surfaces(p % surface) write(ou,*) ' Surface = ' // int_to_str(s % id) else write(ou,*) ' Surface = None' end if - u => universes(p % universe) + u => universes(p % coord % universe) write(ou,*) ' Universe = ' // int_to_str(u % id) write(ou,*) - nullify(c) - nullify(s) - nullify(u) - end subroutine print_particle !=============================================================================== @@ -325,10 +321,6 @@ contains write(ou,*) ' Surface Specification:' // trim(string) write(ou,*) - ! nullify associated pointers - nullify(u) - nullify(m) - end subroutine print_cell !=============================================================================== @@ -353,8 +345,6 @@ contains write(ou,*) ' Cells =' // trim(string) write(ou,*) - nullify(c) - end subroutine print_universe !=============================================================================== @@ -485,8 +475,6 @@ contains end if write(ou,*) - nullify(nuc) - end subroutine print_material !=============================================================================== diff --git a/src/particle_header.f90 b/src/particle_header.f90 index deec93b065..64fbc2e6d5 100644 --- a/src/particle_header.f90 +++ b/src/particle_header.f90 @@ -1,9 +1,33 @@ module particle_header - use constants, only: NEUTRON, ONE + use constants, only: NEUTRON, ONE, NONE + use geometry_header, only: BASE_UNIVERSE implicit none +!=============================================================================== +! LOCALCOORD describes the location of a particle local to a single +! universe. When the geometry consists of nested universes, a particle will have +! a list of coordinates in each level +!=============================================================================== + + type LocalCoord + ! Indices in various arrays for this level + integer :: cell = NONE + integer :: universe = NONE + integer :: lattice = NONE + integer :: lattice_x = NONE + integer :: lattice_y = NONE + + ! Particle position and direction for this level + real(8) :: xyz(3) + real(8) :: uvw(3) + + ! Pointers to next (lower) and previous (higher) universe + type(LocalCoord), pointer :: next => null() + type(LocalCoord), pointer :: prev => null() + end type LocalCoord + !=============================================================================== ! PARTICLE describes the state of a particle being transported through the ! geometry @@ -14,10 +38,12 @@ module particle_header integer(8) :: id ! Unique ID integer :: type ! Particle type (n, p, e, etc) - ! Physical data - real(8) :: xyz(3) ! location - real(8) :: xyz_local(3) ! local location (after transformations) - real(8) :: uvw(3) ! directional cosines + ! Particle coordinates + logical :: in_lower_universe ! is particle in lower universe? + type(LocalCoord), pointer :: coord0 ! coordinates on universe 0 + type(LocalCoord), pointer :: coord ! coordinates on lowest universe + + ! Other physical data real(8) :: wgt ! particle weight real(8) :: E ! energy real(8) :: mu ! angle of scatter @@ -36,15 +62,10 @@ module particle_header real(8) :: interp ! interpolation factor for energy grid ! Indices for various arrays - integer :: cell ! index for current cell + integer :: surface ! index for surface particle is on integer :: cell_born ! index for cell particle was born in - integer :: universe ! index for current universe - integer :: lattice ! index for current lattice - integer :: surface ! index for current surface integer :: material ! index for current material integer :: last_material ! index for last material - integer :: index_x ! lattice index for x direction - integer :: index_y ! lattice index for y direction ! Statistical data integer :: n_collision ! # of collisions @@ -66,22 +87,48 @@ contains ! passed through the fission bank to the source bank, no lookup would be ! needed at the beginning of a cycle - p % type = NEUTRON + p % type = NEUTRON + p % alive = .true. + + ! clear attributes + p % surface = NONE + p % cell_born = NONE + p % material = NONE + p % last_material = NONE p % wgt = ONE p % last_wgt = ONE - p % alive = .true. p % n_bank = 0 - p % cell = 0 - p % cell_born = 0 - p % universe = 0 - p % lattice = 0 - p % surface = 0 - p % material = 0 - p % last_material = 0 - p % index_x = 0 - p % index_y = 0 p % n_collision = 0 + ! remove any original coordinates + call deallocate_coord(p % coord0) + + ! Set up base level coordinates + allocate(p % coord0) + p % coord0 % universe = BASE_UNIVERSE + p % coord => p % coord0 + p % in_lower_universe = .false. + end subroutine initialize_particle +!=============================================================================== +! DEALLOCATE_COORD removes all levels of coordinates below a given level. This +! is used in distance_to_boundary when the particle moves from a lower universe +! to a higher universe since the data for the lower one is not needed anymore. +!=============================================================================== + + recursive subroutine deallocate_coord(coord) + + type(LocalCoord), pointer :: coord + + if (associated(coord)) then + ! recursively deallocate lower coordinates + if (associated(coord % next)) call deallocate_coord(coord%next) + + ! deallocate original coordinate + deallocate(coord) + end if + + end subroutine deallocate_coord + end module particle_header diff --git a/src/physics.f90 b/src/physics.f90 index 420d0f4a2c..2b911329d5 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -11,7 +11,7 @@ module physics use global use interpolation, only: interpolate_tab1 use output, only: write_message, print_particle - use particle_header, only: Particle + use particle_header, only: Particle, LocalCoord use random_lcg, only: prn use search, only: binary_search use string, only: int_to_str @@ -29,29 +29,27 @@ contains type(Particle), pointer :: p - integer :: surf ! surface which particle is on - integer :: last_cell ! most recent cell particle was in - integer :: n_event ! number of collisions/crossings - real(8) :: d_boundary ! distance to nearest boundary - real(8) :: d_collision ! sampled distance to collision - real(8) :: distance ! distance particle travels - logical :: found_cell ! found cell which particle is in? - logical :: in_lattice ! is surface crossing in lattice? - type(Universe), pointer :: univ + integer :: surface_crossed ! surface which particle is on + integer :: last_cell ! most recent cell particle was in + integer :: n_event ! number of collisions/crossings + real(8) :: d_boundary ! distance to nearest boundary + real(8) :: d_collision ! sampled distance to collision + real(8) :: distance ! distance particle travels + logical :: found_cell ! found cell which particle is in? + logical :: lattice_crossed ! is surface crossing in lattice? + type(LocalCoord), pointer :: coord - if (p % cell == 0) then - univ => universes(BASE_UNIVERSE) - call find_cell(univ, p, found_cell) - - ! if particle couldn't be located, print error + if (p % coord % cell == NONE) then + call find_cell(p, found_cell) + ! Particle couldn't be located if (.not. found_cell) then write(message, '(A,3ES11.3)') & - "Could not locate cell for particle at: ", p % xyz + "Could not locate cell for particle at: ", p % coord0 % xyz call fatal_error() end if ! set birth cell attribute - p % cell_born = p % cell + p % cell_born = p % coord % cell end if if (verbosity >= 9) then @@ -60,7 +58,8 @@ contains end if if (verbosity >= 10) then - message = " Born in cell " // trim(int_to_str(cells(p%cell)%id)) + message = " Born in cell " // trim(int_to_str(& + cells(p % coord % cell) % id)) call write_message() end if @@ -74,7 +73,7 @@ contains call calculate_xs(p) ! Find the distance to the nearest boundary - call distance_to_boundary(p, d_boundary, surf, in_lattice) + call distance_to_boundary(p, d_boundary, surface_crossed, lattice_crossed) ! Sample a distance to collision d_collision = -log(prn()) / material_xs % total @@ -83,26 +82,37 @@ contains distance = min(d_boundary, d_collision) ! Advance particle - p % xyz = p % xyz + distance * p % uvw - p % xyz_local = p % xyz_local + distance * p % uvw + coord => p % coord0 + do while (associated(coord)) + coord % xyz = coord % xyz + distance * coord % uvw + coord => coord % next + end do if (d_collision > d_boundary) then - last_cell = p % cell - p % cell = 0 - if (in_lattice) then - p % surface = 0 + last_cell = p % coord % cell + p % coord % cell = NONE + if (lattice_crossed) then + p % surface = NONE call cross_lattice(p) else - p % surface = surf + p % surface = surface_crossed call cross_surface(p, last_cell) end if else ! collision - p % surface = 0 + p % surface = NONE call collision(p) ! Save coordinates at collision for tallying purposes - p % last_xyz = p % xyz + p % last_xyz = p % coord % xyz + + ! Set all uvws to base level -- right now, after a collision, only the + ! base level uvws are changed + coord => p % coord0 + do while(associated(coord)) + coord % uvw = p % coord0 % uvw + coord => coord % next + end do end if ! If particle has too many events, display warning and kill it @@ -712,7 +722,7 @@ contains awr = nuc % awr ! Neutron velocity in LAB - v_n = vel * p % uvw + v_n = vel * p % coord0 % uvw ! Sample velocity of target nucleus call sample_target_velocity(p, nuc, v_t) @@ -750,7 +760,7 @@ contains ! Set energy and direction of particle in LAB frame p % E = E - p % uvw = v_n / vel + p % coord0 % uvw = v_n / vel ! Copy scattering cosine for tallies p % mu = mu @@ -902,13 +912,13 @@ contains end if ! copy directional cosines - u = p % uvw(1) - v = p % uvw(2) - w = p % uvw(3) + u = p % coord0 % uvw(1) + v = p % coord0 % uvw(2) + w = p % coord0 % uvw(3) ! change direction of particle call rotate_angle(u, v, w, mu) - p % uvw = (/ u, v, w /) + p % coord0 % uvw = (/ u, v, w /) ! change energy of particle p % E = E @@ -993,9 +1003,9 @@ contains ! determine direction of target velocity based on the neutron's velocity ! vector and the sampled angle between them - u = p % uvw(1) - v = p % uvw(2) - w = p % uvw(3) + u = p % coord0 % uvw(1) + v = p % coord0 % uvw(2) + w = p % coord0 % uvw(3) call rotate_angle(u, v, w, mu) ! determine speed of target nucleus @@ -1093,7 +1103,7 @@ contains do i = int(n_bank,4) + 1, int(min(n_bank + nu, 3*work),4) ! Bank source neutrons by copying particle data fission_bank(i) % id = p % id - fission_bank(i) % xyz = p % xyz + fission_bank(i) % xyz = p % coord0 % xyz ! sample cosine of angle mu = sample_angle(rxn, E) @@ -1250,13 +1260,13 @@ contains end if ! copy directional cosines - u = p % uvw(1) - v = p % uvw(2) - w = p % uvw(3) + u = p % coord0 % uvw(1) + v = p % coord0 % uvw(2) + w = p % coord0 % uvw(3) ! change direction of particle call rotate_angle(u, v, w, mu) - p % uvw = (/ u, v, w /) + p % coord0 % uvw = (/ u, v, w /) ! change energy of particle p % E = E diff --git a/src/plot.f90 b/src/plot.f90 index 41cc72668a..79a08cd70d 100644 --- a/src/plot.f90 +++ b/src/plot.f90 @@ -6,7 +6,7 @@ module plot cross_lattice, cell_contains use geometry_header, only: Universe, BASE_UNIVERSE use global - use particle_header, only: Particle, initialize_particle + use particle_header, only: Particle, initialize_particle, LocalCoord implicit none @@ -18,20 +18,21 @@ contains subroutine run_plot() - integer :: i ! loop index - integer :: surf ! surface which particle is on - integer :: last_cell ! most recent cell particle was in - real(8) :: coord(3) ! starting coordinates - real(8) :: last_x_coord ! bounding x coordinate - real(8) :: last_y_coord ! bounding y coordinate - real(8) :: d ! distance to boundary - real(8) :: distance ! distance particle travels - logical :: found_cell ! found cell which particle is in? - logical :: in_lattice ! is surface crossing in lattice? + integer :: i ! loop index + integer :: surface_crossed ! surface which particle is on + integer :: last_cell ! most recent cell particle was in + real(8) :: xyz(3) ! starting coordinates + real(8) :: last_x_coord ! bounding x coordinate + real(8) :: last_y_coord ! bounding y coordinate + real(8) :: d ! distance to boundary + real(8) :: distance ! distance particle travels + logical :: found_cell ! found cell which particle is in? + logical :: lattice_crossed ! is surface crossing in lattice? character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file - type(Cell), pointer :: c => null() - type(Universe), pointer :: univ => null() - type(Particle), pointer :: p => null() + type(Cell), pointer :: c => null() + type(Universe), pointer :: univ => null() + type(Particle), pointer :: p => null() + type(LocalCoord), pointer :: coord => null() ! Open plot file for binary writing path_plot = trim(path_input) // "plot.out" @@ -44,9 +45,9 @@ contains write(UNIT=UNIT_PLOT) pixel ! Determine coordinates of the upper-left corner of the plot - coord(1) = plot_origin(1) - plot_width(1) / 2.0 - coord(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0 - coord(3) = plot_origin(3) + xyz(1) = plot_origin(1) - plot_width(1) / 2.0 + xyz(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0 + xyz(3) = plot_origin(3) ! Determine bounding x and y coordinates for plot last_x_coord = plot_origin(1) + plot_width(1) / 2.0 @@ -56,76 +57,75 @@ contains allocate(p) ! loop over horizontal rays - do while(coord(2) > last_y_coord) + do while(xyz(2) > last_y_coord) ! initialize the particle and set starting coordinate and direction call initialize_particle(p) - p % xyz = coord - p % xyz_local = coord - p % uvw = (/ 1, 0, 0 /) + + p % coord % xyz = xyz + p % coord % uvw = (/ 1, 0, 0 /) ! write starting coordinate to file - write(UNIT=UNIT_PLOT) p % xyz + write(UNIT=UNIT_PLOT) p % coord % xyz ! Find cell that particle is currently in - univ => universes(BASE_UNIVERSE) - call find_cell(univ, p, found_cell) + call find_cell(p, found_cell) ! ======================================================================= ! MOVE PARTICLE FORWARD TO NEXT CELL - if (.not. found_cell) then - univ => universes(BASE_UNIVERSE) - do i = 1, univ % n_cells - p % xyz = coord - p % xyz_local = coord - p % cell = univ % cells(i) - - distance = INFINITY - call distance_to_boundary(p, d, surf, in_lattice) - if (d < distance) then - ! Move particle forward to next surface - p % xyz = p % xyz + d * p % uvw - - ! Check to make sure particle is actually going into this cell - ! by moving it slightly forward and seeing if the cell contains - ! that coordinate - - p % xyz = p % xyz + 1e-4 * p % uvw - p % xyz_local = p % xyz - - c => cells(p % cell) - if (.not. cell_contains(c, p)) cycle - - ! Reset coordinate to surface crossing - p % xyz = p % xyz - 1e-4 * p % uvw - p % xyz_local = p % xyz - - ! Set new distance and retain pointer to this cell - distance = d - last_cell = p % cell - end if - end do - - ! No cell was found on this horizontal ray - if (distance == INFINITY) then - p % xyz(1) = last_x_coord - p % cell = 0 - write(UNIT_PLOT) p % xyz, p % cell - - ! Move to next horizontal ray - coord(2) = coord(2) - pixel - cycle - end if - - ! Write coordinate where next cell begins - write(UNIT=UNIT_PLOT) p % xyz, 0 - - ! Process surface crossing for next cell - p % cell = 0 - p % surface = -surf - call cross_surface(p, last_cell) - end if +!!$ if (.not. found_cell) then +!!$ univ => universes(BASE_UNIVERSE) +!!$ do i = 1, univ % n_cells +!!$ p % xyz = coord +!!$ p % xyz_local = coord +!!$ p % cell = univ % cells(i) +!!$ +!!$ distance = INFINITY +!!$ ! call distance_to_boundary(p, d, surf, in_lattice) +!!$ if (d < distance) then +!!$ ! Move particle forward to next surface +!!$ p % xyz = p % xyz + d * p % uvw +!!$ +!!$ ! Check to make sure particle is actually going into this cell +!!$ ! by moving it slightly forward and seeing if the cell contains +!!$ ! that coordinate +!!$ +!!$ p % xyz = p % xyz + 1e-4 * p % uvw +!!$ p % xyz_local = p % xyz +!!$ +!!$ c => cells(p % cell) +!!$ if (.not. cell_contains(c, p)) cycle +!!$ +!!$ ! Reset coordinate to surface crossing +!!$ p % xyz = p % xyz - 1e-4 * p % uvw +!!$ p % xyz_local = p % xyz +!!$ +!!$ ! Set new distance and retain pointer to this cell +!!$ distance = d +!!$ last_cell = p % cell +!!$ end if +!!$ end do +!!$ +!!$ ! No cell was found on this horizontal ray +!!$ if (distance == INFINITY) then +!!$ p % xyz(1) = last_x_coord +!!$ p % cell = 0 +!!$ write(UNIT_PLOT) p % xyz, p % cell +!!$ +!!$ ! Move to next horizontal ray +!!$ xyz(2) = xyz(2) - pixel +!!$ cycle +!!$ end if +!!$ +!!$ ! Write coordinate where next cell begins +!!$ write(UNIT=UNIT_PLOT) p % xyz, 0 +!!$ +!!$ ! Process surface crossing for next cell +!!$ p % cell = 0 +!!$ p % surface = -surf +!!$ call cross_surface(p, last_cell) +!!$ end if ! ======================================================================= ! MOVE PARTICLE ACROSS HORIZONTAL TRACK @@ -133,43 +133,46 @@ contains do while (p % alive) ! Calculate distance to next boundary - call distance_to_boundary(p, distance, surf, in_lattice) + call distance_to_boundary(p, distance, surface_crossed, lattice_crossed) ! Advance particle - p%xyz = p%xyz + distance * p%uvw - p%xyz_local = p%xyz_local + distance * p%uvw + coord => p % coord0 + do while (associated(coord)) + coord % xyz = coord % xyz + distance * coord % uvw + coord => coord % next + end do ! If next boundary crossing is out of range of the plot, only include ! the visible portion and move to next horizontal ray - if (p % xyz(1) >= last_x_coord) then + if (p % coord0 % xyz(1) >= last_x_coord) then p % alive = .false. - p % xyz(1) = last_x_coord + p % coord0 % xyz(1) = last_x_coord ! If there is no cell beyond this boundary, mark it as cell 0 - if (distance == INFINITY) p % cell = 0 + if (distance == INFINITY) p % coord % cell = 0 ! Write ending coordinates to file - write(UNIT=UNIT_PLOT) p % xyz, p % cell + write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell cycle end if ! Write boundary crossing coordinates to file - write(UNIT=UNIT_PLOT) p % xyz, p % cell + write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell - last_cell = p % cell - p % cell = 0 - if (in_lattice) then - p % surface = 0 + last_cell = p % coord % cell + p % coord % cell = 0 + if (lattice_crossed) then + p % surface = NONE call cross_lattice(p) else - p % surface = surf + p % surface = surface_crossed call cross_surface(p, last_cell) ! Since boundary conditions are disabled in plotting mode, we need ! to manually add the last segment - if (surfaces(surf) % bc == BC_VACUUM) then - p % xyz(1) = last_x_coord - write(UNIT=UNIT_PLOT) p % xyz, 0 + if (surfaces(surface_crossed) % bc == BC_VACUUM) then + p % coord0 % xyz(1) = last_x_coord + write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0 exit end if end if @@ -177,7 +180,7 @@ contains end do ! Move y-coordinate to next position - coord(2) = coord(2) - pixel + xyz(2) = xyz(2) - pixel end do ! Close plot file diff --git a/src/source.f90 b/src/source.f90 index f916795930..08d3b6195d 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -79,25 +79,24 @@ contains do j = bank_first, bank_last p => source_bank(j - bank_first + 1) + ! set defaults + call initialize_particle(p) + ! initialize random number seed call set_particle_seed(int(j,8)) ! sample position r = (/ (prn(), k = 1,3) /) p % id = j - p % xyz = p_min + r*(p_max - p_min) - p % xyz_local = p % xyz - p % last_xyz = p % xyz + p % coord0 % xyz = p_min + r*(p_max - p_min) + p % last_xyz = p % coord0 % xyz ! sample angle phi = TWO*PI*prn() mu = TWO*prn() - ONE - p % uvw(1) = mu - p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) - - ! set defaults - call initialize_particle(p) + p % coord0 % uvw(1) = mu + p % coord0 % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) + p % coord0 % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) ! sample energy from Watt fission energy spectrum for U-235 do diff --git a/src/tally.f90 b/src/tally.f90 index ca9f0eb268..253e040d7c 100644 --- a/src/tally.f90 +++ b/src/tally.f90 @@ -343,7 +343,7 @@ contains ! determine next universe bin if (t % n_bins(T_UNIVERSE) > 0) then - bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % universe, i) + bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % coord % universe, i) if (bins(T_UNIVERSE) == NO_BIN_FOUND) cycle else bins(T_UNIVERSE) = 1 @@ -359,7 +359,7 @@ contains ! determine next cell bin if (t % n_bins(T_CELL) > 0) then - bins(T_CELL) = get_next_bin(T_CELL, p % cell, i) + bins(T_CELL) = get_next_bin(T_CELL, p % coord % cell, i) if (bins(T_CELL) == NO_BIN_FOUND) cycle else bins(T_CELL) = 1 @@ -386,7 +386,7 @@ contains m => meshes(t % mesh) ! Determine if we're in the mesh first - call get_mesh_bin(m, p % xyz, mesh_bin, in_mesh) + call get_mesh_bin(m, p % coord0 % xyz, mesh_bin, in_mesh) if (.not. in_mesh) cycle bins(T_MESH) = mesh_bin @@ -610,7 +610,7 @@ contains do i = 1, n_tallies ! Copy starting and ending location of particle xyz0 = p % last_xyz - xyz1 = p % xyz + xyz1 = p % coord0 % xyz ! Get pointer to tally t => tallies(i) @@ -631,7 +631,7 @@ contains if (n_cross == 0) cycle ! Copy particle's direction - uvw = p % uvw + uvw = p % coord0 % uvw ! determine incoming energy bin n = t % n_bins(T_ENERGYIN)