Reflective boundary conditions now work for planes. Closes gh-35.

This commit is contained in:
Paul Romano 2011-09-08 20:28:55 -04:00
parent 70e6eaf458
commit 645149a917
2 changed files with 14 additions and 6 deletions

View file

@ -172,9 +172,10 @@ contains
! CROSS_SURFACE moves a particle into a new cell
!===============================================================================
subroutine cross_surface(p)
subroutine cross_surface(p, last_cell)
type(Particle), pointer :: p
integer, intent(in) :: last_cell ! last cell particle was in
integer :: i ! index of neighbors
integer :: index_cell ! index in cells array
@ -244,6 +245,11 @@ contains
call fatal_error(msg)
end select
! Reassign particle's cell and surface
p % cell = last_cell
p % surface = -p % surface
! Diagnostic message
if (verbosity >= 10) then
msg = " Reflected from surface " // trim(int_to_str(surf%uid))
call message(msg)
@ -522,7 +528,7 @@ contains
surf_p => surfaces(index_surf)
select case (surf_p%type)
case (SURF_PX)
if (u == ZERO) then
if (on_surface .or. u == ZERO) then
d = INFINITY
else
x0 = surf_p % coeffs(1)
@ -531,7 +537,7 @@ contains
end if
case (SURF_PY)
if (v == ZERO) then
if (on_surface .or. v == ZERO) then
d = INFINITY
else
y0 = surf_p % coeffs(1)
@ -540,7 +546,7 @@ contains
end if
case (SURF_PZ)
if (w == ZERO) then
if (on_surface .or. w == ZERO) then
d = INFINITY
else
z0 = surf_p % coeffs(1)
@ -555,7 +561,7 @@ contains
D = surf_p % coeffs(4)
tmp = A*u + B*v + C*w
if (tmp == ZERO) then
if (on_surface .or. tmp == ZERO) then
d = INFINITY
else
d = -(A*x + B*y + C*w - D)/tmp

View file

@ -28,6 +28,7 @@ contains
type(Particle), pointer :: p
integer :: surf ! surface which particle is on
integer :: last_cell ! most recent cell particle was in
integer :: IE ! index on energy grid
real(8) :: d_to_boundary ! distance to nearest boundary
real(8) :: d_to_collision ! sampled distance to collision
@ -87,13 +88,14 @@ contains
p%xyz_local = p%xyz_local + distance * p%uvw
if (d_to_collision > d_to_boundary) then
last_cell = p % cell
p % cell = 0
if (in_lattice) then
p % surface = 0
call cross_lattice(p)
else
p % surface = surf
call cross_surface(p)
call cross_surface(p, last_cell)
end if
else
! collision