[Bug]: Calculating Jacobian failed

[RobinBluehub]

PyBaMM Version

22.6

Python Version

3.8

Describe the bug

When I couple multiple aging mechanisms and conduct multiple experimental cycles to study the battery aging behavior, an error is reported before the cycle period I set is reached. When I increase the number of meshes, I can increase the solution period by a certain amount, but it is far from the period I want to loop, and this increases the program running time. I hope that the Pybamm team can reply in time, which is extremely important to the writing of my paper, thank you very much.

Steps to Reproduce

import pybamm
import os
import numpy as np
import matplotlib.pyplot as plt
os.chdir(pybamm.path[0]+'/..')
pybamm.set_logging_level("NOTICE")

options={
"particle": "Fickian diffusion",
"SEI":"solvent-diffusion limited",
"particle mechanics": "swelling and cracking",
"SEI porosity change":"true",
"loss of active material":"stress-driven",
"lithium plating": "partially reversible"
}
model=pybamm.lithium_ion.DFN(options=options)
params_for_LAM=pybamm.ParameterValues("Ai2020")
param1=pybamm.ParameterValues("OKane2022")

param1.update({'Negative electrode cracking rate':3.9E-20,'Positive electrode cracking rate':3.9E-20}, check_already_exists=False)
param1.update({"Negative electrode Young's modulus [Pa]": 1.5E10}, check_already_exists=False)
param1.update({"Positive electrode Young's modulus [Pa]": 3.75E11}, check_already_exists=False)
param1.update({"Negative electrode Paris' law constant m": 2.2},check_already_exists=False)
param1.update({"Positive electrode Paris' law constant m": 2.2},check_already_exists=False)
param1.update({"Negative electrode initial crack length [m]": 20e-9},check_already_exists=False)
param1.update({"Positive electrode initial crack length [m]": 20e-9},check_already_exists=False)
param1.update({"Negative electrode partial molar volume [m3.mol-1]": 3.1E-6}, check_already_exists=False)
param1.update({"Positive electrode partial molar volume [m3.mol-1]": 1.25E-5}, check_already_exists=False)
param1.update({"Negative electrode Poisson's ratio": 0.3}, check_already_exists=False)
param1.update({"Positive electrode Poisson's ratio": 0.2}, check_already_exists=False)
param1.update({"Negative electrode Paris' law constant b": 1.12},check_already_exists=False)
param1.update({"Positive electrode Paris' law constant b": 1.12},check_already_exists=False)
param1.update({"Negative electrode reference concentration for free of deformation [mol.m-3]": 0}, check_already_exists=False)
param1.update({"Positive electrode reference concentration for free of deformation [mol.m-3]": 0}, check_already_exists=False)
param1.update({"Negative electrode number of cracks per unit area [m-2]": 3.18E15}, check_already_exists=False)
param1.update({"Positive electrode number of cracks per unit area [m-2]": 3.18E15}, check_already_exists=False)
param1.update({"Negative electrode initial crack width [m]": 15e-9}, check_already_exists=False)
param1.update({"Positive electrode initial crack width [m]": 15e-9}, check_already_exists=False)
param1.update({"Cell thermal expansion coefficient [m.K-1]": 1.1E-6}, check_already_exists=False)
param1.update({"Negative electrode volume change": params_for_LAM["Negative electrode volume change"]},check_already_exists=False)
param1.update({"Positive electrode volume change": params_for_LAM["Positive electrode volume change"]},check_already_exists=False)
param1.update({"Negative electrode critical stress [Pa]": 60000000.0},check_already_exists=False)
param1.update({"Positive electrode critical stress [Pa]": 375000000.0},check_already_exists=False)
param1.update({"Negative electrode LAM constant proportional term [s-1]": 2.78e-9}, check_already_exists=False)
param1.update({"Positive electrode LAM constant proportional term [s-1]": 2.78e-9}, check_already_exists=False)
param1.update({"Negative electrode LAM constant exponential term": 2}, check_already_exists=False)
param1.update({"Positive electrode LAM constant exponential term": 2}, check_already_exists=False)
param1.update({"Dead lithium decay rate [s-1]":4E-7})
param1.update({"Lithium plating kinetic rate constant [m.s-1]": 1E-10})
param1.update({"Outer SEI solvent diffusivity [m2.s-1]":2.5E-22})

