Hydrogen transport simulations made easy

 

Webinar | May 15th, 2026

Presented by:
James Dark & Remi Delaporte-Mathurin

Brief history of FESTIM

2019
Start of Development
Predicting tritium retention in 3D ITER monoblocks.

Brief history of FESTIM

2022
Open Source Release
Apache 2.0 License. Increases transparency, collaboration and flexibility.

Brief history of FESTIM

2023

FESTIM v1.0 release
Official release and transition to FESTIM v2.0 development.

Brief history of FESTIM

2023

Start of FESTIM 2 branch

Brief history of FESTIM

2024

FESTIM becomes a NumFOCUS affiliate

Brief history of FESTIM

2026
FESTIM v2.0 release
A complete rethink to eliminate technical debt and improve scaling.

FESTIM2 in a nutshell

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

3,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

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
GitHub Community
31+ contributors

FESTIM Fellowship programme commenced in 2026

"Even though it's quite nice to be the FESTIM v2 expert on site (!!!), I would absolutely recommend the programme to my colleagues."

First FESTIM fellow - Tez Orr (UKAEA)

 

Hands-on-learning - Tez worked directly with the dev team on UKAEA cases and fixed a bug in v2.0

 

Immediate impact  - Became UKAEA's resident FESTIM expert, contributing to both the LIBRTI and STEP programmes

 

FESTIM today

See more FESTIM stats:
festim-dev.github.io/pose

See all FESTIM papers:

festim.readthedocs.io/en/latest/publications.html

FESTIM is used worldwide

Docs users

Trusted by researchers worldwide

Beihang University Logo
Commonwealth Fusion Systems Logo
Cranfield University Logo
ENEA Logo
ENI Logo
Exeter University Logo
Frazer Nash Logo
GenF Systems Logo
IRFM Cadarache Logo
ITER Organization Logo
Institut Jozef Stefan Logo
Institute of plasma physics CAS Logo
Kyoto Fusioneering Logo
LSPM Villetaneuse Logo
MEPhI Moscow Logo
MIT Logo
Marathon Fusion Logo
Oak Ridge National Laboratory Logo
Politecnico di Torino Logo
Thea Energy Logo
Zap Energy Logo
digiLab Logo
Rochester Logo
University San Diego Logo

Are you using FESTIM?

Fill in this form!
github.com/festim-dev/FESTIM/issues

What changed in FESTIM2?

Upgrade from legacy-FEniCS to FEniCSx

  Before

  • Change of variable required - did not scale well
  • Single solver method available
  • No submesh coupling
  • Limited mesh element types

 

 

 Now

  • Mixed-domain: no substitution needed
  • Multiple solver methods (Penalty, Nitches)
  • Submesh coupling: 1D to 2D problems
  • Mixed topology meshes
  • Improved parallel I/O 

 

~11x faster

  Before

  • Monolithic approach
  • Change of variable required
    • Poor scaling
  • Little control over interface equations
  • Adding advection was tricky

  Now

  • Mixed-domain: no substitution needed
  • No change of variable needed
  • Improved performance
  • Full control over interface conditions
  • Very suitable for multi-physics 

 

Multi-material handling

Multi-species and Reactions

  Before

  • Single mobile species only
  • No multi-isotope support
  • No multi-occupancy trapping
  • Rigid trapping reactions

  Now

  • Arbitrary species
    • Multi-isotope
    • Multi-occupancy trapping
    • Defect diffusion
  • Flexible reactions
  • Bulk and surface reactions

 

Examples

FESTIM in action!

Monoblock

CAD based simulations

Read a mesh from GMSH

mesh_data = gmshio.read_from_msh(
    "gmsh/mesh3D.msh", MPI.COMM_WORLD, 0, gdim=3
)
mesh = mesh_data.mesh

facet_tags = mesh_data.facet_tags

cell_tags = mesh_data.cell_tags

Can also import meshes from other sources thanks to io4dolfinx

Heat transfer simulation

Run heat conduction simulation

