Error running PyBaMM EIS with degradation submodels enabled

[Akila1993]

Hello,

I tried using PyBaMM-EIS with the degradation model enabled, However, it is keeping me with the following errors, whereas when using the typical time domain method it is correctly giving the EIS results.

Is it a issue with PyBaMM-EIS.

Also in your examples to import PybaMM EIS you mentioned using import pybammeis but it can only import with import pbeis.

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

# Set random seed
np.random.seed(0)

# Load parameters
pp = pybamm.ParameterValues("OKane2022")

# Define model with configuration
model = pybamm.lithium_ion.SPMe(
    {   
        "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 only",
        "particle mechanics":"swelling and cracking",
        "SEI on cracks": "true",
        #"loss of active material": "stress-driven",
        
    }
)

var_pts = {
    "x_n": 5,
    "x_s": 5,
    "x_p": 5,
    "r_n": 30,
    "r_p": 30,
}

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

# Simulation parameters
cycle_numbers = [1]
frequencies = np.logspace(-4, 3, 30)
I = 1e-3  # Current amplitude for EIS
number_of_periods = 20
samples_per_period = 20

# Store results for Nyquist plot
impedance_data = {}

# Loop through each cycle count
for cycle_number in cycle_numbers:
    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,
    )

    sol = sim.solve()

    # Set initial conditions of aged_model from the solution and EIS Running
    aged_model = model.set_initial_conditions_from(sol)


eis_sim = pbeis.EISSimulation(aged_model)
# calculate impedance at log-spaced frequencies
eis_sim.solve(np.logspace(-2, 4, 30))
eis_sim.nyquist_plot()


---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:559, in ParameterValues.process_symbol(self, symbol)
    558 try:
--> 559     return self._processed_symbols[symbol]
    560 except KeyError:

