mit-crpg/openmc — pull request #168 by bhermanmit (2013-04-10 18:33:58+00:00)
Title: Tracking

Using the awesome new particle restart capability I finally got around to figuring out why a tracking error occurs in my two-region lattice surrounded by reflective boundary conditions. The two region lattice geometry is shown below.

![lattice](https://f.cloud.github.com/assets/1017998/363595/14d0e63c-a208-11e2-83be-8ec11dab4ffa.png)

A tracking error occurs when a neutron reflects off a surface that is close to a corner. An example of this is shown below (not to scale) where a neutron starts in a lattice cell and traverses to a reflective surface that is close to a corner.

![trackingerror](https://f.cloud.github.com/assets/1017998/363579/9fdc6298-a207-11e2-9a4b-229899ef6e27.png)

After the reflection math is performed, the following block of code is called from the cross_surface routine in geometry.F90:
```fortran
if (associated(p % coord0 % next)) then
call deallocate_coord(p % coord0 % next)
call find_cell(found)
if (.not. found) then
call write_particle_restart()
message = "Couldn't find particle after reflecting from surface."
call fatal_error()
end if
end if
```
OpenMC will now try to locate the neutron again after the cell reflection by calling find_cell routine in geometry.F90. In executing find_cell recursively, OpenMC will eventually determine that the neutron is in a lattice. OpenMC always moves the neutron a ```TINY_BIT``` in its direction so that it is not exactly on a surface and should have no problem computing the correct lattice indices ```(i_x, i_y, i_z)```. This is operation is shown in the following code in find_cell:
```fortran
! now move lattice tiny bit along its unit vector
xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
i_x = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
i_y = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
n_x = lat % dimension(1)
n_y = lat % dimension(2)
if (lat % n_dimension == 3) then
i_z = ceiling((xyz(3) - lat % lower_left(3))/lat % width(3))
n_z = lat % dimension(3)
else
i_z = 1
n_z = 1
end if
```
In the case that I am showing above, the neutron is so close to the corner that when it is moved that ```TINY_BIT```, it effectively is now outside the lattice and the index ```i_x``` will be greater than the max x index ```n_x```. When this occurs the code will not find the particle since it is outside of the lattice bounds.

To fix this tracking error I first compute the original indices of the neutron before moving it a ```TINY_BIT```. This assumes that since OpenMC put that particle *exactly* on the surface, it will return indices that are inside the lattice. This is done in the following code:
```fortran
! determine lattice index based on position
xyz = p % coord % xyz
i_x0 = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
i_y0 = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
if (lat % n_dimension == 3) &
i_z0 = ceiling((xyz(3) - lat % lower_left(3))/lat % width(3))
```

Then, I allow the particle to be moved the ```TINY_BIT``` along its unit vector and calculate the new indices. Then I go direction by direction and check and see if by moving this particle, we disturbed its lattice location by moving it out. This is checked with the following code:
```fortran
! Check if neutron out of lattice bounds by moving TINY_BIT.
! This can happen if neutron is within tolerance of TINY_BIT
! to corner, this action will put neutron out of lattice.
if (i_x0 >= 1 .and. i_x < 1) i_x = 1 ! -x direction
if (i_y0 >= 1 .and. i_y < 1) i_y = 1 ! -y direction
if (i_z0 >= 1 .and. i_z < 1) i_z = 1 ! -z direction
if (i_x0 <= n_x .and. i_x > n_x) i_x = n_x ! +x direction
if (i_y0 <= n_x .and. i_y > n_x) i_y = n_y ! +y direction
if (i_z0 <= n_x .and. i_z > n_x) i_z = n_z ! +z direction
```

If any of these checks are satisfied, the lattice indices are either set to the max or the min lattice bounds.