heat_transfer_problem = F.HeatTransferProblem()

...

heat_flux_top = F.HeatFluxBC(subdomain=top_surface, value=10e6)

convective_flux_coolant = F.HeatFluxBC(
    subdomain=cooling_surface, value=lambda T: h_convective * (T_coolant - T)
)

...

heat_transfer_problem.run()
10^6 \ \mathrm{W/m^2}

Convective cooling

H transport problem

Pass the temperature field to the H transport problem

my_model = F.HydrogenTransportProblemDiscontinuous()

my_model.temperature = heat_transfer_problem.u

my_model.method_interface = "penalty"

my_model.subdomains = all_subdomains

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

my_model.mesh = mesh

my_model.surface_to_volume = {
    top_surface: tungsten_volume,
    cooling_surface: cucrzr_volume,
    poloidal_gap_w: tungsten_volume,
    poloidal_gap_cu: copper_volume,
    poloidal_gap_cucrzr: cucrzr_volume,
    toroidal_gap: tungsten_volume,
    bottom: tungsten_volume,
}

penalty_term = 1e20
my_model.interfaces = [
    F.Interface(
        id=16, subdomains=(tungsten_volume, copper_volume), penalty_term=penalty_term
    ),
    F.Interface(
        id=17, subdomains=(copper_volume, cucrzr_volume), penalty_term=penalty_term
    ),
]

import ufl

# Plasma implantation flux BC
phi = 1.61e22
R_p = 9.52e-10
implantation_flux_top = F.FixedConcentrationBC(
    subdomain=top_surface,
    value=lambda T: phi * R_p / (tungsten.D_0 * ufl.exp(-tungsten.E_D / F.k_B / T)),
    species=H,
)

# Instantaneous molecular recombination flux BCs at all other surfaces (fixed concentration of 0)
recombination_fluxes = [
    F.FixedConcentrationBC(subdomain=surf, value=0, species=H)
    for surf in [
        toroidal_gap,
        bottom,
        poloidal_gap_w,
        poloidal_gap_cu,
        poloidal_gap_cucrzr,
        cooling_surface,
    ]
]

my_model.boundary_conditions = [implantation_flux_top] + recombination_fluxes

exports = {
    "poloidal_gap_cu_flux": F.SurfaceFlux(surface=poloidal_gap_cu, field=H),
    "poloidal_gap_cucrzr_flux": F.SurfaceFlux(surface=poloidal_gap_cucrzr, field=H),
    "poloidal_gap_w_flux": F.SurfaceFlux(surface=poloidal_gap_w, field=H),
    "toroidal_gap_flux": F.SurfaceFlux(surface=toroidal_gap, field=H),
    "bottom_flux": F.SurfaceFlux(surface=bottom, field=H),
    "inventory_w": F.TotalVolume(field=H, volume=tungsten_volume),
    "inventory_cu": F.TotalVolume(field=H, volume=copper_volume),
    "inventory_cucrzr": F.TotalVolume(field=H, volume=cucrzr_volume),
}


my_model.exports = list(exports.values())

my_model.settings = F.Settings(
    atol=1e-8,
    rtol=1e-10,
    transient=False,
    max_iterations=10,
)

my_model.initialise()
my_model.run()

Implantation from plasma

c = 0

Examples: post-processing

Complete tutorials available at

festim-workshop.readthedocs.io/en/latest/content/post_process/intro.html

Field Exports

Writes a time-dependent field to a file

  • Solutions: concentrations, temperature
  • Checkpoints: restart a simulation from a previous result
     
  • Reaction rate
  • Custom field export

New in FESTIM2!

Derived Quantities

Compute quantities obtained from the solution

Thermo-Desorption Spectroscopy tutorial
festim-workshop.readthedocs.io/en/latest/content/applications/task02.html

left_flux = F.SurfaceFlux(surface=left_boundary, field=H)
right_flux = F.SurfaceFlux(surface=right_boundary, field=H)
total_trapped_H1 = F.TotalVolume(field=trapped_H1, volume=volume_subdomain)

