Concentration

Concentration#

Initial conditions are an important part of transient simulations for both hydrogen transport and heat-transfer problems. FESTIM defaults to zero values for initial conditions, but can be set using FESTIM’s InitialConditionBase.

This tutorial discusses how to set initial conditions for particle concentrations.

Objectives:

  • Setting the initial concentration for a single species

  • Defining initial concentration conditions for multiple species

Setting the initial concentration for a single species#

The InitialConcentration class can be used for defining initial conditions for the concentrations, which must defined on a volume subdomain in a FESTIM simulation:

import festim as F

material = F.Material(D_0=1, E_D=0)
vol = F.VolumeSubdomain(id=1, material=material)
H = F.Species("H")

IC = F.InitialConcentration(value=2, species=H, volume=vol)

Initial conditions can also be a function of space x and temperature T, such as:

\[\text{IC(x, y, T)} = 2x + 3y + 5T\]
my_function = lambda x,T: 2*x[0] + 3*x[1] + 5*T
IC = F.InitialConcentration(value=my_function, species=H, volume=vol)

Multi-species initial conditions#

The same class InitialConcentration can be used for multiple species, but a separate initial condition is required for each species.

Consider the following 2D example, where we have three species (hydrogen, deuterium, and tritium) with initial and boundary conditions for species’ concentrations:

\[ \text{IC (hydrogen)} = 4 \]
\[ \text{IC (deuterium)} = 5 \]
\[ \text{IC (tritium)} = 6 \]
\[ \text{Concentration BC (left, right, top) = 0} \]

First, we can create a 2D unit square mesh and mark the subdomains:

import festim as F
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import numpy as np

nx, ny = 50, 50
mesh = create_unit_square(MPI.COMM_WORLD, nx, ny)

model = F.HydrogenTransportProblem()
model.mesh = F.Mesh(mesh)

material = F.Material(D_0=1e-2, E_D=0)

top_surface = F.SurfaceSubdomain(id=1, locator = lambda x: np.isclose(x[1], 1.0))
bottom_surface = F.SurfaceSubdomain(id=2, locator = lambda x: np.isclose(x[1], 0.0))
right_surface = F.SurfaceSubdomain(id=3, locator = lambda x: np.isclose(x[0], 1.0))
left_surface = F.SurfaceSubdomain(id=4, locator = lambda x: np.isclose(x[0], 0.0))

vol = F.VolumeSubdomain(id=1, material=material)

model.subdomains = [top_surface, bottom_surface, left_surface, right_surface, vol]

Now, let’s define our species and initial conditions. We must create an InitialConcentration for each species, then we can pass the initial conditions into the initial_conditions attribute:

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

IC_H = F.InitialConcentration(value=4, species=H, volume=vol)
IC_D = F.InitialConcentration(value=5, species=D, volume=vol)
IC_T = F.InitialConcentration(value=6, species=T, volume=vol)

model.species = [H, D, T]
model.initial_conditions = [IC_H, IC_D, IC_T]

Note

To pass the initial conditions into FESTIM, you must use a list!

We also define the zero concentration boundary conditions for each species. Let’s run the simulation for 5 seconds with a stepsize of 1 second, and visualize the final concentration profiles (hydrogen on the left, deuterium in the middle, and tritium on the right):

model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=left_surface, value=0, species=H),
    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=H),
    F.FixedConcentrationBC(subdomain=top_surface, value=0, species=H),
    F.FixedConcentrationBC(subdomain=left_surface, value=0, species=D),
    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=D),
    F.FixedConcentrationBC(subdomain=top_surface, value=0, species=D),
    F.FixedConcentrationBC(subdomain=left_surface, value=0, species=T),
    F.FixedConcentrationBC(subdomain=right_surface, value=0, species=T),
    F.FixedConcentrationBC(subdomain=top_surface, value=0, species=T),
]
model.temperature = 400
model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=50, stepsize=5)

model.initialise()
model.run()

Hide code cell source

 
import pyvista
from dolfinx import plot

def make_ugrid(solution):
    topology, cell_types, geometry = plot.vtk_mesh(solution.function_space)
    u_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    u_grid.point_data["c"] = solution.x.array.real
    u_grid.set_active_scalars("c")
    return u_grid

pyvista.set_jupyter_backend("html")

H_grid = make_ugrid(H.post_processing_solution)
D_grid = make_ugrid(D.post_processing_solution)
T_grid = make_ugrid(T.post_processing_solution)

pl = pyvista.Plotter(shape=(1, 3))
pl.subplot(0, 0)
_ = pl.add_mesh(H_grid,label='h')
pl.view_xy()

pl.subplot(0, 1)
_ = pl.add_mesh(D_grid,label='d')
pl.view_xy()

pl.subplot(0, 2)
_ = pl.add_mesh(T_grid,label='t')
pl.view_xy()

pl.show()

2026-09-05 11:27:34.130 (   0.812s) [    709126556440]vtkXOpenGLRenderWindow.:1416  WARN| bad X server connection. DISPLAY=

Using a dolfinx.fem.Function as initial condition#

It is possible to use a dolfinx.fem.Function() as initial condition. This is particularly useful when using the result of a previous simulation for a restart.

from mpi4py import MPI
import dolfinx


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

V = dolfinx.fem.functionspace(mesh, ("CG", 1))

