History of the Serpent project
The Monte Carlo method has been used for a wide range of radiation
and particle transport applications since the late 1940's. The method is inherently
threedimensional, scalable to an arbitrary level of spatial resolution,
and capable of modeling the transport physics without major approximations
or applicationspecific limitations. The major drawback of the Monte Carlo
method has always been its high computational cost.
Excessive running time may become a limiting factor when the transport simulation is run as part of an iterative loop, or when the calculation needs to be repeated multiple times with variations in the input conditions.
Such applications are common, for example, in nuclear reactor physics.
Up until the end of the 20th century, deterministic transport methods were the
only viable option for burnup calculations, group constant generation
and coupled core physics simulations.
New possibilities began to emerge by the end of the millennium,
along with
the advances
in computer capacity and highperformance parallel computing.
The first Monte Carlo
burnup calculation codes were introduced in the late 1990's,
and production of input parameters for reducedorder
deterministic calculations became a popular research topic.
Work on what later became the Serpent code started at VTT in 2004.
The main motivation was to develop a new continuousenergy
Monte Carlo code, specialized in reactor physics applications, and capable
of producing
homogenized group constants for VTT's inhouse nodal diffusion codes.
Considerable effort was devoted to optimizing the transport routine
for assemblylevel calculations, to overcome the challenges of excessively
high computational cost.
Prehistory
The original working
title for the code was "Probabilistic
Scattering Game", or PSG. All publications
dated before the official prerelease in October
2008 refer to the code using this name.
The first version of PSG was
started in September 2004. The
interaction physics was based on ACE format data libraries, mainly
because the ENDF reaction
laws were reasonably well documented. The same data format was
used by MCNP, which made the verification of the implementation
relatively straightforward.
In an effort to simplify the calculation routines, it was
decided not to use the
continuousenergy cross sections asis, 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 iterations were reduced to minimum,
and materialwise macroscopic cross sections could be precalculated before the transport simulation.
Another feature, originally implemented for the sake of
simplicity, was the Woodcock deltatracking method
used for neutron transport. The standard deltatracking algorithm uses rejection
sampling to avoid calculating the
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 well suited for lattice physics
calculations, in which the neutron meanfreepath is often long compared to the geometry dimensions.
Both the unionized energy grid approach and the deltatracking based transport routine are still used
in Serpent today.
After the first trials
turned out to be successful, the source
code was completely rewritten to better accommodate
features and capabilities needed for lattice physics calculations. One of the
major additions in the new version was the support for
parallel calculation
using the Message Passing Interface (MPI).
Development of the PSG code became a topic for a
doctoral thesis,
which was
completed in 2007.
Serpent 1
The first publicly distributed version,
Serpent 1, began to formulate when the source code was completely rewritten for a
second time by the beginning of 2008. The methodology was then complemented by builtin burnup
calculation capability.
Radioactive decay and fission yield data was read from ENDF format
libraries, and combined with microscopic onegroup transmutation cross sections
obtained from the neutron transport simulation. The formulation of transmutation
and decay chains was completely automated, and setting up the
burnup calculation problem required only minor input from the user.
Builtin burnup calculation
capability also required the development of internal routines
to solve the Bateman depletion equations.
The first implementation was based on the
linear chains method,
which practically implies
the analytical
solution of linearized depletion chains.
The solver routine was relatively easy to implement
using a recursive loop.
The development of an alternative matrix exponential
solution began soon after. Due to the
difficult numerical characteristics of the depletion problem, the computation
of the matrix exponential had not been previously
conceivable for the full nuclide system. However, analysis of the
mathematical properties of burnup matrices revealed that
their eigenvalues are generally confined near the negative
real
axis. This discovery prompted the introduction of a novel matrix exponential
method CRAM (Chebyshev Rational Approximation Method) for solving the
Bateman equations.
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 very short
computation time.
The work was summarized in a doctoral thesis, completed in 2013.
A prerelease version 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 OECD/NEA Data Bank, and released in May 2009. RSICC distribution to North America began one year later.
The code was licensed free of charge for noncommercial research and educational use.
Serpent 2
The methodology in Serpent 1 was implemented over a relatively
short period, without much systematic approach.
New calculation routines were constantly added
on top of existing source code, making the program structure
increasingly complicated and difficult to maintain.
It also became clear that
parallelization based on distributed memory MPI would eventually
become a major limitation, as the total
memory footprint was proportional to the number
of parallel tasks. Excessive memory usage was a problem especially in burnup
calculations. Depletion zone subdivision increased the number of
burnable materials. Each material was associated with a set of pregenerated macroscopic
cross sections, which required significant amounts of storage space.
Optimization tricks implemented to gain efficiency in lattice physics
calculations limited the practical use of Serpent to assemblylevel
burnup problems.
Sooner or later it would have 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 the entire source code again. The
main goal was to extend the applicability
from 2D assemblylevel calculations to 3D fullcore problems
consisting of
hundreds of thousands of depletion zones, without any limitations on
parallelization. This goal was achieved by introducing different
optimization modes for
small and largescale burnup calculation problems, and using sharedmemory techniques (OpenMP)
to distribute the neutron histories between CPU cores without increasing the
overall memory demand.
Serpent 2 was made available to licensed users of Serpent 1, initially for betatesting
purposes, in January 2012.
Multiphysics and new applications
After the public release of Serpent 2, the user community began to grow by some 100 new users each year.
New users, together with increasing computer capacity, lead to new applications and methodologies.
Development of a universal multiphysics interface in 2012 enabled coupling
Serpent to thermal hydraulics and fuel performance codes.
The coupling scheme relies on a superimposed data structure, that enables passing
temperature and density distributions into the transport simulation without any
modifications in the underlying geometry model. The same structures are used for
passing power distributions to the coupled solvers.
The multiphysics interface supports various distribution types, such as regular
mesh, 3D point cloud, userdefined functional dependence and
an unstructured polyhedral mesh in the OpenFOAM mesh format.
The separation of temperature and density data from the
geometry model was accomplished by implementing a deltatracking like
rejection sampling routine, which enables arbitrary,
even continuously varying distributions
for the thermal hydraulic state variables. Variation in material densities
was accomplished using fractional multipliers on the nominal value, but accounting for changes in
temperatures required an onthefly treatment for thermal motion
in neutron collisions. The target motion sampling method (TMS) was the topic of a doctoral thesis, completed in 2015.
Considerable effort was devoted to expanding the capabilities of Serpent from
standalone neutron transport simulations to multiphysics
applications. Methods and procedures were developed for
coupling Serpent to thermal hydraulics and fuel performance codes. New
features, such as a dynamic simulation mode with delayed neutron physics, were
implemented for highfidelity fullcore calculations.
The work carried out between 2012 and 2017 was summarized in another
doctoral thesis.
Development of a collisionbased
domaindecomposition scheme enabled running
even larger burnup calculation problems, and new methods were implemented
for sensitivity and uncertainty analyses.
Mesh routines developed for the OpenFOAM based multiphysics interface
were expanded to mesh and CADbased geometry types.
The capability to model complicated irregular structures using stateoftheart
engineering tools enabled
broadening the scope of Serpent applications to new fields beyond reactor physics.
Complicated geometries are common, for example, in fusion research.
Serpent was coupled to plasma scenario simulations, which provided a source term for the neutron transport simulation.
A
doctoral thesis
on the neutronics calculation chain was completed in 2018.
The work on multiphysics applications also lead to the development of a
photon transport
mode, which was needed to accurately account for the contribution of
gamma radiation in heat deposition
to coolant and structural materials. Expanding the physics routines
from neutrons to photons opened new possibilities in radiation transport.
The builtin burnup routine in Serpent could easily produce source terms for
calculations involving radiation emitted by spent nuclear fuel and materials activated under neutron irradiation.
A radioactive decay source mode combining material compositions to decay spectra
was implemented to automate the procedure.
Radiation transport capability was complemented by a weightwindow based
variance reduction scheme with a builtin response matrix method based
importance solver.
By 2020 the Serpent community had grown to more than
1000 users in 250 organizations in 44 countries around the world.
The most typical Serpent user was 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 750 peerreviewed journal articles and conference papers and 175 Bachelor, Master's and Doctorallevel
theses were published on Serpentrelated topics worldwide between 2005 and 2020. 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.
Current status and future plans
The distribution of Serpent 2 (base version 2.1) was linked to the OECD/NEA Data Bank and RSICC
licensed Serpent 1 until 2021. A new base version 2.2 was released in May 2022, making
the previous versions obsolete.
The use remains free of charge for noncommercial and educational purposes.
Serpent is also distributed under a commercial license by VTT.
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 at the beginning of the Serpent project in 2004. Serpent is
capable of producing all group constants required for steadystate, burnup and dynamic simulations carried out using typical nodal diffusion codes.
The methodology has been tested with various fuel cycle simulator and transient analysis codes developed and used at VTT.
Serpent has also been actively used for generating group constants
for other codes, such as DYN3D and PARCS.
Coupled calculation capability was extensively developed during the EU Horizon 2020 project McSAFE in 20172020.
The work was focused on
steadystate, burnup and transient simulations with twoway coupling to the SUBCHANFLOW thermal hydraulics code and
the TRANSURANUS fuel performance code.
Kraken core physics computational framework
In 2017 it was decided to gradually replace VTT's old legacy codes, such as HEXBU3D, HEXTRAN and TRAB3D with a modern
Serpentbased calculation system. The "Kraken" core physics framework consists of modular neutronics, thermal hydraulics
and fuel behavior solvers, communicating via a versatile multiphysics interface. The neutronics solution in Kraken
is based on Serpent, used either as a highfidelity solver, or as part of a reducedorder calculation sequence producing
group constants for the new nodal neutronics code Ants. Other modular components developed at VTT include
the FINIX fuel behavior code and the Kharon thermal hydraulics solver. Kraken also supports coupling to
thirdparty solvers, such as OpenFOAM. The outer boundary conditions can be further coupled to systemscale simulations
carried out, for example, using TRACE or Apros.
Public release of Kraken
is currently under preparation.
Serpent development related to reactor physics applications is currently carried out within the scope of
the Kraken framework, but the code can also be used without any multiphysics coupling.
Typical standalone applications include criticality safety, modeling of lowpower experimental
reactors and radionuclide inventory calculations.
Applications beyond reactor physics
In recent years, Serpent has also been adopted by the fusion community.
Serpent can be used for various fusionrelated transport
applications, such as tritium breeding, energy deposition, material damage, activation
and shutdown dose rate calculations. The support for CADbased geometries has enabled the modeling of complicated
experimental devices, such as the ITER and JET tokamaks and the Wendelstein 7X stellarator.
Capability to model complicated irregular geometries can also be a major advantage in
radiation transport and shielding calculations. The most natural applications for Serpent
are found at the backend of the fuel cycle, for example, in the transportation, storage
and handling of spent fuel assemblies, or in the decommissioning of old reactors.
Generation of the radiation source
term can be easily included in the calculation sequence. Builtin
global variance reduction scheme enables evaluation of radiation doses in
large and heavily shielded geometries, such as transport casks and hotcell
and spent fuel encapsulation facilities.
The same methods can be used for radiation transport and shielding problems
encountered in the industry, medical physics, accelerator and
space applications.


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. Userdefined detector scores can also be plotted using 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. 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 3. 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 4. Depletion of a BWR fuel assembly with
burnable absorber (7.4
MB gif animation). Infinite 2D lattice model.
Fig 5. A pebblebed type
HTGR fuel pebble containing
15,000 randomly dispersed
TRISO fuel particles. The geometry features the explicit
particle fuel model in Serpent. Infinite 3D cubical array model.
Fig 6. A prismatic HTGR fuel lattice with
control rod inserted in the central fuel block. Infinite 2D lattice model.
Fig 7. Hexagonal SMR core. 3D fullcore calculation.
Fig 8. 1000 MWe PWR core (MIT BEAVRS benchmark), 3D fullcore calculation.
Fig 9. Graphitemoderated HighTemperature Test Reactor (HTTR), 3D fullcore calculation.
Fig 10. INL Advanced Test Reactor (ATR) core, 3D fullcore
calculation.
Fig 11. Jules Horowitz (JHR) material test reactor, currently under construction in Cadarache, France.
3D coupled neutron/photon transport simulation.
Fig 12. Gamma radiation shielding calculation performed for a spent fuel storage cask. The calculation features the radioactive decay source mode, a hybrid CSG/CAD based geometry model and an iterative global variance reduction (GVR) scheme (from Häkkinen, 2018).
Fig 13. Photon transport simulation performed for a spent fuel storage rack (SMRscale). The calculation features the radioactive decay source mode and a Hybrid CSG/CAD based geometry model. The STLformat CAD components are visualized using the Blender software (from Leppänen, 2022).
Fig 14. Fusion neutronics simulation in the Wendelstein 7X stellarator. The geometry is for the most part based on complicated CAD models (from Äkäslompolo et al., 2021).