contribution_trap_1 = -np.diff(total_trapped_H1.data) / np.diff(total_trapped_H1.t)
  • Built-in quantities:
    • Surface flux
    • Average/Total volume
    • Average/Total surface
    • Min/max
  • CustomQuantity

Examples: multi-species framework

Multi-species

H = F.Species("H")

my_model.species = [H]

Defining species

\mathrm{H}

Multi-species

H = F.Species("H")
D = F.Species("D")


my_model.species = [H, D]
\mathrm{H}, \mathrm{D}

Defining species

Multi-species

H = F.Species("H")
empty_traps = F.Species("empty", mobile=False)
trapped_H = F.Species("trapped_H", mobile=False)


my_model.species = [H, empty_traps, trapped_H]

Adding a reaction

\mathrm{H} + [\ ] \xrightleftharpoons[p]{k} [\mathrm{H}]
trapping = F.Reaction(
	reactant=[H, empty_traps],
    product=[trapped_H],
    k_0=...,
    E_k=...,
    p_0=...,
    E_p=...,
)


my_model.reactions = [trapping]

Multi-species

H = F.Species("H")
empty_traps = F.Species("empty", mobile=True)
trapped_H = F.Species("trapped_H", mobile=False)


my_model.species = [H, empty_traps, trapped_H]

Mobile traps?

\mathrm{H} + [\ ] \xrightleftharpoons[p]{k} [\mathrm{H}]
trapping = F.Reaction(
	reactant=[H, empty_traps],
    product=[trapped_H],
    k_0=...,
    E_k=...,
    p_0=...,
    E_p=...,
)


my_model.reactions = [trapping]

Multi-species

H = F.Species("H")
empty_traps = F.Species("empty", mobile=False)
empty_traps2 = F.Species("empty2", mobile=False)
trapped_H = F.Species("trapped_H", mobile=False)
trapped_H2 = F.Species("trapped_H2", mobile=False)

my_model.species = [
  H, empty_traps, trapped_H,
  empty_traps2, trapped_H2
]

2 traps

\mathrm{H} + [\ ] \xrightleftharpoons[p]{k} [\mathrm{H}]
trapping = F.Reaction(
	reactant=[H, empty_traps],
    product=[trapped_H],
	...
)

trapping2 = F.Reaction(
	reactant=[H, empty_traps2],
    product=[trapped_H2],
	...
)


my_model.reactions = [trapping, trapping2]
\mathrm{H} + [\ ]_2 \xrightleftharpoons[p_2]{k_2} [\mathrm{H}]_2

Multi-species

H = F.Species("H")
D = F.Species("D")
empty_traps = F.Species("empty", mobile=False)
trapped_H = F.Species("trapped_H", mobile=False)
trapped_D = F.Species("trapped_D", mobile=False)


my_model.species = [
  H, D, empty_traps,
  trapped_H, trapped_D
]

1 trap, 2 isotopes

\mathrm{H} + [\ ] \xrightleftharpoons[p]{k} [\mathrm{H}]
trapping = F.Reaction(
	reactant=[H, empty_traps],
    product=[trapped_H],
	...
)

trapping2 = F.Reaction(
	reactant=[D, empty_traps],
    product=[trapped_D],
	...
)


my_model.reactions = [trapping, trapping2]
\mathrm{D} + [\ ] \xrightleftharpoons[p]{k} [\mathrm{D}]

Multi-species

H = F.Species("H")
empty_traps = F.Species("empty", mobile=False)
trapped_H = F.Species("trapped_H", mobile=False)
trapped_H2 = F.Species("trapped_H2", mobile=False)

my_model.species = [
  H, empty_traps,
  trapped_H, trapped_H2
]

Multi-occupancy trapping

\mathrm{H} + [\ ] \xrightleftharpoons[p]{k} [\mathrm{H}]
trapping = F.Reaction(
	reactant=[H, empty_traps],
    product=[trapped_H],
	...
)

