[Bug]: Inconsistencies between `pybamm.lithium_ion.SPMe` and Marquis 2019

[ejfdickinson]

PyBaMM Version

24.1

Python Version

23.13

Describe the bug

pybamm.lithium_ion.SPMe claims equivalence with the canonical SPMe model defined in the Marquis 2019 research paper. However, the following inconsistencies have been noted for isothermal SPMe computation using pybamm.lithium_ion.SPMe and Marquis 2019:

• The canonical SPMe in Marquis 2019 enforces that electrolyte diffusivity and electrolyte conductivity are concentration-independent (to first-order as used to define the SPMe, equations [40c] / [40p]), but pybamm.lithium_ion.SPMe uses a lookup to the local concentration, if concentration-dependence is provided in the parameter values.
  • This follows the definition in Brosa Planella 2022 wherein local concentration-dependence is stated as simply optional (although allowing a nonlinear electrolyte diffusion equation will break the leading-order assumptions in Marquis 2019, to the best of my understanding)
  • This contradicts the PyBaMM documentation, which only cites Marquis 2019
• The relation between reported X-averaged reaction overpotential and X-averaged exchange current density does not agree with Marquis 2019 [48b] ([40l] in dimensionless form)
• The relation between reported X-averaged concentration overpotential and X-averaged electrolyte concentrations does not agree with Marquis 2019 [48c] ([40m] in dimensionless form)

By comparison to a benchmark of the Marquis 2019 canonical SPMe implemented in COMSOL Multiphysics, within which the corresponding equalities have been demonstrated, these discrepancies appear to have a measurable consequence for the reported terminal voltage.

In light of the investigative effort that was required to run down these discrepancies, I think there is a case for the following:

