Why the Monte Carlo method?
The continuousenergy 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 threedimensional
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 multiphysics 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 multistage 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 assemblyspecific
multigroup constants, and
2) Core simulation, where the fullscale neutronics solution obtained using diffusion
theory or other reducedorder 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 fullcore 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 threedimensional, 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 generalpurpose 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 continuousenergy 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 prerelease 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 generalpurpose codes available at that time were not
particularly well suited for lattice physics applications. Group constant generation requires the combination of several
userdefined 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.
Prehistory
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
continuousenergy 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
timeconsuming grid search iteration was reduced to minimum.
Another feature, originally implemented for the sake of
simplicity, was the Woodcock deltatracking 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 deltatracking 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 fewgroup 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 rewritten by the end of 2005. The new
version included a boundatom 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 deltatracking in
simple lattice geometries.
New features, such as
temperatureinterpolated material mixtures and a
geometry routine for modeling randomlydispersed 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 inhouse code at VTT.
The project also spawned the development of
MORA ("Monte Carlo Reactor Analysis"),
a simplified fullcore reactor physics code based on a
homogenized multigroup 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 rewritten for a
second time by the beginning of 2008. A persistent
methodological flaw in the freegas 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 deltatracking 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
(RantaAho, 2007). It was soon decided, however, that
the code should also have the capability to run burnup
calculation as a complete standalone application. In order
to accomplish this task without additional user effort, the calculation
of microscopic onegroup transmutation cross sections and the collection of
energydependent 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.
Builtin 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 12001700
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 prerelease 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 "liveandlearn" 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
multicore CPU's and massive parallelization.
In September 2010 it was decided that the best solution to the problems
in Serpent 1 was to rewrite entire source code again. The work
started under working title
"SuperSerpent", and the main goal was to extend the burnup capability
from 2D assemblylevel calculations all the way to fullcore 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 sharedmemory techniques (OpenMP)
for parallelization.
Serpent 2 was made available to licensed users for betatesting
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 20112014, and the Academy of Finland NUMPS project in 20122016
(Leppänen, 2015c).
The Serpent user community consists of almost
800 users in 196 organizations and 40 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 peerreviewed journal articles and conference papers and 100 Bachelor, Master's and Doctorallevel
theses have been published on Serpentrelated 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 mediumsized 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 multiphysics 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 steadystate fuel cycle simulations transient analyses.
The methodology has been tested with various nodal diffusion codes developed and used at VTT, such as HEXBU3D, 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 ongoing and future work is the development of calculation capability for coupled multiphysics applications. This approach represents
a "bruteforce" solution to the coupled core physics problem, and even though fullscale simulations for routine design and safety analyses are beyond the
capabilities of today's computers, the capability to obtain rigorous highfidelity solutions for neutronics and thermal hydraulics are important for several
research applications. The coupling scheme relies on a universal multiphysics 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 multiphysics coupling lead to the development of a photon transport mode and advanced CAD and meshbased 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 shutdown 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 ongoing development topics include:
• 
Coupled neutronphoton transport mode and new physics models for the production of secondary photons

• 
Weighwindow 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.


Gallery
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. VVER440 fuel
lattice surrounding a fluxtrap 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 watercooled
and moderated CANDU fuel cluster. The example demonstrates
the circular array lattice type in the Serpent code. Infinite 2D lattice model.
Fig 6. The VENUS2 reactor
core. The VENUS2 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. VVER440 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
(HighTemperature Test Reactor) 3D fullcore calculation (6
mesh plots of critical and shutdown configurations).
Fig 13. INL Advanced Test Reactor (ATR) core. 3D fullcore
calculation (geometry and 5 mesh plots produced by Serpent 2).
Fig 14. "Stanford Critical Bunny"  reactor physics adaptation of a wellknown 3D computer graphics test model (highenriched 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).
