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.
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
History of the Serpent project
development of Serpent started at VTT in 2004, under the working
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
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
The first version of PSG was
started in September 2004. Already at the beginning it was decided
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.
not widely used for neutron transport applications due to
certain limitations, this method
turned out to be reasonably well suited for lattice physics
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.
results with practical significance were
obtained in early 2005.
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
The first years
After getting the green light from VTT for continuing the development,
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
(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
completed one year later, in June 2007
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
but the project has not been continued since.
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
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
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
axis, which prompted the introduction of a novel matrix exponential
method CRAM (Chebyshev Rational Approximation Method) for solving the
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
The work is summarized in a Doctoral thesis, completed in 2013
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.
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
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)
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
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 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:
||Advanced methods for spatial homogenization
||Coupled multi-physics simulations
||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
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).