# Multi-material problem: initial trap distribution and diffusion barrier

**URL:** <https://festim.discourse.group/t/multi-material-problem-initial-trap-distribution-and-diffusion-barrier/71>\
**Category:** User Support\
**Created:** [August 9, 2024, 2:00pm UTC](https://festim.discourse.group/t/multi-material-problem-initial-trap-distribution-and-diffusion-barrier/71 "2024-08-09T14:00:14Z")\
**Posts on this page:** 4\
**Page:** 1

<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:** [August 9, 2024, 2:00pm UTC](https://festim.discourse.group/t/multi-material-problem-initial-trap-distribution-and-diffusion-barrier/71/1 "2024-08-09T14:00:14Z")

</div>

Hi, all!

Assume a multi-material problem (for example, W-Cu) with traps defined within one material only. There are two points that I’d like to discuss.

1.If the initial condition for traps is provided with a float, then the traps are initialised within the whole geometry rather than within one material. Consider the following MWE:

```auto
import festim as F
import fenics as f
import numpy as np
import matplotlib.pyplot as plt

model = F.Simulation()

vertices = np.linspace(0, 1, num=1000)

model.mesh = F.MeshFromVertices(vertices)

tungsten = F.Material(id=1, D_0=1, E_D=0, S_0=1, E_S=0, borders=[0, 0.5])
copper = F.Material(id=2, D_0=1, E_D=0, S_0=1, E_S=0, borders=[0.5, 1])

model.materials = [tungsten, copper]

trap_1 = F.Trap(
    k_0=1,
    E_k=0,
    p_0=2,
    E_p=0,
    density=1,
    materials=tungsten,
)

model.traps = [trap_1]

model.initial_conditions = [
    F.InitialCondition(field="1", value=5),
]

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

model.T = 300

model.dt = F.Stepsize(
    initial_value=0.1,
    stepsize_change_ratio=1.1,
    max_stepsize=1,
    dt_min=1,
)

model.settings = F.Settings(
    absolute_tolerance=1e10,
    relative_tolerance=1e-10,
    final_time=0,
    chemical_pot=True,
    traps_element_type="DG",
)

model.initialise()

f.plot(model.h_transport_problem.traps[0].previous_solution)

plt.ylabel("Trap concentration")
plt.xlabel("x")
plt.show()

```

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

I suppose a workaround is to use `sp.Piecewise` for the initial condition (or maybe I miss something?). If so, isn’t it confusing?

2.The second point concerns the material interface. If I suppose that hydrogen cannot diffuse from the first material into the second (e.g., there is another thin layer between these materials: W-CuO/CuO2-Cu, that prevents hydrogen migration into copper), what’d be a better option to simulate such a case?

---

<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:** [August 9, 2024, 2:28pm UTC](https://festim.discourse.group/t/multi-material-problem-initial-trap-distribution-and-diffusion-barrier/71/2 "2024-08-09T14:28:36Z")

</div>

Yes this is a current limitation. You could use Piecewise in 1D but it wouldn’t be applicable to complex 3D problems for instance.

Ideally we would have a `volume` argument to `InitialCondition`.  
Btw this would require the traps element type to be DG.

> [@VVKulagin](#):
>
> The second point concerns the material interface. If I suppose that hydrogen cannot diffuse from the first material into the second (e.g., there is another thin layer between these materials: W-CuO/CuO2-Cu, that prevents hydrogen migration into copper), what’d be a better option to simulate such a case?

You could either

- simulate this interlayer explicitely with a very low diffusivity
- if H cannot diffuse through the interface (I guess it’s similar to having a solubility approaching zero? then only simulate the W part and have a non-flux BC on the boundary

---

<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:** [August 9, 2024, 2:36pm UTC](https://festim.discourse.group/t/multi-material-problem-initial-trap-distribution-and-diffusion-barrier/71/3 "2024-08-09T14:36:02Z")

</div>

Thank you for the response.

> You could either
> 
> - simulate this interlayer explicitely with a very low diffusivity
> - if H cannot diffuse through the interface (I guess it’s similar to having a solubility approaching zero? then only simulate the W part and have a non-flux BC on the boundary

If I need both materials to simulate “realistic” temperature evolution in the composite geometry, then I end up with the first approach.

---

<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:** [August 9, 2024, 2:41pm UTC](https://festim.discourse.group/t/multi-material-problem-initial-trap-distribution-and-diffusion-barrier/71/4 "2024-08-09T14:41:22Z")

</div>

This approach isn’t really appropriate for a 2D/3D geometry since you potentially have a very small inter-layer. But ok for 1D I suppose

you could have a heat transfer transient solve on one side and then only read the temperature in the part you want. But that would require a bit of tinkering.

In the `fenicsx` branch, when we [add the interface conditions](https://github.com/festim-dev/FESTIM/issues?page=2&q=is%3Aissue+is%3Aopen) this will not be an issue since it will be possible to define a mobile species in one subdomain only.