import festim as F
import ufl
import dolfinx
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# # ### 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]
dpa_n_i = {
0: [],
0.001: [1e24, 2.5e24, 1e24, 1e24, 2e23],
0.005: [3.5e24, 5e24, 2.5e24, 1.9e24, 1.6e24],
0.023: [2.2e25, 1.5e25, 6.5e24, 2.1e25, 6e24],
0.1: [4.8e25, 3.8e25, 2.6e25, 3.6e25, 1.1e25],
0.23: [5.4e25, 4.4e25, 3.6e25, 3.9e25, 1.4e25],
0.5: [5.5e25, 4.6e25, 4e25, 4.5e25, 1.7e25],
2.5: [5.8e25, 6.5e25, 4.5e25, 5.5e25, 2e25], # re-fit
}
# 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 = 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
center = 0.7e-9
width = 0.5e-9
def festim_sim(densities):
model = F.HydrogenTransportProblem()
vertices = np.concatenate(
[
np.linspace(0, 3e-9, num=200),
np.linspace(3e-9, 8e-6, num=200),
np.linspace(8e-6, 8e-5, num=200),
np.linspace(8e-5, sample_thickness, num=200),
]
)
model.mesh = F.Mesh1D(vertices)
# ### Material ###
damaged_tungsten = F.Material(D_0, E_D)
volume = F.VolumeSubdomain1D(
id=1, borders=[0, sample_thickness], material=damaged_tungsten
)
left_boundary = F.SurfaceSubdomain1D(id=1, x=0)
right_boundary = F.SurfaceSubdomain1D(id=2, x=sample_thickness)
model.subdomains = [volume, left_boundary, right_boundary]
H = F.Species("H")
instrinsic_trapped_H = F.Species(f"Trapped 1", mobile=False)
neutron_induced_trapped_species = [
F.Species(f"Trapped D{i+1}", mobile=False) for i in range(len(densities))
]
damage_dist = lambda x: 1 / (1 + ufl.exp((x[0] - 2.5e-06) / 5e-07))
empty_neutron_induced_traps = []
for i, (density, trapped_spe) in enumerate(
zip(densities, neutron_induced_trapped_species)
):
empty_neutron_induced_traps.append(
F.ImplicitSpecies(
n=lambda x, density=density: density * damage_dist(x),
others=[trapped_spe],
name=f"empty {i+2}",
)
)
empty_intrinsic_traps = F.ImplicitSpecies(
n=2.4e22, others=[instrinsic_trapped_H], name="empty 1"
)
model.species = [H] + [instrinsic_trapped_H] + neutron_induced_trapped_species
# ### Source ###
# Deuterium Beam Profile (S = flux * f(x))
distribution = lambda x: (
1 / (sigma * (2 * ufl.pi) ** 0.5) * ufl.exp(-0.5 * ((x[0] - R_p) / sigma) ** 2)
)
ion_flux = lambda x, t: ufl.conditional(t < t_imp, flux * distribution(x), 0)
source_term = F.ParticleSource(value=ion_flux, volume=volume, species=H)
model.sources = [source_term]
# ### Boundary Conditions ###
model.boundary_conditions = [
F.FixedConcentrationBC(subdomain=left_boundary, value=0, species=H),
F.FixedConcentrationBC(subdomain=right_boundary, value=0, species=H),
]
# ### Temperature ###
def temp_fun(t):
if t < t_imp:
return T_imp
elif t < start_tds:
return T_rest
else:
return min_temp + Beta * (t - start_tds)
model.temperature = temp_fun
# ### Trap Settings ###
k_0 = D_0 / (1.1e-10**2 * 6 * w_atom_density)
trapping_reaction_1 = F.Reaction(
reactant=[H, empty_intrinsic_traps],
product=[instrinsic_trapped_H],
k_0=k_0,
E_k=damaged_tungsten.E_D,
p_0=1e13,
E_p=1.04,
volume=volume,
)
reactions = [trapping_reaction_1]
for i, _ in enumerate(densities):
reactions.append(
F.Reaction(
reactant=[H, empty_neutron_induced_traps[i]],
product=[neutron_induced_trapped_species[i]],
k_0=k_0,
E_k=damaged_tungsten.E_D,
p_0=1e13,
E_p=detrapping_energies[i],
volume=volume,
)
)
model.reactions = reactions
model.settings = F.Settings(
atol=1e10,
rtol=1e-10,
final_time=start_tds + (max_temp - min_temp) / Beta, # time to reach max temp
)
model.settings.stepsize = F.Stepsize(
initial_value=2,
growth_factor=1.1,
cutback_factor=0.9,
target_nb_iterations=4,
# cap dt so the implantation transient and the TDS ramp stay resolved
max_stepsize=lambda t: 50 if t > t_imp + t_rest * 0.5 else 200,
)
derived_quantities = [
F.TotalVolume(field=spe, volume=volume) for spe in model.species[1:]
]
flux_left = F.SurfaceFlux(field=H, surface=left_boundary)
flux_right = F.SurfaceFlux(field=H, surface=right_boundary)
derived_quantities.append(flux_left)
derived_quantities.append(flux_right)
model.exports = derived_quantities
# These are the same solver options FESTIM uses by default, with one
# change: the "bt" line search, without which the solver fails to
# converge once the traps start to fill up. We have to write out the
# whole list because festim 2.0b0 (the version this book uses) drops
# its defaults when given this dictionary. From festim 2.0 onwards,
# the "snes_linesearch_type" line alone would be enough.
model.petsc_options = {
"snes_type": "newtonls",
"snes_linesearch_type": "bt",
"snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
"snes_atol": 1e10,
"snes_rtol": 1e-10,
"snes_max_it": 30,
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
}
model.initialise()
model.run()
return derived_quantities
