Coupling to neutronics

Coupling to neutronics#

This tutorial discusses how to couple FESTIM to neutronics solvers such as OpenMC.

Objectives

  • Understanding the importance of neutronics in hydrogen transport

  • Utilizing external neutronics solvers to generate tritium source terms

  • Reading OpenMC results into FESTIM for accurate particle sources

Where does neutronics couple into hydrogen transport?#

Let us take a look at the governing equation for FESTIM, neglecting reactions and advection:

\[ \frac{\partial c_i}{\partial t} = \nabla\cdot (D_i \nabla c_i) + S_i(\mathbf{x}, t, T) \]

Neutronics solvers can provide a value for the \(S_i\), which represents a volumetric source of mobile hydrogen isotope that can depend on space, time and temperature. In particular, solvers such as OpenMC can provide an accurate source of of tritium generation.

Converting data from OpenMC to FESTIM#

OpenMC can tally spatial distributions of reaction rates and export the results in VTK format. Users can then use the openmc2dolfinx package to convert them into DOLFINx functions.

See also

Check out the OpenMC documentation and openmc2dolfinx repo to learn more.

Solving a hydrogen transport problem in an neutron-irradiated lithium cube#

This section discusses how to setup and solve a coupled problem between OpenMC and FESTIM, as outlined in the FESTIM 2 review paper.

Specifically, we’ll calculate the tritium generation in a neutron-irradiated lithium cube using OpenMC. Then, we will export the source term to FESTIM to solve our diffusion model.

See also

Check out the FESTIM 2 review paper to learn more.

Getting tritium source from OpenMC#

First, run the OpenMC model:

import openmc
import openmc_data_downloader as odd

dim = 60
lithium = openmc.Material(name="lithium")
lithium.set_density("g/cc", 0.534)
lithium.add_element("Li", 1.0)

mats = openmc.Materials([lithium])

odd.download_cross_section_data(
    mats,
    libraries=["FENDL-3.1d"],
    set_OPENMC_CROSS_SECTIONS=True,
    particles=["neutron"],
)

# Geometry
cube_surface = openmc.model.RectangularParallelepiped(
    -dim, dim, -dim, dim, -dim, dim
)
region = -cube_surface
cell = openmc.Cell(region=region, fill=lithium)

vacuum_surf = openmc.Sphere(r=dim * 2, boundary_type="vacuum")
vacuum_region = +cube_surface & -vacuum_surf
vacuum = openmc.Cell(region=vacuum_region, fill=None)

geometry = openmc.Geometry([cell, vacuum])

# Tallies
tally = openmc.Tally(name="tritium_production")
tally.scores = ["(n,Xt)"]
mesh = openmc.RegularMesh()
mesh.dimension = [30, 30, 15]
mesh.lower_left = [-dim, -dim, -dim]
mesh.upper_right = [dim, dim, dim]
tally.filters = [openmc.MeshFilter(mesh)]

tallies = openmc.Tallies([tally])

# Settings
source = openmc.IndependentSource()
source_pos_z = dim + 10
source.space = openmc.stats.Point((0, 0, source_pos_z))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14.1e6], [1.0])

settings = openmc.Settings()
settings.run_mode = "fixed source"
settings.source = source
settings.batches = 10
settings.particles = 10000

my_model = openmc.Model(
    geometry=geometry, settings=settings, materials=mats, tallies=tallies
)

my_model.run(apply_tally_results=True)

mesh.write_data_to_vtk(
    "tritium_production_mesh.vtk", {"tritium_production": tally.mean}
)

This should output a file named tritium_production_mesh.vtk which can then be read into FESTIM.

Tip

We provide the code to run the OpenMC model, but don’t run it this tutorial. To run this code, create a new environment and install OpenMC with the following lines of code:

conda config --add channels conda-forge
conda config --set channel_priority strict
conda create --name openmc-env openmc
conda activate openmc-env

It is also necessary to install vtk and openmc_data_downloader, which can be done with:

conda install -c conda-forge vtk
pip install openmc_data_downloader

Let’s take a look at the OpenMC tritium source output:

Hide code cell source

 
import pyvista as pv

