[Bug]: Degradation submodels fail to run with high dimensional thermal submodels

[yuqin-qy]

PyBaMM Version

25.10.0

Python Version

3.13.9

Describe the bug

Hi all,

I am trying to model degradation for cells with a temperature gradient in the 'z' dimension.

I have used the OKane2022 dataset, and added the relevant thermal parameters as outlined in this page

It seems that in the loss of active material submodel the shape of the initial conditions do not match the rhs (additional points created alongside the z dimension)

Steps to Reproduce

My code is as follow:

import pybamm
param = pybamm.ParameterValues("OKane2022")
L_y = param["Electrode width [m]"]
L_z = param["Electrode height [m]"]
options = {
    "current collector": "potential pair",
    "dimensionality": 1, 
    "thermal": "x-lumped",
    "SEI": "solvent-diffusion limited",
    "SEI porosity change": "true",
    "lithium plating": "partially reversible",
    "lithium plating porosity change": "true",  
    "particle mechanics": ("swelling and cracking", "swelling only"),
    "SEI on cracks": "true",
    "loss of active material": "stress-driven",
}
model = pybamm.lithium_ion.SPMe(options=options)
param.update(
    {
        "Total heat transfer coefficient [W.m-2.K-1]": 100,
        "Negative current collector surface heat transfer coefficient [W.m-2.K-1]": 100,
        "Positive current collector surface heat transfer coefficient [W.m-2.K-1]": 100,
        "Negative tab heat transfer coefficient [W.m-2.K-1]": 0,
        "Positive tab heat transfer coefficient [W.m-2.K-1]": 0,
        "Edge heat transfer coefficient [W.m-2.K-1]": 0,
        "Positive tab width [m]": 0.04,
        "Positive tab centre y-coordinate [m]": L_y * 0.2,
        "Positive tab centre z-coordinate [m]": L_z,
        "Negative tab width [m]": 0.04,
        "Negative tab centre y-coordinate [m]": L_y * 0.8,
        "Negative tab centre z-coordinate [m]": L_z,
    },
    check_already_exists=False,
)
var_pts = {
    "x_n": 5,
    "x_s": 5,  
    "x_p": 5,  
    "r_n": 30,  
    "r_p": 30,  
    "z": 10
}
sim = pybamm.Simulation(
        model,
        parameter_values=param,
        var_pts=var_pts,
        solver=pybamm.IDAKLUSolver()
    )
sol = sim.solve(showprogress=True)  

Relevant log output

---------------------------------------------------------------------------
ModelError                                Traceback (most recent call last)
Cell In[8], line 49
     35 var_pts = {
     36     "x_n": 5,
     37     "x_s": 5,  
   (...)     41     "z": 10
     42 }
     43 sim = pybamm.Simulation(
     44         model,
     45         parameter_values=param,
     46         var_pts=var_pts,
     47         solver=pybamm.IDAKLUSolver()
     48     )
---> 49 sol = sim.solve(showprogress=True)  

File python3.13/site-packages/pybamm/simulation.py:471, in Simulation.solve(self, t_eval, solver, save_at_cycles, calc_esoh, starting_solution, initial_soc, direction, callbacks, showprogress, inputs, t_interp, initial_conditions, **kwargs)
    468 inputs = inputs or {}
    470 if self.operating_mode in ["without experiment", "drive cycle"]:
--> 471     self.build(initial_soc=initial_soc, direction=direction, inputs=inputs)
    472     if save_at_cycles is not None:
    473         raise ValueError(
    474             "'save_at_cycles' option can only be used if simulating an "
    475             "Experiment "
    476         )

File python3.13/site-packages/pybamm/simulation.py:324, in Simulation.build(self, initial_soc, direction, inputs)
    320 self._mesh = pybamm.Mesh(self._geometry, self._submesh_types, self._var_pts)
    321 self._disc = pybamm.Discretisation(
    322     self._mesh, self._spatial_methods, **self._discretisation_kwargs
    323 )
--> 324 self._built_model = self._disc.process_model(
    325     self._model_with_set_params, inplace=False
    326 )
    327 # rebuilt model so clear solver setup
    328 self._solver._model_set_up = {}

File python3.13/site-packages/pybamm/discretisations/discretisation.py:251, in Discretisation.process_model(self, model, inplace)
    249 if self._check_model_flag:
    250     pybamm.logger.verbose(f"Performing model checks for {model.name}")
--> 251     self.check_model(model_disc)
    253 pybamm.logger.info(f"Finish discretising {model.name}")
    255 # Record that the model has been discretised

File python3.13/site-packages/pybamm/discretisations/discretisation.py:1183, in Discretisation.check_model(self, model)
   1181 def check_model(self, model):
   1182     """Perform some basic checks to make sure the discretised model makes sense."""
-> 1183     self.check_initial_conditions(model)

File python3.13/site-packages/pybamm/discretisations/discretisation.py:1211, in Discretisation.check_initial_conditions(self, model)
   1209 for var in model.rhs.keys():
   1210     if model.rhs[var].shape != model.initial_conditions[var].shape:
-> 1211         raise pybamm.ModelError(
   1212             "rhs and initial conditions must have the same shape after "
   1213             "discretisation but rhs.shape = "
   1214             f"{model.rhs[var].shape} and initial_conditions.shape = {model.initial_conditions[var].shape} for variable '{var}'."
   1215         )
   1216 for var in model.algebraic.keys():
   1217     if model.algebraic[var].shape != model.initial_conditions[var].shape:

ModelError: rhs and initial conditions must have the same shape after discretisation but rhs.shape = (10, 1) and initial_conditions.shape = (1, 1) for variable 'Loss of lithium due to loss of active material in negative electrode [mol]'.

[aman-coder03]

Can I work on this?