Issue with F.XDMFExport and F.VTXSpeciesExport

Andrea · November 17, 2025, 4:39pm

Hi all,

I am working with three different materials — tungsten, EUROFER, and FLiBe — in a 2D geometry. I previously used FESTIM 1, but I was getting numerical errors, so I decided to move to FESTIM 2.

When running a transient simulation with:
my_model = F.HydrogenTransportProblem()

I obtain reasonable results for the mobile concentration, but I am not able to export the trapped tritium fields, neither with F.VTXSpeciesExport nor with F.XDMFExport (I am using DG and not CG). I would like to visualise the trapped concentrations to verify whether the numerical issues disappear in FESTIM 2.

I am attaching my mesh files in case you want to reproduce the issue:
mesh_domains.h5 (163.2 KB)
mesh_domains.xdmf (584 Bytes)
mesh_boundaries.h5 (101.5 KB)
mesh_boundaries.xdmf (587 Bytes)

The input data should be correct, as the corresponding 1D model works fine.

Below I include the full script I am using (sorry for the length, but I wanted to provide all details).

import numpy as np
import meshio
#import matplotlib.pyplot as plt
import festim as F
import os
import copy 
from dolfinx import plot

####
flag=1  # 1 per quad, 0 per triangoli
####

folder = "/home/andre/miniconda3/TESI/Pure_Diff/"
results_folder = "/home/andre/miniconda3/TESI/Pure_Diff_traps/"

#my_model =  F.HydrogenTransportProblemDiscontinuous() 
my_model = F.HydrogenTransportProblem()

if flag==1:
    mesh = F.MeshFromXDMF(folder +"quad/mesh_domains.xdmf", folder + "quad/mesh_boundaries.xdmf")
else:
    mesh = F.MeshFromXDMF(folder +"tria/mesh_domains.xdmf", folder + "tria/mesh_boundaries.xdmf")
my_model.mesh = mesh


tungsten_mat = F.Material(
    D_0=4.1e-7/(np.sqrt(3)),
    E_D=0.39,
    K_S_0=1.35e22,
    E_K_S=0.29,

)

euro_fer_mat = F.Material(
    D_0=1.07e-6/(np.sqrt(3)),
    E_D=0.51,
    K_S_0=1.35e22,
    E_K_S=0.16,
)

FLiBe_mat = F.Material(
    D_0=9.3e-7,
    E_D=42.5/96.48,
    K_S_0=(7.9e-2)*(6.022e23),
    E_K_S=35/96.48,
)
tungsten = F.VolumeSubdomain(id=8, material=tungsten_mat)
euro_fer = F.VolumeSubdomain(id=9, material=euro_fer_mat)
FLiBe = F.VolumeSubdomain(id=10, material=FLiBe_mat)




inlet_surf=F.SurfaceSubdomain(id=6, )
dirichelt_surf=F.SurfaceSubdomain(id=7)
neumann_surf=F.SurfaceSubdomain(id=12)

my_model.subdomains = [tungsten, euro_fer, FLiBe, inlet_surf, dirichelt_surf, neumann_surf]

"""
my_model.surface_to_volume = {
    inlet_surf: tungsten,     # o quel volume che tocca l’inlet
    dirichelt_surf: FLiBe,    # quel volume che tocca Dirichlet
    neumann_surf: FLiBe,      # o altro volume che tocca Neumann
}



# Interfaces
my_model.interfaces = [
    F.Interface(id=11, subdomains=[euro_fer, FLiBe], penalty_term=1e10),
    F.Interface(id=13, subdomains=[tungsten, euro_fer], penalty_term=1e10),
]
"""

import ufl

x = ufl.SpatialCoordinate(my_model.mesh.mesh)  # coord x del dominio
m1 = (1240 - 1300) / 0.01     # tra 0 e 0.01
m2 = (850 - 1240) / 0.03     # tra 0.01 e 0.04

T_expr = ufl.conditional(
    ufl.le(x[0], 0.01),
    1300.0 + m1 * x[0],
    ufl.conditional(
        ufl.le(x[0], 0.04),
        1240.1 + m2 * (x[0] - 0.01),
        850,
    ),
)

my_model.temperature = lambda x: T_expr


