Delta- and surface-tracking

This a brief description on the delta-tracking based transport routine used in Serpent. The method is also used by other Monte Carlo codes, most notably in the HOLE geometry package in MONK and MCBEND. The original delta-tracking algorithm was introduced by Woodcock in 1965, and a mathematical verification was derived by Coleman in 1968. Delta-tracking is well described in a text book by Lux and Koblinger, and a description of the methodology used in Serpent is found in an article in Annals of Nuclear Energy from 2010.

The input parameters related to delta-tracking are:

• set dt - sets the probability threshold used for selecting between surface- and delta-tracking
• set forcedt - enforces the use of delta-tracking in a given list of materials
• set blockdt - enforces the use of surface-tracking in a given list of materials
• set minxs - definse the mean-free-path of collisions used to score the collision flux estimator

The output parameters are:

• TODO

Transport algorithms in Monte Carlo simulation

The Monte Carlo simulation consists of a large number particle histories, in which the random walk of an individual particle is followed, or tracked, through the geometry from its birth to eventual absorption or escape. Tracking is most typically carried out in a constructive solid geometry (CSG), composed of homogeneous material cells, which are defined by combinations of elementary and derived surface types. Serpent 2 also has two advanced geometry types, based on STL format CAD models and an unstructured polyhedral mesh. 

The transport simulation follows a random walk from one interaction to the next. The procedure can be described as follows:

1. Sample path length (distance to next collision)
2. Transport particle to the collision point
3. Sample interaction

If the sampled interaction is scattering, the procedure restarts from beginning by sampling the distance to the next collision. The direction and energy are changed in each scattering event.

By definition, the interaction probability per traveled path length is given by the macroscopic total cross section (denoted here as $\Sigma$). If it is assumed that the particle travels through an infinite homogeneous medium characterized by constant total cross section, it can be shown that the free path length follows an exponential distribution. This distribution can be sampled using the inverse method, and the distance to the next collision site is given by:

$l = -\frac{1}{\Sigma}\log\xi$

where $\xi$ is a uniformly distributed random variable on the unit interval.

The prerequisite of using this simple formula for sampling the distance to the next collision site is that the material is infinite and homogeneous. This is the case when the particle remains within a single material region, but not if its track crosses the boundary between two materials.

Surface-tracking

Since the interaction probability per traveled path length does not depend on the events that occured in the particle's history, any point in the path can be considered the stating point of a new path. If the sampled path crosses the boundary between two materials, the track can be stopped at the crossing point and a new path length sampled using the cross section of the next material. This is the general idea in the surface-tracking algorithm.

In order to stop the particle at the boundary between two materials, the algorithm needs to calculate the distance to the nearest boundary surface. The only way to accomplish this is to loop over all candidate surfaces and pick the shortest value.

Surface tracking is considered the standard tracking algorithm and it is used by virtually every Monte Carlo transport calculation code. The method has a few drawbacks related to its efficiency in complex geometries:

• Determining the distance to the nearest boundary can become computationally expensive if the cells are comprised of a large number of surfaces
• The fact that the particle track has to be stopped at each boundary crossing becomes a computational bottleneck when the mean-free-path is long compared to geometry dimensions

Delta-tracking

An alternative to surface-tracking is the Woodcock delta-tracking algorithm, which is based on the rejection sampling of particle path lengths. The procedure relies on the concept of a virtual collision, which is a fictive interaction that preserves the energy and direction of the transported particle.

Since virtual collisions do not change the statistics of the random walk in any way, the material total cross section $\Sigma$ can be adjusted with an arbitrary virtual collision cross section $\Sigma_0$:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.):

without changing the outcome of the simulation. It is therefore possible to adjust the cross sections of all material regions 1, 2, ... in the system such that:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.):

where $\Sigma_\mathrm{m}$ is called the majorant cross section.

In practice, it is not necessary to define the virtual collision cross sections at all if the majorant is simply taken as the maximum of all material totals at each energy point:

Failed to parse (Missing <code>texvc</code> executable. Please see math/README to configure.):

Unlike the physical total cross section, which depends on the material located at the particle position, the majorant cross section is completely independent of the spatial coordinates.

The point of having a macroscopic cross section that is uniform throughout the geometry is that when it is used for sampling path lengths:

$l = -\log(\xi)/\Sigma_\mathrm{m}$

the values are statistically valid regardless of the number of material boundaries crossed.

At the end point of the sampled path the tracking routine performs rejection sampling. The probability to accept the collision is given by ratio of the physical total cross section to the majorant:

