mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Rewrote all geometry routines to allow for nested universes.
This commit is contained in:
parent
06bd0b5626
commit
a843d245b6
10 changed files with 772 additions and 780 deletions
|
|
@ -140,6 +140,7 @@ output.o: string.o
|
|||
output.o: tally_header.o
|
||||
|
||||
particle_header.o: constants.o
|
||||
particle_header.o: geometry_header.o
|
||||
|
||||
physics.o: constants.o
|
||||
physics.o: cross_section_header.o
|
||||
|
|
|
|||
|
|
@ -261,6 +261,9 @@ module constants
|
|||
! ============================================================================
|
||||
! MISCELLANEOUS CONSTANTS
|
||||
|
||||
! indicates that an array index hasn't been set
|
||||
integer, parameter :: NONE = 0
|
||||
|
||||
! Codes for read errors -- better hope these numbers are never used in an
|
||||
! input file!
|
||||
integer, parameter :: ERROR_INT = -huge(0)
|
||||
|
|
|
|||
1092
src/geometry.f90
1092
src/geometry.f90
File diff suppressed because it is too large
Load diff
|
|
@ -347,16 +347,15 @@ contains
|
|||
i_source = i_start + i - 1
|
||||
p => source_bank(i_source)
|
||||
|
||||
p % xyz = temp_bank(i) % xyz
|
||||
p % xyz_local = temp_bank(i) % xyz
|
||||
p % last_xyz = temp_bank(i) % xyz
|
||||
p % uvw = temp_bank(i) % uvw
|
||||
p % E = temp_bank(i) % E
|
||||
p % last_E = p % E
|
||||
|
||||
! set defaults
|
||||
call initialize_particle(p)
|
||||
|
||||
p % coord % xyz = temp_bank(i) % xyz
|
||||
p % coord % uvw = temp_bank(i) % uvw
|
||||
p % last_xyz = temp_bank(i) % xyz
|
||||
p % E = temp_bank(i) % E
|
||||
p % last_E = temp_bank(i) % E
|
||||
|
||||
end do
|
||||
|
||||
end subroutine copy_from_bank
|
||||
|
|
|
|||
|
|
@ -210,44 +210,40 @@ contains
|
|||
case default
|
||||
write(ou,*) 'Unknown Particle ' // int_to_str(p % id)
|
||||
end select
|
||||
write(ou,*) ' x = ' // real_to_str(p % xyz(1))
|
||||
write(ou,*) ' y = ' // real_to_str(p % xyz(2))
|
||||
write(ou,*) ' z = ' // real_to_str(p % xyz(3))
|
||||
write(ou,*) ' x local = ' // real_to_str(p % xyz_local(1))
|
||||
write(ou,*) ' y local = ' // real_to_str(p % xyz_local(2))
|
||||
write(ou,*) ' z local = ' // real_to_str(p % xyz_local(3))
|
||||
write(ou,*) ' u = ' // real_to_str(p % uvw(1))
|
||||
write(ou,*) ' v = ' // real_to_str(p % uvw(2))
|
||||
write(ou,*) ' w = ' // real_to_str(p % uvw(3))
|
||||
write(ou,*) ' x = ' // real_to_str(p % coord0 % xyz(1))
|
||||
write(ou,*) ' y = ' // real_to_str(p % coord0 % xyz(2))
|
||||
write(ou,*) ' z = ' // real_to_str(p % coord0 % xyz(3))
|
||||
write(ou,*) ' x local = ' // real_to_str(p % coord % xyz(1))
|
||||
write(ou,*) ' y local = ' // real_to_str(p % coord % xyz(2))
|
||||
write(ou,*) ' z local = ' // real_to_str(p % coord % xyz(3))
|
||||
write(ou,*) ' u = ' // real_to_str(p % coord0 % uvw(1))
|
||||
write(ou,*) ' v = ' // real_to_str(p % coord0 % uvw(2))
|
||||
write(ou,*) ' w = ' // real_to_str(p % coord0 % uvw(3))
|
||||
write(ou,*) ' Weight = ' // real_to_str(p % wgt)
|
||||
write(ou,*) ' Energy = ' // real_to_str(p % E)
|
||||
write(ou,*) ' x index = ' // int_to_str(p % index_x)
|
||||
write(ou,*) ' y index = ' // int_to_str(p % index_y)
|
||||
write(ou,*) ' x index = ' // int_to_str(p % coord % lattice_x)
|
||||
write(ou,*) ' y index = ' // int_to_str(p % coord % lattice_y)
|
||||
write(ou,*) ' IE = ' // int_to_str(p % IE)
|
||||
write(ou,*) ' Interpolation factor = ' // real_to_str(p % interp)
|
||||
|
||||
if (p % cell > 0) then
|
||||
c => cells(p % cell)
|
||||
if (p % coord % cell /= NONE) then
|
||||
c => cells(p % coord % cell)
|
||||
write(ou,*) ' Cell = ' // int_to_str(c % id)
|
||||
else
|
||||
write(ou,*) ' Cell not determined'
|
||||
end if
|
||||
|
||||
if (p % surface > 0) then
|
||||
if (p % surface /= NONE) then
|
||||
s => surfaces(p % surface)
|
||||
write(ou,*) ' Surface = ' // int_to_str(s % id)
|
||||
else
|
||||
write(ou,*) ' Surface = None'
|
||||
end if
|
||||
|
||||
u => universes(p % universe)
|
||||
u => universes(p % coord % universe)
|
||||
write(ou,*) ' Universe = ' // int_to_str(u % id)
|
||||
write(ou,*)
|
||||
|
||||
nullify(c)
|
||||
nullify(s)
|
||||
nullify(u)
|
||||
|
||||
end subroutine print_particle
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -325,10 +321,6 @@ contains
|
|||
write(ou,*) ' Surface Specification:' // trim(string)
|
||||
write(ou,*)
|
||||
|
||||
! nullify associated pointers
|
||||
nullify(u)
|
||||
nullify(m)
|
||||
|
||||
end subroutine print_cell
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -353,8 +345,6 @@ contains
|
|||
write(ou,*) ' Cells =' // trim(string)
|
||||
write(ou,*)
|
||||
|
||||
nullify(c)
|
||||
|
||||
end subroutine print_universe
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -485,8 +475,6 @@ contains
|
|||
end if
|
||||
write(ou,*)
|
||||
|
||||
nullify(nuc)
|
||||
|
||||
end subroutine print_material
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -1,9 +1,33 @@
|
|||
module particle_header
|
||||
|
||||
use constants, only: NEUTRON, ONE
|
||||
use constants, only: NEUTRON, ONE, NONE
|
||||
use geometry_header, only: BASE_UNIVERSE
|
||||
|
||||
implicit none
|
||||
|
||||
!===============================================================================
|
||||
! LOCALCOORD describes the location of a particle local to a single
|
||||
! universe. When the geometry consists of nested universes, a particle will have
|
||||
! a list of coordinates in each level
|
||||
!===============================================================================
|
||||
|
||||
type LocalCoord
|
||||
! Indices in various arrays for this level
|
||||
integer :: cell = NONE
|
||||
integer :: universe = NONE
|
||||
integer :: lattice = NONE
|
||||
integer :: lattice_x = NONE
|
||||
integer :: lattice_y = NONE
|
||||
|
||||
! Particle position and direction for this level
|
||||
real(8) :: xyz(3)
|
||||
real(8) :: uvw(3)
|
||||
|
||||
! Pointers to next (lower) and previous (higher) universe
|
||||
type(LocalCoord), pointer :: next => null()
|
||||
type(LocalCoord), pointer :: prev => null()
|
||||
end type LocalCoord
|
||||
|
||||
!===============================================================================
|
||||
! PARTICLE describes the state of a particle being transported through the
|
||||
! geometry
|
||||
|
|
@ -14,10 +38,12 @@ module particle_header
|
|||
integer(8) :: id ! Unique ID
|
||||
integer :: type ! Particle type (n, p, e, etc)
|
||||
|
||||
! Physical data
|
||||
real(8) :: xyz(3) ! location
|
||||
real(8) :: xyz_local(3) ! local location (after transformations)
|
||||
real(8) :: uvw(3) ! directional cosines
|
||||
! Particle coordinates
|
||||
logical :: in_lower_universe ! is particle in lower universe?
|
||||
type(LocalCoord), pointer :: coord0 ! coordinates on universe 0
|
||||
type(LocalCoord), pointer :: coord ! coordinates on lowest universe
|
||||
|
||||
! Other physical data
|
||||
real(8) :: wgt ! particle weight
|
||||
real(8) :: E ! energy
|
||||
real(8) :: mu ! angle of scatter
|
||||
|
|
@ -36,15 +62,10 @@ module particle_header
|
|||
real(8) :: interp ! interpolation factor for energy grid
|
||||
|
||||
! Indices for various arrays
|
||||
integer :: cell ! index for current cell
|
||||
integer :: surface ! index for surface particle is on
|
||||
integer :: cell_born ! index for cell particle was born in
|
||||
integer :: universe ! index for current universe
|
||||
integer :: lattice ! index for current lattice
|
||||
integer :: surface ! index for current surface
|
||||
integer :: material ! index for current material
|
||||
integer :: last_material ! index for last material
|
||||
integer :: index_x ! lattice index for x direction
|
||||
integer :: index_y ! lattice index for y direction
|
||||
|
||||
! Statistical data
|
||||
integer :: n_collision ! # of collisions
|
||||
|
|
@ -66,22 +87,48 @@ contains
|
|||
! passed through the fission bank to the source bank, no lookup would be
|
||||
! needed at the beginning of a cycle
|
||||
|
||||
p % type = NEUTRON
|
||||
p % type = NEUTRON
|
||||
p % alive = .true.
|
||||
|
||||
! clear attributes
|
||||
p % surface = NONE
|
||||
p % cell_born = NONE
|
||||
p % material = NONE
|
||||
p % last_material = NONE
|
||||
p % wgt = ONE
|
||||
p % last_wgt = ONE
|
||||
p % alive = .true.
|
||||
p % n_bank = 0
|
||||
p % cell = 0
|
||||
p % cell_born = 0
|
||||
p % universe = 0
|
||||
p % lattice = 0
|
||||
p % surface = 0
|
||||
p % material = 0
|
||||
p % last_material = 0
|
||||
p % index_x = 0
|
||||
p % index_y = 0
|
||||
p % n_collision = 0
|
||||
|
||||
! remove any original coordinates
|
||||
call deallocate_coord(p % coord0)
|
||||
|
||||
! Set up base level coordinates
|
||||
allocate(p % coord0)
|
||||
p % coord0 % universe = BASE_UNIVERSE
|
||||
p % coord => p % coord0
|
||||
p % in_lower_universe = .false.
|
||||
|
||||
end subroutine initialize_particle
|
||||
|
||||
!===============================================================================
|
||||
! DEALLOCATE_COORD removes all levels of coordinates below a given level. This
|
||||
! is used in distance_to_boundary when the particle moves from a lower universe
|
||||
! to a higher universe since the data for the lower one is not needed anymore.
|
||||
!===============================================================================
|
||||
|
||||
recursive subroutine deallocate_coord(coord)
|
||||
|
||||
type(LocalCoord), pointer :: coord
|
||||
|
||||
if (associated(coord)) then
|
||||
! recursively deallocate lower coordinates
|
||||
if (associated(coord % next)) call deallocate_coord(coord%next)
|
||||
|
||||
! deallocate original coordinate
|
||||
deallocate(coord)
|
||||
end if
|
||||
|
||||
end subroutine deallocate_coord
|
||||
|
||||
end module particle_header
|
||||
|
|
|
|||
|
|
@ -11,7 +11,7 @@ module physics
|
|||
use global
|
||||
use interpolation, only: interpolate_tab1
|
||||
use output, only: write_message, print_particle
|
||||
use particle_header, only: Particle
|
||||
use particle_header, only: Particle, LocalCoord
|
||||
use random_lcg, only: prn
|
||||
use search, only: binary_search
|
||||
use string, only: int_to_str
|
||||
|
|
@ -29,29 +29,27 @@ contains
|
|||
|
||||
type(Particle), pointer :: p
|
||||
|
||||
integer :: surf ! surface which particle is on
|
||||
integer :: last_cell ! most recent cell particle was in
|
||||
integer :: n_event ! number of collisions/crossings
|
||||
real(8) :: d_boundary ! distance to nearest boundary
|
||||
real(8) :: d_collision ! sampled distance to collision
|
||||
real(8) :: distance ! distance particle travels
|
||||
logical :: found_cell ! found cell which particle is in?
|
||||
logical :: in_lattice ! is surface crossing in lattice?
|
||||
type(Universe), pointer :: univ
|
||||
integer :: surface_crossed ! surface which particle is on
|
||||
integer :: last_cell ! most recent cell particle was in
|
||||
integer :: n_event ! number of collisions/crossings
|
||||
real(8) :: d_boundary ! distance to nearest boundary
|
||||
real(8) :: d_collision ! sampled distance to collision
|
||||
real(8) :: distance ! distance particle travels
|
||||
logical :: found_cell ! found cell which particle is in?
|
||||
logical :: lattice_crossed ! is surface crossing in lattice?
|
||||
type(LocalCoord), pointer :: coord
|
||||
|
||||
if (p % cell == 0) then
|
||||
univ => universes(BASE_UNIVERSE)
|
||||
call find_cell(univ, p, found_cell)
|
||||
|
||||
! if particle couldn't be located, print error
|
||||
if (p % coord % cell == NONE) then
|
||||
call find_cell(p, found_cell)
|
||||
! Particle couldn't be located
|
||||
if (.not. found_cell) then
|
||||
write(message, '(A,3ES11.3)') &
|
||||
"Could not locate cell for particle at: ", p % xyz
|
||||
"Could not locate cell for particle at: ", p % coord0 % xyz
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! set birth cell attribute
|
||||
p % cell_born = p % cell
|
||||
p % cell_born = p % coord % cell
|
||||
end if
|
||||
|
||||
if (verbosity >= 9) then
|
||||
|
|
@ -60,7 +58,8 @@ contains
|
|||
end if
|
||||
|
||||
if (verbosity >= 10) then
|
||||
message = " Born in cell " // trim(int_to_str(cells(p%cell)%id))
|
||||
message = " Born in cell " // trim(int_to_str(&
|
||||
cells(p % coord % cell) % id))
|
||||
call write_message()
|
||||
end if
|
||||
|
||||
|
|
@ -74,7 +73,7 @@ contains
|
|||
call calculate_xs(p)
|
||||
|
||||
! Find the distance to the nearest boundary
|
||||
call distance_to_boundary(p, d_boundary, surf, in_lattice)
|
||||
call distance_to_boundary(p, d_boundary, surface_crossed, lattice_crossed)
|
||||
|
||||
! Sample a distance to collision
|
||||
d_collision = -log(prn()) / material_xs % total
|
||||
|
|
@ -83,26 +82,37 @@ contains
|
|||
distance = min(d_boundary, d_collision)
|
||||
|
||||
! Advance particle
|
||||
p % xyz = p % xyz + distance * p % uvw
|
||||
p % xyz_local = p % xyz_local + distance * p % uvw
|
||||
coord => p % coord0
|
||||
do while (associated(coord))
|
||||
coord % xyz = coord % xyz + distance * coord % uvw
|
||||
coord => coord % next
|
||||
end do
|
||||
|
||||
if (d_collision > d_boundary) then
|
||||
last_cell = p % cell
|
||||
p % cell = 0
|
||||
if (in_lattice) then
|
||||
p % surface = 0
|
||||
last_cell = p % coord % cell
|
||||
p % coord % cell = NONE
|
||||
if (lattice_crossed) then
|
||||
p % surface = NONE
|
||||
call cross_lattice(p)
|
||||
else
|
||||
p % surface = surf
|
||||
p % surface = surface_crossed
|
||||
call cross_surface(p, last_cell)
|
||||
end if
|
||||
else
|
||||
! collision
|
||||
p % surface = 0
|
||||
p % surface = NONE
|
||||
call collision(p)
|
||||
|
||||
! Save coordinates at collision for tallying purposes
|
||||
p % last_xyz = p % xyz
|
||||
p % last_xyz = p % coord % xyz
|
||||
|
||||
! Set all uvws to base level -- right now, after a collision, only the
|
||||
! base level uvws are changed
|
||||
coord => p % coord0
|
||||
do while(associated(coord))
|
||||
coord % uvw = p % coord0 % uvw
|
||||
coord => coord % next
|
||||
end do
|
||||
end if
|
||||
|
||||
! If particle has too many events, display warning and kill it
|
||||
|
|
@ -712,7 +722,7 @@ contains
|
|||
awr = nuc % awr
|
||||
|
||||
! Neutron velocity in LAB
|
||||
v_n = vel * p % uvw
|
||||
v_n = vel * p % coord0 % uvw
|
||||
|
||||
! Sample velocity of target nucleus
|
||||
call sample_target_velocity(p, nuc, v_t)
|
||||
|
|
@ -750,7 +760,7 @@ contains
|
|||
|
||||
! Set energy and direction of particle in LAB frame
|
||||
p % E = E
|
||||
p % uvw = v_n / vel
|
||||
p % coord0 % uvw = v_n / vel
|
||||
|
||||
! Copy scattering cosine for tallies
|
||||
p % mu = mu
|
||||
|
|
@ -902,13 +912,13 @@ contains
|
|||
end if
|
||||
|
||||
! copy directional cosines
|
||||
u = p % uvw(1)
|
||||
v = p % uvw(2)
|
||||
w = p % uvw(3)
|
||||
u = p % coord0 % uvw(1)
|
||||
v = p % coord0 % uvw(2)
|
||||
w = p % coord0 % uvw(3)
|
||||
|
||||
! change direction of particle
|
||||
call rotate_angle(u, v, w, mu)
|
||||
p % uvw = (/ u, v, w /)
|
||||
p % coord0 % uvw = (/ u, v, w /)
|
||||
|
||||
! change energy of particle
|
||||
p % E = E
|
||||
|
|
@ -993,9 +1003,9 @@ contains
|
|||
|
||||
! determine direction of target velocity based on the neutron's velocity
|
||||
! vector and the sampled angle between them
|
||||
u = p % uvw(1)
|
||||
v = p % uvw(2)
|
||||
w = p % uvw(3)
|
||||
u = p % coord0 % uvw(1)
|
||||
v = p % coord0 % uvw(2)
|
||||
w = p % coord0 % uvw(3)
|
||||
call rotate_angle(u, v, w, mu)
|
||||
|
||||
! determine speed of target nucleus
|
||||
|
|
@ -1093,7 +1103,7 @@ contains
|
|||
do i = int(n_bank,4) + 1, int(min(n_bank + nu, 3*work),4)
|
||||
! Bank source neutrons by copying particle data
|
||||
fission_bank(i) % id = p % id
|
||||
fission_bank(i) % xyz = p % xyz
|
||||
fission_bank(i) % xyz = p % coord0 % xyz
|
||||
|
||||
! sample cosine of angle
|
||||
mu = sample_angle(rxn, E)
|
||||
|
|
@ -1250,13 +1260,13 @@ contains
|
|||
end if
|
||||
|
||||
! copy directional cosines
|
||||
u = p % uvw(1)
|
||||
v = p % uvw(2)
|
||||
w = p % uvw(3)
|
||||
u = p % coord0 % uvw(1)
|
||||
v = p % coord0 % uvw(2)
|
||||
w = p % coord0 % uvw(3)
|
||||
|
||||
! change direction of particle
|
||||
call rotate_angle(u, v, w, mu)
|
||||
p % uvw = (/ u, v, w /)
|
||||
p % coord0 % uvw = (/ u, v, w /)
|
||||
|
||||
! change energy of particle
|
||||
p % E = E
|
||||
|
|
|
|||
189
src/plot.f90
189
src/plot.f90
|
|
@ -6,7 +6,7 @@ module plot
|
|||
cross_lattice, cell_contains
|
||||
use geometry_header, only: Universe, BASE_UNIVERSE
|
||||
use global
|
||||
use particle_header, only: Particle, initialize_particle
|
||||
use particle_header, only: Particle, initialize_particle, LocalCoord
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -18,20 +18,21 @@ contains
|
|||
|
||||
subroutine run_plot()
|
||||
|
||||
integer :: i ! loop index
|
||||
integer :: surf ! surface which particle is on
|
||||
integer :: last_cell ! most recent cell particle was in
|
||||
real(8) :: coord(3) ! starting coordinates
|
||||
real(8) :: last_x_coord ! bounding x coordinate
|
||||
real(8) :: last_y_coord ! bounding y coordinate
|
||||
real(8) :: d ! distance to boundary
|
||||
real(8) :: distance ! distance particle travels
|
||||
logical :: found_cell ! found cell which particle is in?
|
||||
logical :: in_lattice ! is surface crossing in lattice?
|
||||
integer :: i ! loop index
|
||||
integer :: surface_crossed ! surface which particle is on
|
||||
integer :: last_cell ! most recent cell particle was in
|
||||
real(8) :: xyz(3) ! starting coordinates
|
||||
real(8) :: last_x_coord ! bounding x coordinate
|
||||
real(8) :: last_y_coord ! bounding y coordinate
|
||||
real(8) :: d ! distance to boundary
|
||||
real(8) :: distance ! distance particle travels
|
||||
logical :: found_cell ! found cell which particle is in?
|
||||
logical :: lattice_crossed ! is surface crossing in lattice?
|
||||
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
|
||||
type(Cell), pointer :: c => null()
|
||||
type(Universe), pointer :: univ => null()
|
||||
type(Particle), pointer :: p => null()
|
||||
type(Cell), pointer :: c => null()
|
||||
type(Universe), pointer :: univ => null()
|
||||
type(Particle), pointer :: p => null()
|
||||
type(LocalCoord), pointer :: coord => null()
|
||||
|
||||
! Open plot file for binary writing
|
||||
path_plot = trim(path_input) // "plot.out"
|
||||
|
|
@ -44,9 +45,9 @@ contains
|
|||
write(UNIT=UNIT_PLOT) pixel
|
||||
|
||||
! Determine coordinates of the upper-left corner of the plot
|
||||
coord(1) = plot_origin(1) - plot_width(1) / 2.0
|
||||
coord(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
|
||||
coord(3) = plot_origin(3)
|
||||
xyz(1) = plot_origin(1) - plot_width(1) / 2.0
|
||||
xyz(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
|
||||
xyz(3) = plot_origin(3)
|
||||
|
||||
! Determine bounding x and y coordinates for plot
|
||||
last_x_coord = plot_origin(1) + plot_width(1) / 2.0
|
||||
|
|
@ -56,76 +57,75 @@ contains
|
|||
allocate(p)
|
||||
|
||||
! loop over horizontal rays
|
||||
do while(coord(2) > last_y_coord)
|
||||
do while(xyz(2) > last_y_coord)
|
||||
|
||||
! initialize the particle and set starting coordinate and direction
|
||||
call initialize_particle(p)
|
||||
p % xyz = coord
|
||||
p % xyz_local = coord
|
||||
p % uvw = (/ 1, 0, 0 /)
|
||||
|
||||
p % coord % xyz = xyz
|
||||
p % coord % uvw = (/ 1, 0, 0 /)
|
||||
|
||||
! write starting coordinate to file
|
||||
write(UNIT=UNIT_PLOT) p % xyz
|
||||
write(UNIT=UNIT_PLOT) p % coord % xyz
|
||||
|
||||
! Find cell that particle is currently in
|
||||
univ => universes(BASE_UNIVERSE)
|
||||
call find_cell(univ, p, found_cell)
|
||||
call find_cell(p, found_cell)
|
||||
|
||||
! =======================================================================
|
||||
! MOVE PARTICLE FORWARD TO NEXT CELL
|
||||
|
||||
if (.not. found_cell) then
|
||||
univ => universes(BASE_UNIVERSE)
|
||||
do i = 1, univ % n_cells
|
||||
p % xyz = coord
|
||||
p % xyz_local = coord
|
||||
p % cell = univ % cells(i)
|
||||
|
||||
distance = INFINITY
|
||||
call distance_to_boundary(p, d, surf, in_lattice)
|
||||
if (d < distance) then
|
||||
! Move particle forward to next surface
|
||||
p % xyz = p % xyz + d * p % uvw
|
||||
|
||||
! Check to make sure particle is actually going into this cell
|
||||
! by moving it slightly forward and seeing if the cell contains
|
||||
! that coordinate
|
||||
|
||||
p % xyz = p % xyz + 1e-4 * p % uvw
|
||||
p % xyz_local = p % xyz
|
||||
|
||||
c => cells(p % cell)
|
||||
if (.not. cell_contains(c, p)) cycle
|
||||
|
||||
! Reset coordinate to surface crossing
|
||||
p % xyz = p % xyz - 1e-4 * p % uvw
|
||||
p % xyz_local = p % xyz
|
||||
|
||||
! Set new distance and retain pointer to this cell
|
||||
distance = d
|
||||
last_cell = p % cell
|
||||
end if
|
||||
end do
|
||||
|
||||
! No cell was found on this horizontal ray
|
||||
if (distance == INFINITY) then
|
||||
p % xyz(1) = last_x_coord
|
||||
p % cell = 0
|
||||
write(UNIT_PLOT) p % xyz, p % cell
|
||||
|
||||
! Move to next horizontal ray
|
||||
coord(2) = coord(2) - pixel
|
||||
cycle
|
||||
end if
|
||||
|
||||
! Write coordinate where next cell begins
|
||||
write(UNIT=UNIT_PLOT) p % xyz, 0
|
||||
|
||||
! Process surface crossing for next cell
|
||||
p % cell = 0
|
||||
p % surface = -surf
|
||||
call cross_surface(p, last_cell)
|
||||
end if
|
||||
!!$ if (.not. found_cell) then
|
||||
!!$ univ => universes(BASE_UNIVERSE)
|
||||
!!$ do i = 1, univ % n_cells
|
||||
!!$ p % xyz = coord
|
||||
!!$ p % xyz_local = coord
|
||||
!!$ p % cell = univ % cells(i)
|
||||
!!$
|
||||
!!$ distance = INFINITY
|
||||
!!$ ! call distance_to_boundary(p, d, surf, in_lattice)
|
||||
!!$ if (d < distance) then
|
||||
!!$ ! Move particle forward to next surface
|
||||
!!$ p % xyz = p % xyz + d * p % uvw
|
||||
!!$
|
||||
!!$ ! Check to make sure particle is actually going into this cell
|
||||
!!$ ! by moving it slightly forward and seeing if the cell contains
|
||||
!!$ ! that coordinate
|
||||
!!$
|
||||
!!$ p % xyz = p % xyz + 1e-4 * p % uvw
|
||||
!!$ p % xyz_local = p % xyz
|
||||
!!$
|
||||
!!$ c => cells(p % cell)
|
||||
!!$ if (.not. cell_contains(c, p)) cycle
|
||||
!!$
|
||||
!!$ ! Reset coordinate to surface crossing
|
||||
!!$ p % xyz = p % xyz - 1e-4 * p % uvw
|
||||
!!$ p % xyz_local = p % xyz
|
||||
!!$
|
||||
!!$ ! Set new distance and retain pointer to this cell
|
||||
!!$ distance = d
|
||||
!!$ last_cell = p % cell
|
||||
!!$ end if
|
||||
!!$ end do
|
||||
!!$
|
||||
!!$ ! No cell was found on this horizontal ray
|
||||
!!$ if (distance == INFINITY) then
|
||||
!!$ p % xyz(1) = last_x_coord
|
||||
!!$ p % cell = 0
|
||||
!!$ write(UNIT_PLOT) p % xyz, p % cell
|
||||
!!$
|
||||
!!$ ! Move to next horizontal ray
|
||||
!!$ xyz(2) = xyz(2) - pixel
|
||||
!!$ cycle
|
||||
!!$ end if
|
||||
!!$
|
||||
!!$ ! Write coordinate where next cell begins
|
||||
!!$ write(UNIT=UNIT_PLOT) p % xyz, 0
|
||||
!!$
|
||||
!!$ ! Process surface crossing for next cell
|
||||
!!$ p % cell = 0
|
||||
!!$ p % surface = -surf
|
||||
!!$ call cross_surface(p, last_cell)
|
||||
!!$ end if
|
||||
|
||||
! =======================================================================
|
||||
! MOVE PARTICLE ACROSS HORIZONTAL TRACK
|
||||
|
|
@ -133,43 +133,46 @@ contains
|
|||
do while (p % alive)
|
||||
|
||||
! Calculate distance to next boundary
|
||||
call distance_to_boundary(p, distance, surf, in_lattice)
|
||||
call distance_to_boundary(p, distance, surface_crossed, lattice_crossed)
|
||||
|
||||
! Advance particle
|
||||
p%xyz = p%xyz + distance * p%uvw
|
||||
p%xyz_local = p%xyz_local + distance * p%uvw
|
||||
coord => p % coord0
|
||||
do while (associated(coord))
|
||||
coord % xyz = coord % xyz + distance * coord % uvw
|
||||
coord => coord % next
|
||||
end do
|
||||
|
||||
! If next boundary crossing is out of range of the plot, only include
|
||||
! the visible portion and move to next horizontal ray
|
||||
if (p % xyz(1) >= last_x_coord) then
|
||||
if (p % coord0 % xyz(1) >= last_x_coord) then
|
||||
p % alive = .false.
|
||||
p % xyz(1) = last_x_coord
|
||||
p % coord0 % xyz(1) = last_x_coord
|
||||
|
||||
! If there is no cell beyond this boundary, mark it as cell 0
|
||||
if (distance == INFINITY) p % cell = 0
|
||||
if (distance == INFINITY) p % coord % cell = 0
|
||||
|
||||
! Write ending coordinates to file
|
||||
write(UNIT=UNIT_PLOT) p % xyz, p % cell
|
||||
write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell
|
||||
cycle
|
||||
end if
|
||||
|
||||
! Write boundary crossing coordinates to file
|
||||
write(UNIT=UNIT_PLOT) p % xyz, p % cell
|
||||
write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell
|
||||
|
||||
last_cell = p % cell
|
||||
p % cell = 0
|
||||
if (in_lattice) then
|
||||
p % surface = 0
|
||||
last_cell = p % coord % cell
|
||||
p % coord % cell = 0
|
||||
if (lattice_crossed) then
|
||||
p % surface = NONE
|
||||
call cross_lattice(p)
|
||||
else
|
||||
p % surface = surf
|
||||
p % surface = surface_crossed
|
||||
call cross_surface(p, last_cell)
|
||||
|
||||
! Since boundary conditions are disabled in plotting mode, we need
|
||||
! to manually add the last segment
|
||||
if (surfaces(surf) % bc == BC_VACUUM) then
|
||||
p % xyz(1) = last_x_coord
|
||||
write(UNIT=UNIT_PLOT) p % xyz, 0
|
||||
if (surfaces(surface_crossed) % bc == BC_VACUUM) then
|
||||
p % coord0 % xyz(1) = last_x_coord
|
||||
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
|
|
@ -177,7 +180,7 @@ contains
|
|||
end do
|
||||
|
||||
! Move y-coordinate to next position
|
||||
coord(2) = coord(2) - pixel
|
||||
xyz(2) = xyz(2) - pixel
|
||||
end do
|
||||
|
||||
! Close plot file
|
||||
|
|
|
|||
|
|
@ -79,25 +79,24 @@ contains
|
|||
do j = bank_first, bank_last
|
||||
p => source_bank(j - bank_first + 1)
|
||||
|
||||
! set defaults
|
||||
call initialize_particle(p)
|
||||
|
||||
! initialize random number seed
|
||||
call set_particle_seed(int(j,8))
|
||||
|
||||
! sample position
|
||||
r = (/ (prn(), k = 1,3) /)
|
||||
p % id = j
|
||||
p % xyz = p_min + r*(p_max - p_min)
|
||||
p % xyz_local = p % xyz
|
||||
p % last_xyz = p % xyz
|
||||
p % coord0 % xyz = p_min + r*(p_max - p_min)
|
||||
p % last_xyz = p % coord0 % xyz
|
||||
|
||||
! sample angle
|
||||
phi = TWO*PI*prn()
|
||||
mu = TWO*prn() - ONE
|
||||
p % uvw(1) = mu
|
||||
p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
|
||||
p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
|
||||
|
||||
! set defaults
|
||||
call initialize_particle(p)
|
||||
p % coord0 % uvw(1) = mu
|
||||
p % coord0 % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
|
||||
p % coord0 % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
|
||||
|
||||
! sample energy from Watt fission energy spectrum for U-235
|
||||
do
|
||||
|
|
|
|||
|
|
@ -343,7 +343,7 @@ contains
|
|||
|
||||
! determine next universe bin
|
||||
if (t % n_bins(T_UNIVERSE) > 0) then
|
||||
bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % universe, i)
|
||||
bins(T_UNIVERSE) = get_next_bin(T_UNIVERSE, p % coord % universe, i)
|
||||
if (bins(T_UNIVERSE) == NO_BIN_FOUND) cycle
|
||||
else
|
||||
bins(T_UNIVERSE) = 1
|
||||
|
|
@ -359,7 +359,7 @@ contains
|
|||
|
||||
! determine next cell bin
|
||||
if (t % n_bins(T_CELL) > 0) then
|
||||
bins(T_CELL) = get_next_bin(T_CELL, p % cell, i)
|
||||
bins(T_CELL) = get_next_bin(T_CELL, p % coord % cell, i)
|
||||
if (bins(T_CELL) == NO_BIN_FOUND) cycle
|
||||
else
|
||||
bins(T_CELL) = 1
|
||||
|
|
@ -386,7 +386,7 @@ contains
|
|||
m => meshes(t % mesh)
|
||||
|
||||
! Determine if we're in the mesh first
|
||||
call get_mesh_bin(m, p % xyz, mesh_bin, in_mesh)
|
||||
call get_mesh_bin(m, p % coord0 % xyz, mesh_bin, in_mesh)
|
||||
if (.not. in_mesh) cycle
|
||||
|
||||
bins(T_MESH) = mesh_bin
|
||||
|
|
@ -610,7 +610,7 @@ contains
|
|||
do i = 1, n_tallies
|
||||
! Copy starting and ending location of particle
|
||||
xyz0 = p % last_xyz
|
||||
xyz1 = p % xyz
|
||||
xyz1 = p % coord0 % xyz
|
||||
|
||||
! Get pointer to tally
|
||||
t => tallies(i)
|
||||
|
|
@ -631,7 +631,7 @@ contains
|
|||
if (n_cross == 0) cycle
|
||||
|
||||
! Copy particle's direction
|
||||
uvw = p % uvw
|
||||
uvw = p % coord0 % uvw
|
||||
|
||||
! determine incoming energy bin
|
||||
n = t % n_bins(T_ENERGYIN)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue