Why the Monte Carlo method?

The continuous-energy Monte Carlo method has been used for criticality safety analyses, radiation shielding and dose rate calculations, detector modeling and the validation of deterministic transport codes for decades. The main motivator is usually the need to model geometry and interaction physics to within maximum accuracy, often regardless of the computational cost. Monte Carlo codes are well suited for the job, with the capability to handle complicated three-dimensional geometries and to model neutron interactions at the microscopic level without major approximations. In fact, Monte Carlo codes are often used to complement or even replace experimental measurements.

The computational requirements for reactor analysis, in particular the modeling of operating power reactors, are considerably higher. The strong coupling between neutronics, thermal hydraulics and fuel behavior means that the transport problem cannot be solved separately, and even though considerable efforts are devoted to developing various multi-physics calculation systems, the Monte Carlo method cannot yet be considered a practical option for routine design and safety analyses of large LWR cores. The traditional approach to reactor analysis relies instead on a multi-stage calculation scheme, where the complexity of transport physics is gradually simplified, while simultaneously increasing the scale of the modeled system. In practice, the calculation sequence is divided into 1) Spatial homogenization, where the interaction physics at the fuel assembly level is condensed into a set of assembly-specific multi-group constants, and 2) Core simulation, where the full-scale neutronics solution obtained using diffusion theory or other reduced-order transport method is iteratively coupled to thermal hydraulics.

Spatial homogenization (Koebke, 1978; Smith, 1980) is traditionally performed using deterministic lattice transport codes in a 2D geometry, in which a single fuel assembly is separated from its surroundings and modeled as an infinite lattice of identical copies. The calculation is repeated for every assembly type in all reactor operating conditions, and the aim is to preserve the reaction rate balance for the next stage in the calculation chain. The use of Monte Carlo codes for homogenization has gained a lot of interest during the past decade, along with the development of computer capacity and parallel computation. In contrast to the traditional applications, the main advantage of using the method for homogenization is not so much the accuracy, but rather its versatility. The same code and cross section data can be used for modeling any fuel or reactor configuration without compromising the reliability of the calculation scheme. The entire chain from neutron interactions to coupled full-core simulations becomes shorter and more transparent, as Monte Carlo codes are capable of using the best available knowledge on neutron interactions almost directly. The particle transport simulation in is inherently three-dimensional, which makes it possible to capture the effects of axial heterogeneities -- a task not easily accomplished using deterministic 2D lattice codes.

The computational challenges of Monte Carlo based homogenization also differ from its traditional applications. Criticality safety analyses, for example, are typically performed for a relatively limited number of configurations, and the calculations produce a single number, the effective multiplication factor, as the result. Group constant generation, on the other hand, requires covering all assembly types and operating conditions within the reactor core, which translates into thousands of state points and involves simulating fuel burnup. One of the main reasons why Monte Carlo codes have not been widely used for spatial homogenization is the high computational cost, especially for burnup calculation. Another reason is that most general-purpose codes have not been able to simulate fuel burnup until recently, and the generation of certain group constants requires calculation techniques that lie beyond the traditional tally capabilities. One of the original goals of the Serpent project was to show that these limitations can be lifted by developing dedicated Monte Carlo lattice physics codes, and that the continuous-energy Monte Carlo method may become a viable option for spatial homogenization within the near future.

History of the Serpent project

The development of Serpent started at VTT in 2004, under the working title "Probabilistic Scattering Game", or PSG. All publications dated before the official pre-release in October 2008 refer to the code using this name. The name was later changed to Serpent, due to the various ambiguities related to the acronym.

The main reason to start developing a new Monte Carlo neutron transport code from scratch was the fact that the general-purpose codes available at that time were not particularly well suited for lattice physics applications. Group constant generation requires the combination of several user-defined reaction rate tallies and the estimation of combined statistical errors becomes complicated. The overall calculation time increases dramatically along with the number of tally definitions, and the calculation of certain parameters, such as scattering matrices, diffusion coefficients and effective delayed neutron fractions lies completely beyond the standard tally capabilities.


