# Complicated trap distribution

**URL:** <https://festim.discourse.group/t/complicated-trap-distribution/83>\
**Category:** User Support\
**Created:** [November 5, 2024, 1:41pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83 "2024-11-05T13:41:44Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![cruiseraurora](https://avatars.discourse-cdn.com/v4/letter/c/ecd19e/32.png) [@cruiseraurora](https://festim.discourse.group/u/cruiseraurora)\
**Post date:** [November 5, 2024, 1:41pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/1 "2024-11-05T13:41:44Z")

</div>

Hello! I’m not very familiar with python. Using sp.Piecewise, how can I set the trap distribution with these conditions: layer of oxide (10 nm), layer of steel (2 mm), layer of oxide (10 nm), and the trap is located on both sides of the steel with 1 μm trap layer thickness? I think trap distribution should be in the coordinates described by these inequalities:  
(10 nm \< x \< 10 nm + 1 μm) & (10 nm + 1 μm + 2 mm \< x \< 10 nm + 1 μm + 2 mm + 1 μm),  
if in FESTIM the coordinate x zero point is considered to be on the left side of the first (oxide) layer. Unfortunately, I don’t know how set this inequalities in sp.Piecewise  
I know that I can compute only half of the problem due to symmetry but I’m interested in the general principle of setting such complicated trap distribution

---

<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:** [November 5, 2024, 2:28pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/2 "2024-11-05T14:28:47Z")

</div>

Hi @cruiseraurora

If I understand your problem you have 3 materials (2 oxide layers and a steel substrate) and you have a trap that is located only in the steel layer but in the first micron (on each side).

Consider the following steady state example:

```python
import festim as F
import numpy as np
import sympy as sp

my_model = F.Simulation()

L_oxide_1 = 10e-9
L_steel = 2e-3
L_oxide_2 = 10e-9

# Define the mesh with local refinement around the trap layers
my_model.mesh = F.MeshFromVertices(
    np.concatenate(
        [
            np.linspace(0, L_oxide_1, 20),
            np.linspace(L_oxide_1, L_oxide_1 + 1.5e-6, 300),
            np.linspace(L_oxide_1 + 1.5e-6, L_oxide_1 + L_steel - 1.5e-6, 100),
            np.linspace(L_oxide_1 + L_steel - 1.5e-6, L_oxide_1 + L_steel, 300),
            np.linspace(L_oxide_1 + L_steel, L_oxide_1 + L_steel + L_oxide_2, 20),
        ]
    )
)

oxide_1 = F.Material(id=1, D_0=1, E_D=0, borders=[0, L_oxide_1])
steel = F.Material(id=2, D_0=1, E_D=0, borders=[L_oxide_1, L_oxide_1 + L_steel])
oxide_2 = F.Material(
    id=3, D_0=1, E_D=0, borders=[L_oxide_1 + L_steel, L_oxide_1 + L_steel + L_oxide_2]
)

my_model.materials = [oxide_1, steel, oxide_2]

trap_density = 1 # trap/m3
distribution = sp.Piecewise(
    (1, sp.And(L_oxide_1 < F.x, F.x < L_oxide_1 + 1e-6)), # 1 if L_oxide_1 < x < L_oxide_1 + 1e-6
    (1, sp.And(L_oxide_1 + L_steel - 1e-6 < F.x, F.x < L_oxide_1 + L_steel)), # 1 if L_oxide_1 + L_steel - 1e-6 < x < L_oxide_1 + L_steel
    (0, True), # 0 otherwise
)

my_model.traps = F.Trap(
    k_0=1,
    E_k=0,
    p_0=1e-10,
    E_p=0,
    materials=[steel], # trap in steel only
    density=trap_density * distribution,
)

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

my_model.T = 500 # ignored here since activation energies are zero

my_model.settings = F.Settings(
    absolute_tolerance=1e-10,
    relative_tolerance=1e-10,
    traps_element_type="DG",
    transient=False,
)

my_model.exports = [F.XDMFExport(field="1", checkpoint=False)]

my_model.initialise()
my_model.run()

```

Reading the xdmf file in paraview you obtain this on the left side, you have a trapped concentration of 1 in the first micron of the steel:

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

---

<div class="post-metadata">

**Author:** ![cruiseraurora](https://avatars.discourse-cdn.com/v4/letter/c/ecd19e/32.png) [@cruiseraurora](https://festim.discourse.group/u/cruiseraurora)\
**Post date:** [November 8, 2024, 11:22am UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/3 "2024-11-08T11:22:13Z")

</div>

Now I understand how it works! Thank you!

---

<div class="post-metadata">

**Author:** ![cruiseraurora](https://avatars.discourse-cdn.com/v4/letter/c/ecd19e/32.png) [@cruiseraurora](https://festim.discourse.group/u/cruiseraurora)\
**Post date:** [March 3, 2025, 12:12pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/4 "2025-03-03T12:12:44Z")

</div>

Hello again! I tried to run my code with complicated trap distribunion and specific parameters (such as oxide diffusicon energy, steel duffusion energy, etc.) but it seems like the solution doesn’t converge (the concentration oscillations is near to the right side of trap layer)

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

If you try to run code without traps you will get the same TDS-spectrum that I calculate in TMAP7.

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

With near-surface traps it looks like all the deuterium go out immmediately.

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

This is what a TDS spectrum with near-surface traps shoud look like (calculated in TMAP, results in TMAP7 and TMAP8 are the same):

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

Yes, I set `traps_element_type="DG"`.  
I would be grateful if you could point out the mistakes I made in FESTIM.

> **My FESTIM code**
>
> ```auto
> import festim as F
> import sympy as sp
> import numpy as np
> import matplotlib.pyplot as plt
> my_model = F.Simulation()
> 
> oxide_thick = 10e-9
> steel_thick = 2e-3
> 
> oxide_left = F.Material(
> id=1,
> D_0=2e-8,
> E_D=0.63,
> borders=[0, oxide_thick]
> )
> rusfer = F.Material(
> id=2,
> D_0=0.85e-7,
> E_D=0.165,
> borders=[oxide_thick, steel_thick + oxide_thick]
> )
> oxide_right = F.Material(
> id=3,
> D_0=2e-8,
> E_D=0.63,
> borders=[steel_thick + oxide_thick, steel_thick + 2 * oxide_thick]
> )
> my_model.materials = [oxide_left, rusfer, oxide_right]
> 
> my_model.mesh = F.MeshFromVertices(
> np.concatenate(
> [
> np.linspace(0, oxide_thick, 500),
> np.linspace(oxide_thick, oxide_thick + 5e-6, 500),
> np.linspace(oxide_thick + 5e-6, oxide_thick + steel_thick - 5e-6, 500),
> np.linspace(oxide_thick + steel_thick - 5e-6, oxide_thick + steel_thick, 500),
> np.linspace(oxide_thick + steel_thick, oxide_thick + steel_thick + oxide_thick, 500),
> ]
> )
> )
> 
> rusfer_atom_density = 8.21e28 # atom/m3
> trap_density = 1.75e-2
> distribution = sp.Piecewise(
> (1, sp.And(oxide_thick < F.x, F.x < oxide_thick + 1e-6)), # 1 if L_oxide_1 < x < L_oxide_1 + 1e-6
> (1, sp.And(oxide_thick + steel_thick - 1e-6 < F.x, F.x < oxide_thick + steel_thick)), # 1 if L_oxide_1 + L_steel - 1e-6 < x < L_oxide_1 + L_steel
> (0, True), # 0 otherwise
> )
> trap_1 = F.Trap(
> k_0=0.85e-7/(2.87e-10**2),
> E_k=0.165,
> p_0=1e13,
> E_p=0.75,
> # (10 nm < x < 10 nm + 1 mkm) и (10 nm + 1 mkm + 2 mm < x < 10 nm + 1 mkm + 2 mm + 1 mkm)
> density= rusfer_atom_density * trap_density * distribution,
> materials=rusfer
> )
> my_model.traps = trap_1
> 
> fraction_filled = 0.01
> my_model.initial_conditions = [
> F.InitialCondition(field="solute", value=0.85e24),
> F.InitialCondition(field = "1", value = rusfer_atom_density* trap_density * distribution * fraction_filled),
> ]
> 
> my_model.boundary_conditions = [
> F.RecombinationFlux(1e-29, -0.35, 2, surfaces = [1,2])
> ]
> 
> temperature_ramp = 2 # K/s
> my_model.T = 293 + F.t*temperature_ramp
> 
> my_model.dt = F.Stepsize(
> initial_value=1e-3,
> stepsize_change_ratio=1.1,
> max_stepsize = 2,
> dt_min=1e-2,
> milestones = [0, 753]
> )
> my_model.settings = F.Settings(
> absolute_tolerance=1e10,
> relative_tolerance=1e-09,
> final_time=753,
> traps_element_type="DG",
> )
> 
> list_of_derived_quantities = [
> F.TotalVolume("solute", volume=1),
> F.TotalVolume("retention", volume=1),
> F.TotalVolume("1", 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/EK-181.csv", # optional set a filename to export the data to csv
> show_units=True,
> )
> 
> my_model.exports = [derived_quantities]
> #my_model.exports = [F.XDMFExport(field="1", checkpoint=False)]
> my_model.initialise()
> my_model.run()
> 
> t = derived_quantities.t
> flux_left = derived_quantities.filter(fields="solute", surfaces=1).data
> flux_right = derived_quantities.filter(fields="solute", surfaces=2).data
> Temp = derived_quantities.filter(fields="T").data
> flux_total = -np.array(flux_left) - np.array(flux_right)
> 
> plt.plot(Temp, flux_total, '.', linewidth=3)
> plt.xlim(293, 1800)
> # plt.ylim(bottom=-1.25e18, top=0.6e19)
> plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
> plt.xlabel(r"Temperature, K")
> plt.grid()
> plt.show()
> 
> ```

---

<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 3, 2025, 12:42pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/5 "2025-03-03T12:42:52Z")

</div>

Hi @cruiseraurora so if I understand well you have a good agreement between TMAP7, TMAP8 and FESTIM but you want to add additional traps, and that’s not working (ie. oscillations in the concentration profile). Is this correct?

---

<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 3, 2025, 2:06pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/6 "2025-03-03T14:06:49Z")

</div>

I have reproduced the issue locally. I believe this is a parametrisation issue, the trapping rate seems very high. Could you provide the parametrisation for TMAP7 with this near-surface trap please?

You can consider replacing the recombination BC by am homogeneous dirichlet BC (ie. c\_mobile = 0) for simplicity (while debugging).

---

<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 3, 2025, 2:42pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/7 "2025-03-03T14:42:24Z")