• thorough documentation in mathematical equations of the actual contents of different PyBaMM Models, as distinct from just citing references (supports #3392)
• additionally (or perhaps alternatively?) the ability to extract the mathematical content of a model (in legible syntax, not raw source code) from a pybamm.Simulation
  • the high degree of submodel inheritance as well as the parent-child structure of the expression trees in pybamm.BaseModel.rhs and pybamm.BaseModel.algebraic make this hard to tease out in my experience
  • compare for example the COMSOL Multiphysics "Equation View" which is algebraically explicit and, perhaps more importantly, directly syntactically compatible with results evaluation on the corresponding numerical solution objects
  • my apologies if there is already a method in there that does this, of which I'm unaware!

(supersedes #3796 - specific issues have been clarified here and a new report logged)

Steps to Reproduce

import pybamm
import numpy as np
import matplotlib.pyplot as plt

## Solve an SPMe model at fairly high current to provoke transport losses
parameter_values = pybamm.ParameterValues("Chen2020")
model = pybamm.lithium_ion.SPMe()

experiment = pybamm.Experiment(
    [
        "Discharge at 2C for 5 minutes",
    ],
    period="0.1 s",
)

sim = pybamm.Simulation(
    model,
    experiment=experiment,
    parameter_values=parameter_values,
)
sol = sim.solve()

## Diffusivity is evaluated locally
x = sol["Electrolyte concentration [mol.m-3]"].mesh.edges
t_example = 200
cl = sol["Electrolyte concentration [mol.m-3]"](t_example, x=x)
clx = np.gradient(cl) / np.gradient(x)
Nl = sol["Electrolyte diffusion flux [mol.m-2.s-1]"](t_example, x=x)
Beff = sol["Electrolyte transport efficiency"](t_example, x=x)

Dl_implied = -Nl/Beff/clx
Dl_local = parameter_values["Electrolyte diffusivity [m2.s-1]"](cl, parameter_values["Ambient temperature [K]"])

plt.figure()
plt.plot(x, Dl_local, label="Concentration-dependent (local)")
plt.plot(x, Dl_implied, label="Implied from solution data (numerical)")
plt.xlabel("x / m")
plt.ylabel("Electrolyte diffusivity / m2.s-1")
plt.grid()
plt.legend()

## X-averaged overpotentials disagree with Marquis 2019
t_eval = sol["Time [s]"].entries
t = sol["Time [s]"](t=t_eval)

cl_ave_neg = sol["X-averaged negative electrolyte concentration [mol.m-3]"](t=t_eval)
cl_ave_pos = sol["X-averaged positive electrolyte concentration [mol.m-3]"](t=t_eval)
i0_ave_neg = sol["X-averaged negative electrode exchange current density [A.m-2]"](t=t_eval)
i0_ave_pos = sol["X-averaged positive electrode exchange current density [A.m-2]"](t=t_eval)
isurf_neg = sol["X-averaged negative electrode interfacial current density [A.m-2]"](t=t_eval)
isurf_pos = sol["X-averaged positive electrode interfacial current density [A.m-2]"](t=t_eval)

eta_neg = sol["X-averaged negative electrode reaction overpotential [V]"](t=t_eval)
eta_pos = sol["X-averaged positive electrode reaction overpotential [V]"](t=t_eval)
eta_c = sol["X-averaged concentration overpotential [V]"](t=t_eval)

f = (pybamm.constants.F / pybamm.constants.R / parameter_values["Ambient temperature [K]"]).value
cl_0 = parameter_values["Initial concentration in electrolyte [mol.m-3]"]
tplus = parameter_values["Cation transference number"]

# Marquis 2019 [48l]
#   Note: multiple of (1/2) in argument of arcsinh() according to
#   correct exchange current density definition as discussed in erratum to Marquis 2019
eta_neg_manual = (2 / f) * np.arcsinh(isurf_neg / i0_ave_neg / 2)
eta_pos_manual = (2 / f) * np.arcsinh(isurf_pos / i0_ave_pos / 2)
# Marquis 2019 [48m]
eta_c_manual = (2 / f) * (1 - tplus) * (cl_ave_pos - cl_ave_neg) / cl_0

plt.figure()
plt.plot(t, eta_neg, label="Computed")
plt.plot(t, eta_neg_manual, label="Implied (Marquis 2019)")
plt.xlabel("Time / s")
plt.ylabel("Negative electrode overpotential / V")
plt.grid()
plt.legend()

plt.figure()
plt.plot(t, eta_pos, label="Computed")
plt.plot(t, eta_pos_manual, label="Implied (Marquis 2019)")
plt.xlabel("Time / s")
plt.ylabel("Positive electrode overpotential / V")
plt.grid()
plt.legend()

plt.figure()
plt.plot(t, eta_c, label="Computed")
plt.plot(t, eta_c_manual, label="Implied (Marquis 2019)")
plt.xlabel("Time / s")
plt.ylabel("Concentration overpotential / V")
plt.grid()
plt.legend()

Relevant log output

No response

[ejfdickinson]

To add - although the repeat uses default numerical settings, I've done my best to exclude numerical error as a source for the discrepancy, which appears to be roughly mesh-independent within the testing range I've explored.

[valentinsulzer]

Our best effort to automatically display the equations is the model.latexify() function which has an example here: https://github.com/pybamm-team/PyBaMM/blob/develop/docs/source/examples/notebooks/models/latexify.ipynb

It needs more work but it is quite challenging to automatically show the right level of detail with complex expressions like the voltage in the SPMe

[rtimms]

Thanks for all your effort in highlighting these discrepancies. They all come about because in PyBaMM we use a composite expression for all the quantities involving the electrolyte variables whereas in Marquis 2019 everything is explicitly written as leading order plus a correction. You point out one example of this where we use the nonlinear diffusivity.

Essentially anywhere there is something of the form f(c_e) we use the composite expansion f(c_e0 + \delta c_e1) and Marquis 2019 usue f(c_e0) + delta f'(c_e0)c_e1. The documentation needs updating to reflect this, as expecting users to do this re-analysis is unreasonable.

I've attached an updated version of your example where you can see this comparison for the reaction overpotentials. I can complete the other expressions (and make an SPMe notebook to document this) when I have some more free time.

import pybamm
import numpy as np
import matplotlib.pyplot as plt

## Solve an SPMe model at fairly high current to provoke transport losses
parameter_values = pybamm.ParameterValues("Chen2020")
model = pybamm.lithium_ion.SPMe()

experiment = pybamm.Experiment(
    [
        "Discharge at 2C for 5 minutes",
    ],
    period="0.1 s",
)

sim = pybamm.Simulation(
    model,
    experiment=experiment,
    parameter_values=parameter_values,
)
sol = sim.solve()

## Diffusivity is evaluated locally
x = sol["Electrolyte concentration [mol.m-3]"].mesh.edges
t_example = 200
cl = sol["Electrolyte concentration [mol.m-3]"](t_example, x=x)
clx = np.gradient(cl) / np.gradient(x)
Nl = sol["Electrolyte diffusion flux [mol.m-2.s-1]"](t_example, x=x)
Beff = sol["Electrolyte transport efficiency"](t_example, x=x)

Dl_implied = -Nl/Beff/clx
Dl_local = parameter_values["Electrolyte diffusivity [m2.s-1]"](cl, parameter_values["Ambient temperature [K]"])

plt.figure()
plt.plot(x, Dl_local, label="Concentration-dependent (local)")
plt.plot(x, Dl_implied, label="Implied from solution data (numerical)")
plt.xlabel("x / m")
plt.ylabel("Electrolyte diffusivity / m2.s-1")
plt.grid()
plt.legend()

## X-averaged overpotentials disagree with Marquis 2019
t_eval = sol["Time [s]"].entries
t = sol["Time [s]"](t=t_eval)

cl_ave_neg = sol["X-averaged negative electrolyte concentration [mol.m-3]"](t=t_eval)
cl_ave_pos = sol["X-averaged positive electrolyte concentration [mol.m-3]"](t=t_eval)
i0_ave_neg = sol["X-averaged negative electrode exchange current density [A.m-2]"](t=t_eval)
i0_ave_pos = sol["X-averaged positive electrode exchange current density [A.m-2]"](t=t_eval)
i0_neg = sol["Negative electrode exchange current density [A.m-2]"].data
i0_pos = sol["Positive electrode exchange current density [A.m-2]"].data
isurf_neg = sol["X-averaged negative electrode interfacial current density [A.m-2]"](t=t_eval)
isurf_pos = sol["X-averaged positive electrode interfacial current density [A.m-2]"](t=t_eval)

eta_neg = sol["X-averaged negative electrode reaction overpotential [V]"](t=t_eval)
eta_pos = sol["X-averaged positive electrode reaction overpotential [V]"](t=t_eval)
eta_c = sol["X-averaged concentration overpotential [V]"](t=t_eval)

f = (pybamm.constants.F / pybamm.constants.R / parameter_values["Ambient temperature [K]"]).value
cl_0 = parameter_values["Initial concentration in electrolyte [mol.m-3]"]
tplus = parameter_values["Cation transference number"]

# Marquis 2019 [48l]
#   Note: multiple of (1/2) in argument of arcsinh() according to
#   correct exchange current density definition as discussed in erratum to Marquis 2019
eta_neg_manual = (2 / f) * np.arcsinh(isurf_neg / i0_ave_neg / 2)
eta_pos_manual = (2 / f) * np.arcsinh(isurf_pos / i0_ave_pos / 2)
eta_neg_manual_nl = np.mean((2 / f) * np.arcsinh(isurf_neg / i0_neg / 2), axis=0)
eta_pos_manual_nl = np.mean((2 / f) * np.arcsinh(isurf_pos / i0_pos / 2), axis=0)
# Marquis 2019 [48m]
eta_c_manual = (2 / f) * (1 - tplus) * (cl_ave_pos - cl_ave_neg) / cl_0

plt.figure()
plt.plot(t, eta_neg, label="Computed")
plt.plot(t, eta_neg_manual, label="Implied (Marquis 2019)")
plt.plot(t, eta_neg_manual_nl, "x", label="Implied (Marquis 2019, non-linear)")
plt.xlabel("Time / s")
plt.ylabel("Negative electrode overpotential / V")
plt.grid()
plt.legend()

plt.figure()
plt.plot(t, eta_pos, label="Computed")
plt.plot(t, eta_pos_manual, label="Implied (Marquis 2019)")
plt.plot(t, eta_pos_manual_nl, "x", label="Implied (Marquis 2019, non-linear)")
plt.xlabel("Time / s")
plt.ylabel("Positive electrode overpotential / V")
plt.grid()
plt.legend()

plt.figure()
plt.plot(t, eta_c, label="Computed")
plt.plot(t, eta_c_manual, label="Implied (Marquis 2019)")
plt.xlabel("Time / s")
plt.ylabel("Concentration overpotential / V")
plt.grid()
plt.legend()
plt.show()

[ejfdickinson]

@rtimms Thanks for the further analysis!

So, for the record - pybamm.lithium_ion.SPMe does not implement Marquis 2019, and irrespective of anything else, there's a flat documentation bug in the claim that it does?

My understanding is that the composite expansion causes leakage of higher-order terms into the solution where a term has been represented 'functionally', while corresponding higher-order terms in e.g. partial derivatives of dependent variables, are still absent because they have been explicitly excluded from the equation form. This makes it confusing to me to what extent pybamm.lithium_ion.SPMe behaves approximately to the DFN.

• Is the actual equation set documented anywhere in a mathematical syntax, e.g. via some interpretation of Brosa Planella 2022?
• Can you recommend any shortcut workarounds such that pybamm.lithium_ion.SPMe can be obliged to represent Marquis 2019?

[brosaplanella]

Hi! I believe this article I wrote a while ago (https://www.sciencedirect.com/science/article/pii/S0013468621008148) comes quite close to what Rob described above (model presented in Section 3, but note it also includes thermal effects so there might be some extra terms).

[ejfdickinson]

Thanks @brosaplanella - the thermal effects are worth having!

[ejfdickinson]

linking discussion on possible fix in #4274