crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Move cross_surface to tracking (break geometry -> tally dependence)

Browse source
This commit is contained in:
Paul Romano 2017-12-11 14:21:34 +07:00
parent c771d0f0ae
commit fa89104350
2 changed files with 273 additions and 273 deletions
Download patch file Download diff file Expand all files Collapse all files

271
src/geometry.F90
View file

@ -10,7 +10,6 @@ module geometry
use surface_header
use stl_vector, only: VectorInt
use string, only: to_str
use tally, only: score_surface_current
use tally_header
implicit none
@ -379,276 +378,6 @@ contains
end subroutine find_cell
!===============================================================================
! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of
! the geometry, is reflected, or crosses into a new lattice or cell
!===============================================================================
subroutine cross_surface(p)
type(Particle), intent(inout) :: p
real(8) :: u ! x-component of direction
real(8) :: v ! y-component of direction
real(8) :: w ! z-component of direction
real(8) :: norm ! "norm" of surface normal
real(8) :: d ! distance between point and plane
real(8) :: xyz(3) ! Saved global coordinate
integer :: i_surface ! index in surfaces
logical :: rotational ! if rotational periodic BC applied
logical :: found ! particle found in universe?
class(Surface), pointer :: surf
i_surface = abs(p % surface)
surf => surfaces(i_surface)%obj
if (verbosity >= 10 .or. trace) then
call write_message(" Crossing surface " // trim(to_str(surf % id)))
end if
if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then
! =======================================================================
! PARTICLE LEAKS OUT OF PROBLEM
! Kill particle
p % alive = .false.
! Score any surface current tallies -- note that the particle is moved
! forward slightly so that if the mesh boundary is on the surface, it is
! still processed
if (active_current_tallies % size() > 0) then
! TODO: Find a better solution to score surface currents than
! physically moving the particle forward slightly
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
global_tally_leakage = global_tally_leakage + p % wgt
! Display message
if (verbosity >= 10 .or. trace) then
call write_message(" Leaked out of surface " &
&// trim(to_str(surf % id)))
end if
return
elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then
! =======================================================================
! PARTICLE REFLECTS FROM SURFACE
! Do not handle reflective boundary conditions on lower universes
if (p % n_coord /= 1) then
call p % mark_as_lost("Cannot reflect particle " &
// trim(to_str(p % id)) // " off surface in a lower universe.")
return
end if
! Score surface currents since reflection causes the direction of the
! particle to change -- artificially move the particle slightly back in
! case the surface crossing is coincident with a mesh boundary
if (active_current_tallies % size() > 0) then
xyz = p % coord(1) % xyz
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
call score_surface_current(p)
p % coord(1) % xyz = xyz
end if
! Reflect particle off surface
call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw)
! Make sure new particle direction is normalized
u = p%coord(1)%uvw(1)
v = p%coord(1)%uvw(2)
w = p%coord(1)%uvw(3)
norm = sqrt(u*u + v*v + w*w)
p%coord(1)%uvw(:) = [u, v, w] / norm
! Reassign particle's cell and surface
p % coord(1) % cell = p % last_cell(p % last_n_coord)
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.
p % n_coord = 1
call find_cell(p, found)
if (.not. found) then
call p % mark_as_lost("Couldn't find particle after reflecting&
& from surface " // trim(to_str(surf % id)) // ".")
return
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then
call write_message(" Reflected from surface " &
&// trim(to_str(surf%id)))
end if
return
elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then
! =======================================================================
! PERIODIC BOUNDARY
! Do not handle periodic boundary conditions on lower universes
if (p % n_coord /= 1) then
call p % mark_as_lost("Cannot transfer particle " &
// trim(to_str(p % id)) // " across surface in a lower universe.&
& Boundary conditions must be applied to universe 0.")
return
end if
! Score surface currents since reflection causes the direction of the
! particle to change -- artificially move the particle slightly back in
! case the surface crossing is coincident with a mesh boundary
if (active_current_tallies % size() > 0) then
xyz = p % coord(1) % xyz
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
call score_surface_current(p)
p % coord(1) % xyz = xyz
end if
rotational = .false.
select type (surf)
type is (SurfaceXPlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceXPlane)
p % coord(1) % xyz(1) = opposite % x0
type is (SurfaceYPlane)
rotational = .true.
! Rotate direction
u = p % coord(1) % uvw(1)
v = p % coord(1) % uvw(2)
p % coord(1) % uvw(1) = v
p % coord(1) % uvw(2) = -u
! Rotate position
p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0
p % coord(1) % xyz(2) = opposite % y0
end select
type is (SurfaceYPlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceYPlane)
p % coord(1) % xyz(2) = opposite % y0
type is (SurfaceXPlane)
rotational = .true.
! Rotate direction
u = p % coord(1) % uvw(1)
v = p % coord(1) % uvw(2)
p % coord(1) % uvw(1) = -v
p % coord(1) % uvw(2) = u
! Rotate position
p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0
p % coord(1) % xyz(1) = opposite % x0
end select
type is (SurfaceZPlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceZPlane)
p % coord(1) % xyz(3) = opposite % z0
end select
type is (SurfacePlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfacePlane)
! Get surface normal for opposite plane
xyz(:) = opposite % normal(p % coord(1) % xyz)
! Determine distance to plane
norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
d = opposite % evaluate(p % coord(1) % xyz) / norm
! Move particle along normal vector based on distance
p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz
end select
end select
! Reassign particle's surface
if (rotational) then
p % surface = surf % i_periodic
else
p % surface = sign(surf % i_periodic, p % surface)
end if
! Figure out what cell particle is in now
p % n_coord = 1
call find_cell(p, found)
if (.not. found) then
call p % mark_as_lost("Couldn't find particle after hitting &
&periodic boundary on surface " // trim(to_str(surf % id)) // ".")
return
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then
call write_message(" Hit periodic boundary on surface " &
// trim(to_str(surf%id)))
end if
return
end if
! ==========================================================================
! SEARCH NEIGHBOR LISTS FOR NEXT CELL
if (p % surface > 0 .and. allocated(surf%neighbor_pos)) then
! If coming from negative side of surface, search all the neighboring
! cells on the positive side
call find_cell(p, found, surf%neighbor_pos)
if (found) return
elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then
! If coming from positive side of surface, search all the neighboring
! cells on the negative side
call find_cell(p, found, surf%neighbor_neg)
if (found) return
end if
! ==========================================================================
! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
! Remove lower coordinate levels and assignment of surface
p % surface = NONE
p % n_coord = 1
call find_cell(p, found)
if (run_mode /= MODE_PLOTTING .and. (.not. found)) then
! If a cell is still not found, there are two possible causes: 1) there is
! a void in the model, and 2) the particle hit a surface at a tangent. If
! the particle is really traveling tangent to a surface, if we move it
! forward a tiny bit it should fix the problem.
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
! undefined region in the geometry.
if (.not. found) then
call p % mark_as_lost("After particle " // trim(to_str(p % id)) &
// " crossed surface " // trim(to_str(surf % id)) &
// " it could not be located in any cell and it did not leak.")
return
end if
end if
end subroutine cross_surface
!===============================================================================
! CROSS_LATTICE moves a particle into a new lattice element
!===============================================================================

275
src/tracking.F90
View file

@ -4,8 +4,8 @@ module tracking
use cross_section, only: calculate_xs
use error, only: fatal_error, warning
use geometry_header, only: cells
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, check_cell_overlap
use geometry, only: find_cell, distance_to_boundary, cross_lattice, &
check_cell_overlap
use output, only: write_message
use message_passing
use mgxs_header
@ -17,6 +17,7 @@ module tracking
use settings
use simulation_header
use string, only: to_str
use surface_header
use tally_header
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current, &
@ -278,4 +279,274 @@ contains
end subroutine transport
!===============================================================================
! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of
! the geometry, is reflected, or crosses into a new lattice or cell
!===============================================================================
subroutine cross_surface(p)
type(Particle), intent(inout) :: p
real(8) :: u ! x-component of direction
real(8) :: v ! y-component of direction
real(8) :: w ! z-component of direction
real(8) :: norm ! "norm" of surface normal
real(8) :: d ! distance between point and plane
real(8) :: xyz(3) ! Saved global coordinate
integer :: i_surface ! index in surfaces
logical :: rotational ! if rotational periodic BC applied
logical :: found ! particle found in universe?
class(Surface), pointer :: surf
i_surface = abs(p % surface)
surf => surfaces(i_surface)%obj
if (verbosity >= 10 .or. trace) then
call write_message(" Crossing surface " // trim(to_str(surf % id)))
end if
if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then
! =======================================================================
! PARTICLE LEAKS OUT OF PROBLEM
! Kill particle
p % alive = .false.
! Score any surface current tallies -- note that the particle is moved
! forward slightly so that if the mesh boundary is on the surface, it is
! still processed
if (active_current_tallies % size() > 0) then
! TODO: Find a better solution to score surface currents than
! physically moving the particle forward slightly
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
global_tally_leakage = global_tally_leakage + p % wgt
! Display message
if (verbosity >= 10 .or. trace) then
call write_message(" Leaked out of surface " &
&// trim(to_str(surf % id)))
end if
return
elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then
! =======================================================================
! PARTICLE REFLECTS FROM SURFACE
! Do not handle reflective boundary conditions on lower universes
if (p % n_coord /= 1) then
call p % mark_as_lost("Cannot reflect particle " &
// trim(to_str(p % id)) // " off surface in a lower universe.")
return
end if
! Score surface currents since reflection causes the direction of the
! particle to change -- artificially move the particle slightly back in
! case the surface crossing is coincident with a mesh boundary
if (active_current_tallies % size() > 0) then
xyz = p % coord(1) % xyz
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
call score_surface_current(p)
p % coord(1) % xyz = xyz
end if
! Reflect particle off surface
call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw)
! Make sure new particle direction is normalized
u = p%coord(1)%uvw(1)
v = p%coord(1)%uvw(2)
w = p%coord(1)%uvw(3)
norm = sqrt(u*u + v*v + w*w)
p%coord(1)%uvw(:) = [u, v, w] / norm
! Reassign particle's cell and surface
p % coord(1) % cell = p % last_cell(p % last_n_coord)
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.
p % n_coord = 1
call find_cell(p, found)
if (.not. found) then
call p % mark_as_lost("Couldn't find particle after reflecting&
& from surface " // trim(to_str(surf % id)) // ".")
return
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then
call write_message(" Reflected from surface " &
&// trim(to_str(surf%id)))
end if
return
elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then
! =======================================================================
! PERIODIC BOUNDARY
! Do not handle periodic boundary conditions on lower universes
if (p % n_coord /= 1) then
call p % mark_as_lost("Cannot transfer particle " &
// trim(to_str(p % id)) // " across surface in a lower universe.&
& Boundary conditions must be applied to universe 0.")
return
end if
! Score surface currents since reflection causes the direction of the
! particle to change -- artificially move the particle slightly back in
! case the surface crossing is coincident with a mesh boundary
if (active_current_tallies % size() > 0) then
xyz = p % coord(1) % xyz
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
call score_surface_current(p)
p % coord(1) % xyz = xyz
end if
rotational = .false.
select type (surf)
type is (SurfaceXPlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceXPlane)
p % coord(1) % xyz(1) = opposite % x0
type is (SurfaceYPlane)
rotational = .true.
! Rotate direction
u = p % coord(1) % uvw(1)
v = p % coord(1) % uvw(2)
p % coord(1) % uvw(1) = v
p % coord(1) % uvw(2) = -u
! Rotate position
p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0
p % coord(1) % xyz(2) = opposite % y0
end select
type is (SurfaceYPlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceYPlane)
p % coord(1) % xyz(2) = opposite % y0
type is (SurfaceXPlane)
rotational = .true.
! Rotate direction
u = p % coord(1) % uvw(1)
v = p % coord(1) % uvw(2)
p % coord(1) % uvw(1) = -v
p % coord(1) % uvw(2) = u
! Rotate position
p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0
p % coord(1) % xyz(1) = opposite % x0
end select
type is (SurfaceZPlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceZPlane)
p % coord(1) % xyz(3) = opposite % z0
end select
type is (SurfacePlane)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfacePlane)
! Get surface normal for opposite plane
xyz(:) = opposite % normal(p % coord(1) % xyz)
! Determine distance to plane
norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
d = opposite % evaluate(p % coord(1) % xyz) / norm
! Move particle along normal vector based on distance
p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz
end select
end select
! Reassign particle's surface
if (rotational) then
p % surface = surf % i_periodic
else
p % surface = sign(surf % i_periodic, p % surface)
end if
! Figure out what cell particle is in now
p % n_coord = 1
call find_cell(p, found)
if (.not. found) then
call p % mark_as_lost("Couldn't find particle after hitting &
&periodic boundary on surface " // trim(to_str(surf % id)) // ".")
return
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then
call write_message(" Hit periodic boundary on surface " &
// trim(to_str(surf%id)))
end if
return
end if
! ==========================================================================
! SEARCH NEIGHBOR LISTS FOR NEXT CELL
if (p % surface > 0 .and. allocated(surf%neighbor_pos)) then
! If coming from negative side of surface, search all the neighboring
! cells on the positive side
call find_cell(p, found, surf%neighbor_pos)
if (found) return
elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then
! If coming from positive side of surface, search all the neighboring
! cells on the negative side
call find_cell(p, found, surf%neighbor_neg)
if (found) return
end if
! ==========================================================================
! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
! Remove lower coordinate levels and assignment of surface
p % surface = NONE
p % n_coord = 1
call find_cell(p, found)
if (run_mode /= MODE_PLOTTING .and. (.not. found)) then
! If a cell is still not found, there are two possible causes: 1) there is
! a void in the model, and 2) the particle hit a surface at a tangent. If
! the particle is really traveling tangent to a surface, if we move it
! forward a tiny bit it should fix the problem.
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
! undefined region in the geometry.
if (.not. found) then
call p % mark_as_lost("After particle " // trim(to_str(p % id)) &
// " crossed surface " // trim(to_str(surf % id)) &
// " it could not be located in any cell and it did not leak.")
return
end if
end if
end subroutine cross_surface
end module tracking