Using MCBEND for neutron or gamma-ray deterministic calculations

MCBEND 11 is the latest version of the general radiation transport Monte Carlo code from AMEC Foster Wheeler’s ANSWERS ® Software Service. MCBEND is well established in the UK shielding community for radiation shielding and dosimetry assessments. MCBEND supports a number of acceleration techniques, for example the use of an importance map in conjunction with Splitting/Russian Roulette. MCBEND has a well established automated tool to generate this importance map, commonly referred to as the MAGIC module using a diffusion adjoint solution. This method is fully integrated with the MCBEND geometry and material specification, and can easily be run as part of a normal MCBEND calculation. An often overlooked feature of MCBEND is the ability to use this method for forward scoping calculations, which can be run as a very quick deterministic method. Additionally, the development of the Visual Workshop environment for results display provides new capabilities for the use of the forward calculation as a productivity tool. In this paper, we illustrate the use of the combination of the old and new in order to provide an enhanced analysis capability. We also explore the use of more advanced deterministic methods for scoping calculations used in conjunction with MCBEND, with a view to providing a suite of methods to accompany the main Monte Carlo solver.


Introduction
MCBEND 11 [1] is the latest version of the general radiation transport Monte Carlo code from AMEC Foster Wheeler's ANSWERS ® Software Service. MCBEND is well established in the UK shielding community for radiation shielding and dosimetry assessments. The MCBEND package comprises the code, a range of nuclear data libraries, user and application guides, and geometry visualisation and results tools provided by VisualWorkshop.
Some Monte Carlo radiation transport codes have a pedigree of several decades, and the MCBEND code falls into this category. As such, there have been many model developments, including some that have stood the test of time and others that have not. In this paper we explore the history of some of these methods, and show how some of the older methods might have new relevance, in particular when mixed with newer features. In particular we look at two strands of development involving deterministic methods to show how this mix of the old and the new might provide additional tools for shielding applications.
The first strand involves the use of an often overlooked feature in MCBEND that can be used to provide a diffusion solution in a matter of seconds, using the same model input as the full Monte Carlo solution. The use of modern tools for processing and displaying results makes this an attractive option for performing scoping calculations. By revisiting the method to make it compatible with more MCBEND features (eg newer source modules) this becomes relevant to the current user.
The second route we explore is an evolution of the way a tetrahedral based geometry description can be utilised in MCBEND. Originally this feature was developed as a method for representing CAD models in MCBEND, but it has now been extended to provide a toolkit to switch between a modern, fully featured, deterministic methodology based on a tetrahedral mesh, and the MCBEND Monte Carlo method. This is illustrated by use with the FETCH2 code developed at Imperial College, London [2].

Scoping calculations with MCBEND
The Scoping Tool is an often overlooked feature of MCBEND. This mature feature provides a fast scoping capability based on the diffusion solver used by the MCBEND 'MAGIC' option for generating importance maps. Usually this method is used to accelerate the main MCBEND Monte Carlo solution, and is used routinely for this purpose in adjoint mode. However it can also be run in forward mode to provide scoping calculations.
When used to support acceleration, the importance map generated is used in conjunction with Splitting/Russian Roulette. The source for the adjoint is taken from a response function appropriate to the forward calculation. Part of the attraction of the method for acceleration purposes is the integration of MAGIC with the MCBEND geometry and material specification units; this means the method can be run easily as part of a normal MCBEND calculation.
The method is quick and easy to use, with the MCBEND geometry mapped onto either a Cartesian or a Cylindrical 3D mesh automatically.
This ease of use also extends to the forward mode. The integration of the scoping tool with MCBEND means that results can be generated everywhere in the MCBEND model and displayed in the VisualWorkshop environment. An example results display is shown in the Figure below for a PWR case.

