Error running SPM model with Degradation Submodel

[Akila1993]

Hello Everyone,

I tired to use degradation models with SPM as per the code below, when "surface form": "differential" is enabled, I cannot run the cycling experiment. But it works without any issue for DFN and SPMe model.

Can someone comment on this, please?

import pybamm
import matplotlib.pyplot as plt
import pbeis
import numpy as np
import matplotlib.pyplot as plt
import time as timer
from scipy.fft import fft


model = pybamm.lithium_ion.SPM(
    {   
        "thermal": "lumped", #lumped thermal submodel 
        "surface form": "differential", #Adding capacitance to the model 
        "SEI":"solvent-diffusion limited", #"reaction limited",
        "SEI porosity change": "true",
        "lithium plating": "partially reversible",
        "lithium plating porosity change": "true",  # alias for "SEI porosity change"
        "particle mechanics":"swelling and cracking", #"swelling only",
        "SEI on cracks": "true",
        "loss of active material": "stress-driven",
        
    }
)

var_pts = {
    "x_n": 5,  # negative electrode
    "x_s": 5,  # separator
    "x_p": 5,  # positive electrode
    "r_n": 30,  # negative particle
    "r_p": 30,  # positive particle
}

def get_chen2020_parameters():
    parameter_values = pybamm.ParameterValues("Chen2020")
    parameters = parameter_values
    return parameters

def get_okane2022_degradation_parameters():

 degradation_parameters = {
 'Dead lithium decay constant [s-1]': 1e-06,
 'Dead lithium decay rate [s-1]':pp["Dead lithium decay rate [s-1]"],
 'Initial plated lithium concentration [mol.m-3]': 0.0,
 'Lithium metal partial molar volume [m3.mol-1]': 1.3e-05,
 'Lithium plating kinetic rate constant [m.s-1]': 1e-09,
 'Lithium plating transfer coefficient': 0.65,
 'Negative electrode LAM constant exponential term': 2.0,
 'Negative electrode LAM constant proportional term [s-1]': 2.7778e-07,
 "Negative electrode Paris' law constant b": 1.12,
 "Negative electrode Paris' law constant m": 2.2,
 "Negative electrode Poisson's ratio": 0.3,
 "Negative electrode Young's modulus [Pa]": 15000000000.0,
 'Negative electrode cracking rate':pp["Negative electrode cracking rate"],
 'Negative electrode critical stress [Pa]': 60000000.0,
 'Negative electrode initial crack length [m]': 2e-08,
 'Negative electrode initial crack width [m]': 1.5e-08,
 'Negative electrode number of cracks per unit area [m-2]': 3180000000000000.0,
 'Negative electrode partial molar volume [m3.mol-1]': 3.1e-06,
 'Negative electrode reference concentration for free of deformation [mol.m-3]': 0.0,
 'Negative electrode volume change':pp["Negative electrode volume change"],
 'Positive electrode LAM constant exponential term': 2.0,
 'Positive electrode LAM constant proportional term [s-1]': 2.7778e-07,
 'Positive electrode OCP entropic change [V.K-1]': 0.0,
 "Positive electrode Paris' law constant b": 1.12,
 "Positive electrode Paris' law constant m": 2.2,
 "Positive electrode Poisson's ratio": 0.2,
 "Positive electrode Young's modulus [Pa]": 375000000000.0,
 'Positive electrode active material volume fraction': 0.665,
 'Positive electrode cracking rate':pp["Positive electrode cracking rate"],
 'Positive electrode critical stress [Pa]': 375000000.0,
 'Positive electrode initial crack length [m]': 2e-08,
 'Positive electrode initial crack width [m]': 1.5e-08,
 'Positive electrode number of cracks per unit area [m-2]': 3180000000000000.0,
 'Positive electrode partial molar volume [m3.mol-1]': 1.25e-05,
 'Positive electrode reference concentration for free of deformation [mol.m-3]': 0.0,
 'Positive electrode volume change':pp["Positive electrode volume change"],
 'Typical plated lithium concentration [mol.m-3]': 1000.0,
 'Exchange-current density for stripping [A.m-2]':pp["Exchange-current density for stripping [A.m-2]"],
 'Exchange-current density for plating [A.m-2]':pp["Exchange-current density for plating [A.m-2]"],
      
    }
 return degradation_parameters

def update_parameters():
    # Load Chen2020 parameters
    chen2020_params = get_chen2020_parameters()

    # Load OKane2022 degradation parameters
    okane2022_degradation_params = get_okane2022_degradation_parameters()

    # Update Chen2020 parameters with OKane2022 degradation parameters
    chen2020_params.update(okane2022_degradation_params, check_already_exists=False)

    return chen2020_params

updated_parameters = update_parameters()

cycle_number = 5
exp = pybamm.Experiment(
    [
        (
            "Discharge at 1C until 2.5 V",  # ageing cycles
            "Charge at 0.6C until 4.2 V (5 minute period)",
            "Hold at 4.2 V until C/100 (5 minute period)",
        )
    ]
    * cycle_number
)

