# Time step reaching the minimum value during implantation

**URL:** <https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142>\
**Category:** User Support\
**Created:** [September 16, 2025, 2:19pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142 "2025-09-16T14:19:30Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Trinity](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Trinity](https://festim.discourse.group/u/Trinity)\
**Post date:** [September 16, 2025, 2:19pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/1 "2025-09-16T14:19:30Z")

</div>

Hi !Hope you’re having a good day! 😊 I’m trying to use FESTIM to get the time-dependent hydrogen concentration at 1, 5, and 10 microns near the surface during the injection process—it’s one of the final bits of my project, and honestly, my boss is _really_ on my case to wrap this up.

Earlier, I tried coarsening the mesh near the surface, and while the code ran without errors, it just spat out all zeros (bummer!). Now, I’m attempting a finer mesh, but the residuals keep creeping up over time, which eventually resulted in time step reaching the minimum value.

Quick heads-up: Some parts of this were AI-generated, so it might come across as a bit wordy or over-explained—I apologize for that! If you’ve got any tips or could point me , I’d be so grateful. Thanks a ton for your time and help! 🙏

My code :

```auto
import festim as F
my_model = F.Simulation()
my_model.log_level = 20
import numpy as np
import sympy as sp
import os

output_dir = "output"
if not os.path.exists(output_dir):
   os.makedirs(output_dir)
   print(f"Output directory created: {output_dir}")

# Mesh generation with refined regions near the surface

surface_region_1 = np.linspace(0, 10e-9, num=200)

surface_region_2 = np.linspace(10e-9, 50e-9, num=160)

surface_region_3 = np.linspace(50e-9, 200e-9, num=150)

surface_region_4 = np.linspace(200e-9, 1e-6, num=320)

mid_region_1 = np.linspace(1e-6, 10e-6, num=720)

mid_region_2 = np.linspace(10e-6, 100e-6, num=1800)

mid_region_3 = np.linspace(100e-6, 1e-3, num=1800)

far_region = np.linspace(1e-3, 3e-3, num=400)
vertices = np.concatenate([
   surface_region_1,
   surface_region_2,
   surface_region_3,
   surface_region_4,
   mid_region_1,
   mid_region_2,
   mid_region_3,
   far_region
])
my_model.mesh = F.MeshFromVertices(vertices)

# Tungsten material properties
tungsten = F.Material(
   id=1,
   D_0=2.9e-7, # Pre-exponential factor for diffusion (m2/s)
   E_D=0.39, # Activation energy for diffusion (eV)
)
my_model.materials = tungsten

# Total implantation time in seconds
implantation_time = 9000 # s

# Time-dependent ion flux profile for LTD experiment
# The flux changes at different time intervals to simulate experimental conditions
ion_flux = sp.Piecewise(
   (2.3e22 * F.t / 200, (F.t >= 0) & (F.t <= 200)),
   (2.3e22, (F.t > 200) & (F.t <= 1800)),
   (1.9e22, (F.t > 1800) & (F.t <= 3600)),
   (1.03e22, (F.t > 3600) & (F.t <= 5400)),
   (6.23e21, (F.t > 5400) & (F.t <= 7200)),
   (2.5e21, (F.t > 7200) & (F.t <= 9000)),
   (0, True)
)
iflux = sp.simplify(ion_flux)
print(iflux)

# Hydrogen implantation source term
# imp_depth: implantation depth (m)
# width: standard deviation of Gaussian implantation profile (m)
source_term = F.ImplantationFlux(
   flux=ion_flux, # H/m2/s
   imp_depth=4e-9, # m
   width=2e-9, # m
   volume=1
)
my_model.sources = [source_term]

# Tungsten atomic density (atoms/m3)
w_atom_density = 6.3e28 # atom/m3

# Dislocation trap properties
dislocation_trap = F.Trap(
   k_0=2.9e12 / w_atom_density, # Trap formation attempt frequency (s-1)
   E_k=0.29, # Trap formation activation energy (eV)
   p_0=1e13, # Detrap attempt frequency (s-1)
   E_p=1.1, # Detrap activation energy (eV)
   density=4e-3 * w_atom_density, # Trap density (m-3)
   materials=tungsten
)

# Vacancy trap properties
vacancy_trap = F.Trap(
   k_0=2.9e12 / w_atom_density,
   E_k=0.29,
   p_0=1e13,
   E_p=1.5,
   density=2e-4 * w_atom_density,
   materials=tungsten
)
my_model.traps = [dislocation_trap, vacancy_trap]

# Dirichlet boundary conditions (zero concentration at both surfaces)
my_model.boundary_conditions = [
   F.DirichletBC(surfaces=[1, 2], value=0, field=0)
]

# Time-dependent temperature profile for LTD experiment
# Simulates cooling steps during implantation
implantation_temp = sp.Piecewise(
   (570, F.t <= 1800),
   (570 - (570 - 520) / (1850 - 1800) * (F.t - 1800), (F.t > 1800) & (F.t <= 1850)),
   (520, (F.t > 1850) & (F.t <= 3600)),
   (520 - (520 - 470) / (3650 - 3600) * (F.t - 3600), (F.t > 3600) & (F.t <= 3650)),
   (470, (F.t > 3650) & (F.t <= 5400)),
   (470 - (470 - 420) / (5450 - 5400) * (F.t - 5400), (F.t > 5400) & (F.t <= 5450)),
   (420, (F.t > 5450) & (F.t <= 7200)),
   (420 - (420 - 370) / (7250 - 7200) * (F.t - 7200), (F.t > 7200) & (F.t <= 7250)),
   (370, (F.t > 7250) & (F.t <= 9000)),
   (300, True)
)
it = sp.simplify(implantation_temp)
print(it)

# TDS heating rate (K/s)
temperature_ramp = 1 # K/s
# Start time for TDS phase (after implantation)
start_tds = implantation_time + 50 # s

# Create dense time points list for milestones - especially around critical points
temp_change_points = [1800, 1850, 3600, 3650, 5400, 5450, 7200, 7250, 9000, 200]
time_points = list(range(0, 5, 1)) # First 5 seconds
time_points += list(range(5, 20, 2)) # 5-20 seconds
time_points += list(range(20, 100, 5)) # 20-100 seconds
time_points += list(range(100, 1000, 50)) # 100-1000 seconds
time_points += list(range(1000, implantation_time + 1, 200)) # 1000-9000 seconds
time_points += [implantation_time, start_tds]
for point in temp_change_points:
   time_points.extend([point - 10, point - 5, point, point + 5, point + 10])

# Temperature profile definition
# Implantation phase: time-dependent temperature
# Waiting phase: constant 300K
# TDS phase: linear heating at 1K/s
my_model.T = sp.Piecewise(
   (implantation_temp, F.t < 9000),
   (300, (F.t > 9000) & (F.t < start_tds)),
   (300 + temperature_ramp * (F.t - start_tds), True)
)

# Adaptive time stepping strategy
# Smaller steps at critical time points for numerical stability
my_model.dt = F.Stepsize(
   initial_value=0.01,
   stepsize_change_ratio=1,
   max_stepsize=lambda t: 0.001 if t < 10 else (0.01 if t < 100 else (0.1 if t < 1000 else 1)),
   dt_min=1e-10,
   milestones=sorted(set(time_points)),
)

# Simulation settings
my_model.settings = F.Settings(
   absolute_tolerance=1e22,
   relative_tolerance=1e-5,
   final_time=15000,
)

# Define derived quantities to export
list_of_derived_quantities = [
   F.TotalVolume("solute", volume=1),
   F.TotalVolume("retention", volume=1),
   F.TotalVolume("1", volume=1), # Dislocation traps
   F.TotalVolume("2", volume=1), # Vacancy traps
   F.HydrogenFlux(surface=1),
   F.HydrogenFlux(surface=2)
]

# Add point monitoring at specific locations (5um, 10um, 1um)
point_5um_trap1 = F.PointValue(field="1", x=5e-6) # Dislocation traps
point_5um_trap2 = F.PointValue(field="2", x=5e-6) # Vacancy traps
point_10um_trap1 = F.PointValue(field="1", x=10e-6)
point_10um_trap2 = F.PointValue(field="2", x=10e-6)
point_1um_trap1 = F.PointValue(field="1", x=1e-6) # 1 micron
point_1um_trap2 = F.PointValue(field="2", x=1e-6)

# Add point monitoring to export list
list_of_derived_quantities.append(point_5um_trap1)
list_of_derived_quantities.append(point_5um_trap2)
list_of_derived_quantities.append(point_10um_trap1)
list_of_derived_quantities.append(point_10um_trap2)
list_of_derived_quantities.append(point_1um_trap1)
list_of_derived_quantities.append(point_1um_trap2)

# Create exporter for derived quantities
derived_quantities = F.DerivedQuantities(
   list_of_derived_quantities,
   filename=os.path.join(output_dir, "LTD-temp-data.csv"),
   show_units=True,
)
my_model.exports = [derived_quantities]

# Initialize and run the simulation
my_model.initialise()
my_model.run()

# Post-processing: extract and save required data
import pandas as pd

# Read all data
temp_data_path = os.path.join(output_dir, "LTD-temp-data.csv")
df = pd.read_csv(temp_data_path)

# Extract data at 5um, 10um, and 1um
point_columns = [col for col in df.columns if "PointValue" in col]

if len(point_columns) >= 6:
   # Total trapped hydrogen concentration at 5um
   df_5um = df[['t (s)', point_columns[0], point_columns[1]]].copy()
   df_5um['Total trapped'] = df_5um[point_columns[0]] + df_5um[point_columns[1]]
   df_5um = df_5um[['t (s)', 'Total trapped']]
   df_5um.columns = ['Time (s)', 'Concentration (m^-3)']
   
   # Total trapped hydrogen concentration at 10um
   df_10um = df[['t (s)', point_columns[2], point_columns[3]]].copy()
   df_10um['Total trapped'] = df_10um[point_columns[2]] + df_10um[point_columns[3]]
   df_10um = df_10um[['t (s)', 'Total trapped']]
   df_10um.columns = ['Time (s)', 'Concentration (m^-3)']
   
   # Total trapped hydrogen concentration at 1um
   df_1um = df[['t (s)', point_columns[4], point_columns[5]]].copy()
   df_1um['Total trapped'] = df_1um[point_columns[4]] + df_1um[point_columns[5]]
   df_1um = df_1um[['t (s)', 'Total trapped']]
   df_1um.columns = ['Time (s)', 'Concentration (m^-3)']
   
   # Keep only implantation phase data (t <= implantation_time) with 500s interval
   df_5um = df_5um[df_5um['Time (s)'] <= implantation_time]
   df_5um = df_5um[df_5um['Time (s)'] % 500 == 0]
   df_10um = df_10um[df_10um['Time (s)'] <= implantation_time]
   df_10um = df_10um[df_10um['Time (s)'] % 500 == 0]
   df_1um = df_1um[df_1um['Time (s)'] <= implantation_time]
   df_1um = df_1um[df_1um['Time (s)'] % 500 == 0]
   
   # Save data to CSV files
   df_5um.to_csv(os.path.join(output_dir, "LTD-5um.csv"), index=False)
   df_10um.to_csv(os.path.join(output_dir, "LTD-10um.csv"), index=False)
   df_1um.to_csv(os.path.join(output_dir, "LTD-1um.csv"), index=False)
   print(f"Successfully saved 5um, 10um and 1um data to {output_dir} directory")
else:
   print("Warning: Not enough PointValue columns found to extract data")

# Generate TDS curve (desorption phase only)
t = df['t (s)'].values
flux_left = df[[col for col in df.columns if "surface 1" in col and "solute" in col]].values.flatten()
flux_right = df[[col for col in df.columns if "surface 2" in col and "solute" in col]].values.flatten()
flux_total = -flux_left - flux_right

# Only take desorption phase data (t >= implantation_time)
t_desorption = t[t >= implantation_time]
flux_total_desorption = flux_total[t >= implantation_time]

import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.plot(t_desorption, flux_total_desorption, linewidth=3)
plt.ylabel(r"Desorption flux (m$^{-2}$ s$^{-1}$)")
plt.xlabel(r"Time (s)")
plt.title("TDS Curve")
plt.grid(True)
tds_path = os.path.join(output_dir, 'LTD-tds.png')
plt.savefig(tds_path)
plt.close()
print(f"TDS curve successfully saved to {tds_path}")

```

---

<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:** [September 16, 2025, 3:37pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/2 "2025-09-16T15:37:05Z")

</div>

Hi @Trinity thanks for posting here! Can you please format your code with triple quotes ```?

```python
print('this is my python code)

```

---

<div class="post-metadata">

**Author:** ![Trinity](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Trinity](https://festim.discourse.group/u/Trinity)\
**Post date:** [September 16, 2025, 11:30pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/3 "2025-09-16T23:30:56Z")

</div>

Apologize for didn’t format the code,and now it could be fine,thank you!

---

<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:** [September 17, 2025, 12:18pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/4 "2025-09-17T12:18:57Z")

</div>

@Trinity I ran your code and could indeed see that the solution was zero everywhere. As explained in the [troubleshooting section of the documentation](https://festim.readthedocs.io/en/stable/userguide/troubleshooting.html#solution-is-zero-everywhere), this has to do with a high absolute tolerance. Here you set 1E22 which is very high.

I reduced it to 1E15 and it seems to be solving properly (still running though)

---

<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:** [September 17, 2025, 1:09pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/5 "2025-09-17T13:09:39Z")

</div>

![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/fb70ad5a4279b9e645ab15d2b45e8d9b01fc83a9.png)  
I ran the simulation up until about 15% of the total time and when plotting the fields (using `XDMFExport` I noticed some wild oscillations suggesting a mesh issue to investigate

---

<div class="post-metadata">

**Author:** ![Trinity](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Trinity](https://festim.discourse.group/u/Trinity)\
**Post date:** [September 18, 2025, 12:47am UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/6 "2025-09-18T00:47:41Z")

</div>

Got it！I will modify my code according your guidance. Thank you for your help and your time！

---

<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:** [September 18, 2025, 12:08pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/7 "2025-09-18T12:08:30Z")

</div>

Yesterday I found that this mesh was good enough to go through the whole implantation phase:

```python
surface_region_1 = np.linspace(0, 10e-9, num=30)

surface_region_2 = np.linspace(10e-9, 50e-9, num=100)

surface_region_4 = np.linspace(50e-9, 10e-5, num=7000)

far_region = np.linspace(10e-5, 3e-3, num=7000)
vertices = np.concatenate(
    [
        surface_region_1,
        surface_region_2,
        surface_region_4,
        far_region,
    ]
)

```

There is an issue however during hte resting phase when the temperature drops to room temperature, gradients become even stiffer and oscillations appear.

I “fixed” this by setting the resting temperature to 370 K instead of 300.

Maybe another thing to consider is [using discontinuous elements for traps](https://festim.readthedocs.io/en/stable/api/settings.html#settings-api:~:text=traps_element_type%20(str%2C%20optional)%20%E2%80%93%20Finite%20element%20used%20for%20traps.%20If%20traps%20densities%20are%20discontinuous%20(eg.%20different%20materials)%20%E2%80%9CDG%E2%80%9D%20is%20recommended.%20Defaults%20to%20%E2%80%9CCG%E2%80%9D):

---

<div class="post-metadata">

**Author:** ![Trinity](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Trinity](https://festim.discourse.group/u/Trinity)\
**Post date:** [September 20, 2025, 4:05pm UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/8 "2025-09-20T16:05:03Z")

</div>

So sorry I missed your message earlier—my inbox must’ve flagged it as spam. 😅 I’ll **immediately** try adjusting the mesh based on your suggestions and keep testing!

Additionally,since my current work “doesn’t actually need” to focus on the thermal desorption process, I’m thinking of trimming off the last 6000s in the code _just to save time_ —and only keep tracking the hydrogen concentration. Would that make sense?

_(Thanks again for your patience—I really appreciate it!)_ 🙏

---

<div class="post-metadata">

**Author:** ![Trinity](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Trinity](https://festim.discourse.group/u/Trinity)\
**Post date:** [September 21, 2025, 9:47am UTC](https://festim.discourse.group/t/time-step-reaching-the-minimum-value-during-implantation/142/9 "2025-09-21T09:47:58Z")

</div>

Good new! I followed your suggestion, using DG trap\_element\_type,and the code ran smoothly this time, thanks again for your help!