The first version of PSG was started in September 2004. Already at the beginning it was decided to construct the interaction physics based on ACE format data libraries, mainly because the ENDF reaction laws are reasonably well documented, and the same data format is used by MCNP, which made code validation easy and straightforward. In an effort to simplify the calculation routines, it was also decided not to use the continuous-energy cross sections directly, but to reconstruct the data on a master energy grid that was used for all nuclides. Eventually this approach turned out to be very efficient as well, since time-consuming grid search iteration was reduced to minimum.

Another feature, originally implemented for the sake of simplicity, was the Woodcock delta-tracking method (Woodcock, 1965) used for neutron transport. The tracking routine does not involve the calculation of optical distances to material boundaries, which considerably simplified the implementation of the geometry routine. Although not widely used for neutron transport applications due to certain limitations, this method turned out to be reasonably well suited for lattice physics calculations. The main drawback was that the efficiency of the basic delta-tracking method is significantly reduced when localized heavy absorbers, such as control rods or burnable absorber pins are present in the geometry.

The first results with practical significance were obtained in early 2005. The code was able to calculate infinite multiplication factors and homogenized few-group reaction cross sections with promising results. Differences to reference MCNP calculations ranged from 0.05% to about 2%. Some initial results were presented in two international meetings in 2005 (Leppänen, 2005a; 2005b).

The first years

After getting the green light from VTT for continuing the development, the source code was completely re-written by the end of 2005. The new version included a bound-atom scattering model for moderator isotopes, which reduced the differences to reference MCNP calculations to less than 0.5%. Parallel calculation using the Message Passing Interface (MPI) was also implemented. Another major improvement was a new geometry routine that could overcome the localized absorber problems of delta-tracking in simple lattice geometries.

New features, such as temperature-interpolated material mixtures and a geometry routine for modeling randomly-dispersed HTGR particle fuels (Leppänen, 2007b) were added during the next year. (Both of these features were later replaced with better methodologies.) Code development was more or less frozen by the summer of 2006, as the main focus was turned to a doctoral thesis which was completed one year later, in June 2007 (Leppänen, 2007c).

The predecessor of Serpent was used as an in-house code at VTT. The project also spawned the development of MORA ("Monte Carlo Reactor Analysis"), a simplified full-core reactor physics code based on a homogenized multi-group Monte Carlo method. The MORA code was the topic of a conference paper in 2008 (Leppänen, 2008b), but the project has not been continued since.

Serpent 1

The current version of Serpent 1 began to formulate when the source code was completely re-written for a second time by the beginning of 2008. A persistent methodological flaw in the free-gas scattering model was corrected, which finally reduced the differences to MCNP results to the level of statistical accuracy. The old geometry routine was also replaced by a more general approach based on the combination of delta-tracking and the conventional surface tracking method (Leppänen, 2009e).

The development of burnup calculation routines began in early 2008. The first version was based on an external coupling to the ORIGEN2 depletion code using ABURN, a coupling code developed at that time at VTT (Ranta-Aho, 2007). It was soon decided, however, that the code should also have the capability to run burnup calculation as a complete stand-alone application. In order to accomplish this task without additional user effort, the calculation of microscopic one-group transmutation cross sections and the collection of energy-dependent fission product yields and radioactive decay constants had to be completely automated. To avoid becoming dependent on any particular data library or another code, it was decided to read decay and fission yield data directly from unprocessed ENDF format files.

Built-in burnup calculation capability also required the development of internal depletion routines. The first method used for the task was based on Transmutation Trajectory Analysis (TTA) (Cetnar, 2006), which in practice implies the analytical solution of linearized depletion chains. The main advantage of this method is that it can cope with the extensive variation in the transmutation coefficients and time steps. The TTA method was also relatively easy to implement in its basic form using a recursive loop.