sim = pybamm.Simulation(model, parameter_values=updated_parameters,experiment=exp,solver=pybamm.CasadiSolver(mode='safe without grid'),var_pts=var_pts)


solution = sim.solve()


**_The error message is below:_** 

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[80], line 17
     13 sim = pybamm.Simulation(model, parameter_values=updated_parameters,experiment=exp,solver=pybamm.CasadiSolver(mode='safe without grid'),var_pts=var_pts)
     14 #sim = pybamm.Simulation(model, parameter_values=updated_parameters,solver=pybamm.CasadiSolver(mode='safe without grid'),var_pts=var_pts)
---> 17 solution = sim.solve()
     18 #solution = sim.solve([0, 3700])
     19 #sim.plot()

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/simulation.py:724, in Simulation.solve(self, t_eval, solver, check_model, save_at_cycles, calc_esoh, starting_solution, initial_soc, callbacks, **kwargs)
    722 npts = max(int(round(dt / op_conds["period"])) + 1, 2)
    723 try:
--> 724     step_solution = solver.step(
    725         current_solution,
    726         model,
    727         dt,
    728         npts=npts,
    729         save=False,
    730         **kwargs,
    731     )
    732 except pybamm.SolverError as error:
    733     if (
    734         "non-positive at initial conditions" in error.message
    735         and "[experiment]" in error.message
    736     ):

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/solvers/base_solver.py:1177, in BaseSolver.step(self, old_solution, model, dt, npts, external_variables, inputs, save)
   1174 set_up_time = timer.time()
   1176 # (Re-)calculate consistent initial conditions
-> 1177 self._set_initial_conditions(model, ext_and_inputs, update_rhs=False)
   1179 # Non-dimensionalise dt
   1180 dt_dimensionless = dt / model.timescale_eval

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/solvers/base_solver.py:612, in BaseSolver._set_initial_conditions(self, model, inputs_dict, update_rhs)
    608         else:
    609             model.y0 = np.vstack(
    610                 (y0_from_inputs[:len_rhs], y0_from_model[len_rhs:])
    611             )
--> 612     y0 = self.calculate_consistent_state(model, 0, inputs_dict)
    613 # Make y0 a function of inputs if doing symbolic with casadi
    614 model.y0 = y0

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/solvers/base_solver.py:641, in BaseSolver.calculate_consistent_state(self, model, time, inputs)
    639     return model.y0
    640 try:
--> 641     root_sol = self.root_method._integrate(model, np.array([time]), inputs)
    642 except pybamm.SolverError as e:
    643     raise pybamm.SolverError(
    644         "Could not find consistent states: {}".format(e.args[0])
    645     )

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/solvers/casadi_algebraic_solver.py:129, in CasadiAlgebraicSolver._integrate(self, model, t_eval, inputs_dict)
    126 constraints[model.bounds[1] <= 0] = -1
    128 # Set up rootfinder
--> 129 roots = casadi.rootfinder(
    130     "roots",
    131     "newton",
    132     dict(x=y_alg_sym, p=t_sym, g=alg),
    133     {
    134         **self.extra_options,
    135         "abstol": self.tol,
    136         "constraints": list(constraints[len_rhs:]),
    137     },
    138 )
    140 timer = pybamm.Timer()
    141 integration_time = 0

File ~/pybamm-eis/env/lib/python3.9/site-packages/casadi/casadi.py:30310, in rootfinder(*args)
  29600 def rootfinder(*args) -> "casadi::Function":
  29601     """
  29602       [INTERNAL] 
  29603 
   (...)
  30308 
  30309     """
> 30310     return _casadi.rootfinder(*args)

RuntimeError: .../casadi/core/function_internal.cpp:147: Error calling Newton::init for 'roots':
.../casadi/core/rootfinder.cpp:237: Assertion "!sp_jac_.is_singular()" failed:
Rootfinder::init: singularity - the jacobian is structurally rank-deficient. sprank(J)=10 (instead of 11)

[DrSOKane]

Neither Chen2020 or OKane2022 have thermal parameters, so you cannot use the lumped thermal model unless you add thermal parameters from elsewhere.

[Akila1993]

Hi @DrSOKane ,

Thanks for letting me know,

However without the lumped thermal model also, I cannot run this cycling experiment for the SPM model, but it is working for DFN and SPMe. When "surface form": "differential" is disabled I can run it for SPM too. That is my question.

[DrSOKane]

What happens if you run SPMe and set {"x-average side reactions": "true"}? This option also ruins a piece of code I've been working on for at least a year, so maybe we're facing the same problem...

[Akila1993]

Hi @DrSOKane ,

When I run SPMe and set {"x-average side reactions": "true"}, there is no problem. But with SPM still, the issue persists....

Do you have any suggestions?

[DrSOKane]

My suggestion would be to use SPMe. With the surface form, it's almost as accurate than DFN, while being much faster. I realize it's not ideal but it's my honest opinion at this point, especially for high currents where the electrolyte is a major limitation on the amount of power you can draw, and for how long.