[Bug]: PSD + stoichiometry dependant diffusion bug [Half cell]

[MarcAntoine88]

PyBaMM Version

23.4.1

Python Version

23.4.1

Describe the bug

Hello everyone, I would like to run some half cell graphite simulations using particle size distribution and stoichiometry dependant diffusion within the particles. However I'm getting an error related to discretisation. Did someone already face the same issue? Thanks for the help! Cheers, Marc-Antoine

Steps to Reproduce

import pybamm
import matplotlib.pyplot as plt

options_psd = {"working electrode" : "positive", "particle size": "distribution"}
model = pybamm.lithium_ion.DFN(options = options_psd)

parameter_values = pybamm.ParameterValues("OKane2022")

parameter_values = pybamm.get_size_distribution_parameters(parameter_values, sd_n=0.2, sd_p=0.4)

parameter_values["Maximum concentration in positive electrode [mol.m-3]"] = parameter_values["Maximum concentration in negative electrode [mol.m-3]"]
parameter_values["Initial concentration in positive electrode [mol.m-3]"] = 0.001*parameter_values["Initial concentration in negative electrode [mol.m-3]"]
parameter_values["Positive electrode OCP [V]"] = parameter_values["Negative electrode OCP [V]"]
parameter_values["Lower voltage cut-off [V]"] = 0.001
parameter_values["Upper voltage cut-off [V]"] = 1.6
parameter_values["Positive electrode active material volume fraction"] = parameter_values["Negative electrode active material volume fraction"]
parameter_values["Positive electrode porosity"] = parameter_values["Negative electrode porosity"]

def graphite_LGM50_diffusivity_ORegan2022(sto, T):
    """
    LG M50 Graphite diffusivity as a function of stochiometry, in this case the
    diffusivity is taken to be a constant. The value is taken from [1].

    References
    ----------
    .. [1] Kieran O’Regan, Ferran Brosa Planella, W. Dhammika Widanage, and Emma
    Kendrick. "Thermal-electrochemical parameters of a high energy lithium-ion
    cylindrical battery." Electrochimica Acta 425 (2022): 140700

    Parameters
    ----------
    sto: :class:`pybamm.Symbol`
       Electrode stochiometry
    T: :class:`pybamm.Symbol`
       Dimensional temperature

    Returns
    -------
    :class:`pybamm.Symbol`
       Solid diffusivity
    """

    a0 = 11.17
    a1 = -1.553
    a2 = -6.136
    a3 = -9.725
    a4 = 1.85
    b1 = 0.2031
    b2 = 0.5375
    b3 = 0.9144
    b4 = 0.5953
    c0 = -15.11
    c1 = 0.0006091
    c2 = 0.06438
    c3 = 0.0578
    c4 = 0.001356
    d = 2092

    D_ref = (
        10
        ** (
            a0 * sto
            + c0
            + a1 * pybamm.exp(-((sto - b1) ** 2) / c1)
            + a2 * pybamm.exp(-((sto - b2) ** 2) / c2)
            + a3 * pybamm.exp(-((sto - b3) ** 2) / c3)
            + a4 * pybamm.exp(-((sto - b4) ** 2) / c4)
        )
        * 3.0321  # correcting factor (see O'Regan et al 2021)
    )

    E_D_s = d * pybamm.constants.R
    arrhenius = pybamm.exp(E_D_s / pybamm.constants.R * (1 / 298.15 - 1 / T))

    return D_ref * arrhenius

parameter_values["Positive electrode diffusivity [m2.s-1]"] = graphite_LGM50_diffusivity_ORegan2022

experiment = pybamm.Experiment(["Discharge at 1C until 0.01V (5 second period)"])

sim = pybamm.Simulation(model, parameter_values = parameter_values, experiment = experiment, solver = pybamm.CasadiSolver(rtol=1e-3))
sim.solve() 

plt.figure('Half cell model results')
plt.plot(sim.solution['Discharge capacity [A.h]'].entries, sim.solution['Terminal voltage [V]'].entries, label = 'PyBaMM')
plt.xlabel('Capacity [A.h]')
plt.ylabel('Terminal voltage [V]')
plt.legend()
plt.grid(True)

Relevant log output

No response

[brosaplanella]

I suspect the issue is with the stoichiometry dependent diffusion coefficient, because commenting parameter_values["Positive electrode diffusivity [m2.s-1]"] = graphite_LGM50_diffusivity_ORegan2022 works fine. I am not too familiar with how the PSD models work though, so hard for me to tell where the problem is.

[MarcAntoine88]

Hi @brosaplanella , yes exactly, both PSD and stoichiometry dependent diffusion coefficient work well separately, but not together. According to the error, when PSD and stoichiometry dependent diffusion are used together the integration_dimension goes to 'tertiary' when defining the integral matrix in definite_integral_matrix function (pybamm/spatial_methods/finite_volume). I tried to add a tertiary case without great success until now, thanks for the help!

[valentinsulzer]

As a workaround you can delete the variable that is causing problems (assuming you don't need to evaluate it)

model = ...
del model.variables["X-averaged positive particle effective diffusivity distribution [m2.s-1]"]

[MarcAntoine88]

Just getting back to this for my work, thanks a lot for the workaround @tinosulzer , I really appreciate!
Best,
MA