</div>

Hi @cruiseraurora!

Please, consider this code for a simplified problem (a half of the sample, a simpler mesh, and the updated value of the trapping rate).

> **Code**
>
> ```auto
> import festim as F
> import sympy as sp
> import numpy as np
> import matplotlib.pyplot as plt
> 
> my_model = F.Simulation()
> 
> oxide_thick = 10e-9
> steel_thick = 2e-3
> 
> oxide_left = F.Material(id=1, D_0=2e-8, E_D=0.63, borders=[0, oxide_thick])
> rusfer = F.Material(
> id=2, D_0=0.85e-7, E_D=0.165, borders=[oxide_thick, steel_thick / 2 + oxide_thick]
> )
> 
> my_model.materials = [oxide_left, rusfer]
> 
> my_model.mesh = F.MeshFromVertices(
> np.concatenate(
> [
> np.linspace(0, 5 * oxide_thick, 2000),
> np.linspace(5 * oxide_thick, steel_thick / 2 + oxide_thick, num=2000),
> ]
> )
> )
> 
> rusfer_atom_density = 8.21e28 # atom/m3
> trap_density = 1.75e-2
> distribution = sp.Piecewise(
> (
> 1,
> sp.And(oxide_thick < F.x, F.x < oxide_thick + 1e-6),
> ), # 1 if L_oxide_1 < x < L_oxide_1 + 1e-6
> (0, True), # 0 otherwise
> )
> trap_1 = F.Trap(
> k_0=0.85e-7 / (2.87e-10**2 * rusfer_atom_density),
> E_k=0.165,
> p_0=1e13,
> E_p=0.75,
> # (10 nm < x < 10 nm + 1 mkm) и (10 nm + 1 mkm + 2 mm < x < 10 nm + 1 mkm + 2 mm + 1 mkm)
> density=rusfer_atom_density * trap_density * distribution,
> materials=rusfer,
> )
> my_model.traps = trap_1
> 
> fraction_filled = 0.01
> my_model.initial_conditions = [
> F.InitialCondition(field="solute", value=0.85e24),
> F.InitialCondition(
> field="1",
> value=rusfer_atom_density * trap_density * distribution * fraction_filled,
> ),
> ]
> 
> my_model.boundary_conditions = [F.RecombinationFlux(1e-29, -0.35, 2, surfaces=1)]
> 
> temperature_ramp = 2 # K/s
> my_model.T = 293 + F.t * temperature_ramp
> 
> my_model.dt = F.Stepsize(
> initial_value=1e-3,
> stepsize_change_ratio=1.1,
> max_stepsize=5,
> dt_min=1e-2,
> milestones=[0, 753],
> )
> my_model.settings = F.Settings(
> absolute_tolerance=1e12,
> relative_tolerance=1e-8,
> final_time=753,
> traps_element_type="DG",
> )
> 
> list_of_derived_quantities = [
> F.TotalVolume("solute", volume=1),
> F.TotalVolume("retention", volume=1),
> F.TotalVolume("1", 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/EK-181.csv", # optional set a filename to export the data to csv
> show_units=True,
> )
> 
> TXT = [
> F.TXTExport(
> field="1",
> filename="./trap1.txt",
> write_at_last=True,
> times=[753],
> )
> ]
> my_model.exports = [derived_quantities]
> # my_model.exports = [F.XDMFExport(field="1", checkpoint=False)]
> my_model.initialise()
> my_model.run()
> 
> t = derived_quantities.t
> flux_left = derived_quantities.filter(fields="solute", surfaces=1).data
> # flux_right = derived_quantities.filter(fields="solute", surfaces=2).data
> Temp = derived_quantities.filter(fields="T").data
> flux_total = -2 * np.array(flux_left) # - np.array(flux_right)
> 
> plt.plot(Temp, flux_total, linewidth=3)
> plt.xlim(293, 1800)
> # plt.ylim(bottom=-1.25e18, top=0.6e19)
> plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
> plt.xlabel(r"Temperature, K")
> plt.grid()
> plt.show()
> 
> ```

It produces the following TDS spectrum:

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

---

<div class="post-metadata">

**Author:** ![cruiseraurora](https://avatars.discourse-cdn.com/v4/letter/c/ecd19e/32.png) [@cruiseraurora](https://festim.discourse.group/u/cruiseraurora)\
**Post date:** [March 7, 2025, 1:27pm UTC](https://festim.discourse.group/t/complicated-trap-distribution/83/8 "2025-03-07T13:27:42Z")

</div>

Thank you very much! I apologize for not noticing my mistake in the expression for trapping rate