mobile_T = F.Species("T", subdomains=[tungsten, euro_fer, FLiBe], mobile=True)
trapped_T_1 = F.Species("T_trapped_1", mobile=False, subdomains=[tungsten])
trapped_T_2 = F.Species("T_trapped_2", mobile=False, subdomains=[tungsten])
trapped_T_3 = F.Species("T_trapped_3", mobile=False, subdomains=[tungsten])
trapped_T_4= F.Species("T_trapped_4", mobile=False, subdomains=[tungsten])
trapped_T_5 = F.Species("T_trapped_5", mobile=False, subdomains=[tungsten])
trapped_T_6 = F.Species("T_trapped_6", mobile=False, subdomains=[tungsten])
trapped_T_7 = F.Species("T_trapped_7", mobile=False, subdomains=[euro_fer])
trapped_T_8 = F.Species("T_trapped_8", mobile=False, subdomains=[euro_fer])
trapped_T_9 = F.Species("T_trapped_9", mobile=False, subdomains=[euro_fer])
trapped_T_10 = F.Species("T_trapped_10", mobile=False, subdomains=[euro_fer])
trapped_T_11 = F.Species("T_trapped_11", mobile=False, subdomains=[euro_fer])
empty_traps_1 = F.Species("empty_traps_1", mobile=False, subdomains=[tungsten])
empty_traps_2 = F.Species("empty_traps_2", mobile=False, subdomains=[tungsten])
empty_traps_3 = F.Species("empty_traps_3", mobile=False, subdomains=[tungsten])
empty_traps_4 = F.Species("empty_traps_4", mobile=False, subdomains=[tungsten])
empty_traps_5 = F.Species("empty_traps_5", mobile=False, subdomains=[tungsten])
empty_traps_6 = F.Species("empty_traps_6", mobile=False, subdomains=[tungsten])
empty_traps_7 = F.Species("empty_traps_7", mobile=False, subdomains=[euro_fer])
empty_traps_8 = F.Species("empty_traps_8", mobile=False, subdomains=[euro_fer])
empty_traps_9 = F.Species("empty_traps_9", mobile=False, subdomains=[euro_fer])
empty_traps_10 = F.Species("empty_traps_10", mobile=False, subdomains=[euro_fer])
empty_traps_11 = F.Species("empty_traps_11", mobile=False, subdomains=[euro_fer])

#my_model.species = [mobile_T]

