# 2d simulation breaks when introducing traps

**URL:** <https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148>\
**Category:** User Support\
**Created:** [October 21, 2025, 2:03pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148 "2025-10-21T14:03:46Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![miklavrent](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/miklavrent/32/94_2.png) [@miklavrent](https://festim.discourse.group/u/miklavrent)\
**Post date:** [October 21, 2025, 2:03pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/1 "2025-10-21T14:03:46Z")

</div>

Hi,

I am trying to model 2d steel square inside another piece of steel, with traps in the inner piece of steel only. Initially, hydrogen is in inner steel. The code below works ok, but as soon as I uncomment the trap lines, something strange happens: at time = 0, nonzero concentration of solute is at the interface between steels only!

Could you please explain what is wrong with my code?

Best wishes,

Mikhail

> from fenics import \*
> 
> import matplotlib.pyplot as plt
> 
> import festim as F
> 
> from festim import x, y, z
> 
> import numpy as np
> 
> import sympy
> 
> import math
> 
> # GEOMETRY AND MESH
> 
> width = 1
> 
> height = 1
> 
> border\_thickness = 0.1
> 
> nx, ny = 21, 21
> 
> mesh = RectangleMesh(Point(0, 0), Point(width, height), nx, ny)
> 
> # MARK VOLUMES
> 
> volume\_markers = MeshFunction(“size\_t”, mesh, mesh.topology().dim())
> 
> volume\_markers.set\_all(2)
> 
> plot(volume\_markers, title=“Volume Markers (1=steel1, 2=steel2)”)
> 
> plt.show()
> 
> inner\_metal = CompiledSubDomain(
> 
> “x[0] \> x\_min && x[0] \< x\_max && x[1] \> y\_min && x[1] \< y\_max”,
> 
> ```
> x_min=border_thickness,
> 
> x_max=width - border_thickness,
> 
> y_min=border_thickness,
> 
> y_max=height - border_thickness
> 
> ```
> 
> )
> 
> inner\_metal.mark(volume\_markers, 1)
> 
> model = F.Simulation()
> 
> model.mesh = F.Mesh(mesh=mesh, volume\_markers=volume\_markers)
> 
> # MATERIALS
> 
> steel1 = F.Material(
> 
> ```
> id = 1,
> 
> D_0 = 1e-16, #m2/s
> 
> E_D = 0, #eV
> 
> ```
> 
> )
> 
> steel2 = F.Material(
> 
> ```
> id = 2,
> 
> D_0 = 1E-15, #m2/s
> 
> E_D = 0, #eV
> 
> ```
> 
> )
> 
> model.materials = [steel1, steel2]
> 
> # TEMPERATURE
> 
> model.T = 300
> 
> # INITIAL CONDITIONS
> 
> model.initial\_conditions = [
> 
> ```
> F.InitialCondition(
> 
> field="solute",
> 
> value=1e26 \* sympy.Piecewise(
> 
> (1, (x > border_thickness) & (x < width - border_thickness) &
> 
> (y > border_thickness) & (y < height - border_thickness)),
> 
> (0, True)
> 
> )
> 
> )
> 
> ```
> 
> ]
> 
> # SURFACE MARKERS
> 
> surface\_markers = MeshFunction(“size\_t”, mesh, mesh.topology().dim() - 1)
> 
> surface\_markers.set\_all(0)
> 
> steel1\_boundary = CompiledSubDomain(
> 
> “near(x[0], x\_min) || near(x[0], x\_max) || near(x[1], y\_min) || near(x[1], y\_max)”,
> 
> ```
> x_min=border_thickness,
> 
> x_max=width - border_thickness,
> 
> y_min=border_thickness,
> 
> y_max=height - border_thickness
> 
> ```
> 
> )
> 
> steel2\_boundary = CompiledSubDomain(
> 
> “on\_boundary”,
> 
> )
> 
> steel1\_boundary.mark(surface\_markers, 1)
> 
> steel2\_boundary.mark(surface\_markers, 2)
> 
> model.mesh.surface\_markers = surface\_markers
> 
> SS316L\_density = 8.55e28
> 
> #TRAPS
> 
> #SS316L\_trap = F.Trap(
> 
> # k\_0 = 3.82e-7 / ((5.66e-10)\*\*2 \* 6 \* SS316L\_density),
> 
> # E\_k = 0.16,
> 
> # p\_0 = 1e13,
> 
> # E\_p = 0.5,
> 
> # density = 6.08e25,
> 
> # materials = [steel1]
> 
> #)
> 
> #model.traps = [SS316L\_trap]
> 
> # BOUNDARY CONDITIONS
> 
> model.boundary\_conditions = [
> 
> ```
> F.DirichletBC(field = 0,
> 
> value = 0,
> 
> surfaces = \[2\]),
> 
> ```
> 
> ]
> 
> #EXPORTS
> 
> results\_folder = “TWO\_DIMENSION\_PERMEATION”
> 
> model.exports = [
> 
> ```
> F.XDMFExport(
> 
> field="solute",
> 
> filename=results_folder + "/hydrogen_concentration.xdmf",
> 
> checkpoint=False,
> 
> ),
> 
> ```
> 
> ]
> 
> # SETTINGS
> 
> model.settings = F.Settings(
> 
> ```
> absolute_tolerance = 1e16, 
> 
> relative_tolerance = 1e-8,
> 
> final_time= 1e3, 
> 
> )
> 
> ```
> 
> # TIMESTEP
> 
> model.dt = F.Stepsize(
> 
> ```
> initial_value=1,
> 
> stepsize_change_ratio=1.2, 
> 
> dt_min=1e-5,
> 
> max_stepsize=1e3,
> 
> ```
> 
> )
> 
> model.initialise()
> 
> model.run()

---

<div class="post-metadata">

**Author:** ![miklavrent](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/miklavrent/32/94_2.png) [@miklavrent](https://festim.discourse.group/u/miklavrent)\
**Post date:** [October 21, 2025, 2:06pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/2 "2025-10-21T14:06:50Z")

</div>

P.S. This is FESTIM1

---

<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:** [October 21, 2025, 2:28pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/3 "2025-10-21T14:28:40Z")

</div>

Hi @miklavrent could you provide some screenshots to illustrate your issue? The code block isn’t formatted properly so it’s not easy to copy paste

---

<div class="post-metadata">

**Author:** ![miklavrent](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/miklavrent/32/94_2.png) [@miklavrent](https://festim.discourse.group/u/miklavrent)\
**Post date:** [October 21, 2025, 2:46pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/4 "2025-10-21T14:46:55Z")

</div>

Hi Remi,

Thank you for your reply! It looks like there are several unexpected backslashes in the code, perhaps it’ll work if you delete them? Anyway, I attach two screenshots showing initial mobile distribution with and without the trap part.

Best wishes,

Mikhail

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

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

---

<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:** [October 21, 2025, 2:52pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/5 "2025-10-21T14:52:56Z")

</div>

If you look at the colour bar, the first screenshot shows a negative value in the inner square. The red strip is an “overshoot”. This seems like a typical Gibb’s phenomenon issue related to the mesh.

Basically you are trying to approximate very sharp gradients with too few cells. Things to consider:

- refining the mesh
- use discontinuous finite elements (DG) for the trap concentration

To get an idea of the refinement level you need, you may want to setup a 1D simulation an play with the refinement there. By the way, are you sure your problem isn’t a 1D problem?

---

<div class="post-metadata">

**Author:** ![miklavrent](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/miklavrent/32/94_2.png) [@miklavrent](https://festim.discourse.group/u/miklavrent)\
**Post date:** [October 21, 2025, 2:56pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/6 "2025-10-21T14:56:40Z")

</div>

Thank you Remi, I’ll try to refine the mesh and/or use DG. Indeed this looks like a 1D problem, but it’s because I simplified it to bring as close to MWE as possible.

---

<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:** [April 13, 2026, 12:52pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/7 "2026-04-13T12:52:50Z")

</div>

@miklavrent what’s the status of this? can we close this?

---

<div class="post-metadata">

**Author:** ![miklavrent](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/miklavrent/32/94_2.png) [@miklavrent](https://festim.discourse.group/u/miklavrent)\
**Post date:** [April 13, 2026, 1:12pm UTC](https://festim.discourse.group/t/2d-simulation-breaks-when-introducing-traps/148/8 "2026-04-13T13:12:58Z")

</div>

Hi Remi, yes please close! I remember resolving this problem by introducing finer mesh in a border region between two pieces of steel.