Hi,I’m trying to fit 1D TDS simulation of deuterium in tungsten with FESTIM.But since I used sp.Piecewise to set local depth traps,the Newton solver wouldn’t converge.I wonder where my settings went wrong.Thank you！
Here’s my code and raw experimental data:
import festim as F
import numpy as np
import sympy as sp
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)
def TDS(h1,n1, E_p1,h2,n2,E_p2):
"""Runs the simulation with parameters p that represent:
Args:
n (float): concentration of trap 1, at. fr.
E_p (float): detrapping barrier from trap 1, eV
Returns:
F.DerivedQuantities: the derived quantities of the simulation
"""
atom_density = 6.28e28 # atom/m3
# Define Simulation object
synthetic_TDS = F.Simulation()
synthetic_TDS.log_level = 40
# Define a simple mesh
vertices = np.concatenate([
np.linspace(0, 1e-4, num=100),
np.linspace(1e-4, 5e-4, num=200),
np.linspace(5e-4, 1.5e-3, num=200),
])
synthetic_TDS.mesh = F.MeshFromVertices(vertices)
trap_conc1 = sp.Piecewise((n1 * atom_density,F.x< h1),
(0,True))
trap_conc2 = sp.Piecewise((n2 * atom_density,F.x< h2),
(0,True))
# Define material properties
W = F.Material(
id=1,
D_0=2.9e-7,
E_D=0.39,
)
synthetic_TDS.materials = W
# Define traps
trap_1 = F.Trap(
k_0=2.9e12/atom_density,
E_k=0.39,
p_0=1e13,
E_p=E_p1,
density=trap_conc1,
materials=W,
)
trap_2 = F.Trap(
k_0=2.9e12/atom_density,
E_k=0.39,
p_0=1e13,
E_p=E_p2,
density=trap_conc2,
materials=W,
)
synthetic_TDS.traps = [trap_1,trap_2]
# Set initial conditions
synthetic_TDS.initial_conditions = [
F.InitialCondition(field="1", value=trap_conc1),
F.InitialCondition(field="2", value=trap_conc2),
]
# Set boundary conditions
synthetic_TDS.boundary_conditions = [
F.DirichletBC(surfaces=[1, 2], value=0, field=0)
]
# Define the material temperature evolution
ramp = 1 # K/s
synthetic_TDS.T = F.Temperature(value=317.8+ ramp * (F.t))
# Define the simulation settings
synthetic_TDS.dt = F.Stepsize(
initial_value=1,
stepsize_change_ratio=1.2,
dt_min=1e-6,
max_stepsize=2)
synthetic_TDS.settings = F.Settings(
absolute_tolerance=1e12,
relative_tolerance=1e-8,
final_time=1154,
maximum_iterations=100,
)
# Define the exports
derived_quantities = F.DerivedQuantities(
[
F.HydrogenFlux(surface=1),
F.HydrogenFlux(surface=2),
F.AverageVolume(field="T", volume=1),
]
)
synthetic_TDS.exports = [derived_quantities]
synthetic_TDS.initialise()
synthetic_TDS.run()
return derived_quantities
import matplotlib.pyplot as plt
import os
ref = np.genfromtxt('480k-1dpa.csv', delimiter=',')
T = ref[:, 0]
flux = ref[:, 1]
def info(i, p):
"""
Print information during the fitting procedure
"""
print("-" * 40)
print(f"i = {i}")
print("New simulation.")
print(f"Point is: {p}")
from scipy.interpolate import interp1d
prms = []
errors = []
def error_function(prm):
"""
Compute average absolute error between simulation and reference
"""
global i
global prms
global errors
prms.append(prm)
i += 1
info(i, prm)
# Filter the results if a negative value is found
if any([e < 0 for e in prm]):
return 1e30
# Get the simulation result
h1,n1, Ep1,h2,n2,Ep2 = prm
res = TDS(h1,n1, Ep1,h2,n2,Ep2)
T = np.array(res.filter(fields="T").data)
flux = -np.array(res.filter(fields="solute", surfaces=1).data) - np.array(
res.filter(fields="solute", surfaces=2).data)
# Plot the intermediate TDS spectra
if i == 1:
plt.plot(T, flux, color="tab:red", lw=2, label="Initial guess")
else:
plt.plot(T, flux, color="tab:grey", lw=0.5)
interp_tds = interp1d(T, flux, fill_value="extrapolate")
# Compute the mean absolute error between sim and ref
err = np.nanmean(np.abs(interp_tds(ref[:, 0]) - ref[:, 1]).mean())
print(f"Average absolute error is : {err:.2e}")
errors.append(err)
return err
from scipy.optimize import minimize
i = 0 # initialise counter
# Set the tolerances
fatol = 1e16
xatol = 1e-3
initial_guess = [0.0005,0.0005,1,0.0005,0.001,1.26]
# Minimise the error function
res = minimize(
error_function,
np.array(initial_guess),
method="Nelder-Mead",
options={"disp": True, "fatol": fatol, "xatol": xatol},
)
# Process the obtained results
predicted_data = TDS(*res.x)
T = predicted_data.filter(fields="T").data
flux_left = predicted_data.filter(fields="solute", surfaces=1).data
flux_right = predicted_data.filter(fields="solute", surfaces=2).data
flux_total = -(np.array(flux_left) + np.array(flux_right))
# Visualise
plt.plot(ref[:, 0], ref[:, 1], linewidth=2, label="Reference")
plt.plot(T, flux_total, linewidth=2, label="Optimised")
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Temperature (K)")
plt.legend()
plt.show()
plt.savefig('480k-1dpa.png')
reference_prms = [E_p, trap_conc]
plt.scatter(
np.array(prms)[:, 0], np.array(prms)[:, 1], c=np.array(errors), cmap="viridis"
)
plt.plot(np.array(prms)[:, 0], np.array(prms)[:, 1], color="tab:grey", lw=0.5)
plt.scatter(*reference_prms, c="tab:red")
plt.annotate(
"Reference",
xy=reference_prms,
xytext=(reference_prms[0] - 0.003, reference_prms[1] + 0.1),
arrowprops=dict(facecolor="black", arrowstyle="-|>"),
)
plt.annotate(
"Initial guess",
xy=initial_guess,
xytext=(initial_guess[0] - 0.004, initial_guess[1] + 0.05),
arrowprops=dict(facecolor="black", arrowstyle="-|>"),
)
plt.xlabel(r"Trap 1 concentration (at. fr.)")
plt.ylabel(r"Detrapping barrier (eV)")
plt.show()
plt.savefig('scatter.png')
480k-1dpa.zip (4.9 KB)
