# Simulation doesn't converge under local traps condition

**URL:** <https://festim.discourse.group/t/simulation-doesnt-converge-under-local-traps-condition/104>\
**Category:** User Support\
**Created:** [March 9, 2025, 7:00am UTC](https://festim.discourse.group/t/simulation-doesnt-converge-under-local-traps-condition/104 "2025-03-09T07:00:56Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jellycarrot](https://avatars.discourse-cdn.com/v4/letter/j/7ea924/32.png) [@Jellycarrot](https://festim.discourse.group/u/Jellycarrot)\
**Post date:** [March 9, 2025, 7:00am UTC](https://festim.discourse.group/t/simulation-doesnt-converge-under-local-traps-condition/104/1 "2025-03-09T07:00:56Z")

</div>

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:

```auto
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](https://festim.discourse.group/uploads/short-url/a6goki8G2W4q0zEwUjP0TxWXG40.zip) (4.9 KB)

---

<div class="post-metadata">

**Author:** ![remidm](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/remidm/32/4_2.png) [@remidm](https://festim.discourse.group/u/remidm)\
**Post date:** [March 9, 2025, 4:01pm UTC](https://festim.discourse.group/t/simulation-doesnt-converge-under-local-traps-condition/104/2 "2025-03-09T16:01:59Z")

</div>

Hi @Jellycarrot your code is a bit long and hard to parse. Can you make a minimal working example showing your issue?

For example, to show the simulation isn’t converging, you don’t need to have:

- the experimental data
- the optimisation
- the error calculation
- the info function
- the plotting

---

<div class="post-metadata">

**Author:** ![Jellycarrot](https://avatars.discourse-cdn.com/v4/letter/j/7ea924/32.png) [@Jellycarrot](https://festim.discourse.group/u/Jellycarrot)\
**Post date:** [March 10, 2025, 8:47am UTC](https://festim.discourse.group/t/simulation-doesnt-converge-under-local-traps-condition/104/3 "2025-03-10T08:47:27Z")

</div>

OK Professor Remi，thanks for your prompt reply.

Here’s a single simulation:

```auto
#Simulation of a Thermo-Desorption experiment

import festim as F
import os
my_model = F.Simulation()

import numpy as np
import sympy as sp
import warnings

warnings.filterwarnings("ignore", category=DeprecationWarning)
vertices = np.linspace(0,1.5e-3,num =1500)

my_model.mesh = F.MeshFromVertices(vertices)

atom_density =6.28e28 
W = F.Material(
    id=1,
    D_0=2.9e-7, # m2/s
    E_D=0.39, # eV
)

my_model.materials = W

trap_conc1 = sp.Piecewise((5e-4 * atom_density, F.x < 1e-3),
                          (0, True))
trap_conc2 = sp.Piecewise((1e-3 * atom_density, F.x < 1.5e-6),
                          (0, True))

my_model.initial_conditions = [
		F.InitialCondition(field="1", value=trap_conc1),
		F.InitialCondition(field="2", value=trap_conc2)
	]
trap_1 = F.Trap(
        k_0=2.9e12/atom_density,
        E_k=0.39,
        p_0=1e13,
        E_p=1,
        density=trap_conc1,
        materials= W
    )
trap_2 = F.Trap(
        k_0=2.9e12/atom_density,
        E_k=0.39,
        p_0=1e13,
        E_p=1.26,
        density=trap_conc2,
        materials= W
    )
my_model.traps = [trap_1,trap_2]

my_model.boundary_conditions = [
    F.DirichletBC(surfaces=[1, 2], value=0, field=0)
]

ramp = 1 # K/s
my_model.T = F.Temperature(value=304 + ramp * (F.t))

#Stepsize
my_model.dt = F.Stepsize(
    initial_value=0.1,
    stepsize_change_ratio=1.2,
    dt_min=1e-6,
    max_stepsize=1,
)

#Settings
my_model.settings = F.Settings(
    absolute_tolerance=1e10,
    relative_tolerance=1e-10,
    final_time = 1167
)

list_of_derived_quantities = [
        F.TotalVolume("solute", volume=1),
        F.TotalVolume("retention", volume=1),
        F.TotalVolume("1", volume=1),
	F.TotalVolume("2",volume=1),
        F.HydrogenFlux(surface=1),
        F.HydrogenFlux(surface=2),
	F.AverageVolume(field="T", volume=1),
    ]

derived_quantities = F.DerivedQuantities(
    list_of_derived_quantities,
    # filename="tds/derived_quantities.csv" # optional set a filename to export the data to csv
)

my_model.exports = [derived_quantities]

my_model.initialise()
my_model.run()

t = derived_quantities.t
T = derived_quantities.filter(fields="T").data
flux_left = derived_quantities.filter(fields="solute", surfaces=1).data
flux_right = derived_quantities.filter(fields="solute", surfaces=2).data

flux_total = (-np.array(flux_left) - np.array(flux_right))
kk = flux_total.astype(str)
#with open('11.txt', 'w') as file:
    #for row in kk:
        #file.write('\t'.join(map(str, row)) + '\n')

import matplotlib.pyplot as plt
plt.plot(T, flux_total,label="Simulation")
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"T(K)")
plt.title("480k-1dpa")
plt.savefig("480-1.png",dpi=600)

```

---

<div class="post-metadata">

**Author:** ![VVKulagin](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/vvkulagin/32/51_2.png) [@VVKulagin](https://festim.discourse.group/u/VVKulagin)\
**Post date:** [March 10, 2025, 12:50pm UTC](https://festim.discourse.group/t/simulation-doesnt-converge-under-local-traps-condition/104/4 "2025-03-10T12:50:50Z")

</div>

Hi @Jellycarrot

Consider setting `traps_element_type="DG"` in the model settings and optimising your mesh. With the following code I managed to obtain a TDS specturm.

> **Code**
>
> ```auto
> import festim as F
> import numpy as np
> import sympy as sp
> 
> my_model = F.Simulation()
> 
> vertices = np.concatenate(
> [
> np.linspace(0, 1.5e-6, num=250),
> np.linspace(1.5e-6, 1e-3, num=500),
> np.linspace(1e-3, 1.5e-3, num=500),
> ]
> )
> 
> my_model.mesh = F.MeshFromVertices(vertices)
> 
> atom_density = 6.28e28
> W = F.Material(
> id=1,
> D_0=2.9e-7, # m2/s
> E_D=0.39, # eV
> )
> 
> my_model.materials = W
> 
> trap_conc1 = sp.Piecewise((5e-4 * atom_density, F.x < 1e-3), (0, True))
> trap_conc2 = sp.Piecewise((1e-3 * atom_density, F.x < 1.5e-6), (0, True))
> 
> my_model.initial_conditions = [
> F.InitialCondition(field="1", value=trap_conc1),
> F.InitialCondition(field="2", value=trap_conc2),
> ]
> trap_1 = F.Trap(
> k_0=2.9e12 / atom_density,
> E_k=0.39,
> p_0=1e13,
> E_p=1,
> density=trap_conc1,
> materials=W,
> )
> trap_2 = F.Trap(
> k_0=2.9e12 / atom_density,
> E_k=0.39,
> p_0=1e13,
> E_p=1.26,
> density=trap_conc2,
> materials=W,
> )
> my_model.traps = [trap_1, trap_2]
> 
> my_model.boundary_conditions = [F.DirichletBC(surfaces=[1, 2], value=0, field=0)]
> 
> ramp = 1 # K/s
> my_model.T = F.Temperature(value=304 + ramp * (F.t))
> 
> # Stepsize
> my_model.dt = F.Stepsize(
> initial_value=0.01,
> stepsize_change_ratio=1.05,
> dt_min=1e-6,
> max_stepsize=1,
> )
> 
> # Settings
> my_model.settings = F.Settings(
> absolute_tolerance=1e10,
> relative_tolerance=1e-10,
> final_time=1167,
> traps_element_type="DG",
> )
> 
> list_of_derived_quantities = [
> F.HydrogenFlux(surface=1),
> F.HydrogenFlux(surface=2),
> F.AverageVolume(field="T", volume=1),
> ]
> 
> derived_quantities = F.DerivedQuantities(
> list_of_derived_quantities,
> )
> 
> my_model.exports = [derived_quantities]
> 
> my_model.initialise()
> my_model.run()
> 
> t = derived_quantities.t
> T = derived_quantities.filter(fields="T").data
> flux_left = derived_quantities.filter(fields="solute", surfaces=1).data
> flux_right = derived_quantities.filter(fields="solute", surfaces=2).data
> 
> flux_total = -np.array(flux_left) - np.array(flux_right)
> 
> import matplotlib.pyplot as plt
> 
> plt.plot(T, flux_total, label="Simulation")
> plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
> plt.xlabel(r"T(K)")
> plt.title("480k-1dpa")
> plt.show()
> 
> ```

The result:

 ![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/b799060c2571bddcc9b6538e578f099b51e58b90.png)