How the scoping tool works
A diffusion calculation is carried out by the MAGIC module, using the splitting mesh and standard 33, 28 or 22 group schemes, depending on the calculation method. Sources and materials defined in the MCBEND input are mapped to the splitting mesh used by the diffusion calculation.
The results of the diffusion calculation are available in the computational mesh for the diffusion solver (the splitting mesh), but there is also an option to provide results in the material mesh, in which case the results are mapped back to the zones that define the material, facilitating easy direct comparison of results with a Monte-Carlo run. If scoring in the splitting mesh is requested the results are available for VisualWorkshop to display.
The two figures above show the display of dose in a flask model and also the display of an isosurface for flux. We return to this flask model in Section 2.3.

Historical Restrictions
The released version of MCBEND, MCBEND11A Release Update 1, provides for this scoping style calculation for use only in neutron calculations. It also requires the source to be specified in an older, deprecated, set of modules referred to as "Simple Source". This source module, together with an alternative "Complex Source" module have been largely superseded by a "Unified Source" module in MCBEND, which combines the functionality of the two old modules, and is written in modern FORTRAN. There are still however a few features in MCBEND that require the older source modules, and the scoping tool is one of them. Additionally, the source for scoping calculations was restricted to being a histogram source (as a function of energy), with the source separable in space and energy.
These restrictions have now been relaxed, and development versions of MCBEND now allow scoping calculations to be performed using the Unified Source module, and also no longer require the source description to be separable in energy and space. Currently this is done by a Monte Carlo sampling of the source description. Additionally the restriction to neutrons has been removed, so that gamma calculations can also be performed with the scoping tool.

A case study in the use of the scoping tool
In this example, the scoping tool is illustrated by use of a transport flask case study. This example models UO2 fuel rods in a flask, partially flooded. It consists of a set of five fuel baskets in a partially-flooded flask as shown below.

Fig. 4. Cross-section of a flask model.
Initially the splitting mesh from the calculation of the importance map in the original case was used. This was a simple RZ splitting mesh, with no azimuthal discretisation. A few boundaries were added in the Z dimension to coincide with the plane surfaces of external scoring regions. However, the simple radial subdivision within the flask, with no angular subdivision, effectively homogenized the materials within this region, causing the doses in the flask components to be underpredicted by a factor of two or three.
An improved mesh was then specified in the R and θ dimensions, which is illustrated below.  The revised results, give reasonable agreement between the diffusion and Monte Carlo calculations for the flask components. However, the doses at the lid and 1m above the lid were grossly exaggerated. It was judged that this problem was due to a pathological effect -the diffusion equations predict that scattering can occur in a void, permitting neutrons to reach the lid after they have passed through the (thinner) flask wall.
The problem was solved by specifying black absorber round the side of the flask and beyond a radius of 1m, as shown below. In a range of neutron case calculations it has proved possible to achieve an accuracy of about 20 or 30% in most scoring regions, including at the lid and 1m above the lid. However it can be important to understand the limitations of diffusion solutions, as illustrated by this example. In this case the pathological effects, due to neutrons from the side being scattered to the lid could be avoided through the use of suitably positioned black absorbers. The automatic nature of the scoping tool, with minimal user intervention, coupled with the availability of graphical output processing and graphical representations of the geometry make it a very flexible scoping tool. This is further improved by upgrades to provide compatibility with the most up-to-date modules in MCBEND (including the Unified Source module) and a relaxation of historical restrictions on its use.