var_pts={"x_n": 60, "x_s":60, "x_p": 40, "r_n": 60, "r_p":60}
solver=pybamm.CasadiSolver(atol=1e-5,rtol=1e-5,mode="fast with events")

experiment_ini=pybamm.Experiment([
"Hold at 4.2V until 50mA",
"Rest for 2 hours"
])
experiment1=pybamm.Experiment([
"Discharge at 1C until 2.5V",
"Charge at 0.3C until 4.2V",
"Hold at 4.2V until 50mA",
]*1000,termination="80% capacity")

simulation_ini=pybamm.Simulation(model=model,experiment=experiment_ini,parameter_values=param1,solver=solver,var_pts=var_pts)
simulation=pybamm.Simulation(model=model,experiment=experiment1,parameter_values=param1,solver=solver,var_pts=var_pts)
solution_ini=simulation_ini.solve()
solution=simulation.solve(starting_solution=solution_ini)
solution.plot()

Relevant log output

psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
At t = 1.08296 and h = 2.57819e-160, the corrector convergence failed repeatedly or with |h| = hmin.
2022-07-04 20:34:10.573 - [NOTICE] callbacks.on_cycle_end(196): Capacity is now 5.125 Ah (originally 5.132 Ah, will stop at 4.105 Ah)
2022-07-04 20:34:10.573 - [WARNING] callbacks.on_experiment_infeasible(237): 
Experiment is infeasible: 'event: Voltage cut-off [V] [experiment]' was triggered during 'Charge at 0.3C until 4.2V'. The returned solution only contains the first 102 cycles. Try reducing the current, shortening the time interval, or reducing the period.

[brosaplanella]

It seems that the simulation is running fine but it stops early. Does it produce a plot? If so, can you submit a screenshot of the plot at the end time?

[RobinBluehub]

I'll give it a try, thank you very much for taking the time to answer my question.

[RobinBluehub]

Reducing the period of the CV step still doesn't work, the experiment ends early, and the cycle period is even smaller.

[brosaplanella]

Then I am not sure what can help. Can you plot the solution up to that cycle and post it here?

[RobinBluehub]

[RobinBluehub]