$P = \frac{\Sigma(\mathbf{r},E)}{\Sigma_\mathrm{m}(E)}$

If the collision is rejected, a new path length is sampled using \Sigma_\mathrm{m} and the particle is moved to the next tentative collision site.

Since the majorant cross section is always larger than or equal to the total cross section, the path lengths sampled using delta-tracking are, on the average, shorter than those sampled with surface-tracking. The average physical distance between two collisions is preserved, as some paths are extended over multiple virtual collisions.

Advantages and limitations of delta-tracking

The advantage of delta-tracking over the surface-tracking algorithm is that there is no need to calculate the surface distances or stop the particle at the material boundaries. This becomes significant for computational performance in geometries where the particle mean-free-path is long compared to dimensions. The difference between surface- and delta-tracking is emphasized in HTGR particle fuels and advanced CAD and an mesh based geometry types.

Since the majorant cross section used for sampling the path lengths does not depend on the spatial coordinates and the material total is needed only at discrete locations, variations of delta-tracking can be used for transporting particles through inhomogeneous materials. This capability is utilized in the TMS on-the-fly temperature treatment routine and the multi-physics interface in Serpent 2.

Delta-tracking also has its drawbacks. Since the majorant cross section reflects the largest interaction probability within the system, the efficiency of the rejection sampling loop may become poor in the presence of localized heavy absorbers (control rods, burnable absorber pins, etc.) that dominate the majorant cross section, but occupy a relatively small volume of the geometry.

Another drawback is that delta-tracking rules out the use of the track-length estimate (TLE) of particle flux, and integral reaction rate estimates need to be calculated using the potentially less efficient collision estimator (CFE).

-- implements template:multiple image local p = {}

local function isnotempty(s) return s and s:match( '^%s*(.-)%s*$' ) ~= end

local function renderImageCell(image, width, height, link, alt, caption, textalign, istyle) local root = mw.html.create()

local altstr = '|alt=' .. (alt or ) local linkstr = link and ('|link=' .. link) or local widthstr = '|' .. tostring(width) .. 'px'

local imagediv = root:tag('div') imagediv:addClass('thumbimage') imagediv:cssText(istyle) if( height ) then imagediv:css('height', tostring(height) .. 'px') imagediv:css('overflow', 'hidden') end imagediv:wikitext('File:' .. image .. widthstr .. linkstr .. altstr .. '') if isnotempty(caption) then local captiondiv = root:tag('div') captiondiv:addClass('thumbcaption') captiondiv:css('clear', 'left') if isnotempty(textalign) then captiondiv:css('text-align', textalign) end captiondiv:wikitext(caption) end return tostring(root) end

local function getWidth(w1, w2) local w if isnotempty(w1) then w = tonumber(w1) elseif isnotempty(w2) then w = tonumber(w2) end return w or 200 end