Tetrahedral mesh capabilities
MCBEND has the capability to track in a tetrahedral mesh. This capability was originally developed over ten years ago in order to allow import of a CAD mesh. CAD import for ANSWERS codes has subsequently been developed further in a number of ways, see for example [3], but the tetrahedral tracking capability has opened up a route to partnering MCBEND with state-of-the-art deterministic codes, such as Imperial College's FETCH2 code which is currently being evaluated for possible use with MCBEND.
MCBEND provides two implementations of a tetrahedral mesh, the first uses Woodcock tracking in Hole geometries to perform the Monte Carlo calculation. Woodcock tracking is an alternative to conventional tracking and avoids the need to calculate the distance to surfaces by instead using the minimum mean free path of all the materials in the region of interest, and using this to sample for the distance between collisions. To compensate for the modification to the mean free path, the idea of pseudo-collisions is introduced in which a sampling process takes place following a collision to decide whether the collision is a real collision, or a "do nothing" pseudo collision. The sampling process for real collisions is based on the ratio of the minimum mean free path to the actual mean free path. This form of tracking provides an exact representation of the transport equation.
The second method is explicit (conventional) tracking to the tetrahedral surfaces. This is still fast as only plane surfaces are involved. One potential advantage of the explicit tracking approach is that it facilitates the use of track length scoring in the tetrahedral mesh.
Although we now have an implementation of this form of scoring, it may not be such an advantage in practice because it is not necessary for the scoring mesh to be based on the tracking mesh. In MCBEND we have a "Unified Tally" module where scoring bodies can be defined independently of the geometry that defines the materials, with tracking to scoring body surfaces only being necessary at a change of direction of the particle. This may be more useful than scoring in the tetrahedral mesh in practice.
Although originally intended for modelling CAD geometries in MCBEND, the application we are concerned with here is the potential for sharing a mesh/geometry with deterministic codes, including stateof-the-art methodologies. This is explored below using the FETCH2 code.

FETCH2
FETCH2 is a successor to Imperial College's EVENT code which has been widely used for many years for radiation transport calculations. It is a general purpose finite element neutral particle transport solver for neutrons and photons in critical and sub-critical systems. It is written to be highly parallelisable and scales well up to many thousands of processors.
FETCH2 is built around a novel finite element method which combines the advantages of continuous and discontinuous finite element methods. The solver is based on matrix-vector multiplications which ensures that it is future proofed against advances in computer technology, such as use on GPGPUs and Xeon Phi coprocessors. FETCH2 uses domain decomposition to distribute the spatial mesh over multi-processor machines together with a space-angle multigrid scheme to increase the efficiency of the solvers.
FETCH2 can solve the transport equation using the following angular discretisation methods: -Spherical Harmonics (P n ) -Discrete Ordinates (S n ) -Linear Octahedral Wavelets -Haar Wavelets making it suitable for shielding, criticality and reactor physics applications.
To reduce the burden on the user, FETCH2 has the unique ability to anisotropically adapt both the spatial and angular resolutions, either globally or to some goal (such as a detector response).
To illustrate the use of FETCH2 in conjunction with MCBEND, we apply the codes to a waste store model, which can be represented in MCBEND either by using the usual simple body representation, or using a tetrahedral representation.  In order to test the capabilities for sharing a common tetrahedral mesh between FETCH2 and MCBEND, the waste store model shown above was used as a representative shielding problem.

