FESTIM: a framework for hydrogen isotope transport on top of FEniCS

Rémi Delaporte-Mathurin, on behalf of the whole FESTIM community

Simula 25 Years · FEniCS Workshop | Oslo | September 2026

The problem (no nuclear physics)

Hydrogen isotopes in a solid: mobile, drifting, and repeatedly trapped in defects.

Many species at once, and stiff
Concentration jumps at interfaces
1D slabs to CAD-meshed components

\partial_t c_\mathrm{m} = \nabla \cdot \left( D \nabla c_\mathrm{m} - c_\mathrm{m}\, \mathbf{u} \right) + S - \sum_i \partial_t c_{\mathrm{t},i}

↑ drift: thermal, electric, advective

\partial_t c_{\mathrm{t},i} = k_i\, c_\mathrm{m} \left( n_i - c_{\mathrm{t},i} \right) - p_i\, c_{\mathrm{t},i}

one equation per trap population

Why we care: tritium is the fuel, it is scarce, radioactive, and it leaks through everything.

Who we build for

Our users

Know their physics
Can write Python
Have never assembled a bilinear form
Will not debug a solver

So we abstract:

FESTIM owns the vocabulary: materials, traps, species, reactions
FEniCSx owns the mathematics
Users never see a form

We try to encapsulate every FEniCSx feature we adopt.

Two eras

2019 – 2024 · FESTIM 1

Legacy FEniCS
One mobile species + traps
Chemical-potential change of variable
Validated against experiments

2024 – today · FESTIM 2

DOLFINx
Arbitrary species
Submeshes & mixed function spaces
~11× faster, MPI-native

FEniCSx made so many things possible...

The thing that shaped the code

At a material interface the concentration is discontinuous. Two conditions hold across it:

\frac{c_1}{K_1} = \frac{c_2}{K_2}

D_1 \nabla c_1 \cdot n = D_2 \nabla c_2 \cdot n

chemical potential continuous

flux continuous

Legacy workaround

Solve for θ = c/K: one continuous field, one space.

The cost: every source, BC, export and reaction rate rewritten in θ. Nonlinear in T. Didn't scale well.

What we wanted

One field per material, each on its own submesh.

Coupled by a condition we choose, not one baked into the function space.

What DOLFINx gave us

Submeshes, MixedFunctionSpace, assembly over mixed domains. One field per subdomain, one block system:

\sum_i \int_{\Omega_i} \left( \partial_t c_i\, v_i + D_i \nabla c_i \cdot \nabla v_i \right) \mathrm{d}x \;=\; \sum_i \int_{\Omega_i} S_i\, v_i \, \mathrm{d}x

… and the only thing tying the subdomains together is one term on the interface:

\int_\Gamma \gamma \left[\!\left[ \frac{c}{K} \right]\!\right] \left[\!\left[ v \right]\!\right] \mathrm{d}s

γ is a penalty parameter; the Nitsche variant adds consistency terms and drops the tuning
No change of variable required

Our users write none of this.

Weak coupling is a feature

Once the interface condition lives in the form rather than in the function space:

Continuity of chemical potential is just one choice; a flux condition, a kinetic law or a surface reaction work too

Which matters, because equilibrium is sometimes simply wrong.

Delaporte-Mathurin & Dark, When chemical potential continuity fails, arXiv:2608.26231

A subdomain doesn't need to be the same dimension as its parent mesh.

NEW SINCE JULY

Subdomains that are manifolds

A grain boundary is a fast path through a slow bulk. Meshing it as a thin volume is hopeless.

gb = F.VolumeSubdomain(
    id=2,
    material=F.Material(D_0=D_gb, E_D=0.0),
    dim=1,                    # codimension 1 in a 2D mesh
    locator=lambda x: np.isclose(x[0], 0.5),
)

One extra argument: dim. Everything else is the subdomain API users already know.
Under the hood, the same story as before: its own submesh, its own field, same mixed system.
Diffusion along the manifold ~100× faster than through the bulk.

Coupling across dimensions

Bulk on Ω (dim n), manifold on Γ (dim n−1), each carrying its own field:

\partial_t c_\mathrm{b} = \nabla \cdot (D_\mathrm{b} \nabla c_\mathrm{b}) \quad \text{in } \Omega

\partial_t c_\mathrm{gb} = \nabla_\Gamma \cdot (D_\mathrm{gb} \nabla_\Gamma c_\mathrm{gb}) + k\,(c_\mathrm{b} - c_\mathrm{gb}) \quad \text{on } \Gamma

↑ tangential gradient along the manifold

-D_\mathrm{b} \nabla c_\mathrm{b} \cdot n = k\,(c_\mathrm{b} - c_\mathrm{gb})

The flux leaving the bulk is the source entering the manifold: conservative by construction
It nests: a coolant channel is a 1D manifold on a 2D wall, and its inlet and outlet are dim = 0

