Issue with 0 iteration solving

User Support
1

Hello,

I’m reproducing an experiment with FESTIM 2.0 alpha where I expose a tungstene material to D fluxes, with a fixed temperature between the two.
At the third phase the calculation has an issue, probably solving in 0 iterations each step until the end.
I tried to implement milestones in the stepsize objet but I may be missing something.

Thanks in advance for your help !

import festim as F
import ufl
import numpy as np
import numpy.typing as npt
import matplotlib.pyplot as plt

model = F.HydrogenTransportProblem()
L=10e-5

vertice = np.concatenate(
[
    np.linspace(0, 100e-9, num=100, endpoint=False),
    np.linspace(100e-9, 10e-6, num=200, endpoint=False),
    np.linspace(10e-6, L, num=200),
]
)

vertice = np.asarray(vertice)
model.mesh = F.Mesh1D(vertice)

difW = 1.9e-7
E_diff = 0.2
tungst = F.Material(name="W", D_0=difW, E_D=E_diff)
vol = F.VolumeSubdomain1D(id=1, borders=[0,L], material=tungst)
surf_left = F.SurfaceSubdomain1D(id=1, x=0)
surf_right = F.SurfaceSubdomain1D(id=2, x=L)

W_atom_density = 19.25e3 / (183.84 * 1.66e-27)   

model.subdomains= [vol, surf_left, surf_right]

mobile_D = F.Species("D")

trapAD = F.Species("trapAD", mobile = False)

trapA = F.ImplicitSpecies(n=0.0019*W_atom_density, others=[trapAD], name="trapA")   

model.species = [mobile_D, trapAD]#, trapAH, trapBH, trapBD, trapCH, trapCD, mobile_H]

model.reactions = [
    F.Reaction( 
        k_0 = 1e13/W_atom_density,
        E_k = E_diff,
        p_0 = 1e13,
        E_p = 1.65,
        reactant = [mobile_D, trapA],
        product = trapAD,
        volume = vol,
    ),
    
]

D_expo = 48*3600
t_TDS = D_expo + 48*3600
D2_expo = t_TDS + 23.5*3600

def implant_surf_D(t):
    if t < D_expo or t > t_TDS: 
        return 1.41e18
    else:
        return 0

model.boundary_conditions = [
    F.DirichletBC(subdomain=surf_left, value = implant_surf_D, species= mobile_D),
    F.DirichletBC(subdomain=surf_right, value = 0, species=mobile_D),
]

model.temperature = 600

invent_mobD = F.TotalVolume(field=mobile_D, volume=vol)
invent_trapAD = F.TotalVolume(field=trapAD, volume=vol)

model.exports=[
    F.XDMFExport("trap_1_D.xdmf", field=trapAD),
    invent_mobD,
    invent_trapAD,
]

mlstns = [D_expo - 5, D_expo, D_expo + 5, D_expo + 10, D_expo + 15, D_expo + 20, t_TDS,
          t_TDS + 5, t_TDS + 10, t_TDS + 15, t_TDS + 20, t_TDS + 25, t_TDS + 30, t_TDS + 35, t_TDS + 40, t_TDS + 45,
          D2_expo - 45, D2_expo - 40, D2_expo - 35, D2_expo - 30, D2_expo - 25, D2_expo - 20, D2_expo - 15, D2_expo - 10, D2_expo - 5]

model.settings = F.Settings(atol= 1e10, rtol= 1e-10, max_iterations=30, final_time=D2_expo)
model.settings.stepsize = F.Stepsize(initial_value=0.1,growth_factor=1.05,target_nb_iterations=20, cutback_factor=0.95, max_stepsize=100)
model.settings.stepsize.milestones = mlstns


model.initialise()
model.run()

Flat_phase
Flat_phase1680×839 36.2 KB

2

Hey @JonathanGSDUFOUR

Can you try reducing the atol so say 1e8?

Also:

You don’t need this line

3

Thank you for the reply !

Sadly changing atol to 1e8 didn’t solve the issue.
I found that putting 1.45 eV for Ep allowed the simulation to run without issue, but in the end the goal was to add a 2nd and 3rd trap with higher ‘Edt’ than the first at 1.65.
Hope it can lead to a solution.

4

Can you:

  • turn on log level INFO with
import dolfinx
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
  • display the logs of the newton iterations when it’s plateauing like this
5

I have this kinfd of log when it is plateauing :

[2025-03-21 16:12:47.781] [info] Newton iteration 0: r (abs) = 1.4806353513761235e+17 (tol = 100000000), r (rel) = 7.743481576167682e-11 (tol = 1e-10)
[2025-03-21 16:12:47.781] [info] Newton solver finished in 0 iterations and 0 linear solver iterations.
[2025-03-21 16:12:47.781] [info] Adding function to node "/Xdmf/Domain/Grid/Grid"
[2025-03-21 16:12:47.823] [info] Newton iteration 0: r (abs) = 1.4806353513761235e+17 (tol = 100000000), r (rel) = 7.743481576167682e-11 (tol = 1e-10)
[2025-03-21 16:12:47.823] [info] Newton solver finished in 0 iterations and 0 linear solver iterations.
[2025-03-21 16:12:47.823] [info] Adding function to node "/Xdmf/Domain/Grid/Grid"
[2025-03-21 16:12:47.857] [info] Newton iteration 0: r (abs) = 1.4806353513761235e+17 (tol = 100000000), r (rel) = 7.743481576167682e-11 (tol = 1e-10)
[2025-03-21 16:12:47.857] [info] Newton solver finished in 0 iterations and 0 linear solver iterations.
[2025-03-21 16:12:47.857] [info] Adding function to node "/Xdmf/Domain/Grid/Grid"

But I realized that I could solve the issue by increasing linearly the value on the 10 first second instead of just having 0 to 1.41e18 !

6

So I see that the relative error is too high… it’s 1e-10 but on the first newton iteration the error is 1e-11.

I believe this is not a desired behaviour as one would like the first relative residual to be infinity. This has been discussed recently on the FEniCS discourse and it looks like there will be a fix in this PR.

:crystal_ball::magic_wand::woman_elf:Summoning @dokken to see if there’s an obvious fix

Otherwise @JonathanGSDUFOUR you can simply set the rtol to say 1e-11 or 1e-12

7

Not sure what you mean here

8

Thank you it worked with rtol = 1e-11 !

What I meant at the end is that I also managed to make the simulation work by putting a fast ramp (over 10 seconds) then a constant value instead of just having a step for the DirichletBC.

Now I’ll add the other phases I need, but with these two ways it should be good untill the end.

Thanks again !

1 Like