# FESTIM: An open-source code for hydrogen transport simulations

Rémi Delaporte-Mathurin, James Dark, Gabriele Ferrero, Etienne A Hodille, Vladimir Kulagin, Samuele Meschini

► To cite this version:

Rémi Delaporte-Mathurin, James Dark, Gabriele Ferrero, Etienne A Hodille, Vladimir Kulagin, et al. FESTIM: An open-source code for hydrogen transport simulations. International Journal of Hydrogen Energy, 2024, 63, pp.786 - 802. 10.1016/j.ijhydene.2024.03.184. cea-04765555

HAL Id: cea-04765555

https://cea.hal.science/cea-0476555v1

Submitted on 4 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# FESTIM: An open-source code for hydrogen transport simulations

Rémi Delaporte-Mathurin<sup>a,*</sup>, James Dark<sup>b,c</sup>, Gabriele Ferrero<sup>a,d</sup>, Etienne A. Hodille<sup>b</sup>, Vladimir Kulagin<sup>e</sup>, Samuele Meschini<sup>a</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

FESTIM
Hydrogen transport
Finite element
Modelling

## ABSTRACT

FESTIM (Finite Element Simulation of Tritium In Materials), is a versatile open-source finite element code developed in Python for simulating hydrogen transport in materials. FESTIM addresses limitations observed in existing codes by enabling multi-dimensional, multi-material simulations, leveraging the flexible finite element method and the open-source FEniCS library. Use cases illustrating FESTIM's applicability and efficacy are presented: reproduction of thermo-desorption experiments, modelling plasma-facing components (divertor monoblocks), and modelling breeding blankets. A comparative analysis with two other numerical tools (TMAP8 and COMSOL®) is then performed. FESTIM is in very good agreement with these codes and shows similar (or better) computing performances. Finally, the thorough validation and verification (using the method of manufactured solutions) of the code is described, highlighting the code's reliability.

## 1. Introduction

Hydrogen transport modelling is a discipline of great interest in various fields, including materials science [1,2], energy storage [3], nuclear fusion [4,5], and transportation [6]. However, existing tools in these domains have been hampered by inherent limitations, such as focusing mainly on 1D simulations, struggling with handling multiple materials, and being written in outdated programming languages. They also lack openness to collaborative advancement.

For these reasons, FESTIM (Finite Element Simulation of Tritium In Materials) was developed. It was designed to overcome these challenges and provide a simpler, more flexible approach to modelling hydrogen transport. Our primary goal with FESTIM is to make simulating hydrogen migration easier and more user-friendly. Additionally, FESTIM simulates heat transfer for better temperature accuracy.

What makes FESTIM unique is that it uses the flexible finite element method and is built on the FEniCS library [7], which is open source and has a Python interface, making it accessible to many users. The open-source nature of the code means anyone can use and adapt it to their needs. The Apache-2.0 license (under which FESTIM is distributed) allows for flexibility in how the code is used. This also encourages contributions to the code: FESTIM now counts more than 14 contributors from various research institutions and private companies.

Different from other tools like MHIMS [8], TESSIM [9], TMAP7 [10], CRDS [11], mHIT [12], or FACE [13], FESTIM stands out because it is open-source and can handle simulations in multiple dimensions and materials. Its user-friendly design works for simple and complex scenarios, from 1D to 3D, and it covers heat transfer and material interfaces well. Two other open-source tools are also capable of multidimensional and multimaterial simulations: TMAP8 [14] and ACHLYS [15]. These tools are based on the finite element code MOOSE [16]. Efforts have also been made to simulate hydrogen transport using commercial software like Abaqus [17] and COMSOL® [18].

FESTIM is used in academia but also in private companies – going even beyond nuclear fusion – with applications from experimental data analysis to the design and analysis of components interacting with hydrogen.

As FESTIM gains more attention, this paper is a guide to understanding its capabilities and why it is becoming a go-to tool for simulating hydrogen transport.

The paper first explores FESTIM's main features (v1.1.1) and underlying theory (see Section 2) and then delves into its development practices (see Section 3). Following this, some use cases are presented (see Section 4), and a comparative analysis with two other tools is provided (see Section 5). Finally, the paper details the thorough validation and verification processes for FESTIM (see Section 6). The source

<footnote><sup>a</sup> Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote><sup>b</sup> CEA, IRFM/GCFPM, F-13108 Saint-Paul-lez-Durance, France</footnote>

<footnote><sup>c</sup> Université Sorbonne Paris Nord, Laboratoire des Sciences des Procédés et des Matériaux, LSPM, CNRS, UPR 3407, F-93430, Villetaneuse, France</footnote>

<footnote><sup>d</sup> Politecnico di Torino, Dipartimento di Energia “Galileo Ferraris”, Corso Duca degli Abruzzi 24, Turin, Italy</footnote>

<footnote><sup>e</sup> National Research Nuclear University MEPhI, Moscow, 115409, Russian Federation</footnote>

<footnote>* Corresponding author.
E-mail address: remidm@mit.edu (R. Delaporte-Mathurin).</footnote>
**Symbols**

<table><tr><td>D</td><td>Diffusion coefficient (m<sup>2</sup> s<sup>-1</sup>)</td></tr><tr><td>J</td><td>Local particle flux (m<sup>-2</sup> s<sup>-1</sup>)</td></tr><tr><td>K<sub>H</sub></td><td>Henry's constant (solubility) (m<sup>-3</sup> Pa<sup>-1</sup>)</td></tr><tr><td>K<sub>S</sub></td><td>Sievert's constant (solubility) (m<sup>-3</sup> Pa<sup>-1/2</sup>)</td></tr><tr><td>K<sub>d</sub></td><td>Dissociation coefficient (m<sup>-2</sup> s<sup>-1</sup> Pa<sup>-1</sup>)</td></tr><tr><td>K<sub>r</sub></td><td>Recombination coefficient (m s<sup>-1</sup> or m<sup>4</sup> s<sup>-1</sup>)</td></tr><tr><td>P</td><td>Pressure (Pa)</td></tr><tr><td>Q</td><td>Heat volumetric source term (W m<sup>-3</sup>)</td></tr><tr><td>Q*</td><td>Soret coefficient (J mol<sup>-1</sup>)</td></tr><tr><td>S</td><td>Particle source term (m<sup>-3</sup> s<sup>-1</sup>)</td></tr><tr><td>T</td><td>Temperature (K)</td></tr><tr><td>c<sub>m</sub></td><td>Mobile concentration (m<sup>-3</sup>)</td></tr><tr><td>c<sub>t,i</sub></td><td>Particle concentration in the ith trap (m<sup>-3</sup>)</td></tr><tr><td>λ</td><td>Heat conductivity (W m<sup>-1</sup> K<sup>-1</sup>)</td></tr><tr><td>ρ</td><td>Density (kg m<sup>-3</sup>)</td></tr><tr><td>n</td><td>Normal vector</td></tr><tr><td>c<sub>p</sub></td><td>Specific heat capacity (J kg<sup>-1</sup> K<sup>-1</sup>)</td></tr><tr><td>k</td><td>Trapping rate (m<sup>3</sup> s<sup>-1</sup>)</td></tr><tr><td>n</td><td>Trapping site density (m<sup>-3</sup>)</td></tr><tr><td>p</td><td>Detrapping rate (s<sup>-1</sup>)</td></tr><tr><td>t</td><td>Time (s)</td></tr></table>

code of the results presented in this article (scripts, visualisations, etc.) can be found at https://github.com/RemDelaporteMathurin/festim-review [19].

## 2. Key features

### 2.1. Bulk physics

#### 2.1.1. Basic hydrogen transport

The model developed by McNabb & Foster [20] is used to model hydrogen transport in materials in FESTIM. The principle is to separate mobile hydrogen $c_m$ and trapped hydrogen $c_t$ (in dimensions of amount of substance per volume). The diffusion of mobile particles is governed by Fick's law of diffusion where the hydrogen flux is

$$J = -D \nabla c_m \qquad (1)$$

where D is the diffusivity coefficient.

Each trap i is associated with a trapping and a detrapping rate $k_i$ and $p_i$, respectively, and a trap density $n_i$.

The temporal evolution of $c_m$ and $c_t$ is then given by:

$$\frac{\partial c_m}{\partial t} = \nabla \cdot (D \nabla c_m) + S - \sum \frac{\partial c_{t,i}}{\partial t} \qquad (2)$$

$$\frac{\partial c_{t,i}}{\partial t} = k_i c_m (n_i - c_{t,i}) - p_i c_{t,i} \qquad (3)$$

where S is a source of mobile hydrogen. In FESTIM, source terms can be space- and time-dependent. These simulate hydrogen implantation in materials, tritium generation from neutron interactions, etc.

These equations can be solved in Cartesian coordinates, but also in cylindrical and spherical coordinates. This is useful, for instance, when simulating hydrogen transport in a pipe or a pebble. FESTIM can solve steady-state hydrogen transport problems.

An initial condition can be given for transient problems. By default, all concentrations are initialised to zero.

#### 2.1.2. Conservation of chemical potential at interfaces

Continuity of local partial pressure P must be ensured at the interfaces between materials. In the case of a material behaving according to Sievert's law of solubility (metals), the partial pressure is expressed as:

$$P = \left( \frac{c_m}{K_S} \right)^2 \qquad (4)$$

