Issue with XDMF Export for trapped species

User Support
1

Hello there,

I am currently running simulations of Hydrogen transport (with multiple isotopes) in a 1D tungsten slab.
I managed to export and read the concentration profiles of mobile species using XDMFExport. However, I am getting an error every time I try to export the concentration profiles of trapped species using the same XDMFExport function.

I attach here what I tried to set up as MWE:

import numpy as np
import festim as F

my_model = F.HydrogenTransportProblem()

Domain length [m]

L = 12e-3
vertices_uniform = np.linspace(0, L, num=5000)
my_model.mesh = F.Mesh1D(vertices_uniform)

w_atom_density = 6.3e28 # [atoms/m^3] tungsten atomic density (hard-coded)
tungsten = F.Material(D_0=4.1e-7, E_D=0.39, name=“tungsten”)

vol = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten)
surf_left = F.SurfaceSubdomain1D(id=1, x=0)
surf_right = F.SurfaceSubdomain1D(id=2, x=L)
my_model.subdomains = [vol, surf_left, surf_right]

Two mobile species (D, T) and one trap with two occupied states

mobile_D = F.Species(“D”)
mobile_T = F.Species(“T”)
trap1_D = F.Species(“trap1_D”, mobile=False)
trap1_T = F.Species(“trap1_T”, mobile=False)

Empty trap population (implicit) shared by both isotopes

empty_t1 = F.ImplicitSpecies(
name=“empty_trap1”,
n=1.3e-3 * w_atom_density, # [m^-3] trap density
others=[trap1_D, trap1_T],
)

my_model.species = [mobile_D, mobile_T, trap1_D, trap1_T]

R0 = 4.1e-7 / (1.1e-10**2 * 6 * w_atom_density)

reactions = [

D + empty_trap1 → trap1_D

F.Reaction(
k_0=R0, E_k=0.39,
p_0=1e13, E_p=0.87,
reactant=[mobile_D, empty_t1],
product=[trap1_D],
volume=vol,
),

T + empty_trap1 → trap1_T

F.Reaction(
k_0=R0, E_k=0.39,
p_0=1e13, E_p=0.87,
reactant=[mobile_T, empty_t1],
product=[trap1_T],
volume=vol,
),
]
my_model.reactions = reactions

Temp = 500 #T in K
my_model.temperature = Temp

bc_D_value = 2e23
bc_T_value = 6e23

bc_D = F.DirichletBC(subdomain=surf_left, value=bc_D_value, species=“D”)
bc_T = F.DirichletBC(subdomain=surf_left, value=bc_T_value, species=“T”)
my_model.boundary_conditions = [bc_D, bc_T]

folder = “results”

my_model.exports = [
F.XDMFExport(f"{folder}/trap1_D.xdmf", field=trap1_D),
]

my_model.settings = F.Settings(
transient=True,
final_time=1000.0,
atol=1e10,
rtol=1e-10,
max_iterations=50,
stepsize=F.Stepsize(
1e-2,
growth_factor=1.1,
cutback_factor=0.3,
target_nb_iterations=5,
max_stepsize=1,
),
)
my_model.initialise()
my_model.run

Giving the following error:

Solving HydrogenTransportProblem: 0%| | 0.00/1.00k [00:00<?, ?it/s]Traceback (most recent call last):
File “/home/ITER/llealsa/FESTIM-2/discourse_xdmf_export/MWE.py”, line 88, in
my_model.run()
File “/home/ITER/llealsa/miniconda3/envs/festim2a/lib/python3.10/site-packages/festim/problem.py”, line 166, in run
self.iterate()
File “/home/ITER/llealsa/miniconda3/envs/festim2a/lib/python3.10/site-packages/festim/problem.py”, line 197, in iterate
self.post_processing()
File “/home/ITER/llealsa/miniconda3/envs/festim2a/lib/python3.10/site-packages/festim/hydrogen_transport_problem.py”, line 971, in post_processing
export.write(float(self.t))
File “/home/ITER/llealsa/miniconda3/envs/festim2a/lib/python3.10/site-packages/festim/exports/xdmf.py”, line 79, in write
self._writer.write_function(field.post_processing_solution, t)
File “/home/ITER/llealsa/miniconda3/envs/festim2a/lib/python3.10/site-packages/dolfinx/io/utils.py”, line 252, in write_function
super().write_function(getattr(u, “_cpp_object”, u), t, mesh_xpath)
RuntimeError: Function and Mesh dof layouts do not match. Maybe the Function needs to be interpolated?