trapping2 = F.Reaction(
	reactant=[H, empty_traps2],
    product=[trapped_H2],
	...
)


my_model.reactions = [trapping, trapping2]
\mathrm{H} + [\mathrm{H}] \xrightleftharpoons[p]{k} [\mathrm{HH}]

Multi-species

A = F.Species("A")
B = F.Species("B")
C = F.Species("C")
D = F.Species("D")

my_model.species = [A, B, C, D]
\mathrm{A} + \mathrm{B} \xrightleftharpoons[p]{k} \mathrm{C} + \mathrm{D}

"Species" can be anything...

Multi-species

A = F.Species("A")
B = F.Species("B")
C = F.Species("C")
D = F.Species("D")

my_model.species = [A, B, C, D]
\mathrm{A} + \mathrm{B} \xrightleftharpoons[p]{k} \mathrm{C} + \mathrm{D}
reac = F.Reaction(
    reactant=[A, B],
    product=[C, D],
    ...
)

Arbitrary numbers of reactants and products

Multi-species

A = F.Species("A")
B = F.Species("B")


my_model.species = [A, B]
\mathrm{A} + \mathrm{B} \xrightarrow{k} \emptyset
reac = F.Reaction(
    reactant=[A, B],
    product=[],
    ...
)

Arbitrary numbers of reactants and products

Can be used for:

  • Vacancy-Interstitial anhiliation
  • Radioactive decay
  • Defect annealing
  • ...

Multi-species

H = F.Species("H")


my_model.species = [H]

reac = F.SurfaceReactionBC(
    reactant=[H, H],
    pressure=10,
    ...
)

my_model.boundary_conditions = [reac]

Surface reactions

2 \mathrm{H} \xrightleftharpoons[K_d]{K_r} \mathrm{H}_2

Soon: support for dynamically computed gas pressures (enclosures)

Examples: Multiphysics coupling

  

Integrate results from external solvers into FESTIM

Coupling CFD to FESTIM using foam2dolfinx

Using foam2dolfinx we can pass velocity fields calculated from OpenFOAM to FESTIM

Reading OpenFOAM data using foam2dolfinx

from foam2dolfinx import OpenFOAMReader

reader = OpenFOAMReader(filename="my_foam_case.foam", cell_type=10)


# read fields
u = reader.create_dolfinx_function_with_point_data(t=100, name="U")


# read mesh
foam_mesh = reader.dolfinx_meshes_dict["default"]


# read meshtags
facet_mt = reader.create_facet_meshtags()
volume_mt = reader.create_cell_meshtags()

Read a local OpenFOAM case file

 

 

 

 

  • Specify time stamp of data
  • Can extract any field type
  • Extract dolfinx mesh
  • Extract meshtags

Tets

Hexes

Support for:

 

 

 

 

Reading and writing using io4dolfinx

import io4dolfinx

checkpoint_file = Path("sim_checkpoint.bp")

io4dolfinx.write_mesh(checkpoint_file, msh)
io4dolfinx.write_function(checkpoint_file, u, time=100, name="U")
io4dolfinx.write_meshtags(checkpoint_file, msh, facet_mt, meshtag_name="facet_tags")
io4dolfinx.write_meshtags(checkpoint_file, msh, cell_mt, meshtag_name="cell_tags")

..

msh = io4dolfinx.read_mesh(checkpoint_file, COMM)
facet_mt = io4dolfinx.read_meshtags(checkpoint_file, msh, meshtag_name="facet_tags")
cell_mt = io4dolfinx.read_meshtags(checkpoint_file, msh, meshtag_name="cell_tags")

gdim = msh.geometry.dim
el = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1, shape=(gdim,))
w = fem.Function(fem.functionspace(msh, el), name="U")
io4dolfinx.read_function(checkpoint_file, u, time=100, name="U")

Using OpenFOAM data in FESTIM

...


my_model.mesh = foam_mesh