my_model.species = [mobile_T, trapped_T_1,empty_traps_1,trapped_T_2, empty_traps_2, trapped_T_3, empty_traps_3, trapped_T_4, empty_traps_4, trapped_T_5, empty_traps_5, trapped_T_6, empty_traps_6, trapped_T_7, empty_traps_7, trapped_T_8, empty_traps_8, trapped_T_9, empty_traps_9, trapped_T_10, empty_traps_10, trapped_T_11, empty_traps_11]
w_atom_density = 6.3e28  # atom/m3
euro_atom_density = 7.8e3*6.022e23/(57.2*1e-3)   # atom/m3
"""
my_model.initial_conditions = [
    F.InitialConcentration(value=0.0, species=mobile_T, volume=tungsten),
    F.InitialConcentration(value=0.0, species=mobile_T, volume=euro_fer),
    F.InitialConcentration(value=0.0, species=mobile_T, volume=FLiBe),

    F.InitialConcentration(value=2.4e22, species=empty_traps_1, volume=tungsten),
    F.InitialConcentration(value=6.9e25, species=empty_traps_2, volume=tungsten),
    F.InitialConcentration(value=7e25, species=empty_traps_3, volume=tungsten),
    F.InitialConcentration(value=6e25, species=empty_traps_4, volume=tungsten),
    F.InitialConcentration(value=4.7e25, species=empty_traps_5, volume=tungsten),
    F.InitialConcentration(value=2e25, species=empty_traps_6, volume=tungsten),
    F.InitialConcentration(value=6.9e25, species=empty_traps_7, volume=euro_fer),
    F.InitialConcentration(value=7e25, species=empty_traps_8, volume=euro_fer),
    F.InitialConcentration(value=6e25, species=empty_traps_9, volume=euro_fer),
    F.InitialConcentration(value=4.7e25, species=empty_traps_10, volume=euro_fer),
    F.InitialConcentration(value=2e25, species=empty_traps_11, volume=euro_fer),

    F.InitialConcentration(value=0, species=trapped_T_1, volume=tungsten),
    F.InitialConcentration(value=0, species=trapped_T_2, volume=tungsten),
    F.InitialConcentration(value=0, species=trapped_T_3, volume=tungsten),
    F.InitialConcentration(value=0, species=trapped_T_4, volume=tungsten),
    F.InitialConcentration(value=0, species=trapped_T_5, volume=tungsten),
    F.InitialConcentration(value=0, species=trapped_T_6, volume=tungsten),
    F.InitialConcentration(value=0, species=trapped_T_7, volume=euro_fer),
    F.InitialConcentration(value=0, species=trapped_T_8, volume=euro_fer),
    F.InitialConcentration(value=0, species=trapped_T_9, volume=euro_fer),
    F.InitialConcentration(value=0, species=trapped_T_10, volume=euro_fer),
    F.InitialConcentration(value=0, species=trapped_T_11, volume=euro_fer),
]
"""
my_model.reactions = [
    F.Reaction(
        reactant=[mobile_T, empty_traps_1],
        product=[trapped_T_1],
        k_0=4.1e-7/(1.1e-10**2*6*w_atom_density),
        E_k=0.28,
        p_0=1.13,
        E_p=1.04,
        volume=tungsten,
    ),
    F.Reaction(
        reactant=[mobile_T, empty_traps_2],
        product=[trapped_T_2],
        k_0=4.1e-7/(1.1e-10**2*6*w_atom_density),
        E_k=0.39,
        p_0=1.13,
        E_p=1.15,
        volume=tungsten,
    ),

    F.Reaction(
        reactant=[mobile_T, empty_traps_3],
        product=[trapped_T_3],
        k_0=4.1e-7/(1.1e-10**2*6*w_atom_density),
        E_k=0.39,
        p_0=1.13,
        E_p=1.35,
        volume=tungsten,
    ),    
        F.Reaction(
        reactant=[mobile_T, empty_traps_4],
        product=[trapped_T_4],
        k_0=4.1e-7/(1.1e-10**2*6*w_atom_density),
        E_k=0.39,
        p_0=1.13,
        E_p=1.65,
        volume=tungsten,
    ),   
        F.Reaction(
        reactant=[mobile_T, empty_traps_5],
        product=[trapped_T_5],
        k_0=4.1e-7/(1.1e-10**2*6*w_atom_density),
        E_k=0.39,
        p_0=1.13,
        E_p=1.85,
        volume=tungsten,
    ),    F.Reaction(
        reactant=[mobile_T, empty_traps_6],
        product=[trapped_T_6],
        k_0=4.1e-7/(1.1e-10**2*6*w_atom_density),
        E_k=0.39,
        p_0=1.13,
        E_p=2.05,
        volume=tungsten,
    ),    F.Reaction(
        reactant=[mobile_T, empty_traps_7],
        product=[trapped_T_7],
        k_0=4.57e-7/(1.1e-10**2*6*euro_atom_density),
        E_k=0.45,
        p_0=1.13,
        E_p=1.15,
        volume=euro_fer,
    ),    F.Reaction(
        reactant=[mobile_T, empty_traps_8],
        product=[trapped_T_8],
        k_0=4.57e-7/(1.1e-10**2*6*euro_atom_density),
        E_k=0.45,
        p_0=1.13,
        E_p=1.35,
        volume=euro_fer,
    ),    F.Reaction(
        reactant=[mobile_T, empty_traps_9],
        product=[trapped_T_9],
        k_0=4.57e-7/(1.1e-10**2*6*euro_atom_density),
        E_k=0.45,
        p_0=1.13,
        E_p=1.65,
        volume=euro_fer,
    ),    F.Reaction(
        reactant=[mobile_T, empty_traps_10],
        product=[trapped_T_10],
        k_0=4.57e-7/(1.1e-10**2*6*euro_atom_density),
        E_k=0.45,
        p_0=1.13,
        E_p=1.85,
        volume=euro_fer,
    ),    F.Reaction(
        reactant=[mobile_T, empty_traps_11],
        product=[trapped_T_11],
        k_0=4.57e-7/(1.1e-10**2*6*euro_atom_density),
        E_k=0.45,
        p_0=1.13,
        E_p=2.05,
        volume=euro_fer,
    ),
]


my_model._element_for_traps = "DG"

import ufl
phi = 7.15e12
R_p = 9.52e-10
implatantion_flux=F.FixedConcentrationBC(
    subdomain=inlet_surf,
    #value=lambda T_expr: phi * R_p / (tungsten_mat.D_0 * ufl.exp(-tungsten_mat.E_D / F.k_B / T_expr)),
    value=1e14,
    species=mobile_T,
)
Neum_con=F.ParticleFluxBC(subdomain=neumann_surf, value=2.4e+19*2, species=mobile_T)

Diriclet_con=F.FixedConcentrationBC(subdomain=dirichelt_surf, value=4.6e-5*(6.022e23), species=mobile_T)

my_model.boundary_conditions=[implatantion_flux, Diriclet_con, Neum_con]

my_model.sources=[F.ParticleSource(volume=FLiBe, value=1.1e19, species=mobile_T)]

my_model.settings = F.Settings(
    atol=1e-2,
    rtol=1e-2,
    transient=True,
    final_time=1e8,
)

my_model.settings.stepsize = F.Stepsize(
    initial_value=1e7,       # passo iniziale piccolo = stabile
    growth_factor=1.4,        # aumenta lentamente
    cutback_factor=0.5,       # riduce solo quando serve
    target_nb_iterations=20,  # default stabile
)