KeyError: XAverage(0x39c1c2e9f02566c4, x-average, children=["Negative electrode cracking rate * (Negative electrode Young's modulus [Pa] ** Negative electrode Paris' law constant m) * (Negative electrode initial crack length [m] ** ((Negative electrode Paris' law constant m / 2.0) - 1.0)) * ((((Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / 3.0 / (1.0 - Negative electrode Poisson's ratio)) * (broadcast(r-average(X-averaged negative particle concentration)) - broadcast(boundary value(X-averaged negative particle concentration))) * Maximum concentration in negative electrode [mol.m-3] / Negative electrode Young's modulus [Pa]) * Negative electrode Paris' law constant b * sqrt(broadcast(3.141592653589793 * X-averaged negative particle crack length)) * (0.0 <= (Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / 3.0 / (1.0 - Negative electrode Poisson's ratio)) * (broadcast(r-average(X-averaged negative particle concentration)) - broadcast(boundary value(X-averaged negative particle concentration))) * Maximum concentration in negative electrode [mol.m-3] / Negative electrode Young's modulus [Pa])) ** Negative electrode Paris' law constant m) / (3600.0 / (abs(Typical current [A] / Nominal cell capacity [A.h]) * 96485.33212 * Maximum concentration in negative electrode [mol.m-3] * (Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m]) / absolute(Typical current [A] / (Number of electrodes connected in parallel to make a cell * Electrode width [m] * Electrode height [m]))))"], domains={'primary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:559, in ParameterValues.process_symbol(self, symbol)
    558 try:
--> 559     return self._processed_symbols[symbol]
    560 except KeyError:

KeyError: Division(-0x37934c81d85ca00d, /, children=["Negative electrode cracking rate * (Negative electrode Young's modulus [Pa] ** Negative electrode Paris' law constant m) * (Negative electrode initial crack length [m] ** ((Negative electrode Paris' law constant m / 2.0) - 1.0)) * ((((Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / 3.0 / (1.0 - Negative electrode Poisson's ratio)) * (broadcast(r-average(X-averaged negative particle concentration)) - broadcast(boundary value(X-averaged negative particle concentration))) * Maximum concentration in negative electrode [mol.m-3] / Negative electrode Young's modulus [Pa]) * Negative electrode Paris' law constant b * sqrt(broadcast(3.141592653589793 * X-averaged negative particle crack length)) * (0.0 <= (Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / 3.0 / (1.0 - Negative electrode Poisson's ratio)) * (broadcast(r-average(X-averaged negative particle concentration)) - broadcast(boundary value(X-averaged negative particle concentration))) * Maximum concentration in negative electrode [mol.m-3] / Negative electrode Young's modulus [Pa])) ** Negative electrode Paris' law constant m)", '3600.0 / (abs(Typical current [A] / Nominal cell capacity [A.h]) * 96485.33212 * Maximum concentration in negative electrode [mol.m-3] * (Negative electrode thickness [m] + Separator thickness [m] + Positive electrode thickness [m]) / absolute(Typical current [A] / (Number of electrodes connected in parallel to make a cell * Electrode width [m] * Electrode height [m])))'], domains={'primary': ['negative electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:559, in ParameterValues.process_symbol(self, symbol)
    558 try:
--> 559     return self._processed_symbols[symbol]
    560 except KeyError:

KeyError: Multiplication(0xe005a9b6c739eb7, *, children=["Negative electrode cracking rate * (Negative electrode Young's modulus [Pa] ** Negative electrode Paris' law constant m) * (Negative electrode initial crack length [m] ** ((Negative electrode Paris' law constant m / 2.0) - 1.0))", "(((Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / 3.0 / (1.0 - Negative electrode Poisson's ratio)) * (broadcast(r-average(X-averaged negative particle concentration)) - broadcast(boundary value(X-averaged negative particle concentration))) * Maximum concentration in negative electrode [mol.m-3] / Negative electrode Young's modulus [Pa]) * Negative electrode Paris' law constant b * sqrt(broadcast(3.141592653589793 * X-averaged negative particle crack length)) * (0.0 <= (Negative electrode partial molar volume [m3.mol-1] * Negative electrode Young's modulus [Pa] / 3.0 / (1.0 - Negative electrode Poisson's ratio)) * (broadcast(r-average(X-averaged negative particle concentration)) - broadcast(boundary value(X-averaged negative particle concentration))) * Maximum concentration in negative electrode [mol.m-3] / Negative electrode Young's modulus [Pa])) ** Negative electrode Paris' law constant m"], domains={'primary': ['negative electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:559, in ParameterValues.process_symbol(self, symbol)
    558 try:
--> 559     return self._processed_symbols[symbol]
    560 except KeyError:

KeyError: Multiplication(0x1008323fdf220993, *, children=['Negative electrode cracking rate', "(Negative electrode Young's modulus [Pa] ** Negative electrode Paris' law constant m) * (Negative electrode initial crack length [m] ** ((Negative electrode Paris' law constant m / 2.0) - 1.0))"], domains={'primary': ['negative electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:559, in ParameterValues.process_symbol(self, symbol)
    558 try:
--> 559     return self._processed_symbols[symbol]
    560 except KeyError:

KeyError: FunctionParameter(0x311aab223c581fd0, Negative electrode cracking rate, children=['broadcast(broadcast(Volume-averaged cell temperature + Reference temperature [K]))'], domains={'primary': ['negative electrode'], 'secondary': ['current collector']})

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/util.py:59, in FuzzyDict.__getitem__(self, key)
     58 try:
---> 59     return super().__getitem__(key)
     60 except KeyError:

KeyError: 'Negative electrode cracking rate'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
Cell In[12], line 1
----> 1 eis_sim = pbeis.EISSimulation(aged_model)
      2 # calculate impedance at log-spaced frequencies
      3 eis_sim.solve(np.logspace(-2, 4, 30))

File ~/pybamm-eis/pbeis/eis_simulation.py:61, in EISSimulation.__init__(self, model, parameter_values, geometry, submesh_types, var_pts, spatial_methods)
     52 parameter_values["Current function [A]"] = 0
     53 sim = pybamm.Simulation(
     54     self.model,
     55     geometry=geometry,
   (...)
     59     spatial_methods=spatial_methods,
     60 )
---> 61 sim.build()
     62 self.built_model = sim.built_model
     64 # Extract mass matrix and Jacobian

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/simulation.py:436, in Simulation.build(self, check_model, initial_soc)
    434     self._built_model = self.model
    435 else:
--> 436     self.set_parameters()
    437     self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts)
    438     self._disc = pybamm.Discretisation(self._mesh, self._spatial_methods)

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/simulation.py:374, in Simulation.set_parameters(self)
    372     self._model_with_set_params = self._unprocessed_model
    373 else:
--> 374     self._model_with_set_params = self._parameter_values.process_model(
    375         self._unprocessed_model, inplace=False
    376     )
    377     self._parameter_values.process_geometry(self._geometry)
    378 self.model = self._model_with_set_params

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:396, in ParameterValues.process_model(self, unprocessed_model, inplace)
    392 for variable, equation in unprocessed_model.rhs.items():
    393     pybamm.logger.verbose(
    394         "Processing parameters for {!r} (rhs)".format(variable)
    395     )
--> 396     new_rhs[variable] = self.process_symbol(equation)
    397 model.rhs = new_rhs
    399 new_algebraic = {}

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:561, in ParameterValues.process_symbol(self, symbol)
    559     return self._processed_symbols[symbol]
    560 except KeyError:
--> 561     processed_symbol = self._process_symbol(symbol)
    562     self._processed_symbols[symbol] = processed_symbol
    564     return processed_symbol

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:717, in ParameterValues._process_symbol(self, symbol)
    715 # Unary operators
    716 elif isinstance(symbol, pybamm.UnaryOperator):
--> 717     new_child = self.process_symbol(symbol.child)
    718     new_symbol = symbol._unary_new_copy(new_child)
    719     # ensure domain remains the same

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:561, in ParameterValues.process_symbol(self, symbol)
    559     return self._processed_symbols[symbol]
    560 except KeyError:
--> 561     processed_symbol = self._process_symbol(symbol)
    562     self._processed_symbols[symbol] = processed_symbol
    564     return processed_symbol

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:708, in ParameterValues._process_symbol(self, symbol)
    704     return self.process_symbol(function_out)
    706 elif isinstance(symbol, pybamm.BinaryOperator):
    707     # process children
--> 708     new_left = self.process_symbol(symbol.left)
    709     new_right = self.process_symbol(symbol.right)
    710     # make new symbol, ensure domain remains the same

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:561, in ParameterValues.process_symbol(self, symbol)
    559     return self._processed_symbols[symbol]
    560 except KeyError:
--> 561     processed_symbol = self._process_symbol(symbol)
    562     self._processed_symbols[symbol] = processed_symbol
    564     return processed_symbol

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:708, in ParameterValues._process_symbol(self, symbol)
    704     return self.process_symbol(function_out)
    706 elif isinstance(symbol, pybamm.BinaryOperator):
    707     # process children
--> 708     new_left = self.process_symbol(symbol.left)
    709     new_right = self.process_symbol(symbol.right)
    710     # make new symbol, ensure domain remains the same

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:561, in ParameterValues.process_symbol(self, symbol)
    559     return self._processed_symbols[symbol]
    560 except KeyError:
--> 561     processed_symbol = self._process_symbol(symbol)
    562     self._processed_symbols[symbol] = processed_symbol
    564     return processed_symbol

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:708, in ParameterValues._process_symbol(self, symbol)
    704     return self.process_symbol(function_out)
    706 elif isinstance(symbol, pybamm.BinaryOperator):
    707     # process children
--> 708     new_left = self.process_symbol(symbol.left)
    709     new_right = self.process_symbol(symbol.right)
    710     # make new symbol, ensure domain remains the same

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:561, in ParameterValues.process_symbol(self, symbol)
    559     return self._processed_symbols[symbol]
    560 except KeyError:
--> 561     processed_symbol = self._process_symbol(symbol)
    562     self._processed_symbols[symbol] = processed_symbol
    564     return processed_symbol

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:596, in ParameterValues._process_symbol(self, symbol)
    594     else:
    595         new_children.append(self.process_symbol(child))
--> 596 function_name = self[symbol.name]
    598 # Create Function or Interpolant or Scalar object
    599 if isinstance(function_name, tuple):

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/parameters/parameter_values.py:111, in ParameterValues.__getitem__(self, key)
    110 def __getitem__(self, key):
--> 111     return self._dict_items[key]

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/util.py:62, in FuzzyDict.__getitem__(self, key)
     60 except KeyError:
     61     best_matches = self.get_best_matches(key)
---> 62     raise KeyError(f"'{key}' not found. Best matches are {best_matches}")

KeyError: "'Negative electrode cracking rate' not found. Best matches are ['Negative electrode porosity', 'Negative electrode electrons in reaction', 'Negative electrode thickness [m]']"