The development of an alternative solution method began soon after, and the first implementation was completed by late 2008. The focus of the study was to examine if it would be possible to solve a detailed burnup system containing thousands of nuclides by a single matrix exponential method. Due to the difficult numerical characteristics of burnup matrices, the computation of the matrix exponential solution had not been previously conceivable for a full burnup system. However, analyzing the mathematical properties of burnup matrices led to the discovery that their eigenvalues are generally confined to the vicinity of the negative real axis, which prompted the introduction of a novel matrix exponential method CRAM (Chebyshev Rational Approximation Method) for solving the Bateman equations (Pusa, 2010). CRAM can be characterized as the best rational approximation on the negative real axis and it enables simultaneously solving an entire burnup system consisting of 1200-1700 nuclides with an unprecedented accuracy and a very short computation time (Pusa, 2011; Isotalo, 2011a). The work is summarized in a Doctoral thesis, completed in 2013 (Pusa, 2013b).

A pre-release version (1.0.0) of the Serpent code was made available to some universities and research institutes in October 2008. The code was submitted for public distribution to the NEA Data Bank in April 2009 and released one month later. RSICC distribution to North America began in March 2010.

Serpent 2

The calculation methods in Serpent 1 were developed over a relatively short period, without really considering how the solutions should look as a whole. The result of this "live-and-learn" approach was that new calculation routines were constantly being added on top of existing source code, making the program structure increasingly complicated. As the work proceeded, it became more and more difficult to add new features, while keeping the old ones from falling apart. In addition to the problems with maintainability, it became clear that the excessive memory usage and parallelization based on MPI would eventually become major limitations for burnup calculation. Sooner or later it would simply become impossible to keep up with computer development, heading more and more in the direction of multi-core CPU's and massive parallelization.

In September 2010 it was decided that the best solution to the problems in Serpent 1 was to re-write entire source code again. The work started under working title "Super-Serpent", and the main goal was to extend the burnup capability from 2D assembly-level calculations all the way to full-core problems consisting of hundreds of thousands of depletion zones, without any limitations in parallelization. This goal was achieved by introducing different optimization modes for small and large burnup calculation problems (Leppänen, 2012a), and using shared-memory techniques (OpenMP) for parallelization.

Serpent 2 was made available to licensed users for beta-testing purposes in January 2012. Updating and active maintenance of Serpent 1 was discontinued soon after. The public release of Serpent 2 has been postponed several times. The code is no longer considered a "beta version", but the distribution is still handled by a separate software license agreement available to registered users of Serpent 1. A provisional commercial license was made available in 2016.

Current status and future plans

Serpent development has been funded from several national SAFIR research programmes, most recently SAFIR2018. Additional funding was received from the HPMC project in the EU 7th Framework Programme in 2011-2014, and the Academy of Finland NUMPS project in 2012-2016 (Leppänen, 2015c).

The Serpent user community consists of more than 600 users in 172 organizations and 37 countries around the world. The most typical Serpent user is a university student applying the code for academic research and thesis work. This is also seen in the large number of publications and theses on Serpent development and applications. More than 400 peer-reviewed journal articles and conference papers and 100 Bachelor, Master's and Doctoral-level theses have been published on Serpent-related topics worldwide since 2005. Serpent is also used in several research organizations, and in recent years the code has been adopted by the nuclear industry and small and medium-sized companies working on innovative reactor concepts.

The development of Serpent 2 continues, and the work is currently divided between three main topics:

1) Advanced methods for spatial homogenization
2) Coupled multi-physics simulations
3) New applications involving fusion neutronics and radiation transport problems

The development of advanced methods for spatial homogenization continues the work started in the beginning of the Serpent project in 2004. The code is current capable of producing all group constants required for steady-state fuel cycle simulations transient analyses. The methodology has been tested with various nodal diffusion codes developed and used at VTT, such as HEXBU-3D, HEXTRAN, TRAB3D and ARES, as well as the Apros system code. Serpent has also been actively used for generating group constants for codes like DYN3D and PARCS (see various publications by Serpent users).

The second major topic for on-going and future work is the development of calculation capability for coupled multi-physics applications. This approach represents a "brute-force" solution to the coupled core physics problem, and even though full-scale simulations for routine design and safety analyses are beyond the capabilities of today's computers, the capability to obtain rigorous high-fidelity solutions for neutronics and thermal hydraulics are important for several research applications. The coupling scheme relies on a universal multi-physics interface, which handles the transfer of power distributions and state point information between the different solvers with minimal modifications in the original input files. Serpent has been successfully coupled to several thermal hydraulics, CFD and fuel performance codes at VTT and within the user community.

