Error in tracking for tracks nearly tangent to ZCylinders

[geogunow]

Developing my tests for tracking, I have come across some bugs in our current tracking. Specifically, if a create a track that barely intersects a circle, it calculates a segment length of infinity. Whereas the segment lengths should be roughly:

[1.82467365265, 0.00752119300887, 1.82467365265]

OpenMOC calculates [1.82467365265, 0.00752119300887, inf]

The code to produce this result is:

import openmoc
import math

###############################################################################
#                          Main Simulation Parameters
###############################################################################

options = openmoc.options.Options()

num_threads = options.getNumThreads()
track_spacing = options.getTrackSpacing()
num_azim = options.getNumAzimAngles()
tolerance = options.getTolerance()
max_iters = options.getMaxIterations()

openmoc.log.set_log_level('NORMAL')

###############################################################################
#                            Creating Materials
###############################################################################

openmoc.log.py_printf('NORMAL', 'Importing materials data from HDF5...')

materials = openmoc.materialize.load_from_hdf5('c5g7-mgxs.h5', '../')

###############################################################################
#                            Creating Surfaces
###############################################################################

openmoc.log.py_printf('NORMAL', 'Creating surfaces...')

zcylinder = openmoc.ZCylinder(x=0.0, y=0.0, radius=1.0, name='pin')
left = openmoc.XPlane(x=-2.0, name='left')
right = openmoc.XPlane(x=2.0, name='right')
top = openmoc.YPlane(y=2.0, name='top')
bottom = openmoc.YPlane(y=-2.0, name='bottom')

left.setBoundaryType(openmoc.REFLECTIVE)
right.setBoundaryType(openmoc.REFLECTIVE)
top.setBoundaryType(openmoc.REFLECTIVE)
bottom.setBoundaryType(openmoc.REFLECTIVE)

###############################################################################
#                             Creating Cells
###############################################################################

openmoc.log.py_printf('NORMAL', 'Creating cells...')

fuel = openmoc.Cell(name='fuel')
fuel.setFill(materials['UO2'])
fuel.addSurface(halfspace=-1, surface=zcylinder)

moderator = openmoc.Cell(name='moderator')
moderator.setFill(materials['Water'])
moderator.addSurface(halfspace=+1, surface=zcylinder)
moderator.addSurface(halfspace=+1, surface=left)
moderator.addSurface(halfspace=-1, surface=right)
moderator.addSurface(halfspace=+1, surface=bottom)
moderator.addSurface(halfspace=-1, surface=top)

###############################################################################
#                            Creating Universes
###############################################################################

openmoc.log.py_printf('NORMAL', 'Creating universes...')

root_universe = openmoc.Universe(name='root universe')
root_universe.addCell(fuel)
root_universe.addCell(moderator)

###############################################################################
#                         Creating the Geometry
###############################################################################

openmoc.log.py_printf('NORMAL', 'Creating geometry...')

geometry = openmoc.Geometry()
geometry.setRootUniverse(root_universe)

################ Segmentation Tests

# function to segment a track in a geometry
def segment_track(track, geometry):
  geometry.segmentize(track)
  num_segments = track.getNumSegments()
  print "Num segments = ", num_segments
  for i in range(num_segments):
    segment = track.getSegment(i)
    length = segment._length
    fsr_id = segment._region_id
    fwd = segment._cmfd_surface_fwd
    bwd = segment._cmfd_surface_bwd
    mat_id = segment._material.getId()
    print i, length, fwd, bwd, mat_id
  print ''
  track.clearSegments()

nudge_tan_track = openmoc.Track()
offset = math.sqrt(2) - 2
offset -= 1e-5
nudge_tan_track.setValues(offset, -2, 0, 2, -offset, 0, math.atan(1))
segment_track(nudge_tan_track, geometry)

[samuelshaner]

I think this can be fixed by changing the following line in Cell::containsPoint(...):

if (iter->second->_surface->evaluate(point) * iter->second->_halfspace
    < -ON_SURFACE_THRESH)

to:

if (iter->second->_surface->evaluate(point) * iter->second->_halfspace
    < 0.0)

[wbinventor]

Fixed in #262.