Effective diffusivity regime#
This verification case consists of a slab of depth \(l = 1 \times 10^{-3} \ \mathrm{m}\) with one trap under the effective diffusivity regime.
A trapping parameter \(\zeta\) is defined by
\[
\zeta = \frac{\lambda ^ 2 \nu}{D_0 \rho} \exp \left(\frac{E_k - E_p}{k_B T}\right) + \frac{c_m}{\rho}
\]
where
\(\lambda \ \mathrm{(m)}\) is the lattice parameter,
\(\nu \ (\mathrm{s}^{-1})\), the Debye frequency
\(\rho\) is the trapping site fraction,
\(c_m (\text{atom} \ \mathrm{m}^{-3})\) is the mobile atom concentration,
and the effective diffusivity \(\mathrm{D_\text{eff}}\) is defined by
\[
D_\text{eff} = \frac{D}{1 + \frac{1}{\zeta}}
\]
Then with a breakthrough time \(\tau = \frac{l^2}{2\pi^2 D_\text{eff}}\), the exact solution for flux is
\[
J = \frac{}{} \left[ 1 + 2\sum_{m=1}^\infty (-1)^m \exp \left( -m^2 \frac{t}{\tau} \right) \right]
\]
This analytical solution was obtained from TMAP7’s V&V report [3], case Val-1da.
FESTIM Code#
Show code cell content
# referenced https://github.com/gabriele-ferrero/Titans_TT_codecomparison/blob/main/Festim_models/WeakTrap.py
import festim as F
import numpy as np
import matplotlib.pyplot as plt
D_0 = 1.9e-7
N_A = 6.0221408e23
rho_w = 6.3382e28
E_k = 0.2
E_p = 1
T = 1000
D = D_0 * np.exp(-E_k / (F.k_B * T))
S = 2.9e-5 * np.exp(-1 / (F.k_B * T))
sample_depth = 1e-3
c_m = (1e5) ** 0.5 * S * 1.0525e5
zeta = (6 * rho_w * np.exp((E_k - E_p) / (F.k_B * T)) + c_m * N_A) / (rho_w * 1e-3)
D_eff = D / (1 + 1 / zeta)
model = F.Simulation()
model.mesh = F.MeshFromVertices(vertices=np.linspace(0, sample_depth, num=100))
model.materials = F.Material(id=1, D_0=D_0, E_D=E_k)
model.T = F.Temperature(value=T)
model.boundary_conditions = [
F.DirichletBC(surfaces=1, value=c_m * N_A, field=0),
F.DirichletBC(surfaces=2, value=0, field=0),
]
trap = F.Trap(
k_0=1.58e7 / N_A,
E_k=E_k,
p_0=1e13,
E_p=E_p,
density=1e-3 * rho_w,
materials=model.materials[0],
)
model.traps = [trap]
model.settings = F.Settings(
absolute_tolerance=1e10, relative_tolerance=1e-10, final_time=100 # s
)
model.dt = F.Stepsize(0.5)
derived_quantities = F.DerivedQuantities([F.HydrogenFlux(surface=2)])
model.exports = [derived_quantities]
model.initialise()
model.run()
Defining initial values
Defining variational problem
Defining source terms
Defining boundary conditions
Time stepping...
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
0.5 % 5.0e-01 s Elapsed time so far: 1.2 s
1.0 % 1.0e+00 s Elapsed time so far: 1.2 s
1.5 % 1.5e+00 s Elapsed time so far: 1.2 s
2.0 % 2.0e+00 s Elapsed time so far: 1.2 s
2.5 % 2.5e+00 s Elapsed time so far: 1.3 s
3.0 % 3.0e+00 s Elapsed time so far: 1.3 s
3.5 % 3.5e+00 s Elapsed time so far: 1.3 s
4.0 % 4.0e+00 s Elapsed time so far: 1.3 s
4.5 % 4.5e+00 s Elapsed time so far: 1.3 s
5.0 % 5.0e+00 s Elapsed time so far: 1.3 s
5.5 % 5.5e+00 s Elapsed time so far: 1.3 s
6.0 % 6.0e+00 s Elapsed time so far: 1.3 s
6.5 % 6.5e+00 s Elapsed time so far: 1.3 s
7.0 % 7.0e+00 s Elapsed time so far: 1.3 s
7.5 % 7.5e+00 s Elapsed time so far: 1.3 s
8.0 % 8.0e+00 s Elapsed time so far: 1.3 s
8.5 % 8.5e+00 s Elapsed time so far: 1.3 s
9.0 % 9.0e+00 s Elapsed time so far: 1.3 s
9.5 % 9.5e+00 s Elapsed time so far: 1.3 s
10.0 % 1.0e+01 s Elapsed time so far: 1.3 s
10.5 % 1.0e+01 s Elapsed time so far: 1.3 s
11.0 % 1.1e+01 s Elapsed time so far: 1.4 s
11.5 % 1.2e+01 s Elapsed time so far: 1.4 s
12.0 % 1.2e+01 s Elapsed time so far: 1.4 s
12.5 % 1.2e+01 s Elapsed time so far: 1.4 s
13.0 % 1.3e+01 s Elapsed time so far: 1.4 s
13.5 % 1.4e+01 s Elapsed time so far: 1.4 s
14.0 % 1.4e+01 s Elapsed time so far: 1.4 s
14.5 % 1.4e+01 s Elapsed time so far: 1.4 s
15.0 % 1.5e+01 s Elapsed time so far: 1.4 s
15.5 % 1.6e+01 s Elapsed time so far: 1.4 s
16.0 % 1.6e+01 s Elapsed time so far: 1.4 s
16.5 % 1.6e+01 s Elapsed time so far: 1.4 s
17.0 % 1.7e+01 s Elapsed time so far: 1.4 s
17.5 % 1.8e+01 s Elapsed time so far: 1.4 s
18.0 % 1.8e+01 s Elapsed time so far: 1.4 s
18.5 % 1.8e+01 s Elapsed time so far: 1.4 s
19.0 % 1.9e+01 s Elapsed time so far: 1.4 s
19.5 % 2.0e+01 s Elapsed time so far: 1.4 s
20.0 % 2.0e+01 s Elapsed time so far: 1.4 s
/home/docs/checkouts/readthedocs.org/user_builds/festim-vv-report/conda/v1.0/lib/python3.11/site-packages/festim/exports/derived_quantities/derived_quantities.py:129: DeprecationWarning: The current derived_quantities title style will be deprecated in a future release, please use show_units=True instead
warnings.warn(
20.5 % 2.0e+01 s Elapsed time so far: 1.5 s
21.0 % 2.1e+01 s Elapsed time so far: 1.5 s
21.5 % 2.2e+01 s Elapsed time so far: 1.5 s
22.0 % 2.2e+01 s Elapsed time so far: 1.5 s
22.5 % 2.2e+01 s Elapsed time so far: 1.5 s
23.0 % 2.3e+01 s Elapsed time so far: 1.5 s
23.5 % 2.4e+01 s Elapsed time so far: 1.5 s
24.0 % 2.4e+01 s Elapsed time so far: 1.5 s
24.5 % 2.4e+01 s Elapsed time so far: 1.5 s
25.0 % 2.5e+01 s Elapsed time so far: 1.5 s
25.5 % 2.6e+01 s Elapsed time so far: 1.5 s
26.0 % 2.6e+01 s Elapsed time so far: 1.5 s
26.5 % 2.6e+01 s Elapsed time so far: 1.5 s
27.0 % 2.7e+01 s Elapsed time so far: 1.5 s
27.5 % 2.8e+01 s Elapsed time so far: 1.5 s
28.0 % 2.8e+01 s Elapsed time so far: 1.5 s
28.5 % 2.8e+01 s Elapsed time so far: 1.5 s
29.0 % 2.9e+01 s Elapsed time so far: 1.5 s
29.5 % 3.0e+01 s Elapsed time so far: 1.5 s
30.0 % 3.0e+01 s Elapsed time so far: 1.6 s
30.5 % 3.0e+01 s Elapsed time so far: 1.6 s
31.0 % 3.1e+01 s Elapsed time so far: 1.6 s
31.5 % 3.2e+01 s Elapsed time so far: 1.6 s
32.0 % 3.2e+01 s Elapsed time so far: 1.6 s
32.5 % 3.2e+01 s Elapsed time so far: 1.6 s
33.0 % 3.3e+01 s Elapsed time so far: 1.6 s
33.5 % 3.4e+01 s Elapsed time so far: 1.6 s
34.0 % 3.4e+01 s Elapsed time so far: 1.6 s
34.5 % 3.4e+01 s Elapsed time so far: 1.6 s
35.0 % 3.5e+01 s Elapsed time so far: 1.6 s
35.5 % 3.6e+01 s Elapsed time so far: 1.6 s
36.0 % 3.6e+01 s Elapsed time so far: 1.6 s
36.5 % 3.6e+01 s Elapsed time so far: 1.6 s
37.0 % 3.7e+01 s Elapsed time so far: 1.6 s
37.5 % 3.8e+01 s Elapsed time so far: 1.6 s
38.0 % 3.8e+01 s Elapsed time so far: 1.6 s
38.5 % 3.8e+01 s Elapsed time so far: 1.6 s
39.0 % 3.9e+01 s Elapsed time so far: 1.6 s
39.5 % 4.0e+01 s Elapsed time so far: 1.6 s
40.0 % 4.0e+01 s Elapsed time so far: 1.7 s
40.5 % 4.0e+01 s Elapsed time so far: 1.7 s
41.0 % 4.1e+01 s Elapsed time so far: 1.7 s
41.5 % 4.2e+01 s Elapsed time so far: 1.7 s
42.0 % 4.2e+01 s Elapsed time so far: 1.7 s
42.5 % 4.2e+01 s Elapsed time so far: 1.7 s
43.0 % 4.3e+01 s Elapsed time so far: 1.7 s
43.5 % 4.4e+01 s Elapsed time so far: 1.7 s
44.0 % 4.4e+01 s Elapsed time so far: 1.7 s
44.5 % 4.4e+01 s Elapsed time so far: 1.7 s
45.0 % 4.5e+01 s Elapsed time so far: 1.7 s
45.5 % 4.6e+01 s Elapsed time so far: 1.7 s
46.0 % 4.6e+01 s Elapsed time so far: 1.7 s
46.5 % 4.6e+01 s Elapsed time so far: 1.7 s
47.0 % 4.7e+01 s Elapsed time so far: 1.7 s
47.5 % 4.8e+01 s Elapsed time so far: 1.7 s
48.0 % 4.8e+01 s Elapsed time so far: 1.7 s
48.5 % 4.8e+01 s Elapsed time so far: 1.7 s
49.0 % 4.9e+01 s Elapsed time so far: 1.7 s
49.5 % 5.0e+01 s Elapsed time so far: 1.8 s
50.0 % 5.0e+01 s Elapsed time so far: 1.8 s
50.5 % 5.0e+01 s Elapsed time so far: 1.8 s
51.0 % 5.1e+01 s Elapsed time so far: 1.8 s
51.5 % 5.2e+01 s Elapsed time so far: 1.8 s
52.0 % 5.2e+01 s Elapsed time so far: 1.8 s
52.5 % 5.2e+01 s Elapsed time so far: 1.8 s
53.0 % 5.3e+01 s Elapsed time so far: 1.8 s
53.5 % 5.4e+01 s Elapsed time so far: 1.8 s
54.0 % 5.4e+01 s Elapsed time so far: 1.8 s
54.5 % 5.4e+01 s Elapsed time so far: 1.8 s
55.0 % 5.5e+01 s Elapsed time so far: 1.8 s
55.5 % 5.6e+01 s Elapsed time so far: 1.8 s
56.0 % 5.6e+01 s Elapsed time so far: 1.8 s
56.5 % 5.6e+01 s Elapsed time so far: 1.8 s
57.0 % 5.7e+01 s Elapsed time so far: 1.8 s
57.5 % 5.8e+01 s Elapsed time so far: 1.8 s
58.0 % 5.8e+01 s Elapsed time so far: 1.8 s
58.5 % 5.8e+01 s Elapsed time so far: 1.8 s
59.0 % 5.9e+01 s Elapsed time so far: 1.9 s
59.5 % 6.0e+01 s Elapsed time so far: 1.9 s
60.0 % 6.0e+01 s Elapsed time so far: 1.9 s
60.5 % 6.0e+01 s Elapsed time so far: 1.9 s
61.0 % 6.1e+01 s Elapsed time so far: 1.9 s
61.5 % 6.2e+01 s Elapsed time so far: 1.9 s
62.0 % 6.2e+01 s Elapsed time so far: 1.9 s
62.5 % 6.2e+01 s Elapsed time so far: 1.9 s
63.0 % 6.3e+01 s Elapsed time so far: 1.9 s
63.5 % 6.4e+01 s Elapsed time so far: 1.9 s
64.0 % 6.4e+01 s Elapsed time so far: 1.9 s
64.5 % 6.4e+01 s Elapsed time so far: 1.9 s
65.0 % 6.5e+01 s Elapsed time so far: 1.9 s
65.5 % 6.6e+01 s Elapsed time so far: 1.9 s
66.0 % 6.6e+01 s Elapsed time so far: 1.9 s
66.5 % 6.6e+01 s Elapsed time so far: 1.9 s
67.0 % 6.7e+01 s Elapsed time so far: 1.9 s
67.5 % 6.8e+01 s Elapsed time so far: 1.9 s
68.0 % 6.8e+01 s Elapsed time so far: 1.9 s
68.5 % 6.8e+01 s Elapsed time so far: 2.0 s
69.0 % 6.9e+01 s Elapsed time so far: 2.0 s
69.5 % 7.0e+01 s Elapsed time so far: 2.0 s
70.0 % 7.0e+01 s Elapsed time so far: 2.0 s
70.5 % 7.0e+01 s Elapsed time so far: 2.0 s
71.0 % 7.1e+01 s Elapsed time so far: 2.0 s
71.5 % 7.2e+01 s Elapsed time so far: 2.0 s
72.0 % 7.2e+01 s Elapsed time so far: 2.0 s
72.5 % 7.2e+01 s Elapsed time so far: 2.0 s
73.0 % 7.3e+01 s Elapsed time so far: 2.0 s
73.5 % 7.4e+01 s Elapsed time so far: 2.0 s
74.0 % 7.4e+01 s Elapsed time so far: 2.0 s
74.5 % 7.4e+01 s Elapsed time so far: 2.0 s
75.0 % 7.5e+01 s Elapsed time so far: 2.0 s
75.5 % 7.6e+01 s Elapsed time so far: 2.0 s
76.0 % 7.6e+01 s Elapsed time so far: 2.0 s
76.5 % 7.6e+01 s Elapsed time so far: 2.0 s
77.0 % 7.7e+01 s Elapsed time so far: 2.0 s
77.5 % 7.8e+01 s Elapsed time so far: 2.0 s
78.0 % 7.8e+01 s Elapsed time so far: 2.0 s
78.5 % 7.8e+01 s Elapsed time so far: 2.1 s
79.0 % 7.9e+01 s Elapsed time so far: 2.1 s
79.5 % 8.0e+01 s Elapsed time so far: 2.1 s
80.0 % 8.0e+01 s Elapsed time so far: 2.1 s
80.5 % 8.0e+01 s Elapsed time so far: 2.1 s
81.0 % 8.1e+01 s Elapsed time so far: 2.1 s
81.5 % 8.2e+01 s Elapsed time so far: 2.1 s
82.0 % 8.2e+01 s Elapsed time so far: 2.1 s
82.5 % 8.2e+01 s Elapsed time so far: 2.1 s
83.0 % 8.3e+01 s Elapsed time so far: 2.1 s
83.5 % 8.4e+01 s Elapsed time so far: 2.1 s
84.0 % 8.4e+01 s Elapsed time so far: 2.1 s
84.5 % 8.4e+01 s Elapsed time so far: 2.1 s
85.0 % 8.5e+01 s Elapsed time so far: 2.1 s
85.5 % 8.6e+01 s Elapsed time so far: 2.1 s
86.0 % 8.6e+01 s Elapsed time so far: 2.1 s
86.5 % 8.6e+01 s Elapsed time so far: 2.1 s
87.0 % 8.7e+01 s Elapsed time so far: 2.1 s
87.5 % 8.8e+01 s Elapsed time so far: 2.1 s
88.0 % 8.8e+01 s Elapsed time so far: 2.2 s
88.5 % 8.8e+01 s Elapsed time so far: 2.2 s
89.0 % 8.9e+01 s Elapsed time so far: 2.2 s
89.5 % 9.0e+01 s Elapsed time so far: 2.2 s
90.0 % 9.0e+01 s Elapsed time so far: 2.2 s
90.5 % 9.0e+01 s Elapsed time so far: 2.2 s
91.0 % 9.1e+01 s Elapsed time so far: 2.2 s
91.5 % 9.2e+01 s Elapsed time so far: 2.2 s
92.0 % 9.2e+01 s Elapsed time so far: 2.2 s
92.5 % 9.2e+01 s Elapsed time so far: 2.2 s
93.0 % 9.3e+01 s Elapsed time so far: 2.2 s
93.5 % 9.4e+01 s Elapsed time so far: 2.2 s
94.0 % 9.4e+01 s Elapsed time so far: 2.2 s
94.5 % 9.4e+01 s Elapsed time so far: 2.2 s
95.0 % 9.5e+01 s Elapsed time so far: 2.2 s
95.5 % 9.6e+01 s Elapsed time so far: 2.2 s
96.0 % 9.6e+01 s Elapsed time so far: 2.2 s
96.5 % 9.6e+01 s Elapsed time so far: 2.2 s
97.0 % 9.7e+01 s Elapsed time so far: 2.2 s
97.5 % 9.8e+01 s Elapsed time so far: 2.2 s
98.0 % 9.8e+01 s Elapsed time so far: 2.2 s
98.5 % 9.8e+01 s Elapsed time so far: 2.2 s
99.0 % 9.9e+01 s Elapsed time so far: 2.3 s
99.5 % 1.0e+02 s Elapsed time so far: 2.3 s
100.0 % 1.0e+02 s Elapsed time so far: 2.3 s
Comparison with exact solution#
Show code cell source
# plot computed solution
t = np.array(derived_quantities.t)
computed_solution = derived_quantities.filter(surfaces=2).data
plt.plot(t, np.abs(computed_solution) / 2, label="FESTIM", linewidth=3)
# plot exact solution
t = np.array(t)
m = np.arange(1, 10001)
# Calculate the exponential part for all m values at once
tau = sample_depth**2 / (np.pi**2 * D_eff)
exp_part = np.exp(-m ** 2 * t[:, None] / tau)
# Calculate the 'add' part for all m values and sum them up for each t
add = 2 * (-1) ** m * exp_part
exact_solution = 1 + add.sum(axis=1) # Sum along the m dimension and add 1 for the initial ones
exact_solution = N_A * exact_solution * c_m * D / (2 * sample_depth)
plt.plot(t, exact_solution, linestyle="--", color="green", label="exact")
plt.xlabel("Time (s)")
plt.ylabel("Downstream flux (H/m2/s)")
plt.legend()
plt.show()