if flag==1:
    my_model.exports = [
        F.XDMFExport(
            field=mobile_T,
            filename= results_folder +"quad_traps/hydrogen_concentration_festim2_q.xdmf"
        ),
        F.VTXSpeciesExport(filename="quad_traps/Trap_1.bp", field=trapped_T_1),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_2.bp", field=trapped_T_2),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_3.bp", field=trapped_T_3),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_4.bp", field=trapped_T_4),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_5.bp", field=trapped_T_5),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_6.bp", field=trapped_T_6),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_7.bp", field=trapped_T_7),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_8.bp", field=trapped_T_8),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_9.bp", field=trapped_T_9),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_10.bp", field=trapped_T_10),
        #F.VTXSpeciesExport(filename="quad_traps/Trap_11.bp", field=trapped_T_11),
        ]




       # F.XDMFExport(
        #    field=trapped_T,
         #   filename= results_folder +"quad_traps/hydrogen_trapped_T_festim2_t.xdmf"
        #),
        #F.XDMFExport(
         #   field=empty_traps_1,
          #  filename= results_folder +"quad_traps/hydrogen_1_festim2_t.xdmf"
        #),
    
else:
    #profile = ProfileExport(field=mobile_T, volume=[tungsten, euro_fer, FLiBe],)
    F.VTXSpeciesExport(filename= results_folder +"tria_traps/solute.bp", field=trapped_T_1)
    my_model.exports = [
        F.XDMFExport(
           field=mobile_T,
           filename=results_folder +"tria/hydrogen_concentration.xdmf",
                ),
]
        #F.XDMFExport(
         #   field=trapped_T,
          #  filename= results_folder +"tria_traps/hydrogen_trapped_T_festim2_t.xdmf"
        #),
        #F.XDMFExport(
         #   field=empty_traps_1,
          #  filename= results_folder +"tria_traps/hydrogen_1_festim2_t.xdmf"
        #),
    





my_model.initialise()
my_model.run()

import dolfinx 
from dolfinx import fem, mesh 
from ufl import Measure, FacetNormal, jump, grad, Constant, dot 
from mpi4py import MPI 
from petsc4py import PETSc 
domain = my_model.mesh.mesh 
facet_markers = my_model.mesh.define_surface_meshtags() 
dS = Measure("dS", domain=domain, subdomain_data=facet_markers) 

n = FacetNormal(domain) 
c = mobile_T.solution 

D_euro = euro_fer_mat.D_0 
D_flibe = FLiBe_mat.D_0 
D = fem.Constant(domain, PETSc.ScalarType(0.5 * (euro_fer_mat.D_0 + FLiBe_mat.D_0))) 

flux_form = jump(-D * grad(c), n) * dS(11) 
flux_total = fem.assemble_scalar(fem.form(flux_form)) 
flux_total = domain.comm.allreduce(flux_total, op=MPI.SUM) 

area_form = 1.0 * dS(11)
area_dS11 = fem.assemble_scalar(fem.form(area_form)) 
area_dS11 = domain.comm.allreduce(area_dS11, op=MPI.SUM) # Flusso medio 
flux_mean = flux_total / area_dS11 
if domain.comm.rank == 0: 
    print(f"Flusso totale su dS(11): {flux_total:.4e} [atomi/s]") 
    print(f"Area della superficie dS(11): {area_dS11:.4e} [m²]") 
    print(f"Flusso medio su dS(11): {flux_mean:.4e} [atomi/(m²·s)]")

Would F.HydrogenTransportProblemDiscontinuous() be more appropriate for a multi-material 2D problem with several trapping sites?

Thanks in advance for any help or suggestions!

Andrea

remidm · November 17, 2025, 4:59pm

Hi @Andrea

It is very hard for anyone to help you on a several hundred line-long script. Especially when the script isn’t easily reproducible (eg. you have paths to your own file system)
consider writing minimal working examples that reproduce your error (here: not being able to export trapped fields). A MWE for this should be very small: a 1D problem, two materials, 1 trapped concentration, steady state, exports.
F.HydrogenTransportProblem() is suitable for problems with continuity of concentration, which isn’t the case for you. You should consider using F.HydrogenTransportProblemDiscontinuous() to account for interface discontinuities

Topic		Replies	Views
Issue with XDMF Export for trapped species User Support festim2	4	58	October 16, 2025
Problems when reproducing the results of the article User Support	20	208	May 31, 2024
Point is not inside the domain when exporting surface flux User Support	17	93	March 26, 2025
Wrong solution with 2D simulation User Support	15	112	January 29, 2025
Double exports in TXTExport User Support	8	111	December 21, 2024

Issue with F.XDMFExport and F.VTXSpeciesExport

Related topics