where $K_S$ is the material solubility (or Sievert's constant).

In the case of a material behaving according to Henry's law of solubility, the partial pressure is expressed as:

$$P = \frac{c_m}{K_H} \qquad (5)$$

where $K_H$ is the material solubility (or Henry's constant).

Then two different interface cases can occur. At the interface between two Sievert or two Henry materials, the continuity of partial pressure yields:

$$\frac{c_m^-}{K_S^-} = \frac{c_m^+}{K_S^+} \qquad (6)$$

or

$$\frac{c_m^-}{K_H^-} = \frac{c_m^+}{K_H^+} \qquad (7)$$

where the exponents + and – denote both sides of the interface.

At the interface between a Sievert and a Henry material:

$$\left( \frac{c_m^-}{K_S^-} \right)^2 = \frac{c_m^+}{K_H^+} \qquad (8)$$

These equilibrium equations show that a difference in solubilities introduces a concentration jump at interfaces.

In FESTIM, the conservation of chemical potential is obtained by a change of variables [21]:

$$P = \begin{cases} \left( \frac{c_m}{K_S} \right)^2 & \text{in Sievert materials} \\ \frac{c_m}{K_H} & \text{in Henry materials} \end{cases} \qquad (9)$$

P is continuous at interfaces.

Eqs. (2) and (3) are then rewritten and solved for P. Note that the boundary conditions are also rewritten. Once solved, the discontinuous $c_m$ field is obtained from P and the solubilities by solving Eq. (9) for $c_m$.

#### 2.1.3. Soret effect

FESTIM can include the Soret effect [22] (also called thermophoresis, temperature-assisted diffusion, or even thermodiffusion) in hydrogen transport. The flux of hydrogen J is then written as

$$J = -D \nabla c_m - D \frac{Q^* c_m}{R_g T^2} \nabla T \qquad (10)$$

where $Q^*$ is the Soret coefficient (also called heat of transport), and $R_g$ is the gas constant. Few studies have included the Soret effect due to a lack of data for the Soret coefficient [23–25].

#### 2.1.4. Heat transfer

As many processes involved in hydrogen transport are thermally activated, an accurate representation of the temperature field is often required. To this end, FESTIM can solve a heat transfer problem governed by the heat equation:

$$\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (\lambda \nabla T) + Q \qquad (11)$$

where T is the temperature, $c_p$ is the specific heat capacity, $\rho$ is the material's density, $\lambda$ is the thermal conductivity and Q is a volumetric heat source. As for the hydrogen transport problem, the heat equation can be solved in steady state, and initial conditions can be given. In FESTIM, the thermal properties of materials can be arbitrary functions of temperature (typically polynomials).
## 2.2. Surface physics

Boundary conditions are required to fully pose the hydrogen transport problem or the heat transfer problem. Boundary conditions are separated into three categories: (1) enforcing the value of the solution at a boundary (Dirichlet's condition), (2) enforcing the value of the gradient of the solution (Neumann's condition), and (3) enforcing the value of the gradient as a function of the solution itself (Robin's condition).

In FESTIM, users can fix the mobile hydrogen concentration $c_m$ and the temperature $T$ at the boundaries $\partial\Omega$ (Dirichlet):

$$c_m = f(x, y, z, t) \quad \text{on } \partial\Omega \qquad (12)$$

$$T = f(x, y, z, t) \quad \text{on } \partial\Omega \qquad (13)$$

where $f$ is an arbitrary function of spatial coordinates $x, y, z$ and time $t$.

FESTIM has built-in Dirichlet boundary conditions for Sievert's condition and Henry's condition (see Eqs. (14) and (15), respectively).

$$c_m = K_S \sqrt{P} \quad \text{on } \partial\Omega \qquad (14)$$

$$c_m = K_H P \quad \text{on } \partial\Omega \qquad (15)$$

where $P$ is a function of space and time.

Dirichlet boundary conditions can also approximate plasma implantation in near-surface regions to be more computationally efficient [26].

$$c_m = \frac{\varphi_{\text{imp}} R_p}{D} \quad \text{on } \partial\Omega \qquad (16)$$

where $\varphi_{\text{imp}}$ is the implanted flux and $R_p$ is the implantation range. Molecular recombination can be included by setting:

$$c_m = \frac{\varphi_{\text{imp}} R_p}{D} + \sqrt{\frac{\varphi_{\text{imp}}}{K_r}} \quad \text{on } \partial\Omega \qquad (17)$$

where $K_r$ is the recombination coefficient.

Moreover, dissociation can also be accounted for:

$$c_m = \frac{\varphi_{\text{imp}} R_p}{D} + \sqrt{\frac{\varphi_{\text{imp}} + K_d P}{K_r}} \quad \text{on } \partial\Omega \qquad (18)$$

where $K_d$ is the dissociation coefficient and $P$ is the partial pressure of hydrogen.

One can also impose hydrogen or heat fluxes at boundaries (Neumann). In the following, the Soret effect (see Eq. (10)) is not included for simplicity, and $J = -D\nabla c_m$.

$$J \cdot \mathbf{n} = -D\nabla c_m \cdot \mathbf{n} = f(x, y, z, t) \quad \text{on } \partial\Omega \qquad (19)$$

$$-\lambda \nabla T \cdot \mathbf{n} = g(x, y, z, t) \quad \text{on } \partial\Omega \qquad (20)$$

where $\mathbf{n}$ is the normal vector of the boundary.

Recombination and dissociation fluxes can also be applied:

$$J \cdot \mathbf{n} = -D\nabla c_m \cdot \mathbf{n} = K_d P \qquad (21)$$

$$J \cdot \mathbf{n} = -D\nabla c_m \cdot \mathbf{n} = -K_r c_m^{(1,2)} \qquad (22)$$

where $K_d$ is the dissociation coefficient and $K_r$ is the recombination coefficient. In Eq. (22), the exponent of $c_m$ is either 1 or 2, depending on the reaction order. These are Robin boundary conditions since the gradient is imposed as a function of the solution.

Finally, convective heat fluxes can be applied to boundaries:

$$-\lambda \nabla T \cdot \mathbf{n} = h (T - T_{\text{ext}}) \quad \text{on } \partial\Omega \qquad (23)$$

where $h$ is the heat transfer coefficient and $T_{\text{ext}}$ is the fluid's external (bulk) temperature.

## 2.3. Numerical methods

FESTIM is based on FEniCS [7] - a finite element library written in Python and C++. Engineering examples using FEniCS can be found in the literature [27–30]. A thorough description of the finite element method would be beyond the scope of this review but can be found in the FEniCS book [31]. In short, the finite element method relies on writing the strong formulation of a problem (here Eqs. (2), (3), and (11)) in its variational (or weak) formulation. The variational formulation is then discretised by restricting it to mesh elements representing the geometrical domain. In FEniCS, the linear algebra is handled under the hood by PETSc [32].

FESTIM provides some simple meshing capabilities for one-dimensional meshes. For more complex meshes, users should rely on specialised meshing software such as GMSH, Cubit, SALOME, etc. Meshes can then be converted to XDMF format (used by FEniCS and therefore FESTIM) with the package meshio [33]. The code is not restricted to simple geometries and can be used for CAD-based simulations (see Sections 4.2 and 4.3). Nonlinear equations are solved with FEniCS using Newton's method. Users can control the choice of linear solver. Direct (LU, UMFPACK, MUMPS, etc.) and indirect (Conjugate Gradient, GMRES, MINRES, BiCGSTAB, TFQMR, etc.) solvers are available via PETSc can be chosen by the user. Preconditioners for indirect solvers are also available (Algebraic Multigrid, incomplete LU factorisation, incomplete Choleski factorisation, Jacobi, etc.). By default, the Newton solver in FEniCS uses UMFPACK as a linear solver in serial and MUMPS in parallel (which are both sparse direct solvers). FESTIM has other minor features that make it a powerful and versatile tool. One example is the adaptive stepsize, allowing the stepsize to increase when the solver converges rapidly and reduces when the solver struggles to converge.

# 3. Code development

## 3.1. Version control and open source development workflow

The version control of FESTIM and its documentation are handled through the distributed version control system git. The source code for FESTIM is hosted and openly available on GitHub, which facilitates version control and code collaboration. GitHub is also used to track issues, manage contributions via pull requests, and enable a very effective development workflow and collaboration between developers. Thanks to GitHub, anyone can create their own public fork (a copy) of the FESTIM repository. This fork can then be modified by developers without needing write-access to the original repository. If a developer then wants to have their changes added to the official repository, they simply have to open a pull-request and the maintainers will review the changes and decide if they can be merged or not.

Having a code open source offers numerous benefits for code sustainability and community engagement. Firstly, open-source software fosters transparency, allowing users to inspect and understand the underlying code This transparency promotes trust and accountability, essential elements for sustainable development. This aspect is sometimes missing from existing tools for which only a compiled version is distributed. As a result, these tools are used as “black-boxes”. Additionally, open-source projects often attract a diverse community of developers, researchers, and enthusiasts who contribute to its improvement and evolution. This collaborative environment not only leads to more robust and reliable software but also facilitates knowledge sharing and skill development within the community. Moreover, open-source software is more adaptable to changing needs and technological advancements, ensuring its longevity and relevance over time. By empowering users to customise and extend the software to suit their specific requirements, open-source projects foster a sense of ownership and investment in the community, ultimately driving innovation and progress in the field. This open and collaborative approach has led to FESTIM attracting more than 14 contributors from more than ten different institutions.
## 3.2. Documentation

The extensive documentation of FESTIM is available online at https://festim.readthedocs.io/ and contains:

* a **Theory** section containing the physics basis of FESTIM

* a **User Guide** describing the key aspects of FESTIM: the installation instructions, how to create a simulation, the main objects, the post-processing steps, etc.

* a **Contributing Guide**, which describe the development workflow: how to contribute to the source code

* the **API reference** describing in detail the functions and classes in FESTIM

Release notes are also available on the GitHub repository. In addition to the documentation, a workshop is available at https://github.com/festim-dev/FESTIM-workshop. This workshop contains a series of documented tutorials going from very simple simulations to complex problems in FESTIM.

## 3.3. Continuous integration

FESTIM is extensively tested thanks to a test suite verifying various parts of the code, from unit tests to complete system tests. This test suite ensures reliability and is integrated with our Continuous Integration (CI) strategy, as it is automatically run every time a modification is made to the source code. Verification cases are also included in the CI to ensure the correctness of the code (see Section 6).

## 3.4. Distribution and licencing

The software is distributed under the permissive Apache-2.0 licence. FESTIM is packaged and distributed with PyPi, where specific versions of the code can be installed. Development versions of the code can also be installed in order to use the latest unreleased features of FESTIM.

## 3.5. Programming and dependencies

FESTIM is developed in Python, which is one of the most popular programming languages. This guarantees plenty of support from the ever-growing community. Having FESTIM written in Python allows it to easily leverage the power of the Python ecosystem and make use of libraries like NumPy [34], SciPy [35], matplotlib [36], sympy [37], etc. This does not come at the expense of performance since the finite element engine of FESTIM is written in fast C++. Only the user interface is written in user-friendly Python. The main dependency of FESTIM is the finite-element library FEniCS.

# 4. Gallery of use cases

This section will cover three published use cases of FESTIM.

The first use case is the reproduction and analysis of a thermo-desorption spectrum. These are typically 1D simulations of a few mil-limetres commonly with only one material and are used to characterise trapping sites in materials.

The second use case is the estimation of the retention and perme-ation flux in the ITER divertor. This case is more complex as it involves multi-dimensional modelling, multiple materials as well as a coupling with heat transfer. These models are also heavier since the modelled geometry is larger (a few centimetres).

The last case is a model of a tritium breeding blanket of the DEMO fusion reactor. It is also a multi-dimensional, multi-material and multi-physics problem. The added complexity is the introduction of a fluid domain (liquid metal) and integration of a tritium and heat source term from neutronics calculations.

The source code of the results presented in this section (scripts, visualisations, etc.) can be found at https://github.com/RemDelaporte Mathurin/festim-review [19]. Interested readers are invited to refer to the source code or the original papers for a more complete description of these use cases.

## 4.1. TDS analysis

Thermo-desorption experiments consist in loading a sample with hydrogen and then heating it with a well-controlled temperature ramp where the hydrogen desorption flux is measured resulting in a commonly called thermo-desorption spectroscopy (TDS) spectrum. These spectra are useful to characterise trapping sites and hydrogen retention in materials. Analysing these spectra analytically is extremely challenging. The reason for this is that, even if the problem is in most cases 1D and mono-material, solving the transient McNabb & Foster equations with one or more traps is complex even more so considering that the temperature and hydrogen sources are varying in time. It is a highly non-linear problem and therefore, numerical tools are needed. To infer trapping properties (rates, site density, etc.) from TDS data, a modelling tool like FESTIM is required. This section shows two examples of TDS simulations with FESTIM.

The first example is commonly used as a reference in the literature: the tungsten TDS performed by Ogorodnikova et al. [38]. In this experiment, 200 eV deuterium ions were implanted in a polycrystalline tungsten sample at room temperature with an incident flux of $2.5 \times 10^{19} \text{ D m}^{-2} \text{ s}^{-1}$ and a fluence of $10^{22} \text{ D m}^{-2}$. The sample was then heated at a rate of $8 \text{ K s}^{-1}$.

FESTIM reproduced the TDS spectrum with three traps (see Fig. 1). The first two traps are intrinsic traps (i.e. defects initially present in the sample), and the implantation of high-energy deuterium ions induces the third.

The second example is a TDS experiment by Baldwin et al. [39]. Beryllium and deuterium were codeposited on a tungsten sample. A 1 µm thick Be-D layer was formed at room temperature. The sample was then heated at a rate of $0.3 \text{ K s}^{-1}$ and the desorption flux of deuterium was measured. The authors simulated this experiment with TMAP7 and tested three sets of boundary conditions. It was shown that a recombination flux with a modified recombination coefficient was more appropriate for another set of experiments. The contribution of the residual pressure has also been modelled, although it can be shown that it does not influence the results. The spectrum was reproduced with FESTIM with two traps (see Fig. 2). In this model, as proposed by Baldwin et al. the implantation phase is not modelled, but rather an initial trap occupancy is set.

The work of Baldwin et al. [39] proposes a modified recombination coefficient that is dependent on the mobile concentration of hydrogen at the surface. This model is called *dynamically computed surface concentration* (DSC). Even though the DSC model is not natively supported in FESTIM (i.e. not present in the source code), it is very easy for users to quickly implement such a model by taking advantage of the object-oriented structure of the code. In this particular case, we created a custom class inheriting from the `RecombinationFlux` class and very easily implemented the DSC model (see Fig. 3).

FESTIM demonstrated at very early stages that it could easily reproduce TDS experiments [40]. Users can manually change the simulation parameters to fit the experimental data, or FESTIM can be embedded in an automated parametric optimisation routine to facilitate the process [41].

## 4.2. Divertor monoblock

One of the first applications of FESTIM was to simulate the transport of energetic hydrogen isotopes implanted from the ITER plasma in the divertor monoblocks. These components are made of a tungsten substrate of a few centimetres with a CuCrZr cooling pipe running through and a Cu inter-layer (see Fig. 4(b)). Before the development of FESTIM, the way hydrogen transport was simulated in these components was by making a 1D geometry, sometimes only simulating the tungsten domain. The reason for this is likely that the existing tools were limited by the finite difference method. A FESTIM study, leveraging the finite element method, showed that this 1D assumption lead to
<figure><img src="image_4.png" /><figcaption>Fig. 1. FESTIM simulation of a TDS experiment performed on tungsten by Ogorodnikova et al. [38]. The shaded areas are the contribution of each trap.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. FESTIM simulation of TDS experiment performed on Be-D co-deposited layers by Baldwin et al. [39] (recombination flux model). The shaded areas are the contribution of each trap.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Comparison of different surface models in the FESTIM simulation of the TDS experiment performed on Be-D co-deposited layers by Baldwin et al. [39].</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>(a) Toroidal cross-section at symmetry plane</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) 3D geometry showing toroidal and poloidal cross sections.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Poloidal cross-section at symmetry plane</figcaption></figure>

<caption>**Fig. 4.** Steady state temperature field and mobile hydrogen concentration field $c_m$ at $t = 10^6$ s. The streamlines are computed from the gradient of the mobile concentration field $c_m$. The jumps of $c_m$ are due to the interfaces between the different materials, as explained in Section 2.</caption>

wrong estimations of the tritium inventory and permeation fluxes in monoblocks exposed to a high heat and particle fluxes compared to a 2D approximation [40]. This justified the need for tools capable of multi-dimensional, multi-material and multi-physics simulations. The influence of the conservation of chemical potential at interfaces was then investigated [21]. It was shown that the concentration jumps

induced by these interface conditions do not impact the overall inventory as it is dominated by the tungsten. Another difficulty arose when scaling up the thousands of monoblock that make the divertor. To solve this, 2D FESTIM monoblock simulations have been used to build a surrogate model, which was then used to efficiently compute a monoblock inventory at a given time for a given set of exposure
<figure><img src="image_10.png" /><figcaption>Fig. 5. Total retention (mobile + trapped concentrations) in monoblock at t = 10⁶ s.</figcaption></figure>

conditions [42]. This surrogate model was then used to rapidly compute the inventory of the ITER and WEST divertors under various plasma scenarios [43]. An attempt at refining the monoblock model near the plasma-exposed surface has been made by coupling the code with MHIMS and introducing a *supersaturated layer* model [44], which successfully demonstrated that the impact of such a model was negligible as the inventory of the monoblock is dominated by the transport in the bulk.

More recently, the FESTIM monoblock model has been extended to 3D [45]. This study aimed to investigate the influence of molecular recombination on the poloidal gaps of the monoblocks. The geometry is a 6 mm × 23 mm × 25 mm block of tungsten with a copper cooling pipe running through and a CuCrZr inter-layer. The temperature field is computed from a steady state heat transfer problem (see Fig. 4). A heat flux of 10 MW m⁻² is imposed on the plasma-facing surface, and a convective heat transfer boundary condition is applied to the cooling surface. An approximation of the source of implanted particles (Eq. (16)) was set on the plasma-facing surface. Two traps are set in the tungsten subdomain and one in the copper and CuCrZr subdomains. This study showed that mobile particles desorb on this surface when assuming instantaneous recombination on the poloidal gap (see Fig. 4(a)), and the hydrogen retention is drastically reduced. The maximum total hydrogen concentration (including trapped hydrogen) after 10⁶ s of exposure was 1.2 × 10²² H m⁻³ (see Fig. 5).

### 4.3. Breeding blanket

Another application of FESTIM was modelling tritium transport in a Water-Cooled Lithium Lead (WCLL) tritium breeding blanket [46]. The WCLL is a concept to be tested in the ITER TBM programme, and the results will be used to determine which concept will be used in its successor, DEMO. The 2D WCLL model in the study comprised a tungsten first wall, a structural material of EUROFER, and a liquid LiPb tritium breeding material. The study aimed to investigate the influence of tritium trapping on a component scale, which had never been included in breeding blanket studies — likely because suitable tools were not available yet. It also showcases how FESTIM can be coupled with additional physics (see Fig. 6).

In a breeding blanket, there are two sources of hydrogen: the hydrogen coming from the fusion plasma and being implanted in the

tungsten first wall, and the tritium that is created by the interaction of lithium in the liquid metal with the fusion neutrons. Once bred in the liquid metal, tritium diffuses and is advected by the flow. In this case, a fluid dynamics solver (written in FEniCS) solving the Navier–Stokes equations was used to model the flow of the liquid metal LiPb around the breeder. The Navier–Stokes equations were solved in steady state assuming laminar flow for simplification, as the point of the model was to demonstrate the capability to couple FESTIM to fluid dynamics. Natural convection in the LiPb was also taken into account using the Boussinesq density approximation. The FESTIM governing equations had to be modified with an additional advection term $u \cdot \nabla c_m$, where $u$ is the velocity field.

$$ \frac{\partial c_m}{\partial t} = \nabla \cdot (D \nabla c_m) - u \cdot \nabla c_m + S \quad \text{in fluid domains} \quad (24) $$

$$ \frac{\partial c_m}{\partial t} = \nabla \cdot (D \nabla c_m) - \sum \frac{\partial c_{t,i}}{\partial t} \quad \text{in solid domains} \quad (25) $$

$$ \frac{\partial c_{t,i}}{\partial t} = k_i c_m (n_i - c_{t,i}) - p_i c_{t,i} \quad \text{in solid domains} \quad (26) $$

The governing equation for heat transfer was modified accordingly, although in this case the heat transfer Peclet number was very close to zero and therefore heat transfer was dominated by conduction. The reason why FESTIM does not have native support for advection is because extra care needs to be taken when dealing with problems which have both fluid and solid domains.

The tritium and heat source terms ($S$ and $Q$) were taken from OpenMC [47] simulations results that were fitted by a function of the distance from the first wall. While this approximation was satisfactory for this geometry, there are cases where approximating the tritium source term like this will be very challenging or even impossible, like the LIBRA tritium breeding experiment [48]. This is one of the areas of development of FESTIM where a better solution for coupling with OpenMC is needed. An effort towards strongly coupling OpenMC to finite element models has been undergone by the AURORA project [49].

With this implementation, tritium transport models were simulated to evaluate tritium inventories in the structural materials and permeation fluxes to the cooling channels. Additionally, tritium generation in the LiPb was modelled with a volumetric source term characterised from neutronics calculations [50].

This study indicated that the introduction of trapping mechanisms increased inventory values by 15% and delayed permeation to the cooling channels. A parametric study was conducted to explore uncertainties in the literature regarding hydrogen solubility in LiPb (spanning over two orders of magnitude). Varying lithium lead solubility within the literature range was found to cause a 25-fold change in EUROFER tritium inventory. Additionally, permeation fluxes to the cooling channels exhibited a threefold variation.

This work is being developed further by improving the fluid dynamics solver to account for MHD effects on the liquid metal breeder. Furthermore, the influence of neutron damage on trap concentrations in solid domains such as tungsten and EUROFER is being investigated.

## 5. Comparison with other codes

The source code of the results presented in this section can be found at https://github.com/RemDelaporteMathurin/festim-review [19]. Readers are invited to refer to the input scripts for the complete and exact parametrisation of these comparison cases.

### 5.1. TMAP8

TMAP8 is an acronym for Tritium Migration Analysis Program. It is the successor of TMAP7 [10] even though they are completely different codes. It is an open-source finite element code developed at INL that builds on the capabilities of the MOOSE finite element package. Several V&V cases are available on the TMAP8 repository. These cases are
<figure><img src="image_11.png" /><figcaption>Fig. 6. WCLL breeding blanket modelling workflow [46].</figcaption></figure>

inspired by the TMAP4 and TMAP7 V&V and their complete description can be found there [51]. The TMAP8 documentation states that the tool is designed for 0D-1D simulations, which may explain why 2D/3D verification cases are not yet available, even though TMAP8 has been used to simulate 3D problems [52]. For this study, we chose the several cases that were the most relevant: val-2b, ver-1a, ver-1c, ver-1d, and ver-1e. Other cases were available but were very simplistic (e.g. simple diffusion, steady state heat transfer, etc.). A FESTIM model has been created for each case with the same exact same physical parameters and meshes. The TMAP8 input scripts were taken from the TMAP8 code repository with minimal modifications. The TMAP8 version used was commit b84d41b. Both codes were run several times for each case, and the average execution time (on the same computer) was measured. For all cases, FESTIM showed similar or better performance (see Fig. 7). All FESTIM models were run with FEniCS default solver parameters for simplicity (see Section 2.3 for details). A full performances optimisation study of FESTIM would be outside the scope of this paper. Very minimal modifications of the TMAP8 scripts were made compared to what was available on the TMAP8 repository. For example, when needed, the simulation time or stepsize was modified to match FESTIM. Also, the default solver parameters were used for TMAP8: Newton solver with implicit-Euler time integration with ILU as a preconditioner for the Krylov solver.

### 5.1.1. val-2b

The val-2b case is a thermo-desorption experiment on beryllium with a beryllium oxide (BeO) layer. The sample is exposed to a partial pressure of hydrogen and then heated up after rest. To reproduce the experimental desorption curve, the strategy initially proposed by Longhurst and Ambrosek [51] is to modify the BeO diffusivity $D_{BeO}$ from the exposure phase to the desorption phase.

$$D_{BeO} = 1.4 \times 10^{-4} \exp(-2.10/k_B T) \quad \text{during exposure} \quad (27)$$

$$D_{BeO} = 7.0 \times 10^{-5} \exp(-2.33/k_B T) \quad \text{during desorption} \quad (28)$$

The diffusivity and solubility of Be as well as the solubility of BeO are kept constant:

$$D_{Be} = 8.0 \times 10^{-9} \exp(-0.36/k_B T) \quad (29)$$

$$K_{S,Be} = 7.156 \times 10^{27} \exp(-1.00/k_B T) \quad (30)$$

$$K_{S,BeO} = 5.00 \times 10^{20} \exp(-0.81/k_B T) \quad (31)$$

No traps are set in the materials. FESTIM and TMAP8 reproduce the experimental TDS (see Fig. 8).

For this specific case, the script available on the TMAP8 repository had to be modified in order to reproduce the experimental curve. The penalty term was reduced from 0.09 to 0.089, as suggested by the TMAP8 development team. Both codes also show slight differences despite using the same material parameters.

This validation case is far from being the most appropriate, for several reasons. First, the hydrogen exposition phase lasts for 180 015 s, whereas the temperature is only kept at 773 K for 180 000 s, after which it starts decreasing. These 15 s of delay are not very explicitly stated in the TMAP7 or TMAP8 V&V documentation, yet they are the very reason why there is an hydrogen intake in the sample. As the sample starts cooling down while the hydrogen pressure is maintained, the solubility of the BeO increases (because of a negative activation energy) and therefore hydrogen starts to dissolve in the sample. In other words, even though the sample is kept in a hydrogen atmosphere for 50 h, no hydrogen can penetrate before the temperature starts decreasing. This means that if this delay is not 15 s but say 20 s, the hydrogen intake will be much higher and the thermo-desorption spectrum will be very different.

A better TDS validation case would be the one described in Section 4.1 since the peaks observed are due to trapping sites in the material. Unfortunately, it is not possible in TMAP8 – to the best of our knowledge – to have a temperature-dependent trapping rate $k$, making it impossible to model such a case.

### 5.1.2. ver-1a

The ver-1a verification case [10] consists of an enclosure containing an initial quantity of hydrogen gas (at an initial pressure $P_0$). The hydrogen can dissociate on the inner surface of the enclosure and then permeate through the walls. As hydrogen particles escape the enclosure, the pressure decreases and as a result, so does the surface concentration on the inner walls. This 1D problem is, therefore, coupled. At each time step, the flux of particles escaping the enclosure is
<figure><img src="image_12.png" /><figcaption>Fig. 7. Comparison of TMAP8 and FESTIM average execution times.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 8. Comparison of val-2b case results from FESTIM and TMAP8.</figcaption></figure>

computed, and the internal pressure is updated. TMAP8 and FESTIM agree very well (see Fig. 9).

An analytical solution to this problem can be obtained. The analytical expression of the concentration in the wall is given by Longhurst et al. [10]:

$$c(x = 0, t) = 0 \qquad (32)$$

$$c(x = l, t) = K_H P(t) \qquad (33)$$

$$c(x, t) = 2LP_0K_H \sum_{n=1}^{\infty} \frac{e^{-Dt\alpha_n^2} \sin(x\alpha_n)}{(L+l(L^2+\alpha_n^2)) \sin(l\alpha_n)} \qquad (34)$$

where $L = \frac{K_H T A k_B}{V}$, $K_H$ is the solubility of the material, $T$ is the temperature, $A$ is the enclosure surface area, $V$ is the volume of the enclosure, $l$ is the thickness of the wall, and $\alpha_n$ are the roots of:

$$\alpha_n = \frac{L}{\tan(\alpha_n l)} \qquad (35)$$

The flux at the non-exposed surface can be expressed as:

$$\begin{aligned} J(x = 0, t) &= D \nabla c|_{x=0} \\ &= 2LP_0DK_H \sum_{n=1}^{\infty} \frac{e^{-Dt\alpha_n^2} \alpha_n}{(L+l(L^2+\alpha_n^2)) \sin(l\alpha_n)} \end{aligned} \qquad (36)$$

The cumulative release $R$ is then:

$$\begin{aligned} R &= \int_{0}^{t} J(x = 0, t) dt \\ &= 2LP_0 S \sum_{n=1}^{\infty} \frac{1 - e^{-Dt\alpha_n^2}}{(L + l(L^2 + \alpha_n^2)) \alpha_n \sin(l\alpha_n)} \end{aligned} \qquad (37)$$

Multiplying by the area and normalising by the initial quantity in the enclosure, one can obtain the fractional release:

$$\begin{aligned} \text{FR} &= \frac{R A}{P_0 V/(k_B T)} \\ &= \frac{2LSA}{V/(k_B T)} \sum_{n=1}^{\infty} \frac{1 - e^{-Dt\alpha_n^2}}{(L + l(L^2 + \alpha_n^2)) \alpha_n \sin(l\alpha_n)} \end{aligned} \qquad (38)$$

The discrepancies with the analytical solution is due to the coarse mesh (10 vertices). Increasing the mesh density ten-fold reduces the error significantly.

### 5.1.3. ver-1c

The ver-1c verification case is a 1D problem where a 100 m-long slab has an initial concentration of hydrogen of 1 H m⁻³ in the first 10 m only. The diffusivity $D$ is set to 1 m² s⁻¹, no trapping is included, and the surfaces are insulated (zero flux boundary condition). FESTIM and
<figure><img src="image_14.png" /><figcaption>Fig. 9. Comparison of the ver-1a case results from FESTIM and TMAP8 with the analytical solution. Temporal evolution of the cumulative hydrogen release.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 10. Comparison of ver-1c case results from FESTIM and TMAP8. Temporal evolution of the concentration of mobile particles at several locations.</figcaption></figure>

TMAP8 are in perfect agreement (see Fig. 10). FESTIM runs this case more than two times faster (see Fig. 7).

### 5.1.4. ver-1d

The ver-1d case is a 1D problem with one trap. It is the only available case with trapping included. The concentration of mobile particles $c_m$ is kept at $c_0$ on one side ($x = 0$) and set to zero on the other ($x = l$). The case is separated into a diffusion-limited problem and a trapping-limited problem. In the diffusion-limited case, the trap properties (trapping and detrapping rates) and the surface concentration of mobile particles are chosen so that the trap-filling ratio is much lower than unity. In this case, the governing equations can be approximated by a pure diffusion case with no traps but with an effective diffusion coefficient. This regime is also known as the effective diffusivity regime.

The particle flux $J$ at the back surface ($x = l$) is expressed in $H m^{-2} s^{-1}$ by:

$$J_{\text{analytical}}(t) = \frac{c_0 D}{l} \left[ 1 + 2 \sum_{m=1}^{\infty} (-1)^m \exp \left( -m^2 \frac{\pi^2 D_{\text{eff}} l}{l^2} \right) \right] \quad (39)$$

with $D_{\text{eff}}$ ($m^2 s^{-1}$) the effective diffusivity expressed as:

$$D_{\text{eff}} = \frac{D}{1 + \left( \frac{p}{k n} + \frac{c_m}{n} \right)^{-1}} \quad (40)$$

The mean squared error (MSE) with the analytical flux is computed for both codes. They are in extremely good agreement with the analytical solution (see Fig. 11) and the MSE for FESTIM and TMAP8 are small: $7.98 \times 10^{-6}$ and $6.26 \times 10^{-6}$, respectively. The differences may be attributed to different solvers and/or pre-conditioners, which could also explain the differences in execution times (see Fig. 7). It is worth noting that the TMAP8 model for this case scales the trapped concentration to improve convergence whereas FESTIM does not.

The trapping-limited version of the problem is very similar, but the trapping properties are such that the filling ratio is very close to unity. This case was not run, as the solution in TMAP8 is not satisfactory.

### 5.1.5. ver-1e

The ver-1e case is a 1D multi-material problem. The first layer of the composite structure is 33 µm-thick with a diffusivity of $1.274 \times 10^{-7} m^2 s^{-1}$ and the second layer is 66 µm-thick with a diffusivity of $2.622 \times 10^{-11} m^2 s^{-1}$. The conservation of chemical potential is not taken into account in this verification case, meaning that the concentration of mobile particles at the interface is continuous. The concentration of mobile particle is fixed at 50.7079 $H m^{-3}$ on the left surface.

The results produced by both codes are in very good agreement (see Fig. 12).
<figure><img src="image_16.png" /><figcaption>Fig. 11. Comparison of ver-1d case results from FESTIM and TMAP8. Temporal evolution of the particle flux.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 12. Comparison of ver-1e case results from FESTIM and TMAP8. Temporal evolution of the concentration of mobile particles at several locations.</figcaption></figure>

## 5.2. COMSOL multiphysics®

The multidimensional problem showed in Section 4.2 was solved with COMSOL Multiphysics® 6.1 [53] and compared to FESTIM (the complete details of the comparison problem are available as Supplementary data). The exact same meshes and sets of parameters were used.

COMSOL Multiphysics® is a finite element analysis software to analyse different physics (electric, mechanical, fluid, acoustic, and chemical), with an emphasis on how they can interact with each other (multiphysics). COMSOL Multiphysics® can model strongly coupled physics, including fluid dynamics and turbulent flow coupling with hydrogen transport. The software can simulate 1D, 2D, and 3D geometries in steady state and transient. COMSOL Multiphysics® is not open source and a purchased licence is required to run the software. To simulate hydrogen transport with trapping, modules *diffusion of diluted species, heat transfer in solids* and *PDE* (partial differential equations) are used. The diffusion module solves the Fickian diffusion equations. The heat transfer module in the monoblock study is used to generate a temperature distribution, and it is weakly coupled with hydrogen transport. A PDE has been employed to evaluate the trapping and

detrapping terms, which are included as hydrogen generation sources in the transport module.

The hydrogen concentrations and temperature fields produced by COMSOL® and FESTIM are in good agreement (see Fig. 13). Slight differences in the hydrogen concentration were found near the tungsten-copper interface (see Fig. 14) and in the tungsten bulk. These differences are attributed to lower mesh quality in these regions and/or different time-stepping schemes. Indeed, FESTIM was able to converge with a very large initial stepsize of 10⁵ s. On the other hand, COMSOL® required an initial stepsize of 2 × 10⁻⁴ s. For this reason, comparing performances is not trivial. Furthermore, COMSOL® and FESTIM did not use the same linear solver and do not treat interfaces the same way: FESTIM performs a change of variable (see Section 2.1.2 for more details), while COMSOL® has a partition interface condition implemented. In COMSOL®, the linear solver was PARDISO [54] whereas MUMPS was used in FESTIM. The absolute tolerance of the Newton solver was 10⁴ in FESTIM and automatically scaled in COMSOL®. The relative tolerance for both was 10⁻⁵. Complete parametrisation of the COMSOL® model can be found in the .mph file in the paper repository. The execution time of COMSOL® on four cores was about 25 min, whereas the execution time of FESTIM on a single core was less than 200 s. Moreover, since the problem source terms are constant in time,
FESTIM can be run with an initial stepsize of $1 \times 10^6$ s and produce very similar results in less than 20 s. Overall, one should not focus on the computational time of both codes since this is not a direct comparison (different linear solvers, different handling of interfaces discontinuities, different post-processing, etc.). Instead, the main output of this comparison is ensuring both FESTIM and COMSOL® produce the same results with the same physical and geometrical parameters. The various derived quantities, like inventories and particle fluxes at surfaces, were also in very good agreement (see Fig. 15).

## 6. Verification & validation

Before being used for analysis, any code should be verified & validated. The Verification and Validation (V&V) process has two goals: (1) to prove that the governing equations are correctly solved and that the code is error-free and (2) to demonstrate that the governing equations reproduce experimentally observed processes. In other words, verification is answering the question “Are we building the code right?” and validation is answering the question “Are we building the right code?”.

Many hydrogen transport codes in the literature have focused on validation, not verification (or sometimes verification in very simplistic scenarios). The most common validation case is to reproduce a thermo-desorption experiment. While validating a code is necessary, it is not sufficient. For example, let us imagine a code with a bug which adds 0.5 eV to all user-defined detrapping energies ($E_p$). This code could reproduce (fit) a thermo-desorption spectrum just fine: it could be validated. However, the equations solved in the code would not correspond to the real governing equations: it is not verified.

### 6.1. Validation

The McNabb & Foster model [20] has been used for decades to simulate hydrogen transport. It successfully reproduces thermodesorption experiments [8,11,13,38,55-62] and profilometry data [58,60,62-65]. Moreover, permeation experiments have also been reproduced using hydrogen diffusion models [55,56,66-69].

There is no need to repeat here the work that has been done in the past: this mathematical model is validated. It goes without saying that problems outside of this validation domain need additional experimental data.

### 6.2. Verification

The verification process of a code must be up to the code's claims. In other words, for instance, if a code claims to solve the McNabb & Foster model in 3D, in multimaterial cases, and coupled with heat transfer, then restricting the verification to a simple 1D, monomaterial, constant temperature case is not sufficient.

The FESTIM test suite contains dozens of verification cases and more than 280 tests in total (all publicly available on the FESTIM repository). A complete description of all these cases would require its own report; therefore, we will only describe three cases here.

The Method of Manufactured Solutions (MMS) is a mathematical technique used to verify numerical simulation codes [70]. In MMS, an exact analytical solution is intentionally manufactured. The mathematical model (Eqs. (2), (3), or (11)) is then operated on this manufactured solution to obtain an analytic source term. Computer symbolic manipulation can be used to obtain the derivatives of the manufactured solution. Boundary conditions obtained directly from the manufactured solution and the source term then serve as input for the numerical simulation code.

The error between the code output and the manufactured analytical solution is computed and used to assess the code's correctness and accuracy. MMS is a valuable tool for code verification in scientific and engineering computing, ensuring the reliability of computational results.

<figure><img src="image_18.png" /><figcaption>Fig. 13. Monoblock results computed by FESTIM (left) and COMSOL® (right).</figcaption></figure>

This section covers the verification of interface discontinuities, the heat transfer module, and the trapping implementation in FESTIM. For simplicity and clarity, only three distinct verification problems have been defined. All MMS cases in this section are applied on a unit square (for easier visualisation), and the mesh has 100 cells in
<figure><img src="image_19.png" /><figcaption>Fig. 14. Comparison of the mobile particle concentration $c_m$ and total retention profiles between FESTIM and COMSOL®.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 15. Comparison of computed inventories and particle fluxes between FESTIM and COMSOL®.</figcaption></figure>

the x direction and 100 cells in the y direction. Dirichlet boundary conditions enforcing the exact solution are applied on all surfaces. Only steady-state problems for simplicity's sake, but MMS can be – and has been – applied to transient problems. For each case, the L2 error $E$ is calculated in the domain $\Omega$:

$$E = \sqrt{\int_{\Omega} (u_{\text{exact}} - u_{\text{computed}})^2 dx} \qquad (41)$$

where $u_{\text{exact}}$ and $u_{\text{computed}}$ are the exact and computed solutions, respectively.

### 6.2.1. Case 1: discontinuity

The first MMS problem has two materials (denoted, respectively, by left and right). In material left, the solubility is $K_{S,\text{left}} = 3$ and the diffusivity is $D_{\text{left}} = 2$. In material right, the solubility is $K_{S,\text{right}} = 6$ and the diffusivity is $D_{\text{right}} = 5$. Two exact solutions for mobile concentration of hydrogen are manufactured for both subdomains:

$$c_{\text{left,exact}} = 1 + \sin(\pi(2x + 0.5)) + \cos(2\pi y) \qquad (42)$$

$$c_{\text{right,exact}} = \frac{K_{S,\text{right}}}{K_{S,\text{left}}} c_{\text{left,exact}} \qquad (43)$$

Note that the manufactured solutions were chosen so that the particle flux $J = -D\nabla c_m \cdot \mathbf{n}$ is continuous across the materials interface. MMS sources are derived in each material:

$$S_{\text{left}} = 8\pi^2 (\cos(2\pi x) + \cos(2\pi y)) \qquad (44)$$

$$S_{\text{right}} = 40\pi^2 (\cos(2\pi x) + \cos(2\pi y)) \qquad (45)$$

These exact solutions can then determine the MMS fluxes and boundary conditions. The computed and exact solutions agree very well (see Fig. 16(a)). The L2 error is $5.49 \times 10^{-4}$.

### 6.2.2. Case 2: heat transfer

This case verifies the implementation of the heat transfer solver in FESTIM. Two materials with different thermal conductivities are defined: $\lambda_{\text{left}} = 2$ and $\lambda_{\text{right}} = 5$.

The exact solution for temperature is:

$$T_{\text{exact}} = 1 + \sin(\pi(2x + 0.5)) + \cos(2\pi y) \qquad (46)$$

As for Case 1, the manufactured solution is chosen so that the thermal flux $-\lambda\nabla T \cdot \mathbf{n}$ is continuous across the interface.

Different MMS sources are defined (one in each material), and the computed solution agrees very well with the exact solution (see
<figure><img src="image_21.png" /><figcaption>Fig. 16. FESTIM verification results. For each case, line profiles of the 2D solution are plotted.</figcaption></figure>

Fig. 16(b)).

The L2 error is $3.31 \times 10^{-4}$.

$$Q_{\text{left}} = 8\pi^2 (\cos(2\pi x) + \cos(2\pi y))$$

(47) 6.2.3. Case 3: trapping

$$Q_{\text{right}} = 20\pi^2 (\cos(2\pi x) + \cos(2\pi y))$$

(48) Only one trap is considered for simplicity, but the same principle

The last case verifies the implementation of trapping in FESTIM.
applies for more traps. Exact solutions are defined for both the mobile concentration and the trapped concentration.

$$c_{m,exact} = 5 + \sin(2\pi x) + \cos(2\pi y) \qquad (49)$$

$$c_{t,exact} = 5 + \cos(2\pi x) + \sin(2\pi y) \qquad (50)$$

The trap density is $n = 2$ $c_{t,exact}$, the trapping rate is $k = 0.1$, the detrapping rate is $p = 0.2$, and the diffusivity is $D = 5$. MMS sources are obtained for both the mobile concentration and the trapped concentration:

$$S_m = -D\nabla^2 c_{m,exact} + kc_{m,exact}(n - c_{t,exact}) - pc_{t,exact} \qquad (51)$$

$$S_t = -kc_{m,exact}(n - c_{t,exact}) + pc_{t,exact} \qquad (52)$$

Again, the computed solution agrees very well with the exact solution (see Fig. 16(c)). The L2 error is $1.05 \times 10^{-2}$ for the mobile concentration and $7.63 \times 10^{-3}$ for the trapped concentration.

## 7. Future code improvements

There are some limitations to the current state of the tool. For example, FESTIM can only simulate one mobile (diffusing) species. It also does not account for multi-level trapping (i.e. more than one hydrogen atom can be trapped in a trapping site). These features – amongst others – are actively being worked on and will be included in the next major release of FESTIM. This next version will be based on FEniCSx [71], the next generation of the FEniCS library, and will include:

* better handling of materials interfaces using a Discontinuous Galerkin formulation at the interface

* support for mixed-topology meshes

* multiple diffusing species, allowing users to simulate co-diffusion of several isotopologues

* chemical reactions, making it possible to simulate multi-level trapping [44]

* native support for including an advection term

* support for coupling with fluid dynamics codes and neutronics codes (for tritium production and neutron heat generation)

This will also allow the development team to completely rewrite parts of the code and get rid of some technical debt accumulated over the years.

This next version of FESTIM is under active development. It is available on the fenicsx branch of the GitHub repository and can already be installed and tested by users.

## 8. Conclusion

FESTIM stands as a powerful and accessible code that offers robust solutions for simulating hydrogen transport in materials. Its versatility extends across various applications, including experimental analysis, component scale modelling, and safety analysis, making it relevant across a spectrum of research and industrial endeavours. Developed with adherence to modern practices such as GitHub, Continuous Integration, and Continuous Delivery, FESTIM sets a benchmark for efficiency and reliability in numerical tools.

In a comparative assessment against established counterparts like TMAP8 and COMSOL Multiphysics®, FESTIM demonstrated remarkable agreement, affirming its accuracy and efficacy in hydrogen transport simulations. This validation underscores the trustworthiness and suitability of the tool for diverse applications.

Furthermore, FESTIM's verification and validation processes are rigorously executed, ensuring robustness and accuracy in its predictive capabilities. While this paper provides a glimpse into this process, the comprehensive nature of FESTIM's validation instils confidence in its performance across various scenarios and conditions.

In conclusion, FESTIM emerges not only as a cutting-edge tool for hydrogen transport modelling but also as a testament to the efficacy of open-source development and collaborative practices. Its versatility, accuracy, and adherence to modern standards position it as a cornerstone in the pursuit of understanding and harnessing hydrogen dynamics in materials.

## CRediT authorship contribution statement

**Rémi Delaporte-Mathurin:** Conceptualization, Methodology, Software, Writing – original draft, Supervision. **James Dark:** Conceptualization, Software, Writing – original draft. **Gabriele Ferrero:** Investigation, Software, Writing – original draft, Writing – review & editing. **Etienne A. Hodille:** Software, Writing – review & editing. **Vladimir Kulagin:** Software, Writing – review & editing. **Samuele Meschini:** Software, Writing – review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

R. Delaporte-Mathurin would like to acknowledge funding from the Advanced Research Projects Agency-Energy (ARPA-E), USA, U.S. Department of Energy, USA. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

E. Hodille, J. Dark: This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 - EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

The work of G. Ferrero is financially supported by FSE REACT-EU-PON Ricerca e Innovazione 2014–2020, USA, DM 1061/2021.

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ijhydene.2024.03.184.

## References

[1] Behvar A, Haghshenas M, Djukic MB. Hydrogen embrittlement and hydrogen-induced crack initiation in additively manufactured metals: A critical review on mechanical and cyclic loading. Int J Hydrogen Energy 2024;58:1214–39. http://dx.doi.org/10.1016/j.ijhydene.2024.01.232, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0360319924002556.

[2] Choo WY, Lee JY. Thermal analysis of trapped hydrogen in pure iron. Metall Trans A 1982;13(1):135–40. http://dx.doi.org/10.1007/BF02642424.

[3] Wang J, Yang Y, Cai S, Yao J, Xie Q. Pore-scale modelling on hydrogen transport in porous media: Implications for hydrogen storage in saline aquifers. Int J Hydrogen Energy 2023;48(37):13922–33. http://dx.doi.org/10.1016/j.ijhydene.2022.11.299, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0360319922056415.

[4] Pasler V, Arbeiter F, Klein C, Klimenko D, Schlindwein G, von der Weth A. Development and verification of a component-level hydrogen transport model for a DEMO-like HCPB breeder unit with OpenFOAM. Fusion Eng Des 2018;127:249–58. http://dx.doi.org/10.1016/j.fusengdes.2018.01.008, [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0920379618300103.

[5] Baldwin MJ, Doerner RP. Hydrogen isotope transport across tungsten surfaces exposed to a fusion relevant He ion fluence. Nucl Fusion 2017;57(7):076031. http://dx.doi.org/10.1088/1741-4326/aa70b1, Publisher: IOP Publishing.
[6] Baroutaji A, Wilberforce T, Ramadan M, Olabi AG. Comprehensive investigation on hydrogen and fuel cell technology in the aviation and aerospace sectors. Re- new Sustain Energy Rev 2019;106:31–40. http://dx.doi.org/10.1016/j.rser.2019. 02.022, [Online]. Available: https://www.sciencedirect.com/science/article/pii/ S1364032119301157.

[7] Alnæs M, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, Richard- son C, Ring J, Rognes ME, Wells GN. The FEniCS project version 1.5. Arch Numer Softw 2015;3(100). http://dx.doi.org/10.11588/ans.2015.100.20553. [Online]. Available: https://journals.ub.uni-heidelberg.de/index.php/ans/article/ view/20553.

[8] Hodille EA, Bonnin X, Bisson R, Angot T, Becquart CS, Layet JM, Grisolia C. Macroscopic rate equation modeling of trapping/detrapping of hydrogen iso- topes in tungsten materials. J Nucl Mater 2015;467:424–31. http://dx.doi.org/ 10.1016/j.jnucmat.2015.06.041, [Online]. Available: http://www.sciencedirect. com/science/article/pii/S0022311515300660.

[9] Arredondo R, Schmid K, Subba F, Spagnuolo GA. Preliminary estimates of tritium permeation and retention in the first wall of DEMO due to ion bombardment. Nucl Mater Energy 2021;28:101039. http://dx.doi.org/10.1016/j. nme.2021.101039, [Online]. Available: https://www.sciencedirect.com/science/ article/pii/S2352179120101125.

[10] Longhurst GR. TMAP7 user manual. Tech. rep. INEEL/EXT-04-02352, Idaho National Laboratory (INL); 2008, http://dx.doi.org/10.2172/952013, [Online]. Available: https://www.osti.gov/biblio/952013-tmap7-user-manual.

[11] Matveev D, Zlobinski M, Temmerman GD, Unterberg B, Linsmeier C. CRDS modelling of deuterium release from co-deposited beryllium layers in temperature programmed and laser induced desorption experiments. Phys Scr 2020;2020(T171):014053. http://dx.doi.org/10.1088/1402-4896/ab5569, Publisher: IOP Publishing.

[12] Candido L, Alberghi C. Verification and validation of mHIT code over TMAP for hydrogen isotopes transport studies in fusion-relevant environments. Fusion Eng Des 2021;172:112740. http://dx.doi.org/10.1016/j.fusengdes.2021. 112740, [Online]. Available: https://www.sciencedirect.com/science/article/pii/ S0920379621005160.

[13] Smirnov RD, Gutler J, Krasheninnikov SI. Modeling of multispecies dynamics in fusion-related materials with FACE. Fusion Sci Technol 2017;71(1):75–83. http:// dx.doi.org/10.13182/FST16-125, Publisher: Taylor & Francis _eprint: https://doi.org/10.13182/FST16-125.

[14] Lindsay A. TMAP8. 2021, Language: en. [Online]. Available: https://www.osti. gov/doecode/biblio/50762.

[15] Dixon S. Aurora-multisphysics/achlys: Isotope self-diffusion. 2021, [Online]. Available: https://zenodo.org/record/6412090.

[16] Permann CJ, Gaston DR, Andrš D, Carlsen RW, Kong F, Lindsay AD, Miller JM, Peterson JW, Slaughter AE, Stogner RH, Martineau RC. MOOSE: Enabling massively parallel multiphysics simulation. SoftwareX 2020;11:100430. http:// dx.doi.org/10.1016/j.softx.2020.100430, [Online]. Available: https://www. sciencedirect.com/science/article/pii/S2352711019302973.

[17] Benannoue S, Charles Y, Mougenot J, Gaspérini M, Temmerman GD. Multidimensional finite-element simulations of the diffusion and trapping of hydrogen in plasma-facing components including thermal expansion. Phys Scr 2020;T171:014011. http://dx.doi.org/10.1088/1402-4896/ab4335.

[18] Ferrero G, Meschini S, Testoni R. A preliminary CFD and tritium transport analysis for ARC blanket. Fusion Sci Technol 2022;78(8):617–30. http://dx. doi.org/10.1080/15361055.2022.2096365, Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/15361055.2022.2096365.

[19] Delaporte-Mathurin R, Dark J. RemDelaporteMathurin/FESTIM-review: Revision 1. 2024, [Online]. Available: https://zenodo.org/records/10689958.

[20] McNabb A, Foster PK. A new analysis of the diffusion of hydrogen in iron and ferritic steels. Trans Metall Soc AIME 1963;227:618–27.

[21] Delaporte-Mathurin R, Hodille E, Mougenot J, Charles Y, Temmerman GD, Leblond F, Grisolia C. Influence of interface conditions on hydrogen transport studies. Nucl Fusion 2021;61(3):036038. http://dx.doi.org/10.1088/1741- 4326/abd95f, [Online]. Available: http://iopscience.iop.org/article/10.1088/ 1741-4326/abd95f.

[22] Longhurst GR. The soret effect and its implications for fusion reactors. J Nucl Mater 1985;131(1):61–9. http://dx.doi.org/10.1016/0022-3115(85) 90425-8, [Online]. Available: http://www.sciencedirect.com/science/article/pii/ 0022311585904258.

[23] Dasgupta D, Blondel S, Martínez E, Maroudas D, Wirth BD. Impact of Soret effect on hydrogen and helium retention in PFC tungsten under ELM-like conditions. Nucl Fusion 2023;63(7):076029. http://dx.doi.org/10.1088/1741-4326/acda4a, Publisher: IOP Publishing.

[24] Hodille EA, Bernard E, Markelj S, Mougenot J, Becquart CS, Bisson R, Grisolia C. Estimation of the tritium retention in ITER tungsten divertor target using macroscopic rate equations simulations. Phys Scr 2017;T170:014033. http://dx. doi.org/10.1088/1402-4896/aa8787.

[25] Ogorodnikova OV, Raepsaft X, Fütterer MA. Tritium permeation through the first wall of the EU-HCPB blanket. Fusion Eng Des 2000;49–50:921–6. http:// dx.doi.org/10.1016/S0920-3796(00)00339-2, [Online]. Available: https://www. sciencedirect.com/science/article/pii/S092037960003392.

[26] Delaporte-Mathurin R. Hydrogen transport in tokamaks: Estimation of the ITER divertor tritium inventory and influence of helium exposure (These de doctorat), Paris 13; 2022, [Online]. Available: https://www.theses.fr/2022PA131054.

[27] Laughlin JG, Dokken JS, Finsberg HNT, Francis EA, Lee CT, Rognes ME, Rangamani P. Smart: Spatial modeling algorithms for reaction and transport. 2023, arXiv:2306.07368 [q-bio]. [Online]. Available: http://arxiv.org/abs/2306. 07368.

[28] Dokken JS, Johansson A, Massing A, Funke SW. A multimesium finite element method for the Navier-Stokes equations based on projection methods. Comput Methods Appl Mech Engrg 2020;368:113129. http://dx.doi.org/10.1016/j. cma.2020.113129, [Online]. Available: https://www.sciencedirect.com/science/ article/pii/S0045782520303145.

[29] Kamensky D, Bazilevs Y. tIGAr: Automating isogeometric analysis with FEniCS. Comput Methods Appl Mech Engrg 2019;344:477–98. http://dx.doi.org/10. 1016/j.cma.2018.10.002, [Online]. Available: https://www.sciencedirect.com/ science/article/pii/S0045782518304985.

[30] Abali BE. Computational reality. Advanced structured materials, vol. 55, Singapore: Springer; 2017, http://dx.doi.org/10.1007/978-981-10-2444-3, [Online]. Available: http://link.springer.com/10.1007/978-981-10-2444-3.

[31] Logg A, Mardal K-A, Wells G, editors. Automated solution of differential equations by the finite element method: The FEniCS book. Lecture notes in computational science and engineering, vol. 84, Berlin, Heidelberg: Springer; 2012, http://dx.doi.org/10.1007/978-3-642-23099-8, [Online]. Available: https:// link.springer.com/10.1007/978-3-642-23099-8.

[32] Balay S, Abhyankar S, Adams MF, Benson S, Brown J, Brune P, Buschelman K, Constantinescu EM, Dalcin L, Dener A, Eijkhout V, Faibussowsitsch J, Gropp WD, Hapla V, Isaac T, Jolivet P, Karpeev D, Kaushik D, Knepley MG, Kong F, Kruger S, May DA, McInnes LC, Mills RT, Mitchell L, Munson T, Roman JE, Rupp K, Sanan P, Sarich J, Smith BF, Zampini S, Zhang H, Zhang J. PETSc/TAO users manual (rev. 3.19). Tech. rep. ANL-21/39-Rev.3.19, Argonne National Laboratory (ANL), Argonne, IL (United States); 2023, http://dx.doi.org/10.2172/1968587, [Online]. Available: https://www.osti.gov/biblio/1968587.

[33] Schlömer N. meshio: Tools for mesh files. 2022, [Online]. Available: https:// zenodo.org/records/6346837.

[34] Harris CR, Millman KJ, van der Walt SJ, Gommers R, Virtanen P, Cournaape D, Wieser E, Taylor J, Berg S, Smith NJ, Kern R, Picus M, Hoyer S, van Kerkwijk MH, Brett M, Haldane A, del Río JF, Wiebe M, Peterson P, Gérard-Marchant P, Sheppard K, Reddy T, Weckesser W, Abbasi H, Gohlke C, Oliphant TE. Array programming with NumPy. Nature 2020;585(7825):357–62. http://dx.doi.org/10.1038/s41586-020-2649-2, Bandiera-abest: a Cc_license_type: cc_by_Cg_type: Nature Research Journals Number: 7825 Primary_atype: Reviews Publisher: Nature Publishing Group Subject_term: Computational neuroscience;Computational science;Computer science;Software;Solar physics Subject_term_id: computational-neuroscience;computational-science;computer-science;software;solar-physics. [Online]. Available: https://www.nature.com/articles/s41586-020-2649-2.

[35] Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournaape D, Burovski E, Peterson P, Weckesser W, Bright J, van der Walt SJ, Brett M, Wilson J, Millman KJ, Mayorov N, Nelson ARJ, Jones E, Kern R, Larson E, Carey CJ, Polat I, Feng Y, Moore EW, VanderPlas J, Laxalde D, Perktold J, Cimrman R, Henriksen I, Quintero EA, Harris CR, Archibald AM, Ribeiro AH, Pedregosa F, van Mulbregt P. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods 2020;17(3):261–72. http://dx.doi.org/10.1038/s41592- 019-0686-2, [Online]. Available: https://www.nature.com/articles/s41592-019- 0686-2. Number: 3 Publisher: Nature Publishing Group.

[36] Hunter JD. Matplotlib: A 2D graphics environment. Comput Sci Eng 2007;9(3):90–5. http://dx.doi.org/10.1109/MCSE.2007.55, Conference Name: Computing in Science Engineering.

[37] Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka S, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. SymPy: symbolic computing in Python. PeerJ Comput Sci 2017;3:e103. http://dx.doi.org/10.7717/ peerj.cs.103, [Online]. Available: https://peerj.com/articles/cs-103. Publisher: PeerJ Inc.

[38] Ogorodnikova OV, Roth J, Mayer M. Deuterium retention in tungsten in dependence of the surface conditions. J Nucl Mater 2003;313–316:469–77. http://dx.doi.org/10.1016/S0022-3115(02)01375-2, [Online]. Available: http:// www.sciencedirect.com/science/article/pii/S0022311502013752.

[39] Baldwin MJ, Schwarz-Selinger T, Doerner RP. Experimental study and modelling of deuterium thermal release from Be-D co-deposited layers. Nucl Fusion 2014;54(7):073005. http://dx.doi.org/10.1088/0029-5515/54/7/073005, Publisher: IOP Publishing.

[40] Delaporte-Mathurin R, Hodille EA, Mougenot J, Charles Y, Grisolia C. Finite element analysis of hydrogen retention in ITER plasma facing components using FESTIM. Nucl Mater Energy 2019;21:100709. http://dx.doi.org/10.1016/j.nme. 2019.100709, [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/ S2352179119300547.
[41] Delaporte-Mathurin R, Hodille EA, Mougenot J, Charles Y, Grisolia C. Parametric optimisation based on TDS experiments for rapid and efficient identification of hydrogen transport materials properties. Nucl Mater Energy 2021;100984. http://dx.doi.org/10.1016/j.nme.2021.100984, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2352179121000661.

[42] Delaporte-Mathurin R, Hodille E, Mougenot J, De Temmerman G, Charles Y, Grisolia C. Parametric study of hydrogenic inventory in the ITER diver- tor based on machine learning. Sci Rep 2020;10(1):17798. http://dx.doi.org/10.1038/s41598-020-74844-w. [Online]. Available: https://www.nature.com/articles/s41598-020-74844-w. Number: 1 Publisher: Nature Publishing Group.

[43] Delaporte-Mathurin R, Yang H, Denis J, Dark J, Hodille EA, Temmerman GD, Bonnin X, Mougenot J, Charles Y, Bufferand H, Ciraolo G, Grisolia C. Fuel retention in WEST and ITER divertors based on FESTIM monoblock simula- tions. Nucl Fusion 2021;61(12):126001. http://dx.doi.org/10.1088/1741-4326/ac2bbd, Publisher: IOP Publishing.

[44] Hodille EA, Delaporte-Mathurin R, Denis J, Pečovnik M, Bernard E, Ferro Y, Sakamoto R, Charles Y, Mougenot J, De Backer A, Becquart CS, Markelj S, Grisolia C. Modelling of hydrogen isotopes trapping, diffusion and permeation in divertor monoblocks under ITER-like conditions. Nucl Fusion 2021. http://dx.doi.org/10.1088/1741-4326/ac2abc, [Online]. Available: http://iopscience.iop.org/article/10.1088/1741-4326/ac2abc.

[45] Delaporte-Mathurin R, Chochoy R, Mougenot J, Charles Y, Hodille EA, Griso- lia C. 3D effects on hydrogen transport in ITER-like monoblocks. Nucl Fusion 2023. http://dx.doi.org/10.1088/1741-4326/ad1019, [Online]. Available: http://iopscience.iop.org/article/10.1088/1741-4326/ad1019.

[46] Dark J, Delaporte-Mathurin R, Charles Y, Hodille EA, Grisolia C, Mougenot J. Influence of hydrogen trapping on WCLL breeding blanket performances. Nucl Fusion 2021;61(11):116076. http://dx.doi.org/10.1088/1741-4326/ac28b0, Pub- lisher: IOP Publishing.

[47] Romano PK, Horelik NE, Herman BR, Nelson AG, Forget B, Smith K. OpenMC: A state-of-the-art Monte Carlo code for research and development. Ann Nucl Energy 2015;82:90–7. http://dx.doi.org/10.1016/j.anucene.2014.07.048, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S030645491400379X.

[48] Ferry SE, Woller KB, Peterson EE, Sorensen C, Whyte DG. The LIBRA experiment: Investigating robust tritium accountability in molten FLiBe exposed to a D-T fusion neutron spectrum. Fusion Sci Technol 2022;1–23. http://dx.doi.org/10.1080/15361055.2022.2078136, Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/15361055.2022.2078136.

[49] Brooks H, Davis A. Scalable multi-physics for fusion reactors with AURORA. Plasma Phys Control Fusion 2023;65(2):024002. http://dx.doi.org/10.1088/1361-6587/aca998, [Online]. Available: https://iopscience.iop.org/article/10.1088/1361-6587/aca998.

[50] Moro F, Del Nevo A, Flammini D, Martelli E, Mozzillo R, Noce S, Vil- lari R. Neutronic analyses in support of the WCLL DEMO design development. Fusion Eng Des 2018;136:1260–4. http://dx.doi.org/10.1016/j.fusengdes.2018.04.113, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0920379618304071.

[51] Longhurst GR, Ambrosek J. Verification and validation of the tritium trans- port code TMAP7. Fusion Sci Technol 2005;48(1):468–71. http://dx.doi.org/10.13182/FST05-A967.

[52] Simon P-CA, Humrickhouse PW, Lindsay AD. Tritium transport modeling at the pore scale in ceramic breeder materials using TMAP8. IEEE Trans Plasma Sci 2022;50(11):4465–71. http://dx.doi.org/10.1109/TPS.2022.3183525, [Online]. Available: https://ieeexplore.ieee.org/document/9831077/.

[53] COMSOL® software version 6.1 release highlights. [Online]. Available: https://www.comsol.com/release/6.1.

[54] Schenk O, Gärtner K. PARDISO. In: Padua D, editor. Encyclopedia of parallel computing. Boston, MA: Springer US; 2011, p. 1458–64. http://dx.doi.org/10.1007/978-0-387-09766-4_90.

[55] Simoni L, Falcade T, Ferreira DCF, Kwietniewski CEF. An integrated experimental and modeling approach to determine hydrogen diffusion and trapping in a high- strength steel. Int J Hydrogen Energy 2021;46(50):25738–51. http://dx.doi.org/10.1016/j.ijhydene.2021.05.079, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0360319921018334.

[56] Montupet-Leblond F, Hodille EA, Vartanian S, Payet M, Delaporte-Mathurin R, Mougenot J, Charles Y, Bernard E, Grisolia C. Permeation and trapping of hy- drogen in Eurofer97. Nucl Mater Energy 2021;29:101062. http://dx.doi.org/10.1016/j.nme.2021.101062, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2352179121001290.

[57] Simmonds MJ, Baldwin MJ, Temmerman GD, Doerner RP. TMAP modeling of D release from baked multi-layer Be–D co-deposits. Phys Scr 2020;T171:014043. http://dx.doi.org/10.1088/1402-4896/ab4de4, Publisher: IOP Publishing.

[58] Markelj S, Schwarz-Selinger T, Pečovnik M, Chrominski W, Šestan A, Zavašnik J. Deuterium transport and retention in the bulk of tungsten containing he- lium: the effect of helium concentration and microstructure. Nucl Fusion 2020;60(10):106029. http://dx.doi.org/10.1088/1741-4326/abadae, Publisher: IOP Publishing.

[59] Shimada M, Taylor CN. Improved tritium retention modeling with reaction- diffusion code TMAP and bulk depth profiling capability. Nucl Mater Energy 2019;19:273–8. http://dx.doi.org/10.1016/j.nme.2019.03.008, [Online]. Avail- able: https://www.sciencedirect.com/science/article/pii/S2352179118302758.

[60] Pečovnik M, Markelj S, Založnik A, Schwarz-Selinger T. Influence of grain size on deuterium transport and retention in self-damaged tungsten. J Nucl Mater 2019;513:198–208. http://dx.doi.org/10.1016/j.jnucmat.2018.10.026, [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0022311518307256.

[61] Matveev D, Wensing M, Ferry L, Virot F, Barrachin M, Ferro Y, Linsmeier C. Reaction-diffusion modeling of hydrogen transport and surface effects in applica- tion to single-crystalline Be. Nucl Instrum Methods Phys Res B 2018;430:23–30. http://dx.doi.org/10.1016/j.nimb.2018.05.037, [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0168583X18303537.

[62] Hodille EA, Založnik A, Markelj S, Schwarz-Selinger T, Becquart CS, Bisson R, Grisolia C. Simulations of atomic deuterium exposure in self-damaged tungsten. Nucl Fusion 2017;57(5):056002. http://dx.doi.org/10.1088/1741-4326/aa5aa5.

[63] Pečovnik M, Schwarz-Selinger T, Markelj S. Experiments and modelling of multiple sequential MeV ion irradiations and deuterium exposures in tung- sten. J Nucl Mater 2021;152947. http://dx.doi.org/10.1016/j.jnucmat.2021.152947, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0022311521001707.

[64] Simmonds MJ, Yu JH, Wang YQ, Baldwin MJ, Doerner RP, Tynan GR. E x- panding the capability of reaction-diffusion codes using pseudo traps and temperature partitioning: Applied to hydrogen uptake and release from tung- sten. J Nucl Mater 2018;508:472–80. http://dx.doi.org/10.1016/j.jnucmat.2018.05.080, [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0022311517315301.

[65] Hodille EA, Ghiortghiu F, Addab Y, Založnik A, Minissale M, Piazza Z, Martin C, Angot T, Gallais L, Barthe M-F, Becquart CS, Markelj S, Mougenot J, Grisolia C, Bisson R. Retention and release of hydrogen isotopes in tungsten plasma-facing components: the role of grain boundaries and the native oxide layer from a joint experiment-simulation integrated approach. Nucl Fusion 2017;57(7):076019. http://dx.doi.org/10.1088/1741-4326/aa6d24, Publisher: IOP Publishing.

[66] Montupet-Leblond F, Hodille EA, Payet M, Delaporte-Mathurin R, Bernard E, Charles Y, Mougenot J, Vartanian S, Grisolia C. Influence of traps re- versibility on hydrogen permeation and retention in Eurofer97. Nucl Fusion 2022;62(8):086011. http://dx.doi.org/10.1088/1741-4326/ac6e74, Publisher: IOP Publishing.

[67] Calderoni P, Sharpe P, Hara M, Oya Y. Measurement of tritium permeation in flibe (2LIF-BeF2). Fusion Eng Des 2008;83(7–9):1331–4. http://dx.doi.org/10.1016/j.fusengdes.2008.05.016, [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0920379608000926.

[68] Fukada S, Morisaki A. Hydrogen permeability through a mixed molten salt of LiF, NaF and KF (Flinak) as a heat-transfer fluid. J Nucl Mater 2006;358(2):235– 42. http://dx.doi.org/10.1016/j.jnucmat.2006.07.011, [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0022311506004089.

[69] Esteban GA, Peña A, Legarda F, Lindau R. Hydrogen transport and trapping in ODS-EUROFER. Fusion Eng Des 2007;82(15):2634–40. http://dx.doi.org/10.1016/j.fusengdes.2007.02.002, [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0920379607000580.

[70] Roy CJ, Oberkampf WL, editors. Exact solutions. In: Verification and validation in scientific computing. Cambridge: Cambridge University Press; 2010, p. 208– 48. http://dx.doi.org/10.1017/CBO9780511760396.009, [Online]. Available: https://www.cambridge.org/core/books/verification-and-validation-in-scientific-computing/exact-solutions/DFB030CD8A8334FA13DF3B2A627964E4.

[71] Baratta I, Dean J, Dokken JrS, Habera M, Hale J, Richardson C, Rognes M, Scroggs M, Sime N, Wells G. DOLFINx: The next generation FEniCS problem solving environment. 2023.