# Issue with F.XDMFExport and F.VTXSpeciesExport

**URL:** <https://festim.discourse.group/t/issue-with-f-xdmfexport-and-f-vtxspeciesexport/151>\
**Category:** User Support\
**Tags:** festim2\
**Created:** [November 17, 2025, 4:39pm UTC](https://festim.discourse.group/t/issue-with-f-xdmfexport-and-f-vtxspeciesexport/151 "2025-11-17T16:39:34Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Andrea](https://avatars.discourse-cdn.com/v4/letter/a/c37758/32.png) [@Andrea](https://festim.discourse.group/u/Andrea)\
**Post date:** [November 17, 2025, 4:39pm UTC](https://festim.discourse.group/t/issue-with-f-xdmfexport-and-f-vtxspeciesexport/151/1 "2025-11-17T16:39:34Z")

</div>

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](https://festim.discourse.group/uploads/short-url/lUWOJfrqOg3rtratwllPRhPmHAA.h5) (163.2 KB)  
[mesh\_domains.xdmf](https://festim.discourse.group/uploads/short-url/rCc1d8WPme6T81TMNGFjBBcMpKy.xdmf) (584 Bytes)  
[mesh\_boundaries.h5](https://festim.discourse.group/uploads/short-url/wGY8hvEKZ8ilP0HMQhCAlq6whtv.h5) (101.5 KB)  
[mesh\_boundaries.xdmf](https://festim.discourse.group/uploads/short-url/nx35dQuuVJ4prXYBmyRAkAeMSrI.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).

```auto
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)]")

```

1. **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

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [November 17, 2025, 4:59pm UTC](https://festim.discourse.group/t/issue-with-f-xdmfexport-and-f-vtxspeciesexport/151/2 "2025-11-17T16:59:32Z")

</div>

Hi @Andrea

1. 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.
2. `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