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
three-dimensional, scalable to an arbitrary level of spatial resolution,
and capable of modeling the transport physics without major approximations
or application-specific 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,
in computer capacity and high-performance parallel computing.
The first Monte Carlo
burnup calculation codes were introduced in the late 1990's,
and production of input parameters for reduced-order
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 continuous-energy
Monte Carlo code, specialized in reactor physics applications, and capable
homogenized group constants for VTT's in-house nodal diffusion codes.
Considerable effort was devoted to optimizing the transport routine
for assembly-level calculations, to overcome the challenges of excessively
high computational cost.
The original working
title for the code was "Probabilistic
Scattering Game", or PSG. All publications
dated before the official pre-release 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
In an effort to simplify the calculation routines, it was
decided not to use the
continuous-energy cross sections as-is, 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 iterations were reduced to minimum,
and material-wise macroscopic cross sections could be pre-calculated before the transport simulation.
Another feature, originally implemented for the sake of
simplicity, was the Woodcock delta-tracking method
used for neutron transport. The standard delta-tracking 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 mean-free-path is often long compared to the geometry dimensions.
Both the unionized energy grid approach and the delta-tracking based transport routine are still used
in Serpent today.
After the first trials
turned out to be successful, the source
code was completely re-written to better accommodate
features and capabilities needed for lattice physics calculations. One of the
major additions in the new version was the support for
using the Message Passing Interface (MPI).
Development of the PSG code became a topic for a
completed in 2007.
The first publicly distributed version,
Serpent 1, began to formulate when the source code was completely re-written for a
second time by the beginning of 2008. The methodology was then complemented by built-in burnup
Radioactive decay and fission yield data was read from ENDF format
libraries, and combined with microscopic one-group 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.
Built-in 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
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
axis. This discovery 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 very short
The work was summarized in a doctoral thesis, completed in 2013.
A pre-release 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 non-commercial research and educational use.
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 sub-division increased the number of
burnable materials. Each material was associated with a set of pre-generated 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 assembly-level
Sooner or later it would have 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 the entire source code again. The
main goal was to extend the applicability
from 2D assembly-level calculations to 3D full-core problems
hundreds of thousands of depletion zones, without any limitations on
parallelization. This goal was achieved by introducing different
optimization modes for
small and large-scale burnup calculation problems, and using shared-memory 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 beta-testing
purposes, in January 2012.
Multi-physics 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 multi-physics interface in 2012 enabled coupling
Serpent to thermal hydraulics and fuel performance codes.
The coupling scheme relies on a super-imposed 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 multi-physics interface supports various distribution types, such as regular
mesh, 3D point cloud, user-defined 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 delta-tracking 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 on-the-fly 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
stand-alone neutron transport simulations to multi-physics
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 high-fidelity full-core calculations.
The work carried out between 2012 and 2017 was summarized in another
Development of a collision-based
domain-decomposition 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 multi-physics interface
were expanded to mesh- and CAD-based geometry types.
The capability to model complicated irregular structures using state-of-the-art
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.
on the neutronics calculation chain was completed in 2018.
The work on multi-physics applications also lead to the development of a
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 built-in 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 weight-window based
variance reduction scheme with a built-in response matrix method based
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 peer-reviewed journal articles and conference papers and 175 Bachelor, Master's and Doctoral-level
theses were published on Serpent-related 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 medium-sized 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 non-commercial 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:
||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 at the beginning of the Serpent project in 2004. Serpent is
capable of producing all group constants required for steady-state, 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 2017-2020.
The work was focused on
steady-state, burnup and transient simulations with two-way 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 HEXBU-3D, HEXTRAN and TRAB3D with a modern
Serpent-based calculation system. The "Kraken" core physics framework consists of modular neutronics, thermal hydraulics
and fuel behavior solvers, communicating via a versatile multi-physics interface. The neutronics solution in Kraken
is based on Serpent, used either as a high-fidelity solver, or as part of a reduced-order 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
third-party solvers, such as OpenFOAM. The outer boundary conditions can be further coupled to system-scale 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 multi-physics coupling.
Typical stand-alone applications include criticality safety, modeling of low-power 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 fusion-related transport
applications, such as tritium breeding, energy deposition, material damage, activation
and shut-down dose rate calculations. The support for CAD-based geometries has enabled the modeling of complicated
experimental devices, such as the ITER and JET tokamaks and the Wendelstein 7-X 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 back-end 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. Built-in
global variance reduction scheme enables evaluation of radiation doses in
large and heavily shielded geometries, such as transport casks and hot-cell
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
Collaboration with Studsvik Scandpower
In June 2022 VTT and Studsvik Scandpower signed a memorandum of understanding on
the development of a quality-assured version of Serpent, based on the
NQA-1-2015 standard. The aim of the collaboration is
to enable using Serpent for the licensing analyses of nuclear reactors,
inclunding SMRs and other emerging reactor technologies.
Studsvik and VTT also aim to expand the applications of Serpent into new
challenging fields, such as medical physics.
Serpent development continues with two
parallel branches. The NQA-certified version of Serpent will be integrated into
SSP's software family. Distribution of the non-NQA version is continued with the
current licensing options, including public distribution to
non-commercial research and educational use through the NEA Data Bank
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. User-defined 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. 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 3. 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 4. Depletion of a BWR fuel assembly with
burnable absorber (7.4
MB gif animation). Infinite 2D lattice model.
Fig 5. A pebble-bed 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 full-core calculation.
Fig 8. 1000 MWe PWR core (MIT BEAVRS benchmark), 3D full-core calculation.
Fig 9. Graphite-moderated High-Temperature Test Reactor (HTTR), 3D full-core calculation.
Fig 10. INL Advanced Test Reactor (ATR) core, 3D full-core
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 (SMR-scale). The calculation features the radioactive decay source mode and a Hybrid CSG/CAD -based geometry model. The STL-format CAD components are visualized using the Blender software (from Leppänen, 2022).
Fig 14. Fusion neutronics simulation in the Wendelstein 7-X stellarator. The geometry is for the most part based on complicated CAD models (from Äkäslompolo et al., 2021).