Solving HydrogenTransportProblem: 0%| | 0.01/1.00k [00:00<5:37:20, 20.2s/it]

I am using FESTIM’s 2.0a4 version in a conda environment (installed via conda-forge).

I could run the simulations using VTXSpeciesExport for both mobile and trapped species, however I am having a hard time trying to correctly visualize and work with the output data.

Is it possible to export the concentration profiles of trapped species using the XDMFExport function? If so, I would really appreciate the help to solve my current issue.

Thank you very much for your help and best regards.

2

Ok so there are some other issues with the script here, but you have actually uncovered an interesting bug here. Will circle back once I have a better idea of whats going on here.

3

So the issue here in particular is that traps use Discontinuous Galerkin elements, which is why there is then a mismatch between the DOF in the function and in the mesh, as DG elements have more DOF. There is a quick fix for this, which is to override this choice with my_model._element_for_traps = "CG" command.

However, you should be fine with the VTXSpeciesExport, but this is not working correctly at the moment. I’ll make sure to comment back on this post when a PR has been opened.

Here is the full example with that:

import numpy as np
import festim as F

my_model = F.HydrogenTransportProblem()

# Domain length [m]
L = 12e-3
vertices_uniform = np.linspace(0, L, num=5000)
my_model.mesh = F.Mesh1D(vertices_uniform)

w_atom_density = 6.3e28  # [atoms/m^3] tungsten atomic density (hard-coded)
tungsten = F.Material(D_0=4.1e-7, E_D=0.39, name="tungsten")

vol = F.VolumeSubdomain1D(id=1, borders=[0, L], material=tungsten)
surf_left = F.SurfaceSubdomain1D(id=1, x=0)
surf_right = F.SurfaceSubdomain1D(id=2, x=L)
my_model.subdomains = [vol, surf_left, surf_right]

# Two mobile species (D, T) and one trap with two occupied states
mobile_D = F.Species("D")
mobile_T = F.Species("T")
trap1_D = F.Species("trap1_D", mobile=False)
trap1_T = F.Species("trap1_T", mobile=False)

# Empty trap population (implicit) shared by both isotopes
empty_t1 = F.ImplicitSpecies(
    name="empty_trap1",
    n=1.3e-3 * w_atom_density,  # [m^-3] trap density
    others=[trap1_D, trap1_T],
)

my_model.species = [mobile_D, mobile_T, trap1_D, trap1_T]

R0 = 4.1e-7 / (1.1e-10**2 * 6 * w_atom_density)

reactions = [
    # D + empty_trap1 → trap1_D
    F.Reaction(
        k_0=R0,
        E_k=0.39,
        p_0=1e13,
        E_p=0.87,
        reactant=[mobile_D, empty_t1],
        product=[trap1_D],
        volume=vol,
    ),
    # T + empty_trap1 → trap1_T
    F.Reaction(
        k_0=R0,
        E_k=0.39,
        p_0=1e13,
        E_p=0.87,
        reactant=[mobile_T, empty_t1],
        product=[trap1_T],
        volume=vol,
    ),
]
my_model.reactions = reactions

Temp = 500  # T in K
my_model.temperature = Temp

bc_D_value = 2e23
bc_T_value = 6e23

bc_D = F.DirichletBC(subdomain=surf_left, value=bc_D_value, species="D")
bc_T = F.DirichletBC(subdomain=surf_left, value=bc_T_value, species="T")
my_model.boundary_conditions = [bc_D, bc_T]

folder = "results"

my_model.exports = [
    F.XDMFExport(f"{folder}/trap1_D.xdmf", field=trap1_D),
]

my_model.settings = F.Settings(
    transient=True,
    final_time=1000.0,
    atol=1e10,
    rtol=1e-10,
    max_iterations=50,
    stepsize=F.Stepsize(
        1,
        growth_factor=1.1,
        cutback_factor=0.3,
        target_nb_iterations=5,
        max_stepsize=1,
    ),
)

my_model._element_for_traps = "CG"

my_model.initialise()
my_model.run()
4

Thank you very much for the quick response and for providing a clear solution.

Setting my trap elements to “CG" in my simulations did actually solve the issue and now I’m able to analyze and visualize profiles for both mobile and trapped species.

I’ll keep an eye on this thread for your update when the PR is opened, thanks for letting me know!

5

@AdriaLlealS I would recommend using the VTX format for exporting as the XDMF format is not maintained anymore and may soon be deprecated by FEniCS