my_model.surface_meshtags = facet_meshtags
my_model.volume_meshtags = cell_meshtags


...


my_model.advection_terms = [
  F.AdvectionTerm(velocity=u, species=H, subdomain=fluid)
]

Pass mesh and meshtags to FESTIM

 

Pass velocity field using F.AdvectionTerm()

 

 

Adding advection to FESTIM means adding a term in the formulation

\frac{\partial c}{\partial t} = \nabla\cdot(D \nabla c) + S
- \nabla \cdot ( \mathbf{u} c)

Coupling neutronics to FESTIM using openmc2dolfinx

Using openmc2dolfinx we can pass fields calculated from OpenMC to FESTIM

from openmc2dolfinx import UnstructuredMeshReader, StructuredMeshReader

# UnstructuredMeshReader for tets
reader = UnstructuredMeshReader(filename="out.vtk")

# StructuredMeshReader for hexs
reader = StructuredMeshReader(filename="out.vtk")

# read fields
source_from_openmc = reader.create_dolfinx_function(name="mean")

Reading OpenMC data using openmc2dolfinx

Read .vtk file produced by OpenMC

 

 

 

  • Extract any mesh tally:
    • ​Tritium production
    • Heat generation
    • He production ...

Tets

Hexes

Support for:

 

 

 

 

...


my_model.sources = [
  F.ParticleSource(value=source_from_openmc, species=T, subdomain=volume)
]


...

Using OpenMC data in FESTIM

Input into FESTIM as a source term

Examples: FESTIM in the loop

Integrate FESTIM in complex workflows

Machine Learning with FESTIM

1. Parametric FESTIM model

2. Produce training data

3. Train surrogate/emulator

Parametric optimisation with FESTIM

Goal: find the trapping parameters to reproduce an experimental thermo-desorption spectrum

Integrating FESTIM with fuel cycle models

PathView/PathSim: a python package for system modelling

pathsim.org

 

Can be seamlessly coupled with FESTIM for multi-fidelity fuel cycle simulations

How to get started with FESTIM?

FESTIM website


Everything in one place!
festim.readthedocs.io

Fundamentals: Core physics and FESTIM building blocks.

Learning & Tutorials

Visit the tutorial
festim-workshop.readthedocs.io

Misc: troubleshooting, tips & tricks, "build FESTIM from scratch"

Learning & Tutorials

Visit the tutorial
festim-workshop.readthedocs.io

Applications: Full application cases like monoblock, TDS, multi-physics, machine learning...

Learning & Tutorials

Visit the tutorial
festim-workshop.readthedocs.io

FESTIM Fellowship Program

What it is

A selective, project-based
programme for experienced FESTIM users to develop advanced simulation expertise working directly alongside the MIT development team.

How it works

Fall 2026 Intake - Applications open now      |      Deadline: August 1st      Apply: festim.readthedocs.io/en/latest/fellowship.html

What you gain

~10 weeks, flexible, Remote or in person at MIT PSFC. Fellows bring a focused project and get guided technical support, code reviews, and hands-on collaboration.

Become an advanced FESTIM user. Capable of building complex models independently, contributing to the codebase, and supporting colleagues with FESTIM.

FESTIM Community

Slack Channel

Discourse

Monthly Dev meetings

Where we share updates on FESTIM activities and collaborate day-to-day. The best place to get involved and stay connected

Our public forum for questions, troubleshooting and discussion. Post any issue, and the team will repsond

Open Development meetings of the last Thursday of each month.

Next one May 28th!

On the Horizon

Gas Enclosures
Computes gas pressure dynamically.
Hydride Formation
Development already underway.
DG Methods
Discontinuous Galerkin for diffusion-advection.
Codimensional Meshes
Eg. coupling 1D and 2D problems.
Diagnostic Tools
Enhanced troubleshooting.
Interface Conditions
Flux conditions beyond Dirichlet.
FESTIM GUI
New visual companion tool.
Community Events
More workshops and meetups.
Interoperability
Seamless tool integration.