I’m trying to model deuterium retention in tungsten damaged by 20 MeV tungsten ions. Since the D retention is saturated when the damage is above 0.1 dpa, I want to use a following damage depth distribution shown in orange (which is modified from the SRIM calculated depth profile shown in blue).
I tried to modify the damage depth distribution f(x) = \frac{1}{1 + \exp(\frac{x - x_{0}}{\Delta x}) } in Deuterium retention in self damaged tungsten such that x_{0} and \Delta x are now 1.8e-6 and 9e-8, respectively, leading to a sharper drop in damage depth distribution compared to the example. The solver then has trouble converging (i.e., the error message reads “ValueError: stepsize reached minimal value”), even when reducing the minimum timestep from 1e-1 to 1e-5, the relative tolerance from 1e-10 to 1e-5, and mesh near the drop off from ~80 nm to ~0.3 nm. I also tried to replace the damage depth distribution to use a hyperbolic tangent (i.e., \frac{1}{2}(tanh(- \frac{x - x_{0}}{\Delta x}) + 1) with x_{0} and \Delta x set to 1.8e-7 and 1.7e-7, respectively) instead of exponential, but I still can’t get the solver to converge. Any suggestions for how to resolve the convergence issue?
Note that I am using Festim version 1. Below is a copy of my code.
import festim as F
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.use('tkagg')
import os, shutil
file_name = "tds/derived_quantities_example_TDSfromdamagedW.csv"
data_folder_path = '/home/pisces/examples/damagedesorb/tds/'
if os.path.isdir(data_folder_path):
shutil.rmtree(data_folder_path)
os.mkdir(data_folder_path)
# # ### Parameters ###
D_0 = 1.6e-7 # m^2 s^-1
E_D = 0.28 # eV
w_atom_density = 6.3222e28
sample_thickness = 0.8e-3 # m
sample_area = 12e-03 * 15e-03
detrapping_energies = [1.15, 1.35, 1.65, 1.85, 2.05]
num_traps = len(detrapping_energies) + 1 # including intrinsic trap
dpa_n_i = {
0.001: [1e24, 2.5e24, 1e24, 1e24, 2e23],
0.1: [4.8e25, 3.8e25, 2.6e25, 3.6e25, 1.1e25],
}
# Table 2 from Dark et al 10.1088/1741-4326/ad56a0
T_imp = 370 # K
T_rest = 295 # K
R_p = 0.7e-9 # m
sigma = 0.5e-9 # m
t_imp = 72 * 3600 # s
implantation_time = t_imp
t_rest = 60#12 * 3600 # s
resting_time = t_rest
Beta = 3 / 60 # K s^-1
fluence = 1.5e25
flux = fluence / t_imp
start_tds = t_imp + t_rest # s
min_temp, max_temp = 300, 1000
export_times = [implantation_time]
##vertices = sample_thickness/100*np.genfromtxt("mesh.dat",delimiter=",")
def festim_sim(densities):
model = F.Simulation()
vertices = np.concatenate(
[
np.linspace(0, 3e-9, num=100),
np.linspace(3e-9, 3e-6, num = 10000),
np.linspace(3e-6, 8e-6, num = 100),
np.linspace(8e-6, 8e-5, num=100),
np.linspace(8e-5, sample_thickness, num=100),
]
)
model.mesh = F.MeshFromVertices(vertices)
# ### Material ###
damaged_tungsten = F.Material(1, D_0, E_D)
model.materials = damaged_tungsten
# ### Source ###
# Deuterium Beam Profile (S = flux * f(x))
distribution = (
1 / (sigma * (2 * np.pi) ** 0.5) * sp.exp(-0.5 * ((F.x - R_p) / sigma) ** 2)
)
ion_flux = sp.Piecewise((flux * distribution, F.t < t_imp), (0, True))
source_term = F.Source(value=ion_flux, volume=1, field=0)
model.sources = [source_term]
# ### Boundary Conditions ###
model.boundary_conditions = [
F.DirichletBC(surfaces=[1, 2], value=0, field="solute")
]
# ### Temperature ###
model.T = F.Temperature(
value=sp.Piecewise(
(T_imp, F.t < t_imp),
(T_rest, F.t < start_tds),
(min_temp + Beta * (F.t - start_tds), True),
)
)
# ### Trap Settings ###
k_0 = D_0 / (1.1e-10**2 * 6 * w_atom_density)
trap_1 = F.Trap(
E_k=damaged_tungsten.E_D,
k_0=k_0,
E_p=1.04,
p_0=1e13,
density=2.4e22,
materials=damaged_tungsten,
)
neutron_induced_traps = []
damage_dist = 1 / (1 + sp.exp((F.x - 1.8e-6) / 9e-8)) # 0.5 * (sp.tanh(- (F.x - 1.8e-06) / 1.7e-07) + 1)#
for E_p, density in zip(detrapping_energies, densities):
neutron_induced_traps.append(
F.Trap(
k_0=k_0,
E_k=damaged_tungsten.E_D,
p_0=1e13,
E_p=E_p,
density=density * damage_dist,
materials=damaged_tungsten,
)
)
model.traps = [trap_1] + neutron_induced_traps
model.dt = F.Stepsize(
initial_value=1,
stepsize_change_ratio=1.1,
dt_min=1e-5,
t_stop=resting_time - 20,
stepsize_stop_max=50,
milestones=export_times,
)
model.settings = F.Settings(
absolute_tolerance=1e10,
relative_tolerance=1e-5,
final_time=start_tds + (max_temp - min_temp) / Beta, # time to reach max temp
traps_element_type="DG"
)
model.log_level = 20
list_of_derived_quantities = [
F.TotalVolume("solute", volume=1),
F.TotalVolume("retention", volume=1),
F.HydrogenFlux(surface=1),
F.HydrogenFlux(surface=2),
]
list_of_derived_quantities += [
F.TotalVolume(f"{i}", volume=1) for i in range(1, num_traps + 1)
]
derived_quantities = F.DerivedQuantities(list_of_derived_quantities,filename=file_name)
model.exports = [derived_quantities,
F.TXTExport(
field="retention",
times=export_times,
folder="tds",
label="retention",
),
F.TXTExport(
field="solute",
times=export_times,
folder="tds",
label="solute",
),
F.TXTExport(
field="1",
times=export_times,
folder="tds",
label="trap1",
),
F.TXTExport(
field="2",
times=export_times,
folder="tds",
label="trap2",
),
F.TXTExport(
field="3",
times=export_times,
folder="tds",
label="trap3",
),
F.TXTExport(
field="4",
times=export_times,
folder="tds",
label="trap4",
),
F.TXTExport(
field="5",
times=export_times,
folder="tds",
label="trap5",
),
F.TXTExport(
field="6",
times=export_times,
folder="tds",
label="trap6",
),
]
model.initialise()
model.run()
t = np.array(derived_quantities.t)
retention = np.array(derived_quantities.filter(fields="retention").data)
flux_left = np.array(derived_quantities.filter(fields="solute", surfaces=1).data)
flux_right = np.array(derived_quantities.filter(fields="solute", surfaces=2).data)
flux_total = -flux_left - flux_right
temp = min_temp + Beta * (t - start_tds)
idx = np.where(t > start_tds)
plt.figure(1)
plt.ylabel(r"Total Retention (m$^{-1}$)")
plt.xlabel(r"Time (s)")
plt.plot(t,retention, label = f"{dpa}")
plt.legend()
plt.figure(2)
plt.plot(temp[idx], flux_total[idx], linestyle="dashed", color="tab:grey", linewidth=2)
if dpa == 0.1:
plt.figure(3)
plt.title("Damage = 0.1 dpa")
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Temperature (K)")
plt.plot(temp[idx], flux_total[idx], linewidth=3, label="FESTIM")
colors = [(0.9 * (i % 2), 0.2 * (i % 4), 0.4 * (i % 3)) for i in range(num_traps)]
trap_data = [derived_quantities.filter(fields=f"{i}").data for i in range(1, num_traps + 1)]
contributions = [
-np.diff(np.array(trap)[idx]) / np.diff(t[idx]) for trap in trap_data
]
for i, cont in enumerate(contributions):
if i == 0:
label = "Trap 1"
else:
label = f"Trap D{i}"
plt.plot(temp[idx][1:], cont, linestyle="--", color=colors[i], label=label)
plt.fill_between(temp[idx][1:], 0, cont, facecolor="grey", alpha=0.1)
plt.legend()
return derived_quantities
dpa_to_quantities = {}
for dpa, densities in dpa_n_i.items():
dpa_to_quantities[dpa] = festim_sim(densities)
for i in [3, 2]:
plt.figure(i)
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Temperature (K)")
plt.ylim(bottom=0, top=1e17)
ax = plt.gca()
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
from matplotlib import cm, colors
norm = colors.LogNorm(vmin=min(list(dpa_n_i.keys())[1:]), vmax=max(dpa_n_i.keys())) #using [1:] indexing to ignore 0
colorbar = cm.viridis
sm = plt.cm.ScalarMappable(cmap=colorbar, norm=norm)
# Plotting color bar
from mpl_toolkits.axes_grid1 import make_axes_locatable
plt.figure(2)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
plt.colorbar(sm, cax=cax, label="Damage (dpa)")
plt.figure(4)
x = np.genfromtxt("tds/solute_" + str(export_times[0]) + ".0s.txt")[:,0]
solute_x = np.genfromtxt("tds/solute_" + str(export_times[0]) + ".0s.txt")[:,1]
plt.semilogx(x, solute_x, label = "solute")
for i in range(num_traps):
trap_x = np.genfromtxt(f"tds/trap{i+1}_" + str(export_times[0]) + ".0s.txt")[:,1]
plt.semilogx(x, trap_x, label=f"trap {i+1}")
plt.legend()
plt.ylabel(r"Concentration (m$^{-3}$)")
plt.xlabel(r"Depth (m)")
plt.show()
Example_TDSfromdamagedW.zip (3.0 KB)