u_ini = dolfinx.fem.Function(V)
u_ini.interpolate(lambda x: 2e4 * (x[0] * (1 - x[0]) * x[1] * (1 - x[1]))**4)

Hide code cell source

 
import pyvista
from dolfinx import plot

pyvista.set_jupyter_backend("html")

plotter = pyvista.Plotter()
topology, cell_types, geometry = plot.vtk_mesh(u_ini.function_space)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.point_data["u_ini"] = u_ini.x.array
warped = grid.warp_by_scalar()
plotter.add_mesh(warped, show_edges=True, scalars="u_ini", cmap="viridis")
plotter.show()

import festim as F

my_model = F.HydrogenTransportProblem()

my_model.mesh = F.Mesh(mesh)

material = F.Material(D_0=0.1, E_D=0)
vol = F.VolumeSubdomain(id=1, material=material)
my_model.subdomains = [vol]

H = F.Species("H")

my_model.species = [H]
my_model.initial_conditions = [F.InitialConcentration(value=u_ini, species=H, volume=vol)]
my_model.temperature = 300
my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=0.05, final_time=0.3)

my_model.initialise()
my_model.run()

Hide code cell source

 
plotter = pyvista.Plotter()
u_final = H.post_processing_solution
topology, cell_types, geometry = plot.vtk_mesh(u_final.function_space)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.point_data["u_final"] = u_final.x.array
warped = grid.warp_by_scalar()
plotter.add_mesh(warped, show_edges=True, scalars="u_final", cmap="viridis", clim=[0, 0.3])
plotter.show()

Initial condition from a checkpoint file#

A simulation can be restarted from a checkpoint written by an earlier one. Unlike the visualisation formats, a checkpoint stores the field in its own function space, so it comes back exactly as it was written – see Checkpointing.

Here we run a short simulation and checkpoint the hydrogen field with format="checkpoint":

Hide code cell source

 
from mpi4py import MPI
import dolfinx
import numpy as np
import festim as F

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

first_model = F.HydrogenTransportProblem()
first_model.mesh = F.Mesh(mesh)
vol = F.VolumeSubdomain(id=1, material=F.Material(D_0=1e-2, E_D=0))
left = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 0.0))
right = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1.0))
first_model.subdomains = [vol, left, right]

H = F.Species("H")
first_model.species = [H]
first_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=left, value=1, species=H),
    F.FixedConcentrationBC(subdomain=right, value=0, species=H),
]
first_model.temperature = 400
first_model.settings = F.Settings(atol=1e-10, rtol=1e-10, stepsize=1, final_time=10)

first_model.exports = [
    F.SpeciesExport("checkpoint_file.bp", field=H, subdomain=vol, format="checkpoint")
]
first_model.initialise()
first_model.run()

festim.read_function_from_file() reads it back, given the name of the species and the timestamp to read:

u_in = F.read_function_from_file(
    filename="checkpoint_file.bp", name="H", timestamp=10.0
)
print(u_in.x.array[:])

[0.03005561 0.         0.         0.02977179 0.06506472 0.
 0.06470054 0.02960064 0.1104765  0.         0.11010279 0.06442949
 0.02947718 0.17153586 0.         0.17120929 0.10978874 0.06421289
 0.0293761  0.25296247 0.         0.25274836 0.17090498 0.10951651
 0.06402678 0.02928454 0.35826639 0.         0.35823039 0.25250925
 0.17061864 0.10927084 0.0638548  0.02919326 0.48884177 0.
 0.48902436 0.35810715 0.25225337 0.17034557 0.10903844 0.06368324
 0.02909279 0.64311408 0.         0.64349284 0.48904285 0.35792221
 0.25198918 0.17007986 0.10880639 0.06349798 0.02896965 0.81613105
 0.         0.8165577  0.64362717 0.48895119 0.35770136 0.25172209
 0.16981423 0.1085613  0.06328158 0.02879758 1.         0.
 1.         0.81670677 0.64361798 0.48879679 0.35746461 0.25145464
 0.16954093 0.10828856 0.06300844 0.02851048 1.         0.81673189
 0.64353313 0.48861462 0.3572271  0.25118911 0.16925292 0.10797101
 0.06263777 1.         0.81669986 0.64341609 0.48843155 0.35700398
 0.2509304  0.16894419 0.10758865 1.         0.81664501 0.64329842
 0.48827471 0.3568156  0.25068708 0.16860952 1.         0.81658998
 0.64321205 0.48818014 0.35668965 0.25046902 1.         0.8165578
 0.64320175 0.48819822 0.35665636 1.         0.81658388 0.64333797
 0.48838951 1.         0.81673631 0.64372681 1.         0.81716812
 1.        ]

Note

Pass backend="h5py" if the checkpoint was written to a plain .h5 file rather than the ADIOS2 default. Only checkpoints can be read back this way: the visualisation formats store values interpolated onto the mesh nodes rather than the degrees of freedom.

u_in is an ordinary dolfinx.fem.Function, so it can be used as an initial condition for a species exactly as explained in Using a dolfinx.fem.Function as initial condition:

my_model.initial_conditions = [
    F.InitialConcentration(value=u_in, species=H, volume=vol)
]

Warning

The function read back lives on the mesh stored in the checkpoint, which is a different object from any mesh you build in the new script – even one with identical geometry. Build the restarted model on that mesh, reachable as u_in.function_space.mesh, rather than on a freshly created one: interpolating between two distinct meshes is not supported and raises.