Hydrogen takes the fast lane

A real polycrystalline microstructure. Hydrogen enters at the top.

Two fields solved together:
c_b in the grains (dim=2)
c_gb on the boundaries (dim=1)

It reaches depth through the boundary network long before the bulk front arrives.

NEW SINCE JULY

An enclosure is one number

Hydrogen permeates out of a wall into a sealed volume. The pressure obeys an ODE fed by the surface flux, and feeds straight back into the wall’s boundary condition.

\frac{\mathrm{d}P}{\mathrm{d}t} = \frac{R\,T}{V} \int_{\partial \Omega} \varphi \, \mathrm{d}s \qquad P \in \mathbb{R}

One scalar unknown, coupled both ways to a PDE. Our previous answer was an operator splitting and a hand-rolled ODE loop.

We didn’t need one. DOLFINx has real-space elements: P becomes a single global degree of freedom in the same mixed space as the fields, and Newton solves the lot monolithically.

For the user: F.GasSpecies and F.Enclosure. Sieverts’ law then sets c from a pressure with no nodal values, so we impose it weakly, with Nitsche.

An enclosure to simulate a pump

Two membranes with a pumped gap. The 1D mesh is discontiguous: no cell spans the gap, so the gas enclosure is the only path across.

c = c₀

c = 0

layer 1

gap ↓ pump

layer 2

Both faces exchange with one Enclosure via a surface reaction
A Pump is an opening on it; a purge stream is the same object

Pump on: enclosure pressure 25 Pa → 0, and the concentration reaching the second membrane falls by ~5 orders of magnitude.

Permeation flux

Drifts!

Hydrogen does not only follow its own gradient. It drifts: down temperature gradients, under electric fields, and with a moving coolant.

J = -D \nabla c \; - \; \frac{D\, Q^{*} c}{k_B T^{2}} \nabla T

Soret / thermodiffusion; Q_star is the heat of transport

J = -D \nabla c + c\, \mathbf{u}, \qquad \mathbf{u} = \sum_j \mathbf{u}_j

Several drift terms may act on one species; the velocities simply add. SoretTerm, ElectromigrationTerm, AdvectionTerm.

Because the residual is built in UFL, a drift term is literally one more term in the form. No new solver, no new time integrator, and it composes with submeshes and manifolds for free.

NEW SINCE AUGUST

NEW SINCE SEPT

Submeshes broke our I/O

Every feature above multiplies the number of meshes and fields. Users got a directory of .bp files and no way to restart.

F.SpeciesExport("state.bp", field=H, format="checkpoint")

state = F.read_function_from_file("state.bp", name="H", timestamp=10.0)

io4dolfinx (evolution of adios4dolfinx): N-to-M checkpointing, multiple backends

Three codes, one field

OpenMC: tritium source

OpenFOAM: velocity

FESTIM: tritium field

openmc2dolfinx and foam2dolfinx hand over DOLFINx functions, not files, so they compose with submeshes, manifolds and drift terms for free.

ARC breeding blanket · ~243 mg tritium inventory · Dark, Sircar, Bae, Dokken & Delaporte-Mathurin, arXiv:2608.05398

How do we know any of it is right?

Verification

Manufactured solutions: 10 cases, convergence orders checked
Manufactured exact solutions: 6 cases
The TMAP benchmark suite

Validation

8 experimental campaigns
Plasma-driven permeation, deuterium retention in damaged tungsten, Be co-deposition, protium in titanium, EUROFER implantation

f := \partial_t c_\star - \nabla \cdot (D \nabla c_\star) \;\Longrightarrow\; \lVert c_h - c_\star \rVert = \mathcal{O}(h^{p+1})

Every case is a runnable script. The report rebuilds from them, in public.

festim-vv-report.readthedocs.io

Independently cross-compared against MHIMS and mHIT (Nuclear Materials and Energy)

What it looks like in anger

Tritium contamination in a heat exchanger: CFD-driven, multi-material, transient
Divertor monoblock: CAD mesh, coupled heat transfer, three materials, two interfaces

What we want next

On our side

DG for advection-dominated transport
Hydride formation: moving phase boundaries
Better diagnostics when Newton stalls

On our wishlist from you

End-to-end mixed-topology meshes: it seems there currently are limits in UFL
A higher-level mixed-dimensional API
GPU assembly for the big transients

The project

30+

contributors

27+

institutions

8k

forum views / month

Apache 2.0

NumFOCUS affiliated

conda install -c conda-forge festim

Monthly dev meetings, open to anyone. A workshop with runnable tutorials for everything in this talk.

If you maintain a FEniCS-based code and want to compare notes on hiding FEM from users, come find me. 😉

Thank you

Happy 25th, Simula, and thank you for DOLFINx.

github.com/festim-dev/FESTIM

festim.readthedocs.io · festim-workshop.readthedocs.io

Questions?