The work on multi-physics coupling lead to the development of a photon transport mode and advanced CAD and mesh-based geometry types, which has allowed broadening the scope of Serpent applications to new fields beyond reactor physics. New applications are found in particular in fusion research, and Serpent has already been used for modeling complicated experimental devices, such as ITER and JET. Examples of fusion applications include energy deposition, material damage and activation, tritium breeding and shut-down dose rate calculations. Radiation transport and shielding applications are also found in industry, medical physics and space research.

The main part of code development is carried out by the Serpent team at VTT, but several user organizations have also made significant contributions to the development effort. Important on-going development topics include:

Coupled neutron-photon transport mode and new physics models for the production of secondary photons
Weigh-window based variance reduction techniques
Verification of CAD and unstructured mesh based geometry types
Collision history based sensitivity and perturbation calculations
Functional expansion tallies

Feedback and new ideas are greatly appreciated, and all users are encouraged to interactive communication at the Serpent Discussion Forum.


The figures below are produced by the mesh plotter in the Serpent code. The capability can be used in reactor physics applications for visualizing the effect of neutron thermalization on fission rate. The "cold" shades show the relative thermal flux distribution and the "hot" shades the relative fission rate (power) distribution. Bright and dark colors indicate high and low values, respectively. Serpent 2 extends the capability to tally results and multiple color schemes. The output is written in PNG graphics format files. The figures also show some of the geometry modeling capabilities in the Serpent code.

Fig 1. A standard 17 by 17 PWR Fuel assembly. Burnable absorber pins with low fission rate are distinguished in dark red.  Control rods are inserted in the guide tubes, which creates a local depression in the thermal flux. Infinite 2D lattice model.

Fig 2. A BWR fuel assembly with burnable absorber pins. Void fraction in the coolant channel is 60% and the thermal flux is clearly peaked inside the central water cross and outside the channel wall. Infinite 2D lattice model.

Fig 3. VVER-440 fuel lattice surrounding a flux-trap control element. Thermal flux is peaked in the water region where neutrons are thermalized and trapped inside the element. The boron steel absorber is almost black. Infinite 2D lattice model.

Fig 4. A PWR MOX fuel assembly surrounded by UOX assemblies in a reactor lattice. Infinite 2D lattice model.

Fig 5. A heavy water-cooled and moderated CANDU fuel cluster. The example demonstrates the circular array lattice type in the Serpent code. Infinite 2D lattice model.

Fig 6. The VENUS-2 reactor core. The VENUS-2 reactor dosimetry benchmark showed that PSG, the predecessor of the Serpent code can be used for modeling complicated geometries, but the poor efficiency of the collision flux estimator outside the active source limits the applicability in detector calculations (Leppänen, 2006). Full 3D geometry model.

Fig 7. The CROCUS reactor core. Two nested lattices with different pin pitch. Full 3D geometry model.

Fig 8. VVER-440 local power peaking experiment in the LR0 reactor (Leppänen, 2007a). Full 3D geometry model.

Fig 9. Depletion of a BWR fuel assembly with burnable absorber (7.4 MB gif animation). Infinite 2D lattice model.

Fig 10. A total of 15,000 microscopic TRISO fuel particles randomly dispersed inside a PBMR type fuel pebble. Calculated using the explicit particle fuel model in the Serpent code. Infinite 3D cubical array.

Fig 11. A prismatic HTGR fuel lattice with control rod inserted in the central fuel block. Infinite 2D lattice model.

Fig 12. HTTR (High-Temperature Test Reactor) 3D full-core calculation (6 mesh plots of critical and shut-down configurations).

Fig 13. INL Advanced Test Reactor (ATR) core. 3D full-core calculation (geometry and 5 mesh plots produced by Serpent 2).

Fig 14. "Stanford Critical Bunny" -- reactor physics adaptation of a well-known 3D computer graphics test model (high-enriched uranium bunny). Modeled using the new unstructured tetrahedral mesh based geometry type, currently under development in Serpent 2. Geometry model visualized using the utilities in the OpenFOAM toolbox. Fission rate distribution calculated using Serpent 2 (18MB gif animation).