mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Merge pull request #408 from paulromano/improve-threading
Improve OpenMP thread scaling
This commit is contained in:
Will Boyd
commit
0405597701
16 changed files with 466 additions and 387 deletions
|
|
@ -24,12 +24,16 @@ 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)
|
||||
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
|
||||
#===============================================================================
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
||||
|
|
@ -167,6 +174,8 @@ should be used:
|
|||
|
||||
cmake -Ddebug=on /path/to/openmc
|
||||
|
||||
.. _gcov: https://gcc.gnu.org/onlinedocs/gcc/Gcov.html
|
||||
|
||||
Compiling with MPI
|
||||
++++++++++++++++++
|
||||
|
||||
|
|
|
|||
|
|
@ -171,6 +171,26 @@ 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 + global_tally_tracklength
|
||||
global_tallies(K_COLLISION) % value = &
|
||||
global_tallies(K_COLLISION) % value + global_tally_collision
|
||||
global_tallies(LEAKAGE) % value = &
|
||||
global_tallies(LEAKAGE) % value + global_tally_leakage
|
||||
global_tallies(K_ABSORPTION) % value = &
|
||||
global_tallies(K_ABSORPTION) % value + global_tally_absorption
|
||||
!$omp end critical
|
||||
|
||||
! reset private tallies
|
||||
global_tally_tracklength = 0
|
||||
global_tally_collision = 0
|
||||
global_tally_leakage = 0
|
||||
global_tally_absorption = 0
|
||||
!$omp end parallel
|
||||
|
||||
#ifdef _OPENMP
|
||||
! Join the fission bank from each thread into one global fission bank
|
||||
call join_bank_from_threads()
|
||||
|
|
|
|||
373
src/geometry.F90
373
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 ! 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
|
||||
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 ! coordinate level index
|
||||
integer :: i_xyz(3) ! indices in lattice
|
||||
integer :: n ! number of cells to search
|
||||
integer :: index_cell ! index in cells array
|
||||
|
|
@ -138,8 +134,10 @@ 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)
|
||||
do j = p % n_coord + 1, MAX_COORD
|
||||
call p % coord(j) % reset()
|
||||
end do
|
||||
j = p % n_coord
|
||||
|
||||
! set size of list to search
|
||||
if (present(search_cells)) then
|
||||
|
|
@ -147,7 +145,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 +155,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 +167,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
|
||||
|
|
@ -188,28 +186,29 @@ 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
|
||||
! Store lower level coordinates
|
||||
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 +219,41 @@ 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
|
||||
! 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
|
||||
|
||||
! 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,15 +314,13 @@ 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
|
||||
|
||||
! Score to global leakage tally
|
||||
if (tallies_on) then
|
||||
!$omp atomic
|
||||
global_tallies(LEAKAGE) % value = &
|
||||
global_tallies(LEAKAGE) % value + p % wgt
|
||||
global_tally_leakage = global_tally_leakage + p % wgt
|
||||
end if
|
||||
|
||||
! Display message
|
||||
|
|
@ -337,7 +335,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 +346,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 +381,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 +392,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 +403,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 +415,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 +429,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 +443,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 +457,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,28 +475,26 @@ 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)
|
||||
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 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
|
||||
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 +528,7 @@ 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 find_cell(p, found)
|
||||
|
||||
if (run_mode /= MODE_PLOTTING .and. (.not. found)) then
|
||||
|
|
@ -542,9 +537,8 @@ 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
|
||||
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 +562,46 @@ 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 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 +610,14 @@ 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
|
||||
|
||||
! 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 +631,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 +665,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 +691,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 +1114,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) * HALF, u)
|
||||
|
|
@ -1202,14 +1187,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(THREE)/TWO + v/TWO
|
||||
|
|
@ -1225,9 +1206,9 @@ contains
|
|||
! Upper right and lower left sides.
|
||||
edge = -sign(lat % pitch(1)/TWO, beta_dir) ! Oncoming edge
|
||||
if (beta_dir > ZERO) 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(THREE)/TWO + xyz_t(2)/TWO
|
||||
if (abs(beta - edge) < FP_PRECISION) then
|
||||
|
|
@ -1248,9 +1229,9 @@ contains
|
|||
! Lower right and upper left sides.
|
||||
edge = -sign(lat % pitch(1)/TWO, gama_dir) ! Oncoming edge
|
||||
if (gama_dir > ZERO) 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(THREE)/TWO - xyz_t(2)/TWO
|
||||
if (abs(gama - edge) < FP_PRECISION) then
|
||||
|
|
@ -1273,9 +1254,9 @@ contains
|
|||
! Upper and lower sides.
|
||||
edge = -sign(lat % pitch(1)/TWO, v) ! Oncoming edge
|
||||
if (v > ZERO) 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 +1314,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 +1341,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 +1351,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 +1446,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.
|
||||
|
|
@ -1907,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 = max(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
|
||||
|
|
|
|||
|
|
@ -112,11 +112,24 @@ 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(:)
|
||||
|
||||
! 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(:)
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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(&
|
||||
p % coord % lattice_x)) // ',' // trim(to_str(&
|
||||
p % 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
|
||||
|
|
@ -2181,9 +2175,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 +2190,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 +2206,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 +2217,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 +2254,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 +2272,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 +2335,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 +2358,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 +2423,7 @@ contains
|
|||
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end subroutine find_offset
|
||||
|
||||
end module output
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
@ -150,13 +128,30 @@ contains
|
|||
subroutine clear_particle(this)
|
||||
|
||||
class(Particle) :: this
|
||||
integer :: i
|
||||
|
||||
! remove any coordinate levels
|
||||
call deallocate_coord(this % coord0)
|
||||
|
||||
! Make sure coord pointer is nullified
|
||||
nullify(this % coord)
|
||||
do i = 1, MAX_COORD
|
||||
call this % coord(i) % reset()
|
||||
end do
|
||||
|
||||
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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -253,19 +253,19 @@ 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 * &
|
||||
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
global_tally_absorption = global_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
|
||||
!$omp atomic
|
||||
global_tallies(K_ABSORPTION) % value = &
|
||||
global_tallies(K_ABSORPTION) % value + p % wgt * &
|
||||
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
global_tally_absorption = global_tally_absorption + p % wgt * &
|
||||
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
|
||||
end if
|
||||
|
||||
p % alive = .false.
|
||||
p % event = EVENT_ABSORB
|
||||
|
|
@ -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
|
||||
|
|
|
|||
58
src/plot.F90
58
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
|
||||
|
|
@ -57,26 +57,18 @@ contains
|
|||
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()
|
||||
type(Cell), pointer :: c
|
||||
|
||||
call deallocate_coord(p % coord0 % next)
|
||||
p % coord => p % coord0
|
||||
p % n_coord = 1
|
||||
|
||||
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
|
||||
|
|
@ -85,7 +77,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,8 +92,8 @@ 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
|
||||
|
|
@ -160,9 +152,9 @@ contains
|
|||
|
||||
! allocate and initialize particle
|
||||
call p % initialize()
|
||||
p % coord % xyz = xyz
|
||||
p % coord % uvw = [ HALF, HALF, HALF ]
|
||||
p % coord % universe = BASE_UNIVERSE
|
||||
p % coord(1) % xyz = xyz
|
||||
p % coord(1) % uvw = [ HALF, HALF, HALF ]
|
||||
p % coord(1) % universe = BASE_UNIVERSE
|
||||
|
||||
do y = 1, img % height
|
||||
call progress % set_value(dble(y)/dble(img % height)*100)
|
||||
|
|
@ -175,12 +167,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
|
||||
|
|
@ -367,9 +359,9 @@ contains
|
|||
|
||||
! allocate and initialize particle
|
||||
call p % initialize()
|
||||
p % coord0 % xyz = ll
|
||||
p % coord0 % uvw = [ HALF, HALF, HALF ]
|
||||
p % coord0 % universe = BASE_UNIVERSE
|
||||
p % coord(1) % xyz = ll
|
||||
p % coord(1) % uvw = [ HALF, HALF, HALF ]
|
||||
p % coord(1) % universe = BASE_UNIVERSE
|
||||
|
||||
! Open binary plot file for writing
|
||||
open(UNIT=UNIT_PLOT, FILE=pl % path_plot, STATUS='replace', &
|
||||
|
|
@ -393,20 +385,20 @@ 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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -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, &
|
||||
|
|
@ -28,6 +29,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 +39,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 +47,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 +56,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 +89,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 +102,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
|
||||
|
|
@ -111,17 +111,18 @@ 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 * &
|
||||
material_xs % nu_fission
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
global_tally_tracklength = global_tally_tracklength + p % wgt * &
|
||||
distance * material_xs % nu_fission
|
||||
end if
|
||||
|
||||
if (d_collision > d_boundary) then
|
||||
! ====================================================================
|
||||
! 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
|
||||
|
|
@ -138,10 +139,10 @@ contains
|
|||
! PARTICLE HAS COLLISION
|
||||
|
||||
! Score collision estimate of keff
|
||||
!$omp atomic
|
||||
global_tallies(K_COLLISION) % value = &
|
||||
global_tallies(K_COLLISION) % value + p % wgt * &
|
||||
material_xs % nu_fission / material_xs % total
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
global_tally_collision = global_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
|
||||
|
|
@ -168,7 +169,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 +177,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
|
||||
|
||||
|
|
|
|||
|
|
@ -1,11 +1,11 @@
|
|||
k-combined:
|
||||
1.093844E+00 1.626801E-02
|
||||
tally 1:
|
||||
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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue