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, along with the advances 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 of producing 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.

Pre-history

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 relatively straightforward.

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 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 re-written for a second time by the beginning of 2008. The methodology was then complemented by built-in burnup calculation capability. 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 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 1200-1700 nuclides with an unprecedented accuracy and very short computation time. 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.

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 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 burnup problems. 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 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 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 doctoral thesis. 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. A doctoral thesis on the neutronics calculation chain was completed in 2018.

The work on multi-physics 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 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 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 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 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 space applications.

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 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 (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).