pv.set_jupyter_backend("html")

grid = pv.read("openmc_data/tritium_production_mesh.vtk")
scalar_name = "tritium_production"  
grid.set_active_scalars(scalar_name)

plotter = pv.Plotter()
plotter.add_mesh(grid, scalars=scalar_name, cmap="hot", show_edges=False)
plotter.view_xy()
plotter.show()

2026-05-22 15:17:13.678 (   0.712s) [    7CB7E91D3740]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

Read data from OpenMC output#

Now that we have our tritium source data, we can use openmc2dolfinx to create a structured grid reader that outputs a dolfinx.fem.Function to be used in FESTIM:

from openmc2dolfinx import StructuredGridReader

reader = StructuredGridReader("openmc_data/tritium_production_mesh.vtk")
reader.create_dolfinx_mesh()

Setting up our FESTIM model#

Let us set up our FESTIM model now. We will create a refined box mesh and use a dummy material. We will consider insulated boundary conditions (\(c=0\) on all boundaries):

import festim as F
import dolfinx 
from mpi4py import MPI

dim = 60
my_model = F.HydrogenTransportProblem()

tritium = F.Species("T")

refined_mesh = dolfinx.mesh.create_box(
    MPI.COMM_WORLD,
    points=[(-dim, -dim, -dim), (dim, dim, dim)],
    n=(30, 30, 70),
    cell_type=dolfinx.mesh.CellType.tetrahedron,
)

my_model.mesh = F.Mesh(refined_mesh)
boundary = F.SurfaceSubdomain(id=1)

my_mat = F.Material(D_0=1, E_D=0)
volume = F.VolumeSubdomain(id=1, material=my_mat)
my_model.subdomains = [boundary, volume]

my_model.species = [tritium]
my_model.boundary_conditions = [F.FixedConcentrationBC(subdomain=boundary, value=0.0, species=tritium)]
my_model.temperature = 300
my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)

/home/docs/checkouts/readthedocs.org/user_builds/festim-workshop/conda/latest/lib/python3.12/site-packages/festim/coupled_heat_hydrogen_problem.py:1: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  import tqdm.autonotebook

Note

Here we chose to create more refined mesh from the hydrogen transport model than was used in OpenMC. We could, however, use the same mesh from OpenMC by running:

my_model.mesh = F.Mesh(reader.dolfinx_mesh)

Since we chose not to use the OpenMC mesh, we need to interpolate the results onto our new FESTIM mesh:

original_source_term = reader.create_dolfinx_function("tritium_production")
V = dolfinx.fem.functionspace(refined_mesh, ("DG", 0))
source_term = dolfinx.fem.Function(V)
F.helpers.nmm_interpolate(source_term, original_source_term)

Now we can add our interpolated tritium source into our problem:

tritium_source = F.ParticleSource(
    value=source_term, species=tritium, volume=volume
)

my_model.sources = [tritium_source]

Finally, let’s solve and take a look at the results:

my_model.initialise()
my_model.run()

Hide code cell source

 
from dolfinx import plot
import pyvista

c = tritium.post_processing_solution

topology, cell_types, geometry = plot.vtk_mesh(c.function_space)
u_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
u_grid.point_data["c"] = c.x.array.real
u_grid.set_active_scalars("c")
u_plotter = pyvista.Plotter()

u_plotter.add_mesh(u_grid, cmap="viridis", show_edges=False)
u_plotter.view_xy()

if not pyvista.OFF_SCREEN:
    u_plotter.show()
else:
    figure = u_plotter.screenshot("concentration.png")

As we’d expect, there is no diffusion on the outer boundaries of our cube (recall our insulated boundary conditions). Let’s take a look inside the cube to how the tritium source diffuses inside the material:

Hide code cell source

 
bounds = [0, dim, 0, dim, 0, dim]

clipped = u_grid.clip_box(bounds, invert=True)

pl = pv.Plotter()
pl.add_mesh(clipped, cmap="viridis", show_edges=False)
pl.view_xy()
pl.show()

We see a high concentration where the tritium is generated and it diffuses throughout the material. Success!