Corrector convergence failed repeatedly or with |h| = hmin

[CaXYZ]

Hi all,

I would like to model a half-cell experiment, using the DFN halfcell model:

model_options={'working electrode':'positive'}
model=pybamm.lithium_ion.DFN(options=model_options)

The experiment reads:

experimental_prozedure=[("Charge at C/3 until 4.3 V",
                         "Hold at 4.3 V for 60 min",
                         "Discharge at 0.5 C for 120 s",
                          "Discharge at 1 C until 3.2 V",),]

and the solver options are:
solver=pybamm.CasadiSolver(mode="safe",dt_max=50,atol=1e-5,rtol=1e-5,root_tol=1e-5)

But when I run the simulation the warning
At t = 3.09526e-005 and h = 6.71932e-020, the corrector convergence failed repeatedly or with |h| = hmin
is printed several times for various (very small) combinations of t and h, and the simulation takes very long to converge.

Additionally the simulation seems to be very sensitive against small changes in the experiment, e.g., by leaving out the discharge at C/2 step, the simulation does not converge and the error message
SolverError: Maximum number of decreased steps occurred at t=7821.663133788967. Try solving the model up to this time only or reducing dt_max (currently, dt_max=1.5625).
is obtained. I was wondering why the error message states that dt_max is 1.56, as I defined dt_max to be 50 via the solver options (see above) . I tried to decrease dt_max as stated in the error massage. However, I was still not able to obtain a converged solution.

Do you have an idea what might cause this behavior or any suggestions how to best proceed the debugging process?

With best regards
Carsten

[dion-w]

Something that might help in this circumstance is to increase the discretization along the thickness of the battery and in the radial dimension. Also what this error message means: the solver tries to solve the system for the current timestep, and if it can not it will reduce
the timestep. By default the solver will only reduce the timestep 5 times to get a converging solution. You can tell the solver though to reduce the timestep even further before by setting solver = pybamm.CasadiSolver(max_step_decrease_count=15) and then giving this solver setup to the simulation object.
As an example a rather high discretization with the mentioned solver setup would look like this:

    var_pts = {
        "x_n": 50,  # negative electrode
        "x_s": 30,  # separator
        "x_p": 50,  # positive electrode
        "r_n": 150,  # negative particle, n_Shells, this is especially needed if you are using state of charge dependent diffusion coefficients
        "r_p": 150,  # positive particle, n_Shells
    }
solver = pybamm.CasadiSolver(max_step_decrease_count=15)
sim = pybamm.Simulation(model, parameter_values=param, var_pts=var_pts, solver=solver)

[CaXYZ]

Thank you for the detailed explanation and ideas! I will try this.

[CaXYZ]

Hi Dion,
I now intensively tested you suggestions, but, until now, I have not found a way to get a more stable solution. But during the studies I observed, that the convergence of the solver seem to be very sensitive towards the lower discharge limit set in the experiment (3.2V). So when I raised this value, the solving procedure seems to become much more stable.

I now set up a minimum example:

import pybamm
import numpy as np
import matplotlib.pylab as plt


# Experiment
experimental_prozedure=[("Charge at C/2 until 4.3 V",
                         "Hold at 4.3 V for 30 min",
                         "Discharge at 1 C until 2.8 V"),]

experimental_period='1 s'

N_cycles=4
experiment = pybamm.Experiment(
    experimental_prozedure*(N_cycles),period=experimental_period 
)


# Model
model_options={'working electrode':'positive'}
model=pybamm.lithium_ion.DFN(options=model_options)

# Parameters
chemistry=pybamm.parameter_sets.Xu2019
param = pybamm.ParameterValues(chemistry=chemistry)
param['Lower voltage cut-off [V]']=2.0
param['Upper voltage cut-off [V]']=4.4
param['Initial concentration in positive electrode [mol.m-3]']=28938


# Solver
solver=pybamm.CasadiSolver(mode="safe",dt_max=50,atol=1e-5,rtol=1e-5,root_tol=1e-5)

# Mesh
var_pts={'x_n': 20,     # Negative Electrode, default:20
         'x_s': 20,    # Separator, default:20
         'x_p': 20,    # Positive Electrode, default:20
         'r_n': 20,    # Negative Particle Radius, default:20
         'r_p': 20}    # Positive Particle Radius, default:20
    #           'y': 10,      # y-dimension, default:10 -> not used in model, default:10
    #           'z': 10,      # z-dimension, default:10 -> not used in model, default:10
    #           'R_n': 150,   # ?,default:10-> not used in model 
    #           'R_p': 150}   # ?,default:10-> not used in model 
    
#Running simulation
simulation = pybamm.Simulation(model, experiment=experiment, parameter_values=param,solver=solver,var_pts=var_pts)

sim_inputs={}        
solutions=simulation.solve(inputs=sim_inputs)
    
simulation.plot()

When discharging to 3.2 V no numerical problems arise. However, when discharging to 2.8V the warnings mentioned in the previous post are printed.

Do you have any suggestions how to continue the analysis or to increase the robustness of the system? I would really appreciate any suggestions.

With best regards
Carsten

[dion-w]

Looking at the NMC532 curve, what the solver actually does: he tries to solve the system with a certain timestep, if he can't solve he reduces the timestep, and he does up to 5 times by default. This is why also your dt_max is lower at the end then what you set up before. The solver decreases these values itself to increase convergence.
You can imagine it like trying to measure the distance from yourself to a wall which is 1.534 meters away but all you have as measurements tools are sticks that are exactly 1 meter long. You can see that 1 stick would fit between you and the wall, but not two sticks. You cut down the second stick into two halves and see: "1 stick with 1 meter, 1/2 stick with 0.5 meter, and 0.034 meters left". To measure the last bit of distance you have to cut down the third stick a number of times and this of course takes time.
Relating the metaphor with the NMC532 model, the wall is your maximum concentration in the solid, meaning the lithiation of your NMC and the distance to the wall is the timestep of the simulation.

Also currently i'm working on changing the casadi solver so it creates output when i doesn't converge anymore instead of just throwing away all the results it calculated until the "error point". So i also had this problem and am working on it, can't tell you when it's finished though.

[CaXYZ]

Thank you for your fast reply and detailed explanation!