local function getPerRow(pstr, ic) -- split string into array using any non-digit as a dilimiter local pr = mw.text.split(pstr or , '[^%d][^%d]*') -- if split failed, assume a single row if (#pr < 1) then pr = {tostring(ic)} end -- convert the array of strings to an array of numbers, -- adding any implied/missing numbers at the end of the array local r = 1 local thisrow = tonumber(pr[1] or ic) or ic local prownum = {} while( ic > 0 ) do prownum[r] = thisrow ic = ic - thisrow r = r + 1 -- use the previous if the next is missing and -- make sure we don't overstep the number of images thisrow = math.min(tonumber(pr[r] or thisrow) or ic, ic) end return prownum end

local function renderMultipleImages(frame) local pargs = frame:getParent().args local args = frame.args local width = pargs['width'] or local dir = pargs['direction'] or local align = pargs['align'] or args['align'] or local capalign = pargs['caption_align'] or args['caption_align'] or local totalwidth = pargs['total_width'] or args['total_width'] or local imgstyle = pargs['image_style'] or args['image_style'] local header = pargs['header'] or pargs['title'] or local footer = pargs['footer'] or local perrow = nil local thumbclass = { ["left"] = 'tleft', ["none"] = 'tnone', ["center"] = 'tnone', ["centre"] = 'tnone', ["right"] = 'tright' }

-- find all the nonempty images local imagenumbers = {} local imagecount = 0 for k, v in pairs( pargs ) do local i = tonumber(tostring(k):match( '^%s*image([%d]+)%s*$' ) or '0') if( i > 0 and isnotempty(v) ) then table.insert( imagenumbers, i) imagecount = imagecount + 1 end end

-- sort the imagenumbers table.sort(imagenumbers)

-- create an array with the number of images per row perrow = getPerRow(dir == 'vertical' and '1' or pargs['perrow'], imagecount)

-- compute the number of rows local rowcount = #perrow

-- store the image widths and compute row widths and maximum row width local widths = {} local widthmax = 0 local widthsum = {} local k = 0 for r=1,rowcount do widthsum[r] = 0 for c=1,perrow[r] do k = k + 1 if( k <= imagecount ) then local i = imagenumbers[k] widths[k] = getWidth(width, pargs['width' .. i]) widthsum[r] = widthsum[r] + widths[k] end end widthmax = math.max(widthmax, widthsum[r]) end

-- if total_width has been specified, rescale the image widths local heights = {} if( isnotempty(totalwidth) ) then totalwidth = tonumber(totalwidth) widthmax = 0 local k = 0 for r=1,rowcount do local koffset = k local tw = totalwidth - 4 * (perrow[r] - 1) - 12 local ar = {} local arsum = 0 for j=1,perrow[r] do k = k + 1 if( k<= imagecount ) then local i = imagenumbers[k] local h = tonumber( pargs['height' .. i] or ) or 0 if (h > 0) then ar[j] = widths[k]/h heights[k] = h else ar[j] = widths[k]/100 end arsum = arsum + ar[j] end end local ht = tw/arsum local ws = 0 k = koffset for j=1,perrow[r] do k = k + 1 if( k<= imagecount ) then local i = imagenumbers[k] widths[k] = math.floor(ar[j]*ht + 0.5) ws = ws + widths[k] if heights[k] then heights[k] = math.floor(ht) end end end widthsum[r] = ws widthmax = math.max(widthmax, widthsum[r]) end end

-- start building the array of images, if there are images if( imagecount > 0 ) then -- compute width of outer div local bodywidth = 0 for r=1,rowcount do if( widthmax == widthsum[r] ) then bodywidth = widthmax + 4 * (perrow[r] - 1) + 12 end end -- The body has a min-width of 100, which needs to be taken into account on specific widths bodywidth = math.max( 100, bodywidth - 8);

local bg = pargs['background color'] or -- create the array of images local root = mw.html.create('div') root:addClass('thumb') root:addClass('tmulti') root:addClass(thumbclass[align] or 'tright')

if( align == 'center' or align == 'centre' ) then root:addClass('center') end if( pargs['margin_top'] or args['margin_top']) then root:css('margin-top', pargs['margin_top'] or args['margin_top']) end if( pargs['margin_bottom'] or args['margin_bottom']) then root:css('margin-bottom', pargs['margin_bottom'] or args['margin_bottom']) end if( bg ~= ) then root:css('background-color', bg) end

local div = root:tag('div') div:addClass('thumbinner') div:css('width', tostring(bodywidth) .. 'px') :css('max-width', tostring(bodywidth) .. 'px') if( bg ~= ) then div:css('background-color', bg) end -- add the header if( isnotempty(header) ) then div:tag('div') :css('clear', 'both') :css('font-weight', 'bold') :css('text-align', pargs['header_align'] or 'center') :css('background-color', pargs['header_background'] or 'transparent') :wikitext(header) end -- loop through the images local k = 0 for r=1,rowcount do for j=1,perrow[r] do k = k + 1 if( k <= imagecount ) then local imagediv = div:tag('div') imagediv:addClass('tsingle') if dir ~= 'vertical' then imagediv:css('float', 'left') end if bg ~= then imagediv:css('background-color', bg); end imagediv:css('margin', '1px') local i = imagenumbers[k] local img = pargs['image' .. i] local w = widths[k] imagediv:css('width', tostring(2 + w) .. 'px') :css('max-width', tostring(2 + w) .. 'px') imagediv:wikitext(renderImageCell(img, w, heights[k], pargs['link' .. i], pargs['alt' .. i], pargs['caption' .. i], capalign, imgstyle)) end end -- only float content gives a parent height:0, so add a clearing div if dir ~= 'vertical' then div:tag('div') :css('clear', 'left') end end -- add the footer if( isnotempty(footer) ) then div:tag('div') :addClass('thumbcaption') :css('clear', 'left') :css('text-align', pargs['footer_align'] or args['footer_align'] or 'left') :css('background-color', pargs['footer_background'] or 'transparent') :wikitext(footer) end return tostring(root) end return end

function p.render( frame )

   return renderMultipleImages( frame )

end

return p