`import pybamm
import os
import numpy as np
import matplotlib.pyplot as plt
os.chdir(pybamm.path[0]+'/..')
pybamm.set_logging_level("NOTICE")

options={
#"particle": "Fickian diffusion",
"SEI":"solvent-diffusion limited",
"particle mechanics": "swelling and cracking",
"SEI porosity change":"true",
"loss of active material":"stress-driven",
"lithium plating": "partially reversible",
#"SEI film resistance":"distributed",
}
model=pybamm.lithium_ion.DFN(options=options)
params_for_LAM=pybamm.ParameterValues("Ai2020")
param1=pybamm.ParameterValues("OKane2022")

param1.update({'Negative electrode cracking rate':3.9E-21,'Positive electrode cracking rate':3.9E-22}, check_already_exists=False)
param1.update({"Negative electrode Young's modulus [Pa]": 1.5E10}, check_already_exists=False)
param1.update({"Positive electrode Young's modulus [Pa]": 3.75E11}, check_already_exists=False)
param1.update({"Negative electrode Paris' law constant m": 2.2},check_already_exists=False)
param1.update({"Positive electrode Paris' law constant m": 2.2},check_already_exists=False)
param1.update({"Negative electrode initial crack length [m]": 20e-9},check_already_exists=False)
param1.update({"Positive electrode initial crack length [m]": 20e-9},check_already_exists=False)
param1.update({"Negative electrode partial molar volume [m3.mol-1]": 3.1E-6}, check_already_exists=False)
param1.update({"Positive electrode partial molar volume [m3.mol-1]": 1.25E-5}, check_already_exists=False)
param1.update({"Negative electrode Poisson's ratio": 0.3}, check_already_exists=False)
param1.update({"Positive electrode Poisson's ratio": 0.2}, check_already_exists=False)
param1.update({"Negative electrode Paris' law constant b": 1.12},check_already_exists=False)
param1.update({"Positive electrode Paris' law constant b": 1.12},check_already_exists=False)
param1.update({"Negative electrode reference concentration for free of deformation [mol.m-3]": 0}, check_already_exists=False)
param1.update({"Positive electrode reference concentration for free of deformation [mol.m-3]": 0}, check_already_exists=False)
param1.update({"Negative electrode number of cracks per unit area [m-2]": 3.18E15}, check_already_exists=False)
param1.update({"Positive electrode number of cracks per unit area [m-2]": 3.18E15}, check_already_exists=False)
param1.update({"Negative electrode initial crack width [m]": 15e-9}, check_already_exists=False)
param1.update({"Positive electrode initial crack width [m]": 15e-9}, check_already_exists=False)
param1.update({"Cell thermal expansion coefficient [m.K-1]": 1.1E-6}, check_already_exists=False)
param1.update({"Negative electrode volume change": params_for_LAM["Negative electrode volume change"]},check_already_exists=False)
param1.update({"Positive electrode volume change": params_for_LAM["Positive electrode volume change"]},check_already_exists=False)
param1.update({"Negative electrode critical stress [Pa]": 60000000.0},check_already_exists=False)
param1.update({"Positive electrode critical stress [Pa]": 375000000.0},check_already_exists=False)
param1.update({"Negative electrode LAM constant proportional term [s-1]": 2.78e-9}, check_already_exists=False)
param1.update({"Positive electrode LAM constant proportional term [s-1]": 2.78e-9}, check_already_exists=False)
param1.update({"Negative electrode LAM constant exponential term": 2}, check_already_exists=False)
param1.update({"Positive electrode LAM constant exponential term": 2}, check_already_exists=False)
param1.update({"Negative electrode activation energy for cracking rate [J.mol-1]": params_for_LAM["Negative electrode activation energy for cracking rate [J.mol-1]"]},check_already_exists=False)
param1.update({"Positive electrode activation energy for cracking rate [J.mol-1]": params_for_LAM["Positive electrode activation energy for cracking rate [J.mol-1]"]},check_already_exists=False)
param1.update({"Dead lithium decay rate [s-1]":4E-7})
param1.update({"Lithium plating kinetic rate constant [m.s-1]": 1E-10})
param1.update({"Outer SEI solvent diffusivity [m2.s-1]":2.5E-22})
param1.update({"Negative electrolyte volume fraction":0.25},check_already_exists=False)
param1.update({"Positive electrolyte volume fraction":0.335},check_already_exists=False)

var_pts={"x_n": 40, "x_s":40, "x_p": 40, "r_n": 40, "r_p":40}
solver=pybamm.CasadiSolver(mode="fast with events")

experiment_ini=pybamm.Experiment([
"Hold at 4.2V until 50mA ",
"Rest for 2 hours "
],period="10 seconds")
experiment1=pybamm.Experiment([
"Discharge at 1C until 2.5V ",
"Charge at 0.3C until 4.2V ",
"Hold at 4.2V until 50mA (10 seconds period)",
]*500,period="40 seconds",termination="80% capacity")

simulation_ini=pybamm.Simulation(model=model,experiment=experiment_ini,parameter_values=param1,solver=solver,var_pts=var_pts)
simulation=pybamm.Simulation(model=model,experiment=experiment1,parameter_values=param1,solver=solver,var_pts=var_pts)
solution_ini=simulation_ini.solve()
solution=simulation.solve(starting_solution=solution_ini)
solution.plot([
"X-averaged negative particle crack length",
"X-averaged positive particle crack length",
"X-averaged negative particle surface tangential stress",
"X-averaged positive particle surface tangential stress",
"Negative particle crack length",
"Positive particle crack length",
"Negative particle surface tangential stress",
"Positive particle surface tangential stress",
"X-averaged negative particle cracking rate",
"X-averaged positive particle cracking rate",])`

Error：

• [WARNING] callbacks.on_experiment_infeasible(237):

  Experiment is infeasible: 'event: Voltage cut-off [V] [experiment]' was triggered during 'Hold at 4.2V until 50mA (10 seconds period)'. The returned solution only contains the first 16 cycles. Try reducing the current, shortening the time interval, or reducing the period.

[brosaplanella]

Can you plot the porosity and the loss of active material? Note that reducing the period of the "Charge at 0.3C until 4.2V " might also help.

[RobinBluehub]

[RobinBluehub]

Reducing the period of the "Charge at 0.3C until 4.2V" still doesn't work.

[brosaplanella]

I don't know what the issue is... You might need to play with the solver settings to help it converge, as explained in this notebook. The section "Handling instabilities" might be the most relevant one to you.