From f60f27646cc8758c36aaa09281398a480ce42e21 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Jun 2015 18:34:03 +0700 Subject: [PATCH 01/22] Fix bug in print_particle --- src/output.F90 | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 009e5f7e2..bc309556d 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -296,8 +296,8 @@ contains 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(& - p % coord % lattice_y)) // ')' + coord % lattice_x)) // ',' // trim(to_str(& + coord % lattice_y)) // ')' end if ! Print local coordinates @@ -2181,9 +2181,9 @@ contains end select end function get_label - + !=============================================================================== -! FIND_OFFSET uses a given map number, a target cell ID, and a target offset +! FIND_OFFSET uses a given map number, a target cell ID, and a target offset ! to build a string which is the path from the base universe to the target cell ! with the given offset !=============================================================================== @@ -2196,7 +2196,7 @@ contains integer, intent(in) :: final ! Target offset integer, intent(inout) :: offset ! Current offset character(100) :: path ! Path to offset - + integer :: i, j ! Index over cells integer :: k, l, m ! Indices in lattice integer :: old_k, old_l, old_m ! Previous indices in lattice @@ -2212,7 +2212,7 @@ contains class(Lattice), pointer :: lat ! Pointer to current lattice n = univ % n_cells - + ! Write to the geometry stack if (univ%id == 0) then path = trim(path) // to_str(univ%id) @@ -2223,31 +2223,31 @@ contains ! Look through all cells in this universe do i = 1, n - cell_index = univ % cells(i) + cell_index = univ % cells(i) c => cells(cell_index) - + ! If the cell ID matches the goal and the offset matches final, ! write to the geometry stack if (cell_dict % get_key(c % id) == goal .AND. offset == final) then path = trim(path) // "->" // to_str(c%id) return end if - + end do - + ! Find the fill cell or lattice cell that we need to enter do i = 1, n later_cell = .false. - cell_index = univ % cells(i) + cell_index = univ % cells(i) c => cells(cell_index) - this_cell = .false. + this_cell = .false. ! If we got here, we still think the target is in this universe - ! or further down, but it's not this exact cell. - ! Compare offset to next cell to see if we should enter this cell + ! or further down, but it's not this exact cell. + ! Compare offset to next cell to see if we should enter this cell if (i /= n) then do j = i+1, n @@ -2260,8 +2260,8 @@ contains cycle end if - ! Break loop once we've found the next cell with an offset - exit + ! Break loop once we've found the next cell with an offset + exit end do ! Ensure we didn't just end the loop by iteration @@ -2278,13 +2278,13 @@ contains else lat => lattices(c % fill) % obj temp_offset = lat % offset(map, 1, 1, 1) - end if + end if ! If the final offset is in the range of offset - temp_offset+offset ! then the goal is in this cell if (final < temp_offset + offset) then this_cell = .true. - end if + end if end if end if @@ -2341,7 +2341,7 @@ contains ! Loop over lattice coordinates do k = 1, n_x do l = 1, n_y - do m = 1, n_z + do m = 1, n_z if (final >= lat % offset(map, k, l, m) + offset) then if (k == n_x .and. l == n_y .and. m == n_z) then @@ -2364,14 +2364,14 @@ contains ! Target is at this lattice position lat_offset = lat % offset(map, old_k, old_l, old_m) offset = offset + lat_offset - next_univ => universes(lat % universes(old_k, old_l, old_m)) + next_univ => universes(lat % universes(old_k, old_l, old_m)) path = trim(path) // "(" // trim(to_str(old_k)) // & "," // trim(to_str(old_l)) // "," // & trim(to_str(old_m)) // ")" call find_offset(map, goal, next_univ, final, offset, path) return end if - + end do end do end do @@ -2429,7 +2429,7 @@ contains end if end if - end do + end do end subroutine find_offset end module output From 9c536d9fc5716195d7fe73b13113aac966afab70 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Jun 2015 11:25:47 +0700 Subject: [PATCH 02/22] Static allocation for Particle % coord instead of linked list. Not all tests passing yet. --- src/constants.F90 | 9 +- src/geometry.F90 | 282 +++++++++++++++++++-------------------- src/output.F90 | 30 ++--- src/particle_header.F90 | 58 ++++---- src/particle_restart.F90 | 8 +- src/physics.F90 | 12 +- src/plot.F90 | 105 +++++++-------- src/source.F90 | 12 +- src/tally.F90 | 60 ++++----- src/track_output.F90 | 2 +- src/tracking.F90 | 36 +++-- 11 files changed, 284 insertions(+), 330 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 0c9d15f6a..b6def1b46 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -133,6 +133,9 @@ module constants ! Maximum number of lost particles integer, parameter :: MAX_LOST_PARTICLES = 10 + ! Maximum levels of nested universes + integer, parameter :: MAX_COORD = 5 + ! ============================================================================ ! CROSS SECTION RELATED CONSTANTS @@ -317,13 +320,13 @@ module constants OUT_FRONT = 4, & IN_TOP = 5, & OUT_TOP = 6 - + ! Tally trigger types and threshold integer, parameter :: & VARIANCE = 1, & RELATIVE_ERROR = 2, & - STANDARD_DEVIATION = 3 - + STANDARD_DEVIATION = 3 + ! Global tallY parameters integer, parameter :: N_GLOBAL_TALLIES = 4 integer, parameter :: & diff --git a/src/geometry.F90 b/src/geometry.F90 index ba429566e..72293aeff 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -6,7 +6,7 @@ module geometry &RectLattice, HexLattice use global use output, only: write_message - use particle_header, only: LocalCoord, deallocate_coord, Particle + use particle_header, only: LocalCoord, Particle use particle_restart_write, only: write_particle_restart use string, only: to_str use tally, only: score_surface_current @@ -76,20 +76,18 @@ contains type(Particle), intent(inout) :: p integer :: i ! cell loop index on a level + integer :: j + integer :: n_coord integer :: n ! number of cells to search on a level integer :: index_cell ! index in cells array type(Cell), pointer :: c ! pointer to cell type(Universe), pointer :: univ ! universe to search in - type(LocalCoord), pointer :: coord ! particle coordinate to search on - - coord => p % coord0 ! loop through each coordinate level - do while (associated(coord)) - - p % coord => coord - - univ => universes(coord % universe) + n_coord = p % n_coord + do j = 1, n_coord + p % n_coord = j + univ => universes(p % coord(j) % universe) n = univ % n_cells ! loop through each cell on this level @@ -99,10 +97,10 @@ contains if (simple_cell_contains(c, p)) then ! the particle should only be contained in one cell per level - if (index_cell /= coord % cell) then + if (index_cell /= p % coord(j) % cell) then call fatal_error("Overlapping cells detected: " & &// trim(to_str(cells(index_cell) % id)) // ", " & - &// trim(to_str(cells(coord % cell) % id)) & + &// trim(to_str(cells(p % coord(j) % cell) % id)) & &// " on universe " // trim(to_str(univ % id))) end if @@ -111,9 +109,6 @@ contains end if end do - - coord => coord % next - end do end subroutine check_cell_overlap @@ -130,6 +125,7 @@ contains logical, intent(inout) :: found integer, optional :: search_cells(:) integer :: i ! index over cells + integer :: j integer :: i_xyz(3) ! indices in lattice integer :: n ! number of cells to search integer :: index_cell ! index in cells array @@ -138,8 +134,8 @@ contains class(Lattice), pointer :: lat ! pointer to lattice type(Universe), pointer :: univ ! universe to search in - ! Remove coordinates for any lower levels - call deallocate_coord(p % coord % next) + j = p % n_coord + call p % coord(j+1:) % reset() ! set size of list to search if (present(search_cells)) then @@ -147,7 +143,7 @@ contains n = size(search_cells) else use_search_cells = .false. - univ => universes(p % coord % universe) + univ => universes(p % coord(j) % universe) n = univ % n_cells end if @@ -157,7 +153,7 @@ contains if (use_search_cells) then index_cell = search_cells(i) ! check to make sure search cell is in same universe - if (cells(index_cell) % universe /= p % coord % universe) cycle + if (cells(index_cell) % universe /= p % coord(j) % universe) cycle else index_cell = univ % cells(i) end if @@ -169,7 +165,7 @@ contains if (.not. simple_cell_contains(c, p)) cycle ! Set cell on this level - p % coord % cell = index_cell + p % coord(j) % cell = index_cell ! Show cell information on trace if (verbosity >= 10 .or. trace) then @@ -189,27 +185,28 @@ contains ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL ! Create new level of coordinates - allocate(p % coord % next) - p % coord % next % xyz = p % coord % xyz - p % coord % next % uvw = p % coord % uvw + p % coord(j + 1) % xyz = p % coord(j) % xyz + p % coord(j + 1) % uvw = p % coord(j) % uvw ! Move particle to next level and set universe - p % coord => p % coord % next - p % coord % universe = c % fill + j = j + 1 + p % n_coord = j + p % coord(j) % universe = c % fill ! Apply translation if (allocated(c % translation)) then - p % coord % xyz = p % coord % xyz - c % translation + p % coord(j) % xyz = p % coord(j) % xyz - c % translation end if ! Apply rotation if (allocated(c % rotation_matrix)) then - p % coord % xyz = matmul(c % rotation_matrix, p % coord % xyz) - p % coord % uvw = matmul(c % rotation_matrix, p % coord % uvw) - p % coord % rotated = .true. + p % coord(j) % xyz = matmul(c % rotation_matrix, p % coord(j) % xyz) + p % coord(j) % uvw = matmul(c % rotation_matrix, p % coord(j) % uvw) + p % coord(j) % rotated = .true. end if call find_cell(p, found) + j = p % n_coord if (.not. found) exit elseif (c % type == CELL_LATTICE) then CELL_TYPE @@ -220,40 +217,42 @@ contains lat => lattices(c % fill) % obj ! Determine lattice indices - i_xyz = lat % get_indices(p % coord % xyz + TINY_BIT * p % coord % uvw) + i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw) ! Create new level of coordinates - allocate(p % coord % next) - p % coord % next % xyz = lat % get_local_xyz(p % coord % xyz, i_xyz) - p % coord % next % uvw = p % coord % uvw + call p % coord(j + 1) % reset() + p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz) + p % coord(j + 1) % uvw = p % coord(j) % uvw ! set particle lattice indices - p % coord % next% lattice = c % fill - p % coord % next% lattice_x = i_xyz(1) - p % coord % next% lattice_y = i_xyz(2) - p % coord % next% lattice_z = i_xyz(3) + p % coord(j + 1) % lattice = c % fill + p % coord(j + 1) % lattice_x = i_xyz(1) + p % coord(j + 1) % lattice_y = i_xyz(2) + p % coord(j + 1) % lattice_z = i_xyz(3) ! Set the next lowest coordinate level. if (lat % are_valid_indices(i_xyz)) then ! Particle is inside the lattice. - p % coord % next % universe = & - &lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) + p % coord(j + 1) % universe = & + lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) else ! Particle is outside the lattice. if (lat % outer == NO_OUTER_UNIVERSE) then call fatal_error("A particle is outside latttice " & - &// trim(to_str(lat % id)) // " but the lattice has no & + // trim(to_str(lat % id)) // " but the lattice has no & &defined outer universe.") else - p % coord % next % universe = lat % outer + p % coord(j + 1) % universe = lat % outer end if end if ! Move particle to next level and search for the lower cells. - p % coord => p % coord % next + j = j + 1 + p % n_coord = j call find_cell(p, found) + j = p % n_coord if (.not. found) exit end if CELL_TYPE @@ -314,7 +313,7 @@ contains ! TODO: Find a better solution to score surface currents than ! physically moving the particle forward slightly - p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw call score_surface_current(p) end if @@ -337,7 +336,7 @@ contains ! PARTICLE REFLECTS FROM SURFACE ! Do not handle reflective boundary conditions on lower universes - if (.not. associated(p % coord, p % coord0)) then + if (p % n_coord /= 1) then call handle_lost_particle(p, "Cannot reflect particle " & &// trim(to_str(p % id)) // " off surface in a lower universe.") return @@ -348,15 +347,15 @@ contains ! case the surface crossing in coincident with a mesh boundary if (active_current_tallies % size() > 0) then - p % coord0 % xyz = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw + p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw call score_surface_current(p) - p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw end if ! Copy particle's direction cosines - u = p % coord0 % uvw(1) - v = p % coord0 % uvw(2) - w = p % coord0 % uvw(3) + u = p % coord(1) % uvw(1) + v = p % coord(1) % uvw(2) + w = p % coord(1) % uvw(3) select case (surf%type) case (SURF_PX) @@ -383,8 +382,8 @@ contains case (SURF_CYL_X) ! Find y-y0, z-z0 and dot product of direction and surface normal - y = p % coord0 % xyz(2) - surf % coeffs(1) - z = p % coord0 % xyz(3) - surf % coeffs(2) + y = p % coord(1) % xyz(2) - surf % coeffs(1) + z = p % coord(1) % xyz(3) - surf % coeffs(2) R = surf % coeffs(3) dot_prod = v*y + w*z @@ -394,8 +393,8 @@ contains case (SURF_CYL_Y) ! Find x-x0, z-z0 and dot product of direction and surface normal - x = p % coord0 % xyz(1) - surf % coeffs(1) - z = p % coord0 % xyz(3) - surf % coeffs(2) + x = p % coord(1) % xyz(1) - surf % coeffs(1) + z = p % coord(1) % xyz(3) - surf % coeffs(2) R = surf % coeffs(3) dot_prod = u*x + w*z @@ -405,8 +404,8 @@ contains case (SURF_CYL_Z) ! Find x-x0, y-y0 and dot product of direction and surface normal - x = p % coord0 % xyz(1) - surf % coeffs(1) - y = p % coord0 % xyz(2) - surf % coeffs(2) + x = p % coord(1) % xyz(1) - surf % coeffs(1) + y = p % coord(1) % xyz(2) - surf % coeffs(2) R = surf % coeffs(3) dot_prod = u*x + v*y @@ -417,9 +416,9 @@ contains case (SURF_SPHERE) ! Find x-x0, y-y0, z-z0 and dot product of direction and surface ! normal - x = p % coord0 % xyz(1) - surf % coeffs(1) - y = p % coord0 % xyz(2) - surf % coeffs(2) - z = p % coord0 % xyz(3) - surf % coeffs(3) + x = p % coord(1) % xyz(1) - surf % coeffs(1) + y = p % coord(1) % xyz(2) - surf % coeffs(2) + z = p % coord(1) % xyz(3) - surf % coeffs(3) R = surf % coeffs(4) dot_prod = u*x + v*y + w*z @@ -431,9 +430,9 @@ contains case (SURF_CONE_X) ! Find x-x0, y-y0, z-z0 and dot product of direction and surface ! normal - x = p % coord0 % xyz(1) - surf % coeffs(1) - y = p % coord0 % xyz(2) - surf % coeffs(2) - z = p % coord0 % xyz(3) - surf % coeffs(3) + x = p % coord(1) % xyz(1) - surf % coeffs(1) + y = p % coord(1) % xyz(2) - surf % coeffs(2) + z = p % coord(1) % xyz(3) - surf % coeffs(3) R = surf % coeffs(4) dot_prod = (v*y + w*z - R*u*x)/((R + ONE)*R*x*x) @@ -445,9 +444,9 @@ contains case (SURF_CONE_Y) ! Find x-x0, y-y0, z-z0 and dot product of direction and surface ! normal - x = p % coord0 % xyz(1) - surf % coeffs(1) - y = p % coord0 % xyz(2) - surf % coeffs(2) - z = p % coord0 % xyz(3) - surf % coeffs(3) + x = p % coord(1) % xyz(1) - surf % coeffs(1) + y = p % coord(1) % xyz(2) - surf % coeffs(2) + z = p % coord(1) % xyz(3) - surf % coeffs(3) R = surf % coeffs(4) dot_prod = (u*x + w*z - R*v*y)/((R + ONE)*R*y*y) @@ -459,9 +458,9 @@ contains case (SURF_CONE_Z) ! Find x-x0, y-y0, z-z0 and dot product of direction and surface ! normal - x = p % coord0 % xyz(1) - surf % coeffs(1) - y = p % coord0 % xyz(2) - surf % coeffs(2) - z = p % coord0 % xyz(3) - surf % coeffs(3) + x = p % coord(1) % xyz(1) - surf % coeffs(1) + y = p % coord(1) % xyz(2) - surf % coeffs(2) + z = p % coord(1) % xyz(3) - surf % coeffs(3) R = surf % coeffs(4) dot_prod = (u*x + v*y - R*w*z)/((R + ONE)*R*z*z) @@ -477,18 +476,19 @@ contains ! Set new particle direction norm = sqrt(u*u + v*v + w*w) - p % coord0 % uvw = [u, v, w] / norm + p % coord(1) % uvw = [u, v, w] / norm ! Reassign particle's cell and surface - p % coord0 % cell = last_cell + p % coord(1) % cell = last_cell p % surface = -p % surface ! If a reflective surface is coincident with a lattice or universe ! boundary, it is necessary to redetermine the particle's coordinates in ! the lower universes. - if (associated(p % coord0 % next)) then - call deallocate_coord(p % coord0 % next) + if (p % n_coord > 1) then + p % n_coord = 1 + call p % coord(2:) % reset() call find_cell(p, found) if (.not. found) then call handle_lost_particle(p, "Couldn't find particle after reflecting& @@ -498,7 +498,7 @@ contains end if ! Set previous coordinate going slightly past surface crossing - p % last_xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw + p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw ! Diagnostic message if (verbosity >= 10 .or. trace) then @@ -532,8 +532,8 @@ contains ! Remove lower coordinate levels and assignment of surface p % surface = NONE - p % coord => p % coord0 - call deallocate_coord(p % coord % next) + p % n_coord = 1 + call p % coord(2:) % reset() call find_cell(p, found) if (run_mode /= MODE_PLOTTING .and. (.not. found)) then @@ -542,9 +542,9 @@ contains ! the particle is really traveling tangent to a surface, if we move it ! forward a tiny bit it should fix the problem. - p % coord => p % coord0 - call deallocate_coord(p % coord % next) - p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw + p % n_coord = 1 + call p % coord(2:) % reset() + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw call find_cell(p, found) ! Couldn't find next cell anywhere! This probably means there is an actual @@ -568,52 +568,47 @@ contains type(Particle), intent(inout) :: p integer, intent(in) :: lattice_translation(3) + integer :: j integer :: i_xyz(3) ! indices in lattice logical :: found ! particle found in cell? class(Lattice), pointer :: lat - type(LocalCoord), pointer :: parent_coord - lat => lattices(p % coord % lattice) % obj + j = p % n_coord + lat => lattices(p % coord(j) % lattice) % obj if (verbosity >= 10 .or. trace) then call write_message(" Crossing lattice " // trim(to_str(lat % id)) & - &// ". Current position (" // trim(to_str(p % coord % lattice_x)) & - &// "," // trim(to_str(p % coord % lattice_y)) // "," & - &// trim(to_str(p % coord % lattice_z)) // ")") + &// ". Current position (" // trim(to_str(p % coord(j) % lattice_x)) & + &// "," // trim(to_str(p % coord(j) % lattice_y)) // "," & + &// trim(to_str(p % coord(j) % lattice_z)) // ")") end if - ! Find the coordiante level just above the current one. - parent_coord => p % coord0 - do while(.not. associated(parent_coord % next, p % coord)) - parent_coord => parent_coord % next - end do - ! Set the lattice indices. - p % coord % lattice_x = p % coord % lattice_x + lattice_translation(1) - p % coord % lattice_y = p % coord % lattice_y + lattice_translation(2) - p % coord % lattice_z = p % coord % lattice_z + lattice_translation(3) - i_xyz(1) = p % coord % lattice_x - i_xyz(2) = p % coord % lattice_y - i_xyz(3) = p % coord % lattice_z + p % coord(j) % lattice_x = p % coord(j) % lattice_x + lattice_translation(1) + p % coord(j) % lattice_y = p % coord(j) % lattice_y + lattice_translation(2) + p % coord(j) % lattice_z = p % coord(j) % lattice_z + lattice_translation(3) + i_xyz(1) = p % coord(j) % lattice_x + i_xyz(2) = p % coord(j) % lattice_y + i_xyz(3) = p % coord(j) % lattice_z ! Set the new coordinate position. - p % coord % xyz = lat % get_local_xyz(parent_coord % xyz, i_xyz) + p % coord(j) % xyz = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz) OUTSIDE_LAT: if (.not. lat % are_valid_indices(i_xyz)) then - ! The particle is outside the lattice. Search for it from coord0. - call deallocate_coord(p % coord0 % next) - p % coord => p % coord0 + ! The particle is outside the lattice. Search for it from base coord + p % n_coord = 1 + call p % coord(2:) % reset() call find_cell(p, found) if (.not. found) then call handle_lost_particle(p, "Could not locate particle " & - &// trim(to_str(p % id)) // " after crossing a lattice boundary.") + // trim(to_str(p % id)) // " after crossing a lattice boundary.") return end if else OUTSIDE_LAT ! Find cell in next lattice element - p % coord % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) + p % coord(j) % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) call find_cell(p, found) if (.not. found) then @@ -622,15 +617,15 @@ contains ! off all lower-level coordinates and search from universe zero ! Remove lower coordinates - call deallocate_coord(p % coord0 % next) - p % coord => p % coord0 + p % n_coord = 1 + call p % coord(2:) % reset() ! Search for particle call find_cell(p, found) if (.not. found) then call handle_lost_particle(p, "Could not locate particle " & - &// trim(to_str(p % id)) & - &// " after crossing a lattice boundary.") + // trim(to_str(p % id)) & + // " after crossing a lattice boundary.") return end if end if @@ -644,14 +639,17 @@ contains ! that has a parent cell, also include the surfaces of the edge of the universe. !=============================================================================== - subroutine distance_to_boundary(p, dist, surface_crossed, lattice_translation) + subroutine distance_to_boundary(p, dist, surface_crossed, lattice_translation, & + next_level) type(Particle), intent(inout) :: p real(8), intent(out) :: dist integer, intent(out) :: surface_crossed integer, intent(out) :: lattice_translation(3) + integer, intent(out) :: next_level integer :: i ! index for surface in cell + integer :: j integer :: index_surf ! index in surfaces array (with sign) integer :: i_xyz(3) ! lattice indices integer :: level_surf_cross ! surface crossed on current level @@ -675,30 +673,25 @@ contains type(Cell), pointer :: cl type(Surface), pointer :: surf class(Lattice), pointer :: lat - type(LocalCoord), pointer :: coord - type(LocalCoord), pointer :: final_coord - type(LocalCoord), pointer :: parent_coord ! inialize distance to infinity (huge) dist = INFINITY d_lat = INFINITY d_surf = INFINITY lattice_translation(:) = [0, 0, 0] - nullify(final_coord) - ! Get pointer to top-level coordinates - coord => p % coord0 + next_level = 0 ! Loop over each universe level - LEVEL_LOOP: do while(associated(coord)) + LEVEL_LOOP: do j = 1, p % n_coord ! get pointer to cell on this level - cl => cells(coord % cell) + cl => cells(p % coord(j) % cell) ! copy directional cosines - u = coord % uvw(1) - v = coord % uvw(2) - w = coord % uvw(3) + u = p % coord(j) % uvw(1) + v = p % coord(j) % uvw(2) + w = p % coord(j) % uvw(3) ! ======================================================================= ! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL @@ -706,9 +699,9 @@ contains 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) + x = p % coord(j) % xyz(1) + y = p % coord(j) % xyz(2) + z = p % coord(j) % 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 @@ -1129,16 +1122,16 @@ contains ! ======================================================================= ! FIND MINIMUM DISTANCE TO LATTICE SURFACES - LAT_COORD: if (coord % lattice /= NONE) then - lat => lattices(coord % lattice) % obj + LAT_COORD: if (p % coord(j) % lattice /= NONE) then + lat => lattices(p % coord(j) % lattice) % obj LAT_TYPE: select type(lat) type is (RectLattice) ! copy local coordinates - x = coord % xyz(1) - y = coord % xyz(2) - z = coord % xyz(3) + x = p % coord(j) % xyz(1) + y = p % coord(j) % xyz(2) + z = p % coord(j) % xyz(3) ! determine oncoming edge x0 = sign(lat % pitch(1) * 0.5_8, u) @@ -1202,14 +1195,10 @@ contains type is (HexLattice) LAT_TYPE ! Copy local coordinates. - z = coord % xyz(3) - i_xyz(1) = coord % lattice_x - i_xyz(2) = coord % lattice_y - i_xyz(3) = coord % lattice_z - parent_coord => p % coord0 - do while(.not. associated(parent_coord % next, coord)) - parent_coord => parent_coord % next - end do + z = p % coord(j) % xyz(3) + i_xyz(1) = p % coord(j) % lattice_x + i_xyz(2) = p % coord(j) % lattice_y + i_xyz(3) = p % coord(j) % lattice_z ! Compute velocities along the hexagonal axes. beta_dir = u*sqrt(3.0_8)/2.0_8 + v/2.0_8 @@ -1225,9 +1214,9 @@ contains ! Upper right and lower left sides. edge = -sign(lat % pitch(1)/2.0_8, beta_dir) ! Oncoming edge if (beta_dir > 0.0) then - xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, 0, 0]) + xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, 0, 0]) else - xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 0, 0]) + xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 0, 0]) end if beta = xyz_t(1)*sqrt(3.0_8)/2.0_8 + xyz_t(2)/2.0_8 if (abs(beta - edge) < FP_PRECISION) then @@ -1248,9 +1237,9 @@ contains ! Lower right and upper left sides. edge = -sign(lat % pitch(1)/2.0_8, gama_dir) ! Oncoming edge if (gama_dir > 0.0) then - xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[1, -1, 0]) + xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, -1, 0]) else - xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[-1, 1, 0]) + xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 1, 0]) end if gama = xyz_t(1)*sqrt(3.0_8)/2.0_8 - xyz_t(2)/2.0_8 if (abs(gama - edge) < FP_PRECISION) then @@ -1273,9 +1262,9 @@ contains ! Upper and lower sides. edge = -sign(lat % pitch(1)/2.0_8, v) ! Oncoming edge if (v > 0.0) then - xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, 1, 0]) + xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, 1, 0]) else - xyz_t = lat % get_local_xyz(parent_coord % xyz, i_xyz+[0, -1, 0]) + xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, -1, 0]) end if if (abs(xyz_t(2) - edge) < FP_PRECISION) then d = INFINITY @@ -1333,24 +1322,19 @@ contains dist = d_surf surface_crossed = level_surf_cross lattice_translation(:) = [0, 0, 0] - final_coord => coord + next_level = j end if else if ((dist - d_lat)/dist >= FP_REL_PRECISION) then dist = d_lat surface_crossed = None lattice_translation(:) = level_lat_trans - final_coord => coord + next_level = j end if end if - coord => coord % next - end do LEVEL_LOOP - ! Move particle to appropriate coordinate level - if (associated(final_coord)) p % coord => final_coord - end subroutine distance_to_boundary !=============================================================================== @@ -1365,6 +1349,7 @@ contains type(Surface), pointer :: surf ! surface logical :: s ! sense of particle + integer :: j real(8) :: x,y,z ! coordinates of particle real(8) :: func ! surface function evaluated at point real(8) :: A ! coefficient on x for plane @@ -1374,9 +1359,10 @@ contains real(8) :: x0,y0,z0 ! coefficients for quadratic surfaces / box real(8) :: r ! radius for quadratic surfaces - x = p % coord % xyz(1) - y = p % coord % xyz(2) - z = p % coord % xyz(3) + j = p % n_coord + x = p % coord(j) % xyz(1) + y = p % coord(j) % xyz(2) + z = p % coord(j) % xyz(3) select case (surf % type) case (SURF_PX) @@ -1468,7 +1454,7 @@ contains if (abs(func) < FP_COINCIDENT) then ! Particle may be coincident with this surface. Artifically move the ! particle forward a tiny bit. - p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw + p % coord(j) % xyz = p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw s = sense(p, surf) elseif (func > 0) then s = .true. diff --git a/src/output.F90 b/src/output.F90 index bc309556d..7b5e98564 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -258,7 +258,6 @@ contains type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() class(Lattice), pointer :: l => null() - type(LocalCoord), pointer :: coord => null() ! display type of particle select case (p % type) @@ -273,39 +272,34 @@ contains end select ! loop through each level of universes - coord => p % coord0 - i = 0 - do while(associated(coord)) + do i = 1, p % n_coord ! Print level - write(ou,*) ' Level ' // trim(to_str(i)) + write(ou,*) ' Level ' // trim(to_str(i - 1)) ! Print cell for this level - if (coord % cell /= NONE) then - c => cells(coord % cell) + if (p % coord(i) % cell /= NONE) then + c => cells(p % coord(i) % cell) write(ou,*) ' Cell = ' // trim(to_str(c % id)) end if ! Print universe for this level - if (coord % universe /= NONE) then - u => universes(coord % universe) + if (p % coord(i) % universe /= NONE) then + u => universes(p % coord(i) % universe) write(ou,*) ' Universe = ' // trim(to_str(u % id)) end if ! Print information on lattice - if (coord % lattice /= NONE) then - l => lattices(coord % lattice) % obj + if (p % coord(i) % lattice /= NONE) then + l => lattices(p % coord(i) % lattice) % obj write(ou,*) ' Lattice = ' // trim(to_str(l % id)) write(ou,*) ' Lattice position = (' // trim(to_str(& - coord % lattice_x)) // ',' // trim(to_str(& - coord % lattice_y)) // ')' + p % coord(i) % lattice_x)) // ',' // trim(to_str(& + p % coord(i) % lattice_y)) // ')' end if ! Print local coordinates - write(ou,'(1X,A,3ES12.4)') ' xyz = ', coord % xyz - write(ou,'(1X,A,3ES12.4)') ' uvw = ', coord % uvw - - coord => coord % next - i = i + 1 + write(ou,'(1X,A,3ES12.4)') ' xyz = ', p % coord(i) % xyz + write(ou,'(1X,A,3ES12.4)') ' uvw = ', p % coord(i) % uvw end do ! Print surface diff --git a/src/particle_header.F90 b/src/particle_header.F90 index bf03c5635..30737ef52 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,6 +1,6 @@ module particle_header - use constants, only: NEUTRON, ONE, NONE, ZERO + use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_COORD use geometry_header, only: BASE_UNIVERSE implicit none @@ -26,9 +26,8 @@ module particle_header ! Is this level rotated? logical :: rotated = .false. - - ! Pointer to next (more local) set of coordinates - type(LocalCoord), pointer :: next => null() + contains + procedure :: reset => reset_coord end type LocalCoord !=============================================================================== @@ -42,8 +41,8 @@ module particle_header integer :: type ! Particle type (n, p, e, etc) ! Particle coordinates - type(LocalCoord), pointer :: coord0 => null() ! coordinates on universe 0 - type(LocalCoord), pointer :: coord => null() ! coordinates on lowest universe + integer :: n_coord ! number of current coordinates + type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels ! Other physical data real(8) :: wgt ! particle weight @@ -87,26 +86,6 @@ module particle_header contains -!=============================================================================== -! 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 this coord - deallocate(coord) - end if - - end subroutine deallocate_coord - !=============================================================================== ! INITIALIZE_PARTICLE sets default attributes for a particle from the source ! bank @@ -137,9 +116,8 @@ contains this % fission = .false. ! Set up base level coordinates - allocate(this % coord0) - this % coord0 % universe = BASE_UNIVERSE - this % coord => this % coord0 + this % coord(1) % universe = BASE_UNIVERSE + this % n_coord = 1 end subroutine initialize_particle @@ -152,11 +130,25 @@ contains class(Particle) :: this ! remove any coordinate levels - call deallocate_coord(this % coord0) - - ! Make sure coord pointer is nullified - nullify(this % coord) + call this % coord(:) % reset() end subroutine clear_particle +!=============================================================================== +! RESET_COORD +!=============================================================================== + + elemental subroutine reset_coord(this) + class(LocalCoord), intent(inout) :: this + + this % cell = NONE + this % universe = NONE + this % lattice = NONE + this % lattice_x = NONE + this % lattice_y = NONE + this % lattice_z = NONE + this % rotated = .false. + + end subroutine reset_coord + end module particle_header diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 6ff5d63c2..0ed25fdc8 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -89,13 +89,13 @@ contains call pr % read_data(p % id, 'id') call pr % read_data(p % wgt, 'weight') call pr % read_data(p % E, 'energy') - call pr % read_data(p % coord % xyz, 'xyz', length=3) - call pr % read_data(p % coord % uvw, 'uvw', length=3) + call pr % read_data(p % coord(1) % xyz, 'xyz', length=3) + call pr % read_data(p % coord(1) % uvw, 'uvw', length=3) ! Set particle last attributes p % last_wgt = p % wgt - p % last_xyz = p % coord % xyz - p % last_uvw = p % coord % uvw + p % last_xyz = p % coord(1) % xyz + p % last_uvw = p % coord(1) % uvw p % last_E = p % E ! Close hdf5 file diff --git a/src/physics.F90 b/src/physics.F90 index 4e94ea226..1c8120672 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -37,7 +37,7 @@ contains ! Store pre-collision particle properties p % last_wgt = p % wgt p % last_E = p % E - p % last_uvw = p % coord0 % uvw + p % last_uvw = p % coord(1) % uvw ! Add to collision counter for particle p % n_collision = p % n_collision + 1 @@ -332,7 +332,7 @@ contains if (micro_xs(i_nuclide) % index_sab /= NONE) then ! S(a,b) scattering call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, & - p % E, p % coord0 % uvw, p % mu) + p % E, p % coord(1) % uvw, p % mu) else ! get pointer to elastic scattering reaction @@ -340,7 +340,7 @@ contains ! Perform collision physics for elastic scattering call elastic_scatter(i_nuclide, rxn, & - p % E, p % coord0 % uvw, p % mu, p % wgt) + p % E, p % coord(1) % uvw, p % mu, p % wgt) end if p % event_MT = ELASTIC @@ -382,7 +382,7 @@ contains end do ! Perform collision physics for inelastic scattering - call inelastic_scatter(nuc, rxn, p % E, p % coord0 % uvw, & + call inelastic_scatter(nuc, rxn, p % E, p % coord(1) % uvw, & p % mu, p % wgt) p % event_MT = rxn % MT @@ -1074,7 +1074,7 @@ contains if (ufs) then ! Determine indices on ufs mesh for current location - call get_mesh_indices(ufs_mesh, p % coord0 % xyz, ijk, in_mesh) + call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh) if (.not. in_mesh) then call write_particle_restart(p) call fatal_error("Source site outside UFS mesh!") @@ -1112,7 +1112,7 @@ contains p % fission = .true. ! Fission neutrons will be banked do i = int(n_bank,4) + 1, int(min(n_bank + nu, int(size(fission_bank),8)),4) ! Bank source neutrons by copying particle data - fission_bank(i) % xyz = p % coord0 % xyz + fission_bank(i) % xyz = p % coord(1) % xyz ! Set weight of fission bank site fission_bank(i) % wgt = ONE/weight diff --git a/src/plot.F90 b/src/plot.F90 index d9f6b55e2..1a27fff67 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -7,7 +7,7 @@ module plot use global use mesh, only: get_mesh_indices use output, only: write_message - use particle_header, only: deallocate_coord, Particle, LocalCoord + use particle_header, only: Particle, LocalCoord use plot_header use ppmlib, only: Image, init_image, allocate_image, & deallocate_image, set_pixel @@ -46,7 +46,7 @@ contains end subroutine run_plot !=============================================================================== -! POSITION_RGB computes the red/green/blue values for a given plot with the +! POSITION_RGB computes the red/green/blue values for a given plot with the ! current particle's position !=============================================================================== @@ -56,28 +56,21 @@ contains type(ObjectPlot), pointer, intent(in) :: pl integer, intent(out) :: rgb(3) integer, intent(out) :: id - + + integer :: j logical :: found_cell - integer :: level - type(Cell), pointer :: c => null() - type(LocalCoord), pointer :: coord => null() - - call deallocate_coord(p % coord0 % next) - p % coord => p % coord0 + type(Cell), pointer :: c + + p % n_coord = 1 + call p % coord(2:) % reset() call find_cell(p, found_cell) + j = p % n_coord if (check_overlaps) call check_cell_overlap(p) - ! Loop through universes and stop on any specified level - level = 0 - coord => p % coord0 - do - if (level == pl % level) exit - if (.not. associated(coord % next)) exit - coord => coord % next - level = level + 1 - end do - + ! Set coordinate level if specified + if (pl % level >= 0) j = pl % level + 1 + if (.not. found_cell) then ! If no cell, revert to default color rgb = pl % not_found % rgb @@ -85,7 +78,7 @@ contains else if (pl % color_by == PLOT_COLOR_MATS) then ! Assign color based on material - c => cells(coord % cell) + c => cells(p % coord(j) % cell) if (c % material == MATERIAL_VOID) then ! By default, color void cells white rgb = 255 @@ -100,14 +93,14 @@ contains end if else if (pl % color_by == PLOT_COLOR_CELLS) then ! Assign color based on cell - rgb = pl % colors(coord % cell) % rgb - id = cells(coord % cell) % id + rgb = pl % colors(p % coord(j) % cell) % rgb + id = cells(p % coord(j) % cell) % id else rgb = 0 id = -1 end if end if - + end subroutine position_rgb !=============================================================================== @@ -160,9 +153,9 @@ contains ! allocate and initialize particle call p % initialize() - p % coord % xyz = xyz - p % coord % uvw = (/ 0.5, 0.5, 0.5 /) - p % coord % universe = BASE_UNIVERSE + p % coord(1) % xyz = xyz + p % coord(1) % uvw = (/ 0.5, 0.5, 0.5 /) + p % coord(1) % universe = BASE_UNIVERSE do y = 1, img % height call progress % set_value(dble(y)/dble(img % height)*100.) @@ -175,12 +168,12 @@ contains call set_pixel(img, x-1, y-1, rgb(1), rgb(2), rgb(3)) ! Advance pixel in first direction - p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel + p % coord(1) % xyz(in_i) = p % coord(1) % xyz(in_i) + in_pixel end do ! Advance pixel in second direction - p % coord0 % xyz(in_i) = xyz(in_i) - p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel + p % coord(1) % xyz(in_i) = xyz(in_i) + p % coord(1) % xyz(out_i) = p % coord(1) % xyz(out_i) - out_pixel end do ! Draw tally mesh boundaries on the image if requested @@ -201,10 +194,10 @@ contains ! DRAW_MESH_LINES draws mesh line boundaries on an image !=============================================================================== subroutine draw_mesh_lines(pl, img) - + type(ObjectPlot), pointer, intent(in) :: pl type(Image), intent(inout) :: img - + logical :: in_mesh integer :: out_, in_ ! pixel location integer :: r, g, b ! RGB color for meshlines pixels @@ -221,13 +214,13 @@ contains real(8) :: xyz_ll(3) ! lower left xyz real(8) :: xyz_ur(3) ! upper right xyz type(StructuredMesh), pointer :: m => null() - + m => pl % meshlines_mesh - + r = pl % meshlines_color % rgb(1) g = pl % meshlines_color % rgb(2) b = pl % meshlines_color % rgb(3) - + select case (pl % basis) case(PLOT_BASIS_XY) outer = 1 @@ -259,15 +252,15 @@ contains ! check if we're in the mesh for this ijk if (i > 0 .and. i <= m % dimension(outer) .and. & j > 0 .and. j <= m % dimension(inner)) then - + ! get xyz's of lower left and upper right of this mesh cell xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1) xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1) xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j - + ! map the xyz ranges to pixel ranges - + frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer) outrange(1) = int(frac * real(img % width, 8)) frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer) @@ -295,11 +288,11 @@ contains call set_pixel(img, outrange(2) - plus, in_, r, g, b) end do end do - + end if end do end do - + end subroutine draw_mesh_lines !=============================================================================== @@ -350,7 +343,7 @@ contains subroutine create_3d_dump(pl) type(ObjectPlot), pointer :: pl - + integer :: x, y, z ! voxel location indices integer :: rgb(3) ! colors (red, green, blue) from 0-255 integer :: id ! id of cell or material @@ -361,15 +354,15 @@ contains ! compute voxel widths in each direction vox = pl % width/dble(pl % pixels) - + ! initial particle position ll = pl % origin - pl % width / 2.0 ! allocate and initialize particle call p % initialize() - p % coord0 % xyz = ll - p % coord0 % uvw = (/ 0.5, 0.5, 0.5 /) - p % coord0 % universe = BASE_UNIVERSE + p % coord(1) % xyz = ll + p % coord(1) % uvw = (/ 0.5, 0.5, 0.5 /) + p % coord(1) % universe = BASE_UNIVERSE ! Open binary plot file for writing open(UNIT=UNIT_PLOT, FILE=pl % path_plot, STATUS='replace', & @@ -378,7 +371,7 @@ contains ! write plot header info write(UNIT_PLOT) pl % pixels, vox, ll - ! move to center of voxels + ! move to center of voxels ll = ll + vox / 2.0 do x = 1, pl % pixels(1) @@ -393,21 +386,21 @@ contains write(UNIT_PLOT) id ! advance particle in z direction - p % coord0 % xyz(3) = p % coord0 % xyz(3) + vox(3) - + p % coord(1) % xyz(3) = p % coord(1) % xyz(3) + vox(3) + end do - + ! advance particle in y direction - p % coord0 % xyz(2) = p % coord0 % xyz(2) + vox(2) - p % coord0 % xyz(3) = ll(3) - + p % coord(1) % xyz(2) = p % coord(1) % xyz(2) + vox(2) + p % coord(1) % xyz(3) = ll(3) + end do - + ! advance particle in y direction - p % coord0 % xyz(1) = p % coord0 % xyz(1) + vox(1) - p % coord0 % xyz(2) = ll(2) - p % coord0 % xyz(3) = ll(3) - + p % coord(1) % xyz(1) = p % coord(1) % xyz(1) + vox(1) + p % coord(1) % xyz(2) = ll(2) + p % coord(1) % xyz(3) = ll(3) + end do close(UNIT_PLOT) diff --git a/src/source.F90 b/src/source.F90 index 4f931dc90..2bd56af88 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -132,8 +132,8 @@ contains site % xyz = p_min + r*(p_max - p_min) ! Fill p with needed data - p % coord0 % xyz = site % xyz - p % coord0 % uvw = [ ONE, ZERO, ZERO ] + p % coord(1) % xyz = site % xyz + p % coord(1) % uvw = [ ONE, ZERO, ZERO ] ! Now search to see if location exists in geometry call find_cell(p, found) @@ -161,8 +161,8 @@ contains site % xyz = p_min + r*(p_max - p_min) ! Fill p with needed data - p % coord0 % xyz = site % xyz - p % coord0 % uvw = [ ONE, ZERO, ZERO ] + p % coord(1) % xyz = site % xyz + p % coord(1) % uvw = [ ONE, ZERO, ZERO ] ! Now search to see if location exists in geometry call find_cell(p, found) @@ -306,8 +306,8 @@ contains ! copy attributes from source bank site p % wgt = src % wgt p % last_wgt = src % wgt - p % coord % xyz = src % xyz - p % coord % uvw = src % uvw + p % coord(1) % xyz = src % xyz + p % coord(1) % uvw = src % uvw p % last_xyz = src % xyz p % last_uvw = src % uvw p % E = src % E diff --git a/src/tally.F90 b/src/tally.F90 index 52b05d56c..83e41a9d6 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -11,7 +11,7 @@ module tally mesh_intersects_2d, mesh_intersects_3d use mesh_header, only: StructuredMesh use output, only: header - use particle_header, only: LocalCoord, Particle, deallocate_coord + use particle_header, only: LocalCoord, Particle use search, only: binary_search use string, only: to_str use tally_header, only: TallyResult, TallyMapItem, TallyMapElement @@ -487,7 +487,7 @@ contains if (t % estimator == ESTIMATOR_ANALOG) then uvw = p % last_uvw else if (t % estimator == ESTIMATOR_TRACKLENGTH) then - uvw = p % coord0 % uvw + uvw = p % coord(1) % uvw end if ! Find the order for a collection of requested moments ! and store the moment contribution of each @@ -909,7 +909,6 @@ contains type(TallyObject), pointer :: t type(StructuredMesh), pointer :: m type(Material), pointer :: mat - type(LocalCoord), pointer :: coord t => tallies(i_tally) matching_bins(1:t%n_filters) = 1 @@ -918,8 +917,8 @@ contains ! CHECK IF THIS TRACK INTERSECTS THE MESH ! Copy starting and ending location of particle - xyz0 = p % coord0 % xyz - (d_track - TINY_BIT) * p % coord0 % uvw - xyz1 = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw + xyz0 = p % coord(1) % xyz - (d_track - TINY_BIT) * p % coord(1) % uvw + xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw ! Get index for mesh filter i_filter_mesh = t % find_filter(FILTER_MESH) @@ -940,8 +939,8 @@ contains end if ! Reset starting and ending location - xyz0 = p % coord0 % xyz - d_track * p % coord0 % uvw - xyz1 = p % coord0 % xyz + xyz0 = p % coord(1) % xyz - d_track * p % coord(1) % uvw + xyz1 = p % coord(1) % xyz ! ========================================================================= ! CHECK FOR SCORING COMBINATION FOR FILTERS OTHER THAN MESH @@ -953,7 +952,7 @@ contains ! determine next universe bin ! TODO: Account for multiple universes when performing this filter matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord(p % n_coord) % universe, i_tally) case (FILTER_MATERIAL) matching_bins(i) = get_next_bin(FILTER_MATERIAL, & @@ -961,15 +960,12 @@ contains case (FILTER_CELL) ! determine next cell bin - coord => p % coord0 - do while(associated(coord)) + do j = 1, p % n_coord position(FILTER_CELL) = 0 matching_bins(i) = get_next_bin(FILTER_CELL, & - coord % cell, i_tally) + p % coord(j) % cell, i_tally) if (matching_bins(i) /= NO_BIN_FOUND) exit - coord => coord % next end do - nullify(coord) case (FILTER_CELLBORN) ! determine next cellborn bin @@ -1009,7 +1005,7 @@ contains n_cross = sum(abs(ijk1(:m % n_dimension) - ijk0(:m % n_dimension))) + 1 ! Copy particle's direction - uvw = p % coord0 % uvw + uvw = p % coord(1) % uvw ! Bounding coordinates do j = 1, m % n_dimension @@ -1135,12 +1131,12 @@ contains logical, intent(out) :: found_bin integer :: i ! loop index for filters + integer :: j integer :: n ! number of bins for single filter integer :: offset ! offset for distribcell real(8) :: E ! particle energy type(TallyObject), pointer :: t type(StructuredMesh), pointer :: m - type(LocalCoord), pointer :: coord found_bin = .true. t => tallies(i_tally) @@ -1154,13 +1150,13 @@ contains m => meshes(t % filters(i) % int_bins(1)) ! Determine if we're in the mesh first - call get_mesh_bin(m, p % coord0 % xyz, matching_bins(i)) + call get_mesh_bin(m, p % coord(1) % xyz, matching_bins(i)) case (FILTER_UNIVERSE) ! determine next universe bin ! TODO: Account for multiple universes when performing this filter matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord(p % n_coord) % universe, i_tally) case (FILTER_MATERIAL) if (p % material /= MATERIAL_VOID) then @@ -1170,37 +1166,31 @@ contains case (FILTER_CELL) ! determine next cell bin - coord => p % coord0 - do while(associated(coord)) + do j = 1, p % n_coord position(FILTER_CELL) = 0 matching_bins(i) = get_next_bin(FILTER_CELL, & - coord % cell, i_tally) + p % coord(j) % cell, i_tally) if (matching_bins(i) /= NO_BIN_FOUND) exit - coord => coord % next end do - nullify(coord) case (FILTER_DISTRIBCELL) ! determine next distribcell bin matching_bins(i) = NO_BIN_FOUND - coord => p % coord0 offset = 0 - do while(associated(coord)) - if (cells(coord % cell) % type == CELL_FILL) then - offset = offset + cells(coord % cell) % & + do j = 1, p % n_coord + if (cells(p % coord(j) % cell) % type == CELL_FILL) then + offset = offset + cells(p % coord(j) % cell) % & offset(t % filters(i) % offset) - elseif(cells(coord % cell) % type == CELL_LATTICE) then - offset = offset + lattices(coord % next % lattice) % obj % & - offset(t % filters(i) % offset, coord % next % lattice_x, & - coord % next % lattice_y, coord % next % lattice_z) + elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then + offset = offset + lattices(p % coord(j + 1) % lattice) % obj % & + offset(t % filters(i) % offset, p % coord(j + 1) % lattice_x, & + p % coord(j + 1) % lattice_y, p % coord(j + 1) % lattice_z) end if - if (t % filters(i) % int_bins(1) == coord % cell) then + if (t % filters(i) % int_bins(1) == p % coord(j) % cell) then matching_bins(i) = offset + 1 exit end if - coord => coord % next end do - nullify(coord) case (FILTER_CELLBORN) ! determine next cellborn bin @@ -1296,7 +1286,7 @@ contains TALLY_LOOP: do i = 1, active_current_tallies % size() ! Copy starting and ending location of particle xyz0 = p % last_xyz - xyz1 = p % coord0 % xyz + xyz1 = p % coord(1) % xyz ! Get pointer to tally i_tally = active_current_tallies % get_item(i) @@ -1328,7 +1318,7 @@ contains end if ! Copy particle's direction - uvw = p % coord0 % uvw + uvw = p % coord(1) % uvw ! determine incoming energy bin j = t % find_filter(FILTER_ENERGYIN) diff --git a/src/track_output.F90 b/src/track_output.F90 index 60ab7e968..ad3a8d1ee 100644 --- a/src/track_output.F90 +++ b/src/track_output.F90 @@ -46,7 +46,7 @@ contains end if ! Write current coordinates into the newest column. - coords(:, n_tracks) = p % coord0 % xyz + coords(:, n_tracks) = p % coord(1) % xyz end subroutine write_particle_track !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index c14d85b5d..f7c58bffb 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -28,6 +28,8 @@ contains type(Particle), intent(inout) :: p + integer :: j ! coordinate level + integer :: next_level ! next coordinate level to check integer :: surface_crossed ! surface which particle is on integer :: lattice_translation(3) ! in-lattice translation vector integer :: last_cell ! most recent cell particle was in @@ -36,7 +38,6 @@ contains real(8) :: d_collision ! sampled distance to collision real(8) :: distance ! distance particle travels logical :: found_cell ! found cell which particle is in? - type(LocalCoord), pointer :: coord ! Display message if high verbosity or trace is on if (verbosity >= 9 .or. trace) then @@ -45,7 +46,7 @@ contains ! If the cell hasn't been determined based on the particle's location, ! initiate a search for the current cell - if (p % coord % cell == NONE) then + if (p % coord(p % n_coord) % cell == NONE) then call find_cell(p, found_cell) ! Particle couldn't be located @@ -54,7 +55,7 @@ contains end if ! set birth cell attribute - p % cell_born = p % coord % cell + p % cell_born = p % coord(p % n_coord) % cell end if ! Initialize number of events to zero @@ -87,7 +88,7 @@ contains ! Find the distance to the nearest boundary call distance_to_boundary(p, d_boundary, surface_crossed, & - &lattice_translation) + lattice_translation, next_level) ! Sample a distance to collision if (material_xs % total == ZERO) then @@ -100,10 +101,8 @@ contains distance = min(d_boundary, d_collision) ! Advance particle - coord => p % coord0 - do while (associated(coord)) - coord % xyz = coord % xyz + distance * coord % uvw - coord => coord % next + do j = 1, p % n_coord + p % coord(j) % xyz = p % coord(j) % xyz + distance * p % coord(j) % uvw end do ! Score track-length tallies @@ -120,8 +119,9 @@ contains ! ==================================================================== ! PARTICLE CROSSES SURFACE - last_cell = p % coord % cell - p % coord % cell = NONE + if (next_level > 0) p % n_coord = next_level + last_cell = p % coord(p % n_coord) % cell + p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then ! Particle crosses lattice boundary p % surface = NONE @@ -168,7 +168,7 @@ contains p % fission = .false. ! Save coordinates for tallying purposes - p % last_xyz = p % coord0 % xyz + p % last_xyz = p % coord(1) % xyz ! Set last material to none since cross sections will need to be ! re-evaluated @@ -176,19 +176,15 @@ contains ! 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 % next)) - if (coord % next % rotated) then + do j = 1, p % n_coord - 1 + if (p % coord(j + 1) % rotated) then ! If next level is rotated, apply rotation matrix - coord % next % uvw = matmul(cells(coord % cell) % & - rotation_matrix, coord % uvw) + p % coord(j + 1) % uvw = matmul(cells(p % coord(j) % cell) % & + rotation_matrix, p % coord(j) % uvw) else ! Otherwise, copy this level's direction - coord % next % uvw = coord % uvw + p % coord(j + 1) % uvw = p % coord(j) % uvw end if - - ! Advance coordinate level - coord => coord % next end do end if From deb7ce9ce79230103fd7564ca073491e925cee3b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Jun 2015 20:02:06 +0700 Subject: [PATCH 03/22] Fix handling of reflective boundary conditions --- src/geometry.F90 | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 72293aeff..dc8c4adc3 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -220,7 +220,6 @@ contains i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw) ! Create new level of coordinates - call p % coord(j + 1) % reset() p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz) p % coord(j + 1) % uvw = p % coord(j) % uvw @@ -486,15 +485,13 @@ contains ! boundary, it is necessary to redetermine the particle's coordinates in ! the lower universes. - if (p % n_coord > 1) then - p % n_coord = 1 - call p % coord(2:) % reset() - call find_cell(p, found) - if (.not. found) then - call handle_lost_particle(p, "Couldn't find particle after reflecting& - & from surface.") - return - end if + p % n_coord = 1 + call p % coord(2:) % reset() + call find_cell(p, found) + if (.not. found) then + call handle_lost_particle(p, "Couldn't find particle after reflecting& + & from surface.") + return end if ! Set previous coordinate going slightly past surface crossing From 22f04ce1ce59e0caf43b1287689a8cb9b47db55c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Jun 2015 20:03:42 +0700 Subject: [PATCH 04/22] Update results for test_filter_universe --- tests/test_filter_universe/results_true.dat | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_filter_universe/results_true.dat b/tests/test_filter_universe/results_true.dat index e840d6aac..904b053ca 100644 --- a/tests/test_filter_universe/results_true.dat +++ b/tests/test_filter_universe/results_true.dat @@ -1,11 +1,11 @@ k-combined: 1.093844E+00 1.626801E-02 tallies: -6.144371E+01 -7.795909E+02 -7.330533E+00 -1.137816E+01 -4.850651E+01 -4.933509E+02 +6.369043E+01 +8.385422E+02 +7.330762E+00 +1.137874E+01 +5.011385E+01 +5.251276E+02 5.132453E+00 5.801019E+00 From 825f037dbad82f04fb17c01d066a91828f43efb9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Jun 2015 15:51:18 +0700 Subject: [PATCH 05/22] Call LocalCoord % reset() in gfortran-4.6 compatible manner --- src/geometry.F90 | 9 +++------ src/particle_header.F90 | 5 ++++- src/plot.F90 | 1 - 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index dc8c4adc3..6f6b5f20b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -134,8 +134,10 @@ contains class(Lattice), pointer :: lat ! pointer to lattice type(Universe), pointer :: univ ! universe to search in + do j = p % n_coord + 1, MAX_COORD + call p % coord(j) % reset() + end do j = p % n_coord - call p % coord(j+1:) % reset() ! set size of list to search if (present(search_cells)) then @@ -486,7 +488,6 @@ contains ! the lower universes. p % n_coord = 1 - call p % coord(2:) % reset() call find_cell(p, found) if (.not. found) then call handle_lost_particle(p, "Couldn't find particle after reflecting& @@ -530,7 +531,6 @@ contains ! Remove lower coordinate levels and assignment of surface p % surface = NONE p % n_coord = 1 - call p % coord(2:) % reset() call find_cell(p, found) if (run_mode /= MODE_PLOTTING .and. (.not. found)) then @@ -540,7 +540,6 @@ contains ! forward a tiny bit it should fix the problem. p % n_coord = 1 - call p % coord(2:) % reset() p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw call find_cell(p, found) @@ -594,7 +593,6 @@ contains OUTSIDE_LAT: if (.not. lat % are_valid_indices(i_xyz)) then ! The particle is outside the lattice. Search for it from base coord p % n_coord = 1 - call p % coord(2:) % reset() call find_cell(p, found) if (.not. found) then call handle_lost_particle(p, "Could not locate particle " & @@ -615,7 +613,6 @@ contains ! Remove lower coordinates p % n_coord = 1 - call p % coord(2:) % reset() ! Search for particle call find_cell(p, found) diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 30737ef52..1aabc2570 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -128,9 +128,12 @@ contains subroutine clear_particle(this) class(Particle) :: this + integer :: i ! remove any coordinate levels - call this % coord(:) % reset() + do i = 1, MAX_COORD + call this % coord(i) % reset() + end do end subroutine clear_particle diff --git a/src/plot.F90 b/src/plot.F90 index 1a27fff67..f761b4040 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -62,7 +62,6 @@ contains type(Cell), pointer :: c p % n_coord = 1 - call p % coord(2:) % reset() call find_cell(p, found_cell) j = p % n_coord From 52b31d8ac9ac4b3f6127c07c5dd55dcf2aab7e98 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 3 Jun 2015 16:42:55 +0700 Subject: [PATCH 06/22] Use omp reduction instead of atomic on global tallies --- src/eigenvalue.F90 | 26 +++++++++++++++++++++++++- src/geometry.F90 | 4 +--- src/global.F90 | 4 ++++ src/physics.F90 | 8 ++------ src/tracking.F90 | 8 ++------ 5 files changed, 34 insertions(+), 16 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 3b7f149a8..edf3e2985 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -75,7 +75,11 @@ contains ! ==================================================================== ! LOOP OVER PARTICLES -!$omp parallel do schedule(static) firstprivate(p) +!$omp parallel do schedule(static) firstprivate(p) & +!$omp reduction(+:tally_tracklength) & +!$omp reduction(+:tally_collision) & +!$omp reduction(+:tally_leakage) & +!$omp reduction(+:tally_absorption) PARTICLE_LOOP: do i_work = 1, work current_work = i_work @@ -171,6 +175,22 @@ contains subroutine finalize_generation() + ! Update global tallies with the omp private accumulation variables + global_tallies(K_TRACKLENGTH) % value = & + global_tallies(K_TRACKLENGTH) % value + tally_tracklength + global_tallies(K_COLLISION) % value = & + global_tallies(K_COLLISION) % value + tally_collision + global_tallies(LEAKAGE) % value = & + global_tallies(LEAKAGE) % value + tally_leakage + global_tallies(K_ABSORPTION) % value = & + global_tallies(K_ABSORPTION) % value + tally_absorption + + ! reset private tallies + tally_tracklength = 0 + tally_collision = 0 + tally_leakage = 0 + tally_absorption = 0 + #ifdef _OPENMP ! Join the fission bank from each thread into one global fission bank call join_bank_from_threads() @@ -210,6 +230,10 @@ contains if (.not. active_batches) then call reset_result(global_tallies) n_realizations = 0 + tally_tracklength = 0 + tally_collision = 0 + tally_leakage = 0 + tally_absorption = 0 end if ! Perform CMFD calculation if on diff --git a/src/geometry.F90 b/src/geometry.F90 index 6f6b5f20b..faf4cd1d0 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -320,9 +320,7 @@ contains ! Score to global leakage tally if (tallies_on) then -!$omp atomic - global_tallies(LEAKAGE) % value = & - global_tallies(LEAKAGE) % value + p % wgt + tally_leakage = tally_leakage + p % wgt end if ! Display message diff --git a/src/global.F90 b/src/global.F90 index 2ec111dc0..86e638e62 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -116,6 +116,10 @@ module global ! 3) leakage fraction type(TallyResult), allocatable, target :: global_tallies(:) + real(8) :: tally_tracklength = 0 + real(8) :: tally_collision = 0 + real(8) :: tally_leakage = 0 + real(8) :: tally_absorption = 0 ! Tally map structure type(TallyMap), allocatable :: tally_maps(:) diff --git a/src/physics.F90 b/src/physics.F90 index 1c8120672..c85f102b6 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -253,18 +253,14 @@ contains p % last_wgt = p % wgt ! Score implicit absorption estimate of keff -!$omp atomic - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + p % absorb_wgt * & + tally_absorption = tally_absorption + p % absorb_wgt * & micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption else ! See if disappearance reaction happens if (micro_xs(i_nuclide) % absorption > & prn() * micro_xs(i_nuclide) % total) then ! Score absorption estimate of keff -!$omp atomic - global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + p % wgt * & + tally_absorption = tally_absorption + p % wgt * & micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption p % alive = .false. diff --git a/src/tracking.F90 b/src/tracking.F90 index f7c58bffb..59c31c6b1 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -110,9 +110,7 @@ contains call score_tracklength_tally(p, distance) ! Score track-length estimate of k-eff -!$omp atomic - global_tallies(K_TRACKLENGTH) % value = & - global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * & + tally_tracklength = tally_tracklength + p % wgt * distance * & material_xs % nu_fission if (d_collision > d_boundary) then @@ -138,9 +136,7 @@ contains ! PARTICLE HAS COLLISION ! Score collision estimate of keff -!$omp atomic - global_tallies(K_COLLISION) % value = & - global_tallies(K_COLLISION) % value + p % wgt * & + tally_collision = tally_collision + p % wgt * & material_xs % nu_fission / material_xs % total ! score surface current tallies -- this has to be done before the collision From 64d3ccba9d566e6d2d3c5f8ec6a34bcedec37a51 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 12 Jun 2015 00:15:02 +0700 Subject: [PATCH 07/22] Rather than use reduction clause in !$omp do (which doesn't work in gfortran or ifort), just accumulate global tallies by hand in a critical section. --- src/eigenvalue.F90 | 20 ++++++++------------ src/global.F90 | 2 ++ 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index edf3e2985..ceb65f724 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -75,11 +75,7 @@ contains ! ==================================================================== ! LOOP OVER PARTICLES -!$omp parallel do schedule(static) firstprivate(p) & -!$omp reduction(+:tally_tracklength) & -!$omp reduction(+:tally_collision) & -!$omp reduction(+:tally_leakage) & -!$omp reduction(+:tally_absorption) +!$omp parallel do schedule(static) firstprivate(p) PARTICLE_LOOP: do i_work = 1, work current_work = i_work @@ -100,7 +96,7 @@ contains end do GENERATION_LOOP call finalize_batch() - + if (satisfy_triggers) exit BATCH_LOOP end do BATCH_LOOP @@ -176,6 +172,8 @@ contains subroutine finalize_generation() ! Update global tallies with the omp private accumulation variables +!$omp parallel +!$omp critical global_tallies(K_TRACKLENGTH) % value = & global_tallies(K_TRACKLENGTH) % value + tally_tracklength global_tallies(K_COLLISION) % value = & @@ -184,12 +182,14 @@ contains global_tallies(LEAKAGE) % value + tally_leakage global_tallies(K_ABSORPTION) % value = & global_tallies(K_ABSORPTION) % value + tally_absorption +!$omp end critical ! reset private tallies tally_tracklength = 0 tally_collision = 0 tally_leakage = 0 tally_absorption = 0 +!$omp end parallel #ifdef _OPENMP ! Join the fission bank from each thread into one global fission bank @@ -230,10 +230,6 @@ contains if (.not. active_batches) then call reset_result(global_tallies) n_realizations = 0 - tally_tracklength = 0 - tally_collision = 0 - tally_leakage = 0 - tally_absorption = 0 end if ! Perform CMFD calculation if on @@ -244,12 +240,12 @@ contains ! Calculate combined estimate of k-effective if (master) call calculate_combined_keff() - + ! Check_triggers if (master) call check_triggers() #ifdef MPI call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, & - MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) #endif if (satisfy_triggers .or. & (trigger_on .and. current_batch == n_max_batches)) then diff --git a/src/global.F90 b/src/global.F90 index 86e638e62..4f57d3006 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -120,6 +120,8 @@ module global real(8) :: tally_collision = 0 real(8) :: tally_leakage = 0 real(8) :: tally_absorption = 0 +!$omp threadprivate(tally_tracklength, tally_collision, tally_leakage, & +!$omp& tally_absorption) ! Tally map structure type(TallyMap), allocatable :: tally_maps(:) From 2bdd579ba535cf134e44845d38e6a8191509f7f2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 12 Jun 2015 00:37:39 +0700 Subject: [PATCH 08/22] Only accumulate k-effective global tallies if eigenvalue calculation --- src/physics.F90 | 12 ++++++++---- src/tracking.F90 | 13 +++++++++---- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index c85f102b6..ae0a5fdd1 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -253,15 +253,19 @@ contains p % last_wgt = p % wgt ! Score implicit absorption estimate of keff - tally_absorption = tally_absorption + p % absorb_wgt * & - micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption + if (run_mode == MODE_EIGENVALUE) then + tally_absorption = tally_absorption + p % absorb_wgt * & + micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption + end if else ! See if disappearance reaction happens if (micro_xs(i_nuclide) % absorption > & prn() * micro_xs(i_nuclide) % total) then ! Score absorption estimate of keff - tally_absorption = tally_absorption + p % wgt * & - micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption + if (run_mode == MODE_EIGENVALUE) then + tally_absorption = tally_absorption + p % wgt * & + micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption + end if p % alive = .false. p % event = EVENT_ABSORB diff --git a/src/tracking.F90 b/src/tracking.F90 index 59c31c6b1..b40faa314 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -1,5 +1,6 @@ module tracking + use constants, only: MODE_EIGENVALUE use cross_section, only: calculate_xs use error, only: fatal_error, warning use geometry, only: find_cell, distance_to_boundary, cross_surface, & @@ -110,8 +111,10 @@ contains call score_tracklength_tally(p, distance) ! Score track-length estimate of k-eff - tally_tracklength = tally_tracklength + p % wgt * distance * & - material_xs % nu_fission + if (run_mode == MODE_EIGENVALUE) then + tally_tracklength = tally_tracklength + p % wgt * distance * & + material_xs % nu_fission + end if if (d_collision > d_boundary) then ! ==================================================================== @@ -136,8 +139,10 @@ contains ! PARTICLE HAS COLLISION ! Score collision estimate of keff - tally_collision = tally_collision + p % wgt * & - material_xs % nu_fission / material_xs % total + if (run_mode == MODE_EIGENVALUE) then + tally_collision = tally_collision + p % wgt * & + material_xs % nu_fission / material_xs % total + end if ! score surface current tallies -- this has to be done before the collision ! since the direction of the particle will change and we need to use the From f698e9d31ca4b6bb61d65028d291952f87f2d121 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Jul 2015 20:27:09 +0700 Subject: [PATCH 09/22] Rename threadprivate global tally variables and add a few comments --- src/eigenvalue.F90 | 16 ++++++++-------- src/geometry.F90 | 2 +- src/global.F90 | 23 +++++++++++++++-------- src/physics.F90 | 4 ++-- src/tracking.F90 | 6 +++--- 5 files changed, 29 insertions(+), 22 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index ceb65f724..97198d5f4 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -175,20 +175,20 @@ contains !$omp parallel !$omp critical global_tallies(K_TRACKLENGTH) % value = & - global_tallies(K_TRACKLENGTH) % value + tally_tracklength + global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength global_tallies(K_COLLISION) % value = & - global_tallies(K_COLLISION) % value + tally_collision + global_tallies(K_COLLISION) % value + global_tally_collision global_tallies(LEAKAGE) % value = & - global_tallies(LEAKAGE) % value + tally_leakage + global_tallies(LEAKAGE) % value + global_tally_leakage global_tallies(K_ABSORPTION) % value = & - global_tallies(K_ABSORPTION) % value + tally_absorption + global_tallies(K_ABSORPTION) % value + global_tally_absorption !$omp end critical ! reset private tallies - tally_tracklength = 0 - tally_collision = 0 - tally_leakage = 0 - tally_absorption = 0 + global_tally_tracklength = 0 + global_tally_collision = 0 + global_tally_leakage = 0 + global_tally_absorption = 0 !$omp end parallel #ifdef _OPENMP diff --git a/src/geometry.F90 b/src/geometry.F90 index faf4cd1d0..8f73314cb 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -320,7 +320,7 @@ contains ! Score to global leakage tally if (tallies_on) then - tally_leakage = tally_leakage + p % wgt + global_tally_leakage = global_tally_leakage + p % wgt end if ! Display message diff --git a/src/global.F90 b/src/global.F90 index 4f57d3006..ef183e620 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -112,16 +112,23 @@ module global ! Global tallies ! 1) collision estimate of k-eff - ! 2) track-length estimate of k-eff - ! 3) leakage fraction + ! 2) absorption estimate of k-eff + ! 3) track-length estimate of k-eff + ! 4) leakage fraction type(TallyResult), allocatable, target :: global_tallies(:) - real(8) :: tally_tracklength = 0 - real(8) :: tally_collision = 0 - real(8) :: tally_leakage = 0 - real(8) :: tally_absorption = 0 -!$omp threadprivate(tally_tracklength, tally_collision, tally_leakage, & -!$omp& tally_absorption) + + ! It is possible to protect accumulate operations on global tallies by using + ! an atomic update. However, when multiple threads accumulate to the same + ! global tally, it can cause a higher cache miss rate due to + ! invalidation. Thus, we use threadprivate variables to accumulate global + ! tallies and then reduce at the end of a generation. + real(8) :: global_tally_collision = ZERO + real(8) :: global_tally_absorption = ZERO + real(8) :: global_tally_tracklength = ZERO + real(8) :: global_tally_leakage = ZERO +!$omp threadprivate(global_tally_collision, global_tally_absorption, & +!$omp& global_tally_tracklength, global_tally_leakage) ! Tally map structure type(TallyMap), allocatable :: tally_maps(:) diff --git a/src/physics.F90 b/src/physics.F90 index ae0a5fdd1..706ef0954 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -254,7 +254,7 @@ contains ! Score implicit absorption estimate of keff if (run_mode == MODE_EIGENVALUE) then - tally_absorption = tally_absorption + p % absorb_wgt * & + global_tally_absorption = global_tally_absorption + p % absorb_wgt * & micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption end if else @@ -263,7 +263,7 @@ contains prn() * micro_xs(i_nuclide) % total) then ! Score absorption estimate of keff if (run_mode == MODE_EIGENVALUE) then - tally_absorption = tally_absorption + p % wgt * & + global_tally_absorption = global_tally_absorption + p % wgt * & micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption end if diff --git a/src/tracking.F90 b/src/tracking.F90 index b40faa314..e5faeb209 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -112,8 +112,8 @@ contains ! Score track-length estimate of k-eff if (run_mode == MODE_EIGENVALUE) then - tally_tracklength = tally_tracklength + p % wgt * distance * & - material_xs % nu_fission + global_tally_tracklength = global_tally_tracklength + p % wgt * & + distance * material_xs % nu_fission end if if (d_collision > d_boundary) then @@ -140,7 +140,7 @@ contains ! Score collision estimate of keff if (run_mode == MODE_EIGENVALUE) then - tally_collision = tally_collision + p % wgt * & + global_tally_collision = global_tally_collision + p % wgt * & material_xs % nu_fission / material_xs % total end if From d0694bb865481845d43b033c42157a6dbddf3c1d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 5 Jul 2015 08:37:46 +0700 Subject: [PATCH 10/22] Increase maximum number of nested coordinate levels to 10 --- src/constants.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constants.F90 b/src/constants.F90 index b6def1b46..38a0d23e1 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -134,7 +134,7 @@ module constants integer, parameter :: MAX_LOST_PARTICLES = 10 ! Maximum levels of nested universes - integer, parameter :: MAX_COORD = 5 + integer, parameter :: MAX_COORD = 10 ! ============================================================================ ! CROSS SECTION RELATED CONSTANTS From c321b135ca149db3c46042c92563e278981254bd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 Jul 2015 07:45:21 +0700 Subject: [PATCH 11/22] Set maximum number of coordinates from CMakeLists.txt --- CMakeLists.txt | 4 ++++ src/constants.F90 | 3 --- src/particle_header.F90 | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e335dfd57..f2f04452b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,6 +30,10 @@ if (verbose) set(CMAKE_VERBOSE_MAKEFILE on) endif() +# Maximum number of nested coordinates levels +set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels") +add_definitions(-DMAX_COORD=${maxcoord}) + #=============================================================================== # MPI for distributed-memory parallelism / HDF5 for binary output #=============================================================================== diff --git a/src/constants.F90 b/src/constants.F90 index 38a0d23e1..b62cf08b0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -133,9 +133,6 @@ module constants ! Maximum number of lost particles integer, parameter :: MAX_LOST_PARTICLES = 10 - ! Maximum levels of nested universes - integer, parameter :: MAX_COORD = 10 - ! ============================================================================ ! CROSS SECTION RELATED CONSTANTS diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 1aabc2570..e78a89513 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,6 +1,6 @@ module particle_header - use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_COORD + use constants, only: NEUTRON, ONE, NONE, ZERO use geometry_header, only: BASE_UNIVERSE implicit none From 9060c42e5a9960cb373cb51179fee55688221a92 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 5 Jul 2015 18:52:13 -0600 Subject: [PATCH 12/22] Allow >5 spaces continuation lines in style guide --- docs/source/devguide/styleguide.rst | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index 19e44ec7f..25c58ae77 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -125,8 +125,16 @@ program, subroutine, function, if, associate, etc. Emacs users should set the variables f90-if-indent, f90-do-indent, f90-continuation-indent, f90-type-indent, f90-associate-indent, and f90-program indent to 2. -Continuation lines should be indented by an extra 5 spaces. This is the default -value of f90-continuation-indent in Emacs. +Continuation lines should be indented by at least 5 spaces. They may be indented +more in order to make the content match the context. For example, either of +these are valid continuation indentations: + +.. code-block:: fortran + + local_xyz(1) = xyz(1) - (this % lower_left(1) + & + (i_xyz(1) - HALF)*this % pitch(1)) + call which_data(scatt_type, get_scatt, get_nuscatt, get_chi_t, get_chi_p, & + get_chi_d, scatt_order) Whitespace in Expressions ------------------------- From 8c54fe6f334a83e18fa262d642459ed7304c68e7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 Jul 2015 08:12:34 +0700 Subject: [PATCH 13/22] List coverage and maxcoord CMake options in documentation. --- CMakeLists.txt | 2 +- docs/source/usersguide/install.rst | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f2f04452b..84aeb2e3f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,7 +24,7 @@ option(profile "Compile with profiling flags" OFF) option(debug "Compile with debug flags" OFF) option(optimize "Turn on all compiler optimization flags" OFF) option(verbose "Create verbose Makefiles" OFF) -option(coverage "Compile with flags" OFF) +option(coverage "Compile with coverage analysis flags" OFF) if (verbose) set(CMAKE_VERBOSE_MAKEFILE on) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 848c6ecaa..8d98a5c72 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -160,6 +160,13 @@ openmp Enables shared-memory parallelism using the OpenMP API. The Fortran compiler being used must support OpenMP. +coverage + Compile and link code instrumented for coverage analysis. This is typically + used in conjunction with gcov. + +maxcoord + Maximum number of nested coordinate levels in geometry. Defaults to 10. + To set any of these options (e.g. turning on debug mode), the following form should be used: From 1a0969e9ca7fe9d3b3a259b296cb4c4e707e2b99 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 Jul 2015 08:53:49 +0700 Subject: [PATCH 14/22] Check for too many nested coordinate levels. --- src/geometry.F90 | 77 ++++++++++++++++++++++++++++++++++++++++++++++ src/initialize.F90 | 66 ++++++++++++++++++++++----------------- 2 files changed, 115 insertions(+), 28 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 433747a0e..79f1dd3ee 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -1885,5 +1885,82 @@ contains end subroutine count_instance +!=============================================================================== +! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in +! the geometry +!=============================================================================== + + recursive function maximum_levels(univ) result(levels) + + type(Universe), intent(in) :: univ ! universe to search through + integer :: levels ! maximum number of levels for this universe + + integer :: i ! index over cells + integer :: j, k, m ! indices in lattice + integer :: levels_below ! max levels below this universe + type(Cell), pointer :: c ! pointer to current cell + type(Universe), pointer :: next_univ ! next universe to loop through + class(Lattice), pointer :: lat ! pointer to current lattice + + levels_below = 0 + do i = 1, univ % n_cells + c => cells(univ % cells(i)) + + ! ==================================================================== + ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL + if (c % type == CELL_FILL) then + + next_univ => universes(c % fill) + levels_below = maximum_levels(next_univ) + + ! ==================================================================== + ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL + elseif (c % type == CELL_LATTICE) then + + ! Set current lattice + lat => lattices(c % fill) % obj + + select type (lat) + + type is (RectLattice) + + ! Loop over lattice coordinates + do j = 1, lat % n_cells(1) + do k = 1, lat % n_cells(2) + do m = 1, lat % n_cells(3) + next_univ => universes(lat % universes(j, k, m)) + levels_below = max(levels_below, maximum_levels(next_univ)) + end do + end do + end do + + type is (HexLattice) + + ! Loop over lattice coordinates + do m = 1, lat % n_axial + do k = 1, 2*lat % n_rings - 1 + do j = 1, 2*lat % n_rings - 1 + ! This array location is never used + if (j + k < lat % n_rings + 1) then + cycle + ! This array location is never used + else if (j + k > 3*lat % n_rings - 1) then + cycle + else + next_univ => universes(lat % universes(j, k, m)) + levels_below = max(levels_below, maximum_levels(next_univ)) + end if + end do + end do + end do + + end select + + end if + end do + + levels = 1 + levels_below + + end function maximum_levels end module geometry diff --git a/src/initialize.F90 b/src/initialize.F90 index 4c3cc8147..518a285a4 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -7,7 +7,8 @@ module initialize use set_header, only: SetInt use energy_grid, only: logarithmic_grid, grid_method, unionized_grid use error, only: fatal_error, warning - use geometry, only: neighbor_lists, count_instance, calc_offsets + use geometry, only: neighbor_lists, count_instance, calc_offsets, & + maximum_levels use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& &BASE_UNIVERSE use global @@ -95,11 +96,20 @@ contains ! Initialize distribcell_filters call prepare_distribcell() - + ! After reading input and basic geometry setup is complete, build lists of ! neighboring cells for efficient tracking call neighbor_lists() + ! Check to make sure there are not too many nested coordinate levels in the + ! geometry since the coordinate list is statically allocated for performance + ! reasons + if (maximum_levels(universes(BASE_UNIVERSE)) > MAX_COORD) then + call fatal_error("Too many nested coordinate levels in the geometry. & + &Try increasing the maximum number of coordinate levels by & + &providing the CMake -Dmaxcoord= option.") + end if + if (run_mode /= MODE_PLOTTING) then ! With the AWRs from the xs_listings, change all material specifications ! so that they contain atom percents summing to 1 @@ -936,7 +946,7 @@ contains count_all = .false. - ! Loop over tallies + ! Loop over tallies do i = 1, n_tallies ! Get pointer to tally @@ -954,25 +964,25 @@ contains if (size(tally % filters(j) % int_bins) > 1) then call fatal_error("A distribcell filter was specified with & &multiple bins. This feature is not supported.") - end if + end if end if end do end do - + if (count_all) then - + univ => universes(BASE_UNIVERSE) ! sum the number of occurrences of all cells call count_instance(univ) - ! Loop over tallies - do i = 1, n_tallies + ! Loop over tallies + do i = 1, n_tallies ! Get pointer to tally - tally => tallies(i) + tally => tallies(i) ! Initialize the filters do j = 1, tally % n_filters @@ -993,7 +1003,7 @@ contains ! Calculate offsets for each target distribcell do i = 1, n_maps - do j = 1, n_universes + do j = 1, n_universes univ => universes(j) call calc_offsets(univ_list(i), i, univ, counts, found) end do @@ -1003,7 +1013,7 @@ contains deallocate(counts) deallocate(found) deallocate(univ_list) - + end subroutine prepare_distribcell !=============================================================================== @@ -1018,31 +1028,31 @@ contains logical, intent(out), allocatable :: found(:,:) ! Target found integer :: i, j, k, l, m ! Loop counters - type(SetInt) :: cell_list ! distribells to track + type(SetInt) :: cell_list ! distribells to track type(Universe), pointer :: univ ! pointer to universe class(Lattice), pointer :: lat ! pointer to lattice type(TallyObject), pointer :: tally ! pointer to tally type(TallyFilter), pointer :: filter ! pointer to filter - + ! Begin gathering list of cells in distribcell tallies n_maps = 0 - + ! Populate list of distribcells to track do i = 1, n_tallies tally => tallies(i) do j = 1, tally % n_filters - filter => tally % filters(j) + filter => tally % filters(j) if (filter % type == FILTER_DISTRIBCELL) then if (.not. cell_list % contains(filter % int_bins(1))) then call cell_list % add(filter % int_bins(1)) end if - end if + end if end do end do - + ! Compute the number of unique universes containing these distribcells ! to determine the number of offset tables to allocate do i = 1, n_universes @@ -1053,7 +1063,7 @@ contains end if end do end do - + ! Allocate the list of offset tables for each unique universe allocate(univ_list(n_maps)) @@ -1071,34 +1081,34 @@ contains univ => universes(i) do j = 1, univ % n_cells - + if (cell_list % contains(univ % cells(j))) then - - ! Loop over all tallies + + ! Loop over all tallies do l = 1, n_tallies tally => tallies(l) - + do m = 1, tally % n_filters filter => tally % filters(m) - + ! Loop over only distribcell filters ! If filter points to cell we just found, set offset index - if (filter % type == FILTER_DISTRIBCELL) then + if (filter % type == FILTER_DISTRIBCELL) then if (filter % int_bins(1) == univ % cells(j)) then filter % offset = k end if end if end do - end do - + end do + univ_list(k) = univ % id k = k + 1 end if end do end do - - ! Allocate the offset tables for lattices + + ! Allocate the offset tables for lattices do i = 1, n_lattices lat => lattices(i) % obj From a7815ad913b2e3c8ae2e654daf9f07b1e5820245 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 Jul 2015 10:28:03 +0700 Subject: [PATCH 15/22] Fix bug in maximum_levels() pointed out by @wbinventor and respond to his other comments as well. --- docs/source/usersguide/install.rst | 4 +++- src/geometry.F90 | 12 ++++++------ src/initialize.F90 | 1 + src/particle_header.F90 | 2 +- 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 8d98a5c72..7d3cda173 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -162,7 +162,7 @@ openmp coverage Compile and link code instrumented for coverage analysis. This is typically - used in conjunction with gcov. + used in conjunction with gcov_. maxcoord Maximum number of nested coordinate levels in geometry. Defaults to 10. @@ -174,6 +174,8 @@ should be used: cmake -Ddebug=on /path/to/openmc +.. _gcov: https://gcc.gnu.org/onlinedocs/gcc/Gcov.html + Compiling with MPI ++++++++++++++++++ diff --git a/src/geometry.F90 b/src/geometry.F90 index 79f1dd3ee..39d95817b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -76,8 +76,8 @@ contains type(Particle), intent(inout) :: p integer :: i ! cell loop index on a level - integer :: j - integer :: n_coord + integer :: j ! coordinate level index + integer :: n_coord ! saved number of coordinate levels integer :: n ! number of cells to search on a level integer :: index_cell ! index in cells array type(Cell), pointer :: c ! pointer to cell @@ -125,7 +125,7 @@ contains logical, intent(inout) :: found integer, optional :: search_cells(:) integer :: i ! index over cells - integer :: j + integer :: j ! coordinate level index integer :: i_xyz(3) ! indices in lattice integer :: n ! number of cells to search integer :: index_cell ! index in cells array @@ -186,7 +186,7 @@ contains ! ====================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - ! Create new level of coordinates + ! Store lower level coordinates p % coord(j + 1) % xyz = p % coord(j) % xyz p % coord(j + 1) % uvw = p % coord(j) % uvw @@ -221,7 +221,7 @@ contains ! Determine lattice indices i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw) - ! Create new level of coordinates + ! Store lower level coordinates p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz) p % coord(j + 1) % uvw = p % coord(j) % uvw @@ -1911,7 +1911,7 @@ contains if (c % type == CELL_FILL) then next_univ => universes(c % fill) - levels_below = maximum_levels(next_univ) + levels_below = max(levels_below, maximum_levels(next_univ)) ! ==================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL diff --git a/src/initialize.F90 b/src/initialize.F90 index 518a285a4..599923713 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -104,6 +104,7 @@ contains ! Check to make sure there are not too many nested coordinate levels in the ! geometry since the coordinate list is statically allocated for performance ! reasons + print *, 'Max levels = ', maximum_levels(universes(BASE_UNIVERSE)) if (maximum_levels(universes(BASE_UNIVERSE)) > MAX_COORD) then call fatal_error("Too many nested coordinate levels in the geometry. & &Try increasing the maximum number of coordinate levels by & diff --git a/src/particle_header.F90 b/src/particle_header.F90 index e78a89513..7c2586c3e 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -41,7 +41,7 @@ module particle_header integer :: type ! Particle type (n, p, e, etc) ! Particle coordinates - integer :: n_coord ! number of current coordinates + integer :: n_coord ! number of current coordinates type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels ! Other physical data From 46625bbc5dbfad68a21aaa5b4f1dab54fd0b2142 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 Jul 2015 10:41:01 +0700 Subject: [PATCH 16/22] Remove print statement that was added accidentally --- src/initialize.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 599923713..518a285a4 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -104,7 +104,6 @@ contains ! Check to make sure there are not too many nested coordinate levels in the ! geometry since the coordinate list is statically allocated for performance ! reasons - print *, 'Max levels = ', maximum_levels(universes(BASE_UNIVERSE)) if (maximum_levels(universes(BASE_UNIVERSE)) > MAX_COORD) then call fatal_error("Too many nested coordinate levels in the geometry. & &Try increasing the maximum number of coordinate levels by & From 57a82f80dc7e882373479b22c039bd9803565900 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 8 Jul 2015 21:52:58 -0600 Subject: [PATCH 17/22] Fix Python API length checking for hex lat pitch --- openmc/geometry.py | 2 +- openmc/universe.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/geometry.py b/openmc/geometry.py index f1c038403..2ea9705fa 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -205,7 +205,7 @@ class GeometryFile(object): sort_xml_elements(self._geometry_file) clean_xml_indentation(self._geometry_file) - # Write the XML Tree to the materials.xml file + # Write the XML Tree to the geometry.xml file tree = ET.ElementTree(self._geometry_file) tree.write("geometry.xml", xml_declaration=True, encoding='utf-8', method="xml") diff --git a/openmc/universe.py b/openmc/universe.py index 62055049d..6d9bbccb6 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -1059,7 +1059,7 @@ class HexLattice(Lattice): @Lattice.pitch.setter def pitch(self, pitch): check_type('lattice pitch', pitch, Iterable, Real) - check_length('lattice pitch', pitch, 2, 3) + check_length('lattice pitch', pitch, 1, 2) for dim in pitch: check_greater_than('lattice pitch', dim, 0) self._pitch = pitch From 10f4615d3018cb809fc7abb8799d103eb6a26c89 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 8 Jul 2015 23:32:54 -0600 Subject: [PATCH 18/22] Make PyAPI export lattice outer universes --- openmc/universe.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/universe.py b/openmc/universe.py index 6d9bbccb6..bd3b4726c 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -935,6 +935,7 @@ class RectLattice(Lattice): if self._outer is not None: outer = ET.SubElement(lattice_subelement, "outer") outer.text = '{0}'.format(self._outer._id) + self._outer.create_xml_subelement(xml_element) # Export Lattice cell dimensions dimension = ET.SubElement(lattice_subelement, "dimension") @@ -1197,6 +1198,7 @@ class HexLattice(Lattice): if self._outer is not None: outer = ET.SubElement(lattice_subelement, "outer") outer.text = '{0}'.format(self._outer._id) + self._outer.create_xml_subelement(xml_element) lattice_subelement.set("n_rings", str(self._num_rings)) From af3ff7b9e371adb7062df4bcec779258d12283c8 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 13 Jul 2015 20:43:24 -0600 Subject: [PATCH 19/22] Fix indexing bug for HexLattice in Python API --- openmc/universe.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index bd3b4726c..84ce3696e 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -1238,8 +1238,8 @@ class HexLattice(Lattice): universe.create_xml_subelement(xml_element) # Initialize the remaining universes. - for r in range(self._num_rings-1): - for theta in range(2*(self._num_rings - r)): + for r in range(self._num_rings - 1): + for theta in range(6*(self._num_rings - 1 - r)): universe = self._universes[r][theta] universe.create_xml_subelement(xml_element) From 42886182fd929afaf371431c9fe42945d213c883 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 Jul 2015 10:47:46 +0200 Subject: [PATCH 20/22] Transition to container-based Travis CI infrastructure --- .travis.yml | 21 ++++++++++----- tests/travis_install.sh | 58 +++++++++++++++++++++++------------------ 2 files changed, 47 insertions(+), 32 deletions(-) diff --git a/.travis.yml b/.travis.yml index 21d110b40..7af617f85 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,11 +1,21 @@ +sudo: false language: python python: - "2.7" - "3.4" +addons: + apt: + packages: + - gfortran + - g++ +cache: + directories: + - $HOME/mpich_install + - $HOME/hdf5_install + - $HOME/phdf5_install before_install: # ============== Handle Python third-party packages ============== - - sudo apt-get update - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then wget https://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh -O miniconda.sh; else @@ -20,13 +30,12 @@ before_install: - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy h5py - source activate test-environment - # ============== Install GCC, MPICH, HDF5, PHDF5 ============== - - sudo apt-get install -qq -y gfortran g++ + # Install GCC, MPICH, HDF5, PHDF5 - ./tests/travis_install.sh - export FC=gfortran - - export MPI_DIR=$PWD/mpich_install - - export PHDF5_DIR=$PWD/phdf5_install - - export HDF5_DIR=$PWD/hdf5_install + - export MPI_DIR=$HOME/mpich_install + - export PHDF5_DIR=$HOME/phdf5_install + - export HDF5_DIR=$HOME/hdf5_install install: true diff --git a/tests/travis_install.sh b/tests/travis_install.sh index 4976ad092..2d21532de 100755 --- a/tests/travis_install.sh +++ b/tests/travis_install.sh @@ -1,38 +1,44 @@ -#!/bin/sh +#!/bin/bash set -ev # Build MPICH and HDF5 for rest of debug tests -if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then +if [[ "$TRAVIS_PULL_REQUEST" != "false" ]]; then # Build MPICH - wget -q http://www.mpich.org/static/downloads/3.1.3/mpich-3.1.3.tar.gz - tar -xzvf mpich-3.1.3.tar.gz >/dev/null 2>&1 - cd mpich-3.1.3 - ./configure --prefix=$PWD/../mpich_install -q - make -j >/dev/null 2>&1 - make install >/dev/null 2>&1 - cd .. + if [[ ! -d $HOME/mpich_install ]]; then + wget -q http://www.mpich.org/static/downloads/3.1.3/mpich-3.1.3.tar.gz + tar -xzvf mpich-3.1.3.tar.gz >/dev/null 2>&1 + cd mpich-3.1.3 + ./configure --prefix=$HOME/mpich_install -q + make >/dev/null 2>&1 + make install >/dev/null 2>&1 + cd .. + fi # Build PHDF5 - wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz - tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1 - mv hdf5-1.8.15 phdf5-1.8.15; cd phdf5-1.8.15 - CC=$PWD/../mpich_install/bin/mpicc FC=$PWD/../mpich_install/bin/mpif90 \ - ./configure \ - --prefix=$PWD/../phdf5_install -q --enable-fortran \ - --enable-fortran2003 --enable-parallel - make -j >/dev/null 2>&1 - make install >/dev/null 2>&1 - cd .. + if [[ ! -d $HOME/phdf5_install ]]; then + wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz + tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1 + mv hdf5-1.8.15 phdf5-1.8.15; cd phdf5-1.8.15 + CC=$HOME/mpich_install/bin/mpicc FC=$HOME/mpich_install/bin/mpif90 \ + ./configure \ + --prefix=$HOME/phdf5_install -q --enable-fortran \ + --enable-fortran2003 --enable-parallel + make >/dev/null 2>&1 + make install >/dev/null 2>&1 + cd .. + fi # Build HDF5 - tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1 - cd hdf5-1.8.15 - CC=gcc FC=gfortran ./configure --prefix=$PWD/../hdf5_install -q \ - --enable-fortran --enable-fortran2003 - make -j >/dev/null 2>&1 - make install >/dev/null 2>&1 - cd .. + if [[ ! -d $HOME/hdf5_install ]]; then + tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1 + cd hdf5-1.8.15 + CC=gcc FC=gfortran ./configure --prefix=$HOME/hdf5_install -q \ + --enable-fortran --enable-fortran2003 + make -j >/dev/null 2>&1 + make install >/dev/null 2>&1 + cd .. + fi fi From 6a8db6e458f94ebd7bea91a23fca27fd3a8f1581 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 14 Jul 2015 22:47:12 -0600 Subject: [PATCH 21/22] Fix distribcell lattice outer error --- src/tally.F90 | 14 ++++++++++--- .../case-1/geometry.xml | 10 ++++++---- .../case-1/results_true.dat | 20 +++++++++---------- 3 files changed, 27 insertions(+), 17 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 83e41a9d6..4f0692877 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1182,9 +1182,17 @@ contains offset = offset + cells(p % coord(j) % cell) % & offset(t % filters(i) % offset) elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then - offset = offset + lattices(p % coord(j + 1) % lattice) % obj % & - offset(t % filters(i) % offset, p % coord(j + 1) % lattice_x, & - p % coord(j + 1) % lattice_y, p % coord(j + 1) % lattice_z) + if (lattices(p % coord(j + 1) % lattice) % obj & + % are_valid_indices([& + p % coord(j + 1) % lattice_x, & + p % coord(j + 1) % lattice_y, & + p % coord(j + 1) % lattice_z])) then + offset = offset + lattices(p % coord(j + 1) % lattice) % obj % & + offset(t % filters(i) % offset, & + p % coord(j + 1) % lattice_x, & + p % coord(j + 1) % lattice_y, & + p % coord(j + 1) % lattice_z) + end if end if if (t % filters(i) % int_bins(1) == p % coord(j) % cell) then matching_bins(i) = offset + 1 diff --git a/tests/test_filter_distribcell/case-1/geometry.xml b/tests/test_filter_distribcell/case-1/geometry.xml index ddbd23afe..507559615 100644 --- a/tests/test_filter_distribcell/case-1/geometry.xml +++ b/tests/test_filter_distribcell/case-1/geometry.xml @@ -4,21 +4,23 @@ + 2 2 -2.0 -2.0 2.0 2.0 + 3 2 2 2 2 - - - - + + + + diff --git a/tests/test_filter_distribcell/case-1/results_true.dat b/tests/test_filter_distribcell/case-1/results_true.dat index fb68a38f9..49bf3ed4e 100644 --- a/tests/test_filter_distribcell/case-1/results_true.dat +++ b/tests/test_filter_distribcell/case-1/results_true.dat @@ -1,14 +1,14 @@ k-combined: 0.000000E+00 0.000000E+00 tally 1: -8.565980E-03 -7.337601E-05 -8.045879E-03 -6.473617E-05 -9.027116E-03 -8.148882E-05 -7.326285E-03 -5.367445E-05 +1.440759E-02 +2.075788E-04 +1.222930E-02 +1.495558E-04 +1.407292E-02 +1.980471E-04 +1.034365E-02 +1.069911E-04 tally 2: -3.296526E-02 -1.086708E-03 +5.105347E-02 +2.606457E-03 From 9205e5f526c2232a9119d841cec7d700e8345b13 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 15 Jul 2015 07:29:17 +0200 Subject: [PATCH 22/22] Build/load dependencies for pushes and pull requests --- tests/travis_install.sh | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/tests/travis_install.sh b/tests/travis_install.sh index 2d21532de..c01a980c0 100755 --- a/tests/travis_install.sh +++ b/tests/travis_install.sh @@ -2,11 +2,8 @@ set -ev -# Build MPICH and HDF5 for rest of debug tests -if [[ "$TRAVIS_PULL_REQUEST" != "false" ]]; then - - # Build MPICH - if [[ ! -d $HOME/mpich_install ]]; then +# Build MPICH +if [[ ! -e $HOME/mpich_install/bin/mpiexec ]]; then wget -q http://www.mpich.org/static/downloads/3.1.3/mpich-3.1.3.tar.gz tar -xzvf mpich-3.1.3.tar.gz >/dev/null 2>&1 cd mpich-3.1.3 @@ -14,10 +11,10 @@ if [[ "$TRAVIS_PULL_REQUEST" != "false" ]]; then make >/dev/null 2>&1 make install >/dev/null 2>&1 cd .. - fi +fi - # Build PHDF5 - if [[ ! -d $HOME/phdf5_install ]]; then +# Build PHDF5 +if [[ ! -e $HOME/phdf5_install/bin/h5pfc ]]; then wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1 mv hdf5-1.8.15 phdf5-1.8.15; cd phdf5-1.8.15 @@ -28,10 +25,10 @@ if [[ "$TRAVIS_PULL_REQUEST" != "false" ]]; then make >/dev/null 2>&1 make install >/dev/null 2>&1 cd .. - fi +fi - # Build HDF5 - if [[ ! -d $HOME/hdf5_install ]]; then +# Build HDF5 +if [[ ! -e $HOME/hdf5_install/bin/h5fc ]]; then tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1 cd hdf5-1.8.15 CC=gcc FC=gfortran ./configure --prefix=$HOME/hdf5_install -q \ @@ -39,6 +36,4 @@ if [[ "$TRAVIS_PULL_REQUEST" != "false" ]]; then make -j >/dev/null 2>&1 make install >/dev/null 2>&1 cd .. - fi - fi