diff --git a/src/geometry.f90 b/src/geometry.f90 index 4a1051ce8c..2c3baa78e1 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -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 diff --git a/src/physics.f90 b/src/physics.f90 index 028bf19201..2df86795d3 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -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