FESTIM v2.0: Upgraded framework for multi-species hydrogen transport

Remi Delaporte-Mathurin, James Dark, Jorgen Dokken, Huihua Yang, Kaelyn Dunnell, Chirag Khurana, Ethan Peterson

FEniCS 2026 | June, 2026 | Paris, France

Nuclear Fusion

Deuterium

Tritium

Helium

neutron

Fusion

The tritium issue

The maths behind FESTIM

\frac{\partial c_\mathrm{m}}{\partial t} = \nabla \cdot (D \nabla c_\mathrm{m}) + S - \sum \frac{\partial c_\mathrm{t,i}}{\partial t}

\frac{\partial c_\mathrm{t, i}}{\partial t} = k_i \ c_\mathrm{m} \ (n_\mathrm{trap,i} - c_\mathrm{t}) - p_i c_\mathrm{t}

Diffusion-reaction

- D_1 \nabla c_\mathrm{m,1} = - D_2 \nabla c_\mathrm{m,2}

\frac{c_1}{K_{1}} =\frac{c_2}{K_{2}}

\frac{c_1}{K_{1}} =\left( \frac{c_2}{K_{2}} \right) ^ {2}

At interfaces

Hard

Very hard!

Flux continuity

Conservation of chemical potential

FESTIM Evolution Timeline

2019 - 2024

FESTIM ≤ 1

Based on Legacy FEniCS
Initial development and validation
Established core H-transport models

2024 - Today

FESTIM 2

Migration to DOLFINx
New multi-species paradigm
Improved performance & scalability
Modern FEniCS ecosystem integration

github.com/festim-dev/festim

FESTIM2 in a nutshell github.com/festim-dev/festim

User friendly - Intuitive, Python API easy to use
Apache 2.0 license - Suitable for academic and commercial
Easy install - conda install -c conda-forge festim
Online tutorials - Dozens of tutorials
Online documentation - Auto-generated

8,000+ forum views/month - Public, searchable Q&A
27+ institutions - Adopted by universities, labs and industry
30+ contributors - Open to new developers and users
Monthly dev meetings - Community-driven roadmap
130+ Slack users - Active developer and user support

Built on FEniCSx - Proven FEM backend, 1.5M+ downloads
5 team members at MIT - more around the world
On GitHub - Transparent, version-controlled development
NumFOCUS affiliation - Backed by a sustainability non-profit
CI/CD workflows - Rigorous testing and streamlined releases

Performance

FEniCSx - JIT to fast C++ rapid assembly
High-performance solvers - PETSc back-end
Built for HPC - Native MPI, runs on 1000’s of cores
Tuned control - Custom meshes, solvers, finite elements
Fast I/O - Parallel mesh loading, low memory footprint

Usability

Sustainability

Community

Examples

The FESTIM team at MIT and beyond

James Dark

Postdoctoral associate

Huihua Yang

Postdoctoral associate

Rémi Delaporte-Mathurin

Research Scientist

Kaelyn Dunnell

Graduate Student

Chirag Khurana

Graduate Student

MIT Team

Joe Dean

Research Associate
University of Cambridge

Jørgen Dokken

Senior Research Scientist
Simula Research Laboratory

FEniCS collaborators

FESTIM Community

31+ contributors

FESTIM is used worldwide

Docs users

Fusion machine developers

Research organisations

Average person in this room

knows maths
knows how to code

can't read greek
Not FEM experts

\nu \nabla \tau \eta \xi \delta

import festim as F
import numpy as np

my_model = F.HydrogenTransportProblem()

my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=1001))

mat = F.Material(D_0=1, E_D=0.0)

volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat)
boundary_left = F.SurfaceSubdomain1D(id=1, x=0)
boundary_right = F.SurfaceSubdomain1D(id=2, x=1)
my_model.subdomains = [volume_subdomain, boundary_left, boundary_right]

H = F.Species("H")
my_model.species = [H]

my_model.temperature = 300

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=boundary_left, value=1, species=H),
    F.FixedConcentrationBC(subdomain=boundary_right, value=0, species=H),
]

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=2)

my_model.settings.stepsize = F.Stepsize(0.05, milestones=[0.05, 0.1, 0.2, 0.5, 1])  # s

profile = F.Profile1DExport(field=H, subdomain=volume_subdomain, times=my_model.settings.stepsize.milestones)

my_model.exports = [
    F.VTXSpeciesExport(
        field=H,
        filename="hydrogen_concentration.bp",
    ),
    profile,
]

my_model.initialise()

my_model.run()

Before FEniCSx....

\frac{c_1}{K_{1}} =\frac{c_2}{K_{2}}

\theta =\frac{c}{K}

\frac{\partial c_\mathrm{m}}{\partial t} = \nabla \cdot (D \nabla c_\mathrm{m}) + S - \sum \frac{\partial c_\mathrm{t,i}}{\partial t}

\frac{\partial c_\mathrm{t, i}}{\partial t} = k_i \ c_\mathrm{m} \ (n_\mathrm{trap,i} - c_\mathrm{t}) - p_i c_\mathrm{t}

\theta_1 =\theta_2

\frac{\partial (\theta K)}{\partial t} = \nabla \cdot (D \nabla (\theta K) )+ S - \sum \frac{\partial c_\mathrm{t,i}}{\partial t}

\frac{\partial c_\mathrm{t, i}}{\partial t} = k_i \ (\theta K) \ (n_\mathrm{trap,i} - c_\mathrm{t}) - p_i c_\mathrm{t}

Change of variable method

With FEniCSx

Mixed formulation

With FEniCSx

Mixed formulation

Submesh 1

Submesh 2

With FEniCSx

Mixed formulation

\partial_t c_1 = \nabla \cdot (D \nabla c_1) + S

\partial_t c_2 = \nabla \cdot (D \nabla c_2) + S

\frac{c_1}{K_{1}} =\frac{c_2}{K_{2}}

1. Penalty method

- \int_{\Gamma} \langle D \nabla c^n\rangle \cdot \mathbf{n} \llbracket v \rrbracket - \langle D \nabla v \rangle\cdot\mathbf{n} \llbracket \frac{c^n}{K}\rrbracket + \frac{\gamma}{\langle h\rangle}\llbracket \frac{c^n}{K} \rrbracket\llbracket v \rrbracket ~\mathrm{d}s

2. Nitsche method

- \int_{\Gamma} \gamma \left( \frac{c_1}{K_1} - \frac{c_2}{K_2} \right) \llbracket v \rrbracket ~\mathrm{d}s

Dark, James, Rémi Delaporte-Mathurin, Jørgen S. Dokken, et al. “FESTIM v2.0: Upgraded Framework for Multi-Species Hydrogen Transport and Enhanced Performance.” International Journal of Hydrogen Energy 220 (March 2026): 153987. https://doi.org/10.1016/j.ijhydene.2026.153987.

Please redirect complex maths questions to this guy!

With FEniCSx

Mixed formulation

c_1

c_2

With FEniCSx

Mixed formulation

No need for change of variable
but requires a penalty term
Provides flexibility on interface physics
- Flux interface conditions
- Multi-species reactions