Pre-loaded semi-infinite slab#
This verification case from case Val-1c of TMAP7’s V&V report [3] consists of a semi-infinite slab with no traps under a constant concentration \(c_0\) on the first \(10 \ \mathrm{m}\) of the slab. The concentration at the boundaries is set to zero.
FESTIM Code#
Show code cell content
import festim as F
import numpy as np
import sympy as sp
from matplotlib import pyplot as plt
preloaded_length = 10 # m
C_0 = 1 # atom m^-3
D = 1 # 1 m^2 s^-1
model = F.Simulation()
### Mesh Settings ###
vertices = np.concatenate(
[
np.linspace(0, 10, 400),
np.linspace(10, 100, 1000),
]
)
model.mesh = F.MeshFromVertices(vertices)
model.boundary_conditions = [F.DirichletBC(surfaces=[1], value=0, field="solute")]
initial_concentration = sp.Piecewise((C_0, F.x <= preloaded_length), (0, True))
model.initial_conditions = [
F.InitialCondition(field="solute", value=initial_concentration)
]
model.materials = F.Material(id=1, D_0=D, E_D=0)
model.T = F.Temperature(500) # ignored in this problem
model.dt = F.Stepsize(
initial_value=0.01,
stepsize_change_ratio=1.1,
)
test_points = [0.5, preloaded_length, 12] # m
final_times = [100, 100, 50]
profile_times = [0.1] + np.linspace(0, 100, num=10).tolist()[1:]
derived_quantities = F.DerivedQuantities(
[F.PointValue("solute", x=v) for v in test_points]
)
model.exports = [
derived_quantities,
F.TXTExport(
field="solute", filename="./tmap_1c_data/c_profiles.txt", times=profile_times
),
]
model.settings = F.Settings(
absolute_tolerance=1e-10, relative_tolerance=1e-10, final_time=max(final_times)
)
model.initialise()
model.run()
Defining initial values
Defining variational problem
Defining source terms
Defining boundary conditions
Time stepping...
0.0 % 1.0e-02 s Elapsed time so far: 0.0 s
0.0 % 2.1e-02 s Elapsed time so far: 0.0 s
0.0 % 3.3e-02 s Elapsed time so far: 0.0 s
0.1 % 4.6e-02 s Elapsed time so far: 0.0 s
0.1 % 6.1e-02 s Elapsed time so far: 0.0 s
0.1 % 7.7e-02 s Elapsed time so far: 0.0 s
0.1 % 9.5e-02 s Elapsed time so far: 0.1 s
0.1 % 1.0e-01 s Elapsed time so far: 0.1 s
0.1 % 1.1e-01 s Elapsed time so far: 0.1 s
0.1 % 1.1e-01 s Elapsed time so far: 0.1 s
0.1 % 1.2e-01 s Elapsed time so far: 0.1 s
0.1 % 1.3e-01 s Elapsed time so far: 0.1 s
0.1 % 1.3e-01 s Elapsed time so far: 0.1 s
0.1 % 1.4e-01 s Elapsed time so far: 0.1 s
0.1 % 1.5e-01 s Elapsed time so far: 0.1 s
0.2 % 1.6e-01 s Elapsed time so far: 0.1 s
0.2 % 1.8e-01 s Elapsed time so far: 0.1 s
0.2 % 1.9e-01 s Elapsed time so far: 0.1 s
0.2 % 2.0e-01 s Elapsed time so far: 0.1 s
0.2 % 2.2e-01 s Elapsed time so far: 0.1 s
0.2 % 2.4e-01 s Elapsed time so far: 0.2 s
0.3 % 2.6e-01 s Elapsed time so far: 0.2 s
0.3 % 2.8e-01 s Elapsed time so far: 0.2 s
0.3 % 3.0e-01 s Elapsed time so far: 0.2 s
0.3 % 3.3e-01 s Elapsed time so far: 0.2 s
0.4 % 3.6e-01 s Elapsed time so far: 0.2 s
0.4 % 3.9e-01 s Elapsed time so far: 0.2 s
0.4 % 4.2e-01 s Elapsed time so far: 0.2 s
/home/docs/checkouts/readthedocs.org/user_builds/festim-vv-report/conda/v1.1/lib/python3.11/site-packages/festim/generic_simulation.py:439: UserWarning: To ensure that TXTExport exports data at the desired times TXTExport.times are added to milestones
warnings.warn(msg)
0.5 % 4.6e-01 s Elapsed time so far: 0.2 s
0.5 % 5.0e-01 s Elapsed time so far: 0.2 s
0.6 % 5.5e-01 s Elapsed time so far: 0.2 s
0.6 % 6.0e-01 s Elapsed time so far: 0.2 s
0.7 % 6.5e-01 s Elapsed time so far: 0.2 s
0.7 % 7.2e-01 s Elapsed time so far: 0.2 s
0.8 % 7.8e-01 s Elapsed time so far: 0.2 s
0.9 % 8.6e-01 s Elapsed time so far: 0.2 s
0.9 % 9.4e-01 s Elapsed time so far: 0.2 s
1.0 % 1.0e+00 s Elapsed time so far: 0.3 s
1.1 % 1.1e+00 s Elapsed time so far: 0.3 s
1.2 % 1.2e+00 s Elapsed time so far: 0.3 s
1.4 % 1.4e+00 s Elapsed time so far: 0.3 s
1.5 % 1.5e+00 s Elapsed time so far: 0.3 s
1.6 % 1.6e+00 s Elapsed time so far: 0.3 s
1.8 % 1.8e+00 s Elapsed time so far: 0.3 s
2.0 % 2.0e+00 s Elapsed time so far: 0.3 s
2.1 % 2.2e+00 s Elapsed time so far: 0.3 s
2.4 % 2.4e+00 s Elapsed time so far: 0.3 s
2.6 % 2.6e+00 s Elapsed time so far: 0.3 s
2.9 % 2.9e+00 s Elapsed time so far: 0.3 s
3.1 % 3.1e+00 s Elapsed time so far: 0.3 s
3.4 % 3.4e+00 s Elapsed time so far: 0.3 s
3.8 % 3.8e+00 s Elapsed time so far: 0.4 s
4.2 % 4.2e+00 s Elapsed time so far: 0.4 s
4.6 % 4.6e+00 s Elapsed time so far: 0.4 s
5.0 % 5.0e+00 s Elapsed time so far: 0.4 s
5.5 % 5.5e+00 s Elapsed time so far: 0.4 s
6.1 % 6.1e+00 s Elapsed time so far: 0.4 s
6.7 % 6.7e+00 s Elapsed time so far: 0.4 s
7.3 % 7.3e+00 s Elapsed time so far: 0.4 s
8.1 % 8.1e+00 s Elapsed time so far: 0.4 s
8.9 % 8.9e+00 s Elapsed time so far: 0.4 s
9.7 % 9.7e+00 s Elapsed time so far: 0.4 s
10.7 % 1.1e+01 s Elapsed time so far: 0.4 s
11.1 % 1.1e+01 s Elapsed time so far: 0.4 s
11.6 % 1.2e+01 s Elapsed time so far: 0.5 s
12.0 % 1.2e+01 s Elapsed time so far: 0.5 s
12.6 % 1.3e+01 s Elapsed time so far: 0.5 s
13.2 % 1.3e+01 s Elapsed time so far: 0.5 s
13.8 % 1.4e+01 s Elapsed time so far: 0.5 s
14.5 % 1.5e+01 s Elapsed time so far: 0.5 s
15.3 % 1.5e+01 s Elapsed time so far: 0.5 s
16.2 % 1.6e+01 s Elapsed time so far: 0.5 s
17.1 % 1.7e+01 s Elapsed time so far: 0.5 s
18.2 % 1.8e+01 s Elapsed time so far: 0.5 s
19.3 % 1.9e+01 s Elapsed time so far: 0.5 s
20.6 % 2.1e+01 s Elapsed time so far: 0.5 s
22.0 % 2.2e+01 s Elapsed time so far: 0.5 s
22.2 % 2.2e+01 s Elapsed time so far: 0.5 s
22.5 % 2.3e+01 s Elapsed time so far: 0.6 s
22.8 % 2.3e+01 s Elapsed time so far: 0.6 s
23.1 % 2.3e+01 s Elapsed time so far: 0.6 s
23.5 % 2.4e+01 s Elapsed time so far: 0.6 s
23.9 % 2.4e+01 s Elapsed time so far: 0.6 s
24.4 % 2.4e+01 s Elapsed time so far: 0.6 s
24.9 % 2.5e+01 s Elapsed time so far: 0.6 s
25.4 % 2.5e+01 s Elapsed time so far: 0.6 s
26.0 % 2.6e+01 s Elapsed time so far: 0.6 s
26.7 % 2.7e+01 s Elapsed time so far: 0.6 s
27.4 % 2.7e+01 s Elapsed time so far: 0.6 s
28.2 % 2.8e+01 s Elapsed time so far: 0.6 s
29.1 % 2.9e+01 s Elapsed time so far: 0.6 s
30.1 % 3.0e+01 s Elapsed time so far: 0.6 s
31.1 % 3.1e+01 s Elapsed time so far: 0.6 s
32.3 % 3.2e+01 s Elapsed time so far: 0.7 s
33.3 % 3.3e+01 s Elapsed time so far: 0.7 s
34.5 % 3.4e+01 s Elapsed time so far: 0.7 s
35.8 % 3.6e+01 s Elapsed time so far: 0.7 s
37.1 % 3.7e+01 s Elapsed time so far: 0.7 s
38.7 % 3.9e+01 s Elapsed time so far: 0.7 s
40.4 % 4.0e+01 s Elapsed time so far: 0.7 s
42.2 % 4.2e+01 s Elapsed time so far: 0.7 s
44.3 % 4.4e+01 s Elapsed time so far: 0.7 s
44.4 % 4.4e+01 s Elapsed time so far: 0.7 s
44.6 % 4.5e+01 s Elapsed time so far: 0.8 s
44.8 % 4.5e+01 s Elapsed time so far: 0.8 s
45.1 % 4.5e+01 s Elapsed time so far: 0.8 s
45.3 % 4.5e+01 s Elapsed time so far: 0.8 s
45.6 % 4.6e+01 s Elapsed time so far: 0.8 s
45.9 % 4.6e+01 s Elapsed time so far: 0.8 s
46.2 % 4.6e+01 s Elapsed time so far: 0.8 s
46.6 % 4.7e+01 s Elapsed time so far: 0.8 s
47.0 % 4.7e+01 s Elapsed time so far: 0.8 s
47.4 % 4.7e+01 s Elapsed time so far: 0.8 s
47.9 % 4.8e+01 s Elapsed time so far: 0.8 s
48.4 % 4.8e+01 s Elapsed time so far: 0.8 s
49.0 % 4.9e+01 s Elapsed time so far: 0.8 s
49.6 % 5.0e+01 s Elapsed time so far: 0.8 s
50.3 % 5.0e+01 s Elapsed time so far: 0.8 s
51.1 % 5.1e+01 s Elapsed time so far: 0.8 s
52.0 % 5.2e+01 s Elapsed time so far: 0.9 s
52.9 % 5.3e+01 s Elapsed time so far: 0.9 s
53.9 % 5.4e+01 s Elapsed time so far: 0.9 s
55.0 % 5.5e+01 s Elapsed time so far: 0.9 s
55.6 % 5.6e+01 s Elapsed time so far: 0.9 s
56.1 % 5.6e+01 s Elapsed time so far: 0.9 s
56.7 % 5.7e+01 s Elapsed time so far: 0.9 s
57.4 % 5.7e+01 s Elapsed time so far: 0.9 s
58.1 % 5.8e+01 s Elapsed time so far: 0.9 s
58.9 % 5.9e+01 s Elapsed time so far: 0.9 s
59.8 % 6.0e+01 s Elapsed time so far: 0.9 s
60.8 % 6.1e+01 s Elapsed time so far: 0.9 s
61.9 % 6.2e+01 s Elapsed time so far: 1.0 s
63.0 % 6.3e+01 s Elapsed time so far: 1.0 s
64.3 % 6.4e+01 s Elapsed time so far: 1.0 s
65.8 % 6.6e+01 s Elapsed time so far: 1.0 s
66.7 % 6.7e+01 s Elapsed time so far: 1.0 s
67.7 % 6.8e+01 s Elapsed time so far: 1.0 s
68.8 % 6.9e+01 s Elapsed time so far: 1.0 s
70.0 % 7.0e+01 s Elapsed time so far: 1.0 s
71.3 % 7.1e+01 s Elapsed time so far: 1.0 s
72.8 % 7.3e+01 s Elapsed time so far: 1.0 s
74.4 % 7.4e+01 s Elapsed time so far: 1.0 s
76.1 % 7.6e+01 s Elapsed time so far: 1.0 s
77.8 % 7.8e+01 s Elapsed time so far: 1.1 s
79.6 % 8.0e+01 s Elapsed time so far: 1.1 s
81.6 % 8.2e+01 s Elapsed time so far: 1.1 s
83.8 % 8.4e+01 s Elapsed time so far: 1.1 s
86.2 % 8.6e+01 s Elapsed time so far: 1.1 s
88.8 % 8.9e+01 s Elapsed time so far: 1.1 s
88.9 % 8.9e+01 s Elapsed time so far: 1.1 s
88.9 % 8.9e+01 s Elapsed time so far: 1.1 s
89.0 % 8.9e+01 s Elapsed time so far: 1.1 s
89.0 % 8.9e+01 s Elapsed time so far: 1.2 s
89.1 % 8.9e+01 s Elapsed time so far: 1.2 s
89.2 % 8.9e+01 s Elapsed time so far: 1.2 s
89.2 % 8.9e+01 s Elapsed time so far: 1.2 s
89.3 % 8.9e+01 s Elapsed time so far: 1.2 s
89.4 % 8.9e+01 s Elapsed time so far: 1.2 s
89.5 % 9.0e+01 s Elapsed time so far: 1.2 s
89.6 % 9.0e+01 s Elapsed time so far: 1.2 s
89.8 % 9.0e+01 s Elapsed time so far: 1.2 s
89.9 % 9.0e+01 s Elapsed time so far: 1.2 s
90.0 % 9.0e+01 s Elapsed time so far: 1.2 s
90.2 % 9.0e+01 s Elapsed time so far: 1.2 s
90.4 % 9.0e+01 s Elapsed time so far: 1.2 s
90.6 % 9.1e+01 s Elapsed time so far: 1.2 s
90.8 % 9.1e+01 s Elapsed time so far: 1.2 s
91.0 % 9.1e+01 s Elapsed time so far: 1.2 s
91.3 % 9.1e+01 s Elapsed time so far: 1.2 s
91.5 % 9.2e+01 s Elapsed time so far: 1.3 s
91.9 % 9.2e+01 s Elapsed time so far: 1.3 s
92.2 % 9.2e+01 s Elapsed time so far: 1.3 s
92.6 % 9.3e+01 s Elapsed time so far: 1.3 s
93.0 % 9.3e+01 s Elapsed time so far: 1.3 s
93.5 % 9.3e+01 s Elapsed time so far: 1.3 s
94.0 % 9.4e+01 s Elapsed time so far: 1.3 s
94.5 % 9.5e+01 s Elapsed time so far: 1.3 s
95.1 % 9.5e+01 s Elapsed time so far: 1.3 s
95.8 % 9.6e+01 s Elapsed time so far: 1.3 s
96.5 % 9.7e+01 s Elapsed time so far: 1.3 s
97.3 % 9.7e+01 s Elapsed time so far: 1.3 s
98.2 % 9.8e+01 s Elapsed time so far: 1.3 s
99.2 % 9.9e+01 s Elapsed time so far: 1.3 s
100.0 % 1.0e+02 s Elapsed time so far: 1.4 s
Comparison with exact solution#
This is a comparison of the computed concentration profiles at different times with the exact analytical solution (shown in dashed lines).
The analytical solution is given by
\[
c(x, t) = \frac{c_0}{2}\left[
2\mathrm{erf}\left(\frac{x}{2 \sqrt{Dt}}\right)
- \mathrm{erf}\left(\frac{x - h}{2 \sqrt{Dt}}\right)
- \mathrm{erf}\left(\frac{x + h}{2 \sqrt{Dt}}\right)
\right]
\]
where \(h\) is the thickness of the pre-loaded region.
Show code cell source
from matplotlib import cm
from matplotlib.colors import Normalize
from scipy.special import erf
norm = Normalize(vmin=0, vmax=max(profile_times))
cmap = cm.viridis
def exact_solution(x, t):
sqrt_term = np.sqrt(4 * D * t)
return (
C_0
/ 2
* (
2 * erf(x / sqrt_term)
- erf((x - preloaded_length) / sqrt_term)
- erf((x + preloaded_length) / sqrt_term)
)
)
plt.figure()
filename = model.exports[1].filename
data = np.genfromtxt(filename, delimiter=",", names=True)
for i, t in enumerate(profile_times):
label = "exact" if i == 0 else ""
x = data["x"]
y_name = f"t{t:.2e}s".replace("+", "").replace("-", "").replace(".", "")
y = data[y_name]
x, y = zip(*sorted(zip(x, y)))
exact_y = exact_solution(np.array(x), t)
plt.plot(x, exact_y, linestyle="dashed", color="tab:grey", linewidth=3, label=label)
plt.plot(x, y, color=cmap(norm(t)))
plt.xlabel("x")
plt.ylabel("Concentration (atom / m^3)")
# Add colorbar
sm = cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([]) # You can also set this to the range of your data
plt.colorbar(sm, label="Time (s)", ax=plt.gca())
plt.legend()
plt.show()
The results can also be compared with the results obtained by TMAP7.
Show code cell source
fig, axs = plt.subplots(1, len(test_points), sharey=True)
for i, x in enumerate(test_points):
plt.sca(axs[i])
# plotting computed data
computed_solution = derived_quantities[i].data
t = np.array(derived_quantities[i].t)
plt.plot(t, computed_solution, label="FESTIM", linewidth=3)
# plotting exact solution
plt.plot(t, exact_solution(x, t), label="Exact", color="green", linestyle="--")
# plotting TMAP data
tmap_data = np.genfromtxt(
f"./tmap_1c_data/tmap_point_data_{i}.txt", delimiter=" ", names=True
)
tmap_t = tmap_data["t"]
tmap_solution = tmap_data["tmap"]
plt.scatter(tmap_t, tmap_solution, label="TMAP7", color="purple")
plt.title(f"x={x} m")
if i == 0:
plt.ylabel("Concentration (atom / m$^3$)")
plt.xlabel("t (s)")
fig.tight_layout()
plt.legend()
plt.show()
