Simulation stopping on simple permeation experiment

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.

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%.

Hey there Ignacio!

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


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:

Hi!

Try setting the log_level to 20 with

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 of the documentation for details on the Stepsize parameters.

Here’s the complete code:

# 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()

@jhdark maybe the log_level should default to 30 instead of 40?

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

1 Like

I think we broke the Olympic Record for response time. :medal_sports:

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

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

If your solubility is in

s_0 = 400  # mol/m3*MPa^0.5

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

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

This is why I recommend using HTM for handling properties as it takes care of unit conversion.

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

        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.

Also your mm to m conversion is wrong

1 mm = 1E-3 m not

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

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

I also overlapped the analytical steady state permeation flux

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()