# Simulation stopping on simple permeation experiment

**URL:** <https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69>\
**Category:** User Support\
**Created:** [August 1, 2024, 1:12pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69 "2024-08-01T13:12:40Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ignacio-Marcos](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/ignacio-marcos/32/72_2.png) [@Ignacio-Marcos](https://festim.discourse.group/u/Ignacio-Marcos)\
**Post date:** [August 1, 2024, 1:12pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/1 "2024-08-01T13:12:40Z")

</div>

Hi everyone,

I’m doing a simple series of permeation experiments with a set of different values for each parameter (diffusivity, Temperature and such), specifically, a full factorial DOE. I’m having trouble with the simulation stopping without an error message when I increase the time or just when running the script, such as with the case that I am presenting.

The script is based on task 3, with a little modification to make it suitable for the full DOE. I don’t have much background with Finite Element Analysis or solvers, but I believe the issue might be there. Hope someone can point me in the right direction. I’ll paste the code below.

```auto
import festim as F
import numpy as np

Temp=298 #K
Pres=1 #bar
d_0=1e-9 #m2/s
e_D=5 #KJ/mol
s_0=10 #mol/m3*MPa^0.5
e_S=5 #KJ/mol
Thick=0.1 #mm

Presc=Pres*1e5 #Pa
e_Dc=e_D*0.010364 #eV
s_0c=s_0*6.022e23 #H/m3*0.5Pa
e_Sc=e_S*0.010364 #eV
Thickc=0.1*1e-4 

my_model = F.Simulation()
my_model.mesh = F.MeshFromVertices(
    vertices=np.linspace(0, Thickc , num=1001)
)

my_model.materials = F.Material(
    id=1, 
    D_0=d_0, #m2/s
    E_D=e_Dc #eV
    )

my_model.T = F.Temperature(
    value=Temp
    )

my_model.boundary_conditions = [
    F.SievertsBC(
        surfaces=1, 
        S_0=s_0c,
        E_S=e_Sc,
        pressure=Presc
        ),
    F.DirichletBC(
        surfaces=2, 
        value=0, 
        field=0
        )
]

my_model.settings = F.Settings(
    absolute_tolerance=1e-1,
    relative_tolerance=1e-1,
    final_time=500 # s
    )

my_model.dt = F.Stepsize(
    initial_value=1/20,
     dt_min=1e-2
     )

derived_quantities = F.DerivedQuantities([F.HydrogenFlux(surface=2)])

my_model.exports = [derived_quantities]

my_model.initialise()

my_model.run()

times = derived_quantities.t
computed_flux = derived_quantities.filter(surfaces=2).data

import matplotlib.pyplot as plt
plt.scatter(times, np.abs(computed_flux), alpha=0.1, label="computed",linewidths=0.1)
plt.ylim(bottom=0)
plt.xlabel("Time (s)")
plt.ylabel("Downstream flux (H/m2/s)")
# Add a footnote below and to the right side of the chart
plt.figtext(
    0.6, 0.2, 
    f"Variables:\n- Temperature: {Temp} ºK\n- Pressure: {Pres} bar, {Presc} Pa\n- D_0: {d_0} m2/s \n- E_D: {e_D} KJ/mol, {e_Dc} eV\n- S_0: {s_0} mol/m3*MPa^0.5, {s_0c} H/m3*Pa^0.5\n- E_S: {e_S} KJ/mol, E_S: {e_Sc} eV", 
    ha="left", 
    fontsize=10, 
)
plt.show()

```

Thanks to everyone for reading and let me know if I can give any more information. The error I get, if you call it that, is the following, the simulation just stops at a certain number, in this case 0.5%.

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

---

<div class="post-metadata">

**Author:** ![jhdark](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/jhdark/32/13_2.png) [@jhdark](https://festim.discourse.group/u/jhdark)\
**Post date:** [August 1, 2024, 1:12pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/2 "2024-08-01T13:12:48Z")

</div>

Hey there Ignacio!

I think the issue is regarding the absolute tolerance, if you increase it up to 1e10 it converges much better:

```auto

my_model.settings = F.Settings(

absolute_tolerance=1e10, relative_tolerance=1e-10, final_time=500 # s

)

my_model.dt = F.Stepsize(initial_value=1e-3, stepsize_change_ratio=1.05)

```

Plus seems to converge within a second:

 ![image](https://global.discourse-cdn.com/free1/uploads/festim/original/1X/e2e391d85f12f3e0eb0b75d61892713cac45e5fe.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:** [August 1, 2024, 1:12pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/3 "2024-08-01T13:12:49Z")

</div>

Hi!

Try [setting the `log_level` to 20](https://festim.readthedocs.io/en/stable/userguide/troubleshooting.html#solver-doesn-t-converge) with

```auto
my_model.log_level = 20 # 20 = INFO

```

You will see that the solver struggles to converge.

This is a tolerance issue. If you increase the absolute tolerance to 1E10 and set the relative tolerance to 1E-10, then the solver converges: it converges in “zero iteration” meaning it reached steady state.

If you set the final time to 1 second and the initial timestep to 1E-3 s with a `stepsize_change_ratio` of 1.1 then you obtain the following graph.  
See the [this section](https://festim.readthedocs.io/en/stable/userguide/stepsize.html#stepsize) of the documentation for details on the `Stepsize` parameters.

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

Here’s the complete code:

```python
# this is the script on txt. just paste it on a python file.

import festim as F
import numpy as np

Temp = 298 # K
Pres = 1 # bar
d_0 = 1e-9 # m2/s
e_D = 5 # KJ/mol
s_0 = 10 # mol/m3*MPa^0.5
e_S = 5 # KJ/mol
Thick = 0.1 # mm

my_model = F.Simulation()
my_model.mesh = F.MeshFromVertices(vertices=np.linspace(0, Thick * 1e-4, num=1001))

my_model.materials = F.Material(
    id=1,
    D_0=d_0, # m2/s
    E_D=e_D * 0.010364, # eV
)

my_model.T = F.Temperature(value=Temp)

my_model.boundary_conditions = [
    F.SievertsBC(
        surfaces=1,
        S_0=s_0 * 6.022e23, # H/m3*0.5Pa
        E_S=e_S * 0.010364, # eV
        pressure=Pres * 1e5,
    ),
    F.DirichletBC(surfaces=2, value=0, field=0),
]

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

my_model.dt = F.Stepsize(initial_value=1e-3, stepsize_change_ratio=1.1)

derived_quantities = F.DerivedQuantities([F.HydrogenFlux(surface=2)])

my_model.exports = [derived_quantities]

# my_model.log_level = 20
my_model.initialise()

my_model.run()

times = derived_quantities.t
computed_flux = derived_quantities.filter(surfaces=2).data

import matplotlib.pyplot as plt

plt.scatter(times, np.abs(computed_flux), alpha=0.8, label="computed", linewidths=0.1)
plt.ylim(bottom=0)
plt.xlabel("Time (s)")
plt.ylabel("Downstream flux (H/m2/s)")
# Add a footnote below and to the right side of the chart
plt.figtext(
    0.6,
    0.2,
    f"Variables:\n- Temperature: {Temp} ºK\n- Pressure: {Pres} bar\n- D_0: {d_0} m2/s\n- E_D: {e_D} KJ/mol\n- S_0: {s_0} mol/m3*MPa^0.5\n- E_S: {e_S} KJ/mol",
    ha="left",
    fontsize=10,
)
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:** [August 1, 2024, 1:18pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/4 "2024-08-01T13:18:47Z")

</div>

@jhdark maybe the log\_level should default to 30 instead of 40?

---

<div class="post-metadata">

**Author:** ![Ignacio-Marcos](https://yyz2.discourse-cdn.com/free1/user_avatar/festim.discourse.group/ignacio-marcos/32/72_2.png) [@Ignacio-Marcos](https://festim.discourse.group/u/Ignacio-Marcos)\
**Post date:** [August 1, 2024, 1:26pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/5 "2024-08-01T13:26:18Z")

</div>

Hi, thanks for the amazingly quick response.

This does solve the issue with this simulation, but I’m having trouble with some others combination of values. particullarly on the next value, I run into an stepsize problem. why might this be?  
Temp Pres d\_0 e\_D s\_0 e\_S Thick  
77 1 0,0001 5 400 5 1

---

<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 1, 2024, 1:28pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/6 "2024-08-01T13:28:56Z")

</div>

> [@Ignacio-Marcos](#):
>
> amazingly quick response.

I think we broke the Olympic Record for response time. 🏅

Can you format your code please? as well as providing exactly what kind of “stepsize problem” you run into

---

<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 1, 2024, 1:39pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/7 "2024-08-01T13:39:08Z")

</div>

Btw @Ignacio-Marcos you’ve got a conversion error in your script.

If your solubility is in

```python
s_0 = 400 # mol/m3*MPa^0.5

```

Then simply multiplying by Avogadro converts it to H/m3_0.5MPa and not H/m3_0.5Pa.

```python
        S_0=s_0 * 6.022e23, # H/m3*0.5Pa

```

This is why I recommend using [HTM](https://h-transport-materials.readthedocs.io/en/latest/user/defining_property.html) for handling properties as it takes care of unit conversion.

```python
import h_transport_materials as htm

solubility = htm.Solubility(
    S_0=400 * htm.ureg.mol / htm.ureg.m **3 * htm.ureg.MPa** -0.5,
    E_S=5 * htm.ureg.kJ / htm.ureg.mol,
)

print(solubility)

print(s_0 * 6.022e23)

```

Produces

```bash
        Author: 
        Material: 
        Year: None
        Isotope: None
        Pre-exponential factor: 2.41×10²³ particle/m³/Pa⁰⋅⁵
        Activation energy: 5.18×10⁻² eV/particle
        
2.4088e+26

```

You see you are overestimating your solubility by three orders of magnitude.

---

<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 1, 2024, 1:41pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/8 "2024-08-01T13:41:23Z")

</div>

Also your mm to m conversion is wrong

1 mm = 1E-3 m not

```auto
Thick = 1 # mm
my_model.mesh = F.MeshFromVertices(vertices=np.linspace(0, Thick * 1e-4, num=1001))

```

---

<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 1, 2024, 1:46pm UTC](https://festim.discourse.group/t/simulation-stopping-on-simple-permeation-experiment/69/9 "2024-08-01T13:46:41Z")

</div>

@Ignacio-Marcos here’s a working script with the correct thickness in m.  
(haven’t fixed the other conversion mistakes though)

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

I also overlapped the analytical steady state permeation flux

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

Temp = 77 # K
Pres = 1 # bar
d_0 = 0.0001 # m2/s
e_D = 5 # KJ/mol
s_0 = 400 # mol/m3*MPa^0.5
e_S = 5 # KJ/mol
Thick = 1 # mm

my_model = F.Simulation()
my_model.mesh = F.MeshFromVertices(vertices=np.linspace(0, Thick * 1e-3, num=1001))

my_model.materials = F.Material(
    id=1,
    D_0=d_0, # m2/s
    E_D=e_D * 0.010364, # eV
)

my_model.T = F.Temperature(value=Temp)

my_model.boundary_conditions = [
    F.SievertsBC(
        surfaces=1,
        S_0=s_0 * 6.022e23, # H/m3*0.5Pa
        E_S=e_S * 0.010364, # eV
        pressure=Pres * 1e5,
    ),
    F.DirichletBC(surfaces=2, value=0, field=0),
]

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

my_model.dt = F.Stepsize(initial_value=1e-3, stepsize_change_ratio=1.1, dt_min=1e-5)

derived_quantities = F.DerivedQuantities([F.HydrogenFlux(surface=2)])

my_model.exports = [derived_quantities]

my_model.log_level = 20
my_model.initialise()

my_model.run()

times = derived_quantities.t
computed_flux = derived_quantities.filter(surfaces=2).data

import matplotlib.pyplot as plt

plt.scatter(times, np.abs(computed_flux), alpha=0.8, label="computed", linewidths=0.1)
plt.ylim(bottom=0)
plt.xlabel("Time (s)")
plt.ylabel("Downstream flux (H/m2/s)")
# Add a footnote below and to the right side of the chart
plt.figtext(
    0.6,
    0.2,
    f"Variables:\n- Temperature: {Temp} ºK\n- Pressure: {Pres} bar\n- D_0: {d_0} m2/s\n- E_D: {e_D} KJ/mol\n- S_0: {s_0} mol/m3*MPa^0.5\n- E_S: {e_S} KJ/mol",
    ha="left",
    fontsize=10,
)

# compute analytical steady state permeation flux
S = my_model.boundary_conditions[0].S_0 * np.exp(
    -my_model.boundary_conditions[0].E_S / F.k_B / Temp
)
c_0 = S * my_model.boundary_conditions[0].pressure ** 0.5

D = my_model.materials[0].D_0 * np.exp(-my_model.materials[0].E_D / F.k_B / Temp)

analytical_flux = c_0 * D / (Thick * 1e-3)
plt.axhline(y=analytical_flux, color="r", linestyle="--", label="analytical")
print(f"Analytical flux: {analytical_flux:.2e} H/m2/s")
plt.show()

```