Hi there,
I am attempting to replicate some results from a recent paper entitled entitled “Modeling hydrogen isotope transport behavior in the first wall: role of the bonding region in W-RAFM steel with an improved evaluation of transport coefficients” (see also attached) and I was curious to better understand the results of simulated retention in the first wall. I have attempted to replicate the retention results shown in Figure 8b, but my current outputs (left) are about twice as high as the values shown there (right):
As you can see, each piece is consistently larger in my plot, with the W and bonding regions having ~2.5x increase in retention, and the RAFM having a ~1.3x increase. There are a few areas I think my work could be going astray:
The paper mentions they use an implantation flux of 1.5e21/m^2/s to set the concentration on the plasma-facing edge of the tungsten, according to the expression
c_m = \frac{\phi_{imp} R_p}{D}
where Rp = 5.6e-9 m is the implantation depth and D = D(T) is the arrhenius diffusion coefficient. As described here (see Eq. 16) I am using this to set a Dirichlet boundary condition.
Since the paper doesn’t provide jump frequency or detrapping pre-exponentials k_0 or p_0, respectively, I’m also currently assuming that
for each domain, with lambda = 1.1e-10 and N_surf = 6 across all 3 domains. For p_0, I’m using 1.0e13 across all of them. Therefore, k0 and p0 are both rather rough approximations at the moment, and may be the cause of the discrepancy.
For reference, my full code for the plot above is attached. create_multilayer.py constructs a quasi-1D geometry in two dimensions using the gmsh package:
import gmsh
from dolfinx.io import gmsh as gmshio
from mpi4py import MPI
from dataclasses import dataclass
@dataclass
class MeshTags:
"""Maps logical names to Gmsh physical group IDs."""
top_surface: int = 1
bottom_surface: int = 2
left_surface: int = 3
right_surface: int = 4
left_volume: int = 5
middle_volume: int = 6
right_volume: int = 7
interface_left: int = 8 # between left and middle volumes
interface_right: int = 9 # between middle and right volumes
def create_3layer_mesh(width_left=0.4, width_middle=0.1,
width_right=0.5, height=1.0,
mesh_size=0.02):
"""
Create a 2D mesh of a rectangle split vertically into 3 layers
(left | middle | right), using Gmsh.
Returns (dolfinx mesh_data, MeshTags) ready for FESTIM.
"""
tags = MeshTags()
total_width = width_left + width_middle + width_right
x1 = width_left # left/middle interface x-coord
x2 = width_left + width_middle # middle/right interface x-coord
gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 1)
gmsh.model.add("three_layer")
# Corner + interface points
# Bottom row (y=0), left to right
p_bl = gmsh.model.occ.addPoint(0, 0, 0, meshSize=mesh_size)
p_b1 = gmsh.model.occ.addPoint(x1, 0, 0, meshSize=mesh_size)
p_b2 = gmsh.model.occ.addPoint(x2, 0, 0, meshSize=mesh_size)
p_br = gmsh.model.occ.addPoint(total_width, 0, 0, meshSize=mesh_size)
# Top row (y=height), left to right
p_tl = gmsh.model.occ.addPoint(0, height, 0, meshSize=mesh_size)
p_t1 = gmsh.model.occ.addPoint(x1, height, 0, meshSize=mesh_size)
p_t2 = gmsh.model.occ.addPoint(x2, height, 0, meshSize=mesh_size)
p_tr = gmsh.model.occ.addPoint(total_width, height, 0, meshSize=mesh_size)
# Horizontal edges, split per layer so we can tag top/bottom
bottom_left = gmsh.model.occ.addLine(p_bl, p_b1)
bottom_middle = gmsh.model.occ.addLine(p_b1, p_b2)
bottom_right = gmsh.model.occ.addLine(p_b2, p_br)
top_left = gmsh.model.occ.addLine(p_tl, p_t1)
top_middle = gmsh.model.occ.addLine(p_t1, p_t2)
top_right = gmsh.model.occ.addLine(p_t2, p_tr)
# Vertical edges (outer walls + internal interfaces)
left_edge = gmsh.model.occ.addLine(p_bl, p_tl)
interface_left = gmsh.model.occ.addLine(p_b1, p_t1)
interface_right = gmsh.model.occ.addLine(p_b2, p_t2)
right_edge = gmsh.model.occ.addLine(p_br, p_tr)
# Curve loops and surfaces (all CCW)
left_loop = gmsh.model.occ.addCurveLoop(
[bottom_left, interface_left, -top_left, -left_edge]
)
left_surf = gmsh.model.occ.addPlaneSurface([left_loop])
middle_loop = gmsh.model.occ.addCurveLoop(
[bottom_middle, interface_right, -top_middle, -interface_left]
)
middle_surf = gmsh.model.occ.addPlaneSurface([middle_loop])
right_loop = gmsh.model.occ.addCurveLoop(
[bottom_right, right_edge, -top_right, -interface_right]
)
right_surf = gmsh.model.occ.addPlaneSurface([right_loop])
gmsh.model.occ.synchronize()
# Physical groups
# Volumes (2D surfaces)
gmsh.model.addPhysicalGroup(2, [left_surf], tag=tags.left_volume)
gmsh.model.setPhysicalName(2, tags.left_volume, "left_volume")
gmsh.model.addPhysicalGroup(2, [middle_surf], tag=tags.middle_volume)
gmsh.model.setPhysicalName(2, tags.middle_volume, "middle_volume")
gmsh.model.addPhysicalGroup(2, [right_surf], tag=tags.right_volume)
gmsh.model.setPhysicalName(2, tags.right_volume, "right_volume")
# Outer boundaries (grouped across layers where appropriate)
gmsh.model.addPhysicalGroup(
1, [top_left, top_middle, top_right], tag=tags.top_surface
)
gmsh.model.setPhysicalName(1, tags.top_surface, "top_surface")
gmsh.model.addPhysicalGroup(
1, [bottom_left, bottom_middle, bottom_right], tag=tags.bottom_surface
)
gmsh.model.setPhysicalName(1, tags.bottom_surface, "bottom_surface")
gmsh.model.addPhysicalGroup(1, [left_edge], tag=tags.left_surface)
gmsh.model.setPhysicalName(1, tags.left_surface, "left_surface")
gmsh.model.addPhysicalGroup(1, [right_edge], tag=tags.right_surface)
gmsh.model.setPhysicalName(1, tags.right_surface, "right_surface")
# Internal interfaces
gmsh.model.addPhysicalGroup(1, [interface_left], tag=tags.interface_left)
gmsh.model.setPhysicalName(1, tags.interface_left, "interface_left")
gmsh.model.addPhysicalGroup(1, [interface_right], tag=tags.interface_right)
gmsh.model.setPhysicalName(1, tags.interface_right, "interface_right")
# Generate mesh & import into dolfinx
model_rank = 0
gmsh.model.mesh.generate(2)
mesh_data = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank)
gmsh.finalize()
return mesh_data, tags
Shen2026Replication.py was made with FESTIM 2.0, specifically version 2.0b2.post2 (which I believe is the most current). It is pasted below:
import festim as F
import numpy as np
import matplotlib.pyplot as plt
import dolfinx, ufl
from mpi4py import MPI
from scipy.constants import k, e
from FESTIM.create_multilayer import create_3layer_mesh
# Physical constants / global configuration
k_B = k / e # Boltzmann constant [eV/K]
T_sim = 823.0 # FW temperature [K] (section 4.2)
# Implantation (eq. 17, section 2.3)
phi_imp = 1.5e21 # [m^-2 s^-1]
R_p = 5.6e-9 # [m]
# W diffusivity (intrinsic, values from section 4 of paper)
D0_W, ED_W = 2.40e-7, 0.20 # [m^2/s], [eV]
# Surface concentration approximation for implantation (eq. 17): c_m = phi*R_p/D
D_W_at_T = D0_W * np.exp(-ED_W / (k_B * T_sim))
c_surface_W = phi_imp * R_p / D_W_at_T
# Layer thicknesses for the FW (1D model, section 4.2): W | bond | RAFM
L_W = 2.0e-3 # 2 mm W armour
L_BOND = 0.1e-3 # 0.1 mm bonding region
L_RAFM = 20.0e-3 # 20 mm RAFM structural steel
HEIGHT = 0.1e-3 # height of slab => TotalVolume ~ atoms/m^2 of FW
# Atomic site densities (used to build k_0 = D_0 / (lambda^2 * n_surf * N_L))
LAMBDA = 1.1e-10 # jump distance [m]
N_SURF = 6.0 # coordination number
N_L_W = 6.3e28 # W atomic density [m^-3] (paper uses 6.3e28)
N_L_RAFM = 8.4e28 # Fe (bcc) atomic density [m^-3]
N_L_BOND = 7.0e28 # intermediate, used for the bonding region
# Detrapping pre-exponential--assumed for now
P0_TRAP = 1.0e13 # [s^-1]
# Trap parameters from the paper
# Table 4 -- pure W
TRAPS_W = [
{"E_p": 1.5, "n_t": 9.0e22},
{"E_p": 1.7, "n_t": 1.5e23},
{"E_p": 2.9, "n_t": 1.0e23},
]
# Table 2 -- RAFM steel
TRAPS_RAFM = [
{"E_p": 1.4, "n_t": 6.0e22},
{"E_p": 2.2, "n_t": 8.0e21},
{"E_p": 2.8, "n_t": 3.6e22},
]
# Table 5 -- bonding regions (one entry per joint type)
TRAPS_BOND_WRAFM = [
{"E_p": 1.0, "n_t": 1.0e24},
{"E_p": 1.4, "n_t": 5.0e23},
{"E_p": 2.9, "n_t": 3.0e22},
]
TRAPS_BOND_WNiRAFM = [
{"E_p": 1.0, "n_t": 5.5e24},
{"E_p": 1.5, "n_t": 1.0e23},
{"E_p": 2.4, "n_t": 5.0e24},
]
TRAPS_BOND_WCuRAFM = [
{"E_p": 0.9, "n_t": 2.0e24},
{"E_p": 1.5, "n_t": 1.0e23},
{"E_p": 2.4, "n_t": 5.0e24},
]
# Material definitions (diffusivity + Sievert solubility).
#
# Diffusivities are the quasi-intrinsic values from the paper (eqs. 24-27).
# In festim 2.0, HydrogenTransportProblemDiscontinuous needs K_S to evaluate
# the chemical-potential continuity condition at the interfaces (paper eq. 16):
# c_m^- / K_S^- = c_m^+ / K_S^+
#
# Eq. 29 of the paper uses the series-resistor formula L/D = sum(L_i / D_i),
# which assumes concentration continuity at the interfaces. If my understanding
# is correct, this is equivalent to using a common K_S for all three layers.
#
# We choose an arbitrary value for uniform K_s. It shouldn't affect the result
# but FESTIM needs a concrete number.
K_S_0_COMMON = 1.0e24 # at/m^3/Pa^0.5 (arbitrary, uniform across layers)
E_K_S_COMMON = 0.0 # eV (arbitrary, uniform across layers)
_sol = dict(K_S_0=K_S_0_COMMON, E_K_S=E_K_S_COMMON, solubility_law="sievert")
# Diffusivity parameters (taken from paper)
W_MAT_PARAMS = dict(D_0=2.40e-7, E_D=0.20, **_sol)
RAFM_MAT_PARAMS = dict(D_0=1.95e-7, E_D=0.12, **_sol)
# Bonding regions (eqs. 25-27 for D; same K_S as everything else)
BOND_PARAMS_WRAFM = dict(D_0=2.40e-7, E_D=0.20, **_sol)
BOND_PARAMS_WNi = dict(D_0=3.80e-8, E_D=0.28, **_sol)
BOND_PARAMS_WCu = dict(D_0=1.45e-9, E_D=0.13, **_sol)
def k0_from_D0(D0, N_L):
"""
Trapping pre-exponential from the standard McNabb-Foster form:
k = D0 / (lambda^2 * n_surf * N_L) * exp(-E_D / kT)
"""
return D0 / ((LAMBDA**2) * N_SURF * N_L)
# Integrate a dolfinx Function over its (sub)mesh to get total retention
def _integrate_solution(dolfinx_fn):
mesh = dolfinx_fn.function_space.mesh
form = dolfinx.fem.form(dolfinx_fn * ufl.dx)
local = dolfinx.fem.assemble_scalar(form)
return mesh.comm.allreduce(local, op=MPI.SUM)
# Core simulation routine
def run_case(case_name, bond_params=None, bond_traps=None, N_L_bond=N_L_BOND):
"""
Build and solve a FW diffusion-trapping model for one configuration.
Parameters:
case_name : str
Label used in plot / printouts.
bond_params : dict or None
D_0/E_D for the middle (bonding) layer. If None, the middle layer
gets the W properties and NO traps.
bond_traps : list of dicts or None
Trap parameters for the bonding region. Ignored if bond_params is None.
N_L_bond : float
Lattice site density used to compute k_0 for the bonding region.
Returns:
dict with per-layer retention (atoms / m^2 of FW surface).
"""
print(f"\n{'-'*15}\nRunning case: {case_name}\n{'-'*15}")
# Use a quasi-1D two-dimensional simulation of the layer geometry
mesh_size = L_BOND / 10.0 # resolves bond layer without too much overhead.
mesh_data, tags = create_3layer_mesh(
width_left=L_W,
width_middle=L_BOND,
width_right=L_RAFM,
height=HEIGHT,
mesh_size=mesh_size,
)
festim_mesh = F.Mesh(mesh_data.mesh)
# Model + mesh tags (these come directly from gmsh via create_3layer_mesh)
model = F.HydrogenTransportProblemDiscontinuous()
model.mesh = festim_mesh
model.facet_meshtags = mesh_data.facet_tags
model.volume_meshtags = mesh_data.cell_tags
# Materials & volume/surface subdomains
tungsten = F.Material(name="tungsten", **W_MAT_PARAMS)
# For the "no bonding region" baseline we give the middle layer W
# properties (same as its left neighbour) and no traps, so it's physically
# invisible but still lets us reuse the same 3-layer mesh for all cases.
if bond_params is None:
bond_mat_params = W_MAT_PARAMS
bond_traps_effective = []
else:
bond_mat_params = bond_params
bond_traps_effective = bond_traps or []
bond_mat = F.Material(name="bonding_region", **bond_mat_params)
rafm = F.Material(name="rafm_steel", **RAFM_MAT_PARAMS)
vol_W = F.VolumeSubdomain(id=tags.left_volume, material=tungsten)
vol_bond = F.VolumeSubdomain(id=tags.middle_volume, material=bond_mat)
vol_RAFM = F.VolumeSubdomain(id=tags.right_volume, material=rafm)
# Plasma-facing side = left edge of the W layer.
# Coolant / vacuum side = right edge of the RAFM layer.
surf_left = F.SurfaceSubdomain(id=tags.left_surface) # W surface
surf_right = F.SurfaceSubdomain(id=tags.right_surface) # RAFM surface
model.subdomains = [vol_W, vol_bond, vol_RAFM, surf_left, surf_right]
# Interfaces between the three layers (chemical potential continuity).
# A large penalty keeps the jump condition well-enforced; see here:
# https://festim.discourse.group/t/interface-condition-issue-in-multi-material-diffusion-with-mixed-solubility-laws/158
model.interfaces = [
F.Interface(
id=tags.interface_left,
subdomains=(vol_W, vol_bond),
penalty_term=1e20,
),
F.Interface(
id=tags.interface_right,
subdomains=(vol_bond, vol_RAFM),
penalty_term=1e20,
),
]
model.surface_to_volume = {surf_left: vol_W, surf_right: vol_RAFM}
# Species: 1 mobile T + one trapped-T immobile species per (layer, trap)
H = F.Species("T")
H.subdomains = [vol_W, vol_bond, vol_RAFM]
all_species = [H]
all_reactions = []
# Per-layer bookkeeping of trapped species for the retention integrals.
layer_to_trapped = {vol_W: [], vol_bond: [], vol_RAFM: []}
def add_traps(volume, traps_list, D0_mat, ED_mat, N_L_mat, tag):
"""Register trap reactions for one layer."""
k0 = k0_from_D0(D0_mat, N_L_mat)
for i, trap in enumerate(traps_list):
trapped = F.Species(f"trap_{tag}_{i}", mobile=False)
trapped.subdomains = [volume]
empty = F.ImplicitSpecies(n=trap["n_t"], others=[trapped])
all_species.append(trapped)
layer_to_trapped[volume].append(trapped)
all_reactions.append(
F.Reaction(
reactant=[H, empty],
product=[trapped],
k_0=k0,
E_k=ED_mat,
p_0=P0_TRAP,
E_p=trap["E_p"],
volume=volume,
)
)
add_traps(vol_W, TRAPS_W, W_MAT_PARAMS["D_0"], W_MAT_PARAMS["E_D"], N_L_W, "W")
add_traps(vol_bond, bond_traps_effective, bond_mat_params["D_0"], bond_mat_params["E_D"], N_L_bond, "bond")
add_traps(vol_RAFM, TRAPS_RAFM, RAFM_MAT_PARAMS["D_0"], RAFM_MAT_PARAMS["E_D"], N_L_RAFM, "RAFM")
model.species = all_species
model.reactions = all_reactions
# Boundary conditions & temperature
model.temperature = T_sim
model.boundary_conditions = [
F.FixedConcentrationBC(subdomain=surf_left, species=H, value=c_surface_W),
F.FixedConcentrationBC(subdomain=surf_right, species=H, value=0.0),
]
# Solver settings --> steady state. At steady state the trapping reactions
# equilibrate to the Oriani relation (trap occupancy determined locally by c_m)
model.settings = F.Settings(
atol=1e-10,
rtol=1e-10,
transient=False,
max_iterations=1000,
)
# Run
model.initialise()
model.run()
# Post-process: per-layer retention = integral of (c_m + sum c_trapped)/area
results = {}
for volume, label in [(vol_W, "W"), (vol_bond, "bond"), (vol_RAFM, "RAFM")]:
mobile_sol = H.subdomain_to_post_processing_solution[volume]
mob_int = _integrate_solution(mobile_sol)
trap_int = 0.0
for tsp in layer_to_trapped[volume]:
tsp_sol = tsp.subdomain_to_post_processing_solution[volume]
trap_int += _integrate_solution(tsp_sol)
# Divide by slab height -> atoms / m^2 of FW surface
results[label] = (mob_int + trap_int) / HEIGHT
print(f" {label:4s}: mobile = {mob_int/HEIGHT:.3e},"
f" trapped = {trap_int/HEIGHT:.3e} -> "
f"total = {results[label]:.3e} atoms/m^2")
results["total"] = sum(results[l] for l in ("W", "bond", "RAFM"))
return results
# Run the four FW configurations of figure 8(b)
if __name__ == "__main__":
cases = {}
cases["W-RAFM"] = run_case("W-RAFM (no bonding region)", bond_params=None)
cases["W-bond-RAFM"] = run_case("W-bonding-RAFM",
bond_params=BOND_PARAMS_WRAFM,
bond_traps=TRAPS_BOND_WRAFM)
cases["W-Ni-RAFM"] = run_case("W-Ni-RAFM",
bond_params=BOND_PARAMS_WNi,
bond_traps=TRAPS_BOND_WNiRAFM)
cases["W-Cu-RAFM"] = run_case("W-Cu-RAFM",
bond_params=BOND_PARAMS_WCu,
bond_traps=TRAPS_BOND_WCuRAFM)
# Plot figure 8(b) --> stacked bar w/ W (bottom) + bonding (middle) + RAFM (top)
labels = list(cases.keys())
W_vals = np.array([cases[c]["W"] for c in labels]) / 1e20
bond_vals = np.array([cases[c]["bond"] for c in labels]) / 1e20
rafm_vals = np.array([cases[c]["RAFM"] for c in labels]) / 1e20
fig, ax = plt.subplots(figsize=(7, 5))
x = np.arange(len(labels))
width = 0.55
ax.bar(x, W_vals, width, label="W", color="#d9bc2c")
ax.bar(x, bond_vals, width, bottom = W_vals, label="Bond region", color="#c79a3a")
ax.bar(x, rafm_vals, width, label="RAFM",
bottom = W_vals+ bond_vals, color="#6e4a1a")
ax.set_xticks(x)
ax.set_xticklabels(labels)
ax.set_ylabel(r"Retention (10$^{20}$ T m$^{-2}$)")
ax.set_title("Reproduction of Shen et al. (2026), Fig. 8(b)")
ax.legend(loc="upper left")
ax.grid(axis="y", linestyle=":", alpha=0.5)
totals = W_vals + bond_vals + rafm_vals
for xi, tot in zip(x, totals):
ax.text(xi, tot + 0.2, f"{tot:.1f}", ha="center", va="bottom", fontsize=9)
plt.tight_layout()
plt.show()
# Summary table
print("\n\nSummary (retention in 10^20 atoms/m^2):")
print(f"{'case':15s} {'W':>8s} {'bond':>10s} {'RAFM':>8s} {'total':>8s}")
for c in labels:
r = cases[c]
print(f"{c:15s} {r['W']/1e20:8.2f} {r['bond']/1e20:10.2e} "
f"{r['RAFM']/1e20:8.2f} {r['total']/1e20:8.2f}")
Any thoughts would be greatly appreciated.