Waste Store Model
Here the material in green is concrete, and the duct has dimensions 40 cm in the x direction, 60 cm in the y direction. We are interested in the gamma-ray dose-rate both at the exit of the duct, but also in the variation along the duct. The dominant contribution to the dose-rate arises from particles which enter the mouth of the duct and scatter at the corners before reaching the duct exit. The first FETCH2 model used for the waste store example had a fixed 3D mesh (ie with no mesh adaptivity). The open source GMSH tool [4] was used to generate the mesh which incorporated approximately 423000 elements, with the tetrahedral element size mostly of order 65 cm, but refined to about 5 cm in the vicinity of the duct. Eleven gamma groups were used from the Bugle96 library. The methods used were linear continuous plus piecewise constant Subgrid Scale. The scatter model was Spherical Harmonics anisotropic scatter (using P 0 up to P 5 ). Matrix free spatial multigrid solvers were applied using PETSc [5] FGMRES smoothing. For the 3D model, the calculations were run on 64 parallel processors. Doses were scored at duct centre points (linear interpolation within bounding element). Doses were calculated using ICRP-51 (AP and LAT). The angular discretisation was Equal Weight S n quadrature. A comparison of the FETCH2 calculations with MCBEND using the tetrahedral mesh and explicit tracking is given in the table below. For the P 5 case, the FETCH2 calculations are generally within a factor two of the MCBEND calculations, and so FETCH2 is in reasonable agreement with MCBEND. The MCBEND calculations used the Forced Flight acceleration method in addition to other acceleration methods. Forced Flight causes a particle to be transported directly to a defined geometrical interface at every collision, with the particle weight suitably adjusted to avoid any bias. A feature of forced flight that was new in MCBEND 11 was the ability to have multiple forced flight interfaces acting in series, and this facility was used here.
There are still some known model differences between the FETCH2 calculations and the MCBEND calculations in the area of source definition, and nuclear data. Additionally there is ongoing analysis to assess a range of methods available in FETCH2, including: x Wavelets for angular discretisation x The available Multi-Grid Solvers x The Adaptive resolution capability FETCH2's powerful adaptive resolution capability is briefly illustrated below, as it incorporates features of interest in their own right.
Two approaches to adaptive resolution are available in FETCH2: Regular adaptivity which aims to reduce error everywhere; and Goal based adaptivity which aims to reduce the error for a specified goal. The Goal based method relies on solving the forward and adjoint transport equation -the latter providing the "importance" similar to MAGIC.
The adaptive resolution is fully anisotropic (anisotropic in space and angle). However the illustration below uses spatial adaptivity alone for a two dimensional slice of the waste store model, using: x S 4 angular quadrature x Solution adapted 10 times x Two adaptivity tolerances The modelling choices allow fast turnaround of calculations with modest computing resources.
The Figure below illustrates the use of regular spatial adaptivity. The plot on the right has the finer adaptive tolerance. The figure illustrates: x Resolution highlighting ray effects (S4) x Little resolution in the duct (flux varies from order 10 to 10E-9 making defining a regular adaptive tolerance challenging) x Decreasing the tolerance begins to add more resolution to the duct as well as substantially more resolution elsewhere. The next figure now uses Goal based spatial adaptivity, with the goal set as average flux in a thin region near the duct exit. The goal based adaptivity clearly out-performs the regular adaptivity in this case, in the sense of achieving much more resolution in and around the duct, and also in the dominant pathway to the duct. The figure below shows a close up of the duct region, where the effects of the spatial adaptivity can be seen more clearly. The darker regions where the mesh is of finer resolution will likely be a mix of ray effects and real physical effects.

Conclusions
This paper describes how a number of features in MCBEND of various vintages have been combined to allow additional tools for shielding applications. The adjoint diffusion capability used for generating an importance map for the main Monte Carlo calculation also provides for a very quick running forward scoping tool for neutron calculations. This has now been extended to include gamma calculations, and also to be compatible with the newer source module, Unified Source, which was released with MCBEND 10. Additionally the improved facilities for results display in Visual Workshop, where results can be overlaid on the geometry gives this old feature a new lease of life which can allow calculations and results to be displayed in seconds.
At the other end of the deterministic scale is the ability to track and score on a tetrahedral mesh, which allows MCBEND to partner state-of-the-art deterministic codes if these codes use a tetrahedral spatial discretisation. One such code is the FETCH2 code developed at Imperial College, London. This makes it possible to run MCBEND calculations on the same mesh as is used for the deterministic calculation, with benefits for porting models between codes. Although MCBEND can score in the tetrahedral mesh, it can also score in bodies dedicated to scoring, such as subdivided cylinders and cuboids, using the latest scoring module, Unified Tally, which was released with MCBEND 11. This module can be used in conjunction with a model based on a tetrahedral mesh for tracking, providing benefits for sharing model between codes.
The dusting off and extension of older methods, combined with some more recent capabilities has allowed the old and the new to come together to create new tools for the shielding practitioner, and illustrates a benefit of keeping older methods current, and also as testing methods advance to much more automation, there is a benefit to ensuring that older features are comprehensively covered in any test suite as far as practical.