Repository navigation
[Bug]: Lithium is not conserved with concentration-dependent particle diffusivity on a non-uniform particle mesh #5754
Description
PyBaMM Version: 26.8.0.1.dev11+g861f9315c
Python Version: 3.12.11
Describe the bug
When the particle submesh is non-uniform and the particle diffusivity depends on
concentration, the Fickian particle imposes the wrong surface flux, so lithium is
silently created or destroyed. On ORegan2022 with an Exponential1DSubMesh and
r_n = r_p = 16, the cell gains 26 % of its particle lithium over three 1C
cycles. Nothing in the model reports a loss or gain — Total lithium in particles [mol] simply drifts.
The error is a systematic multiplicative bias on the surface flux, applied in the
same direction on every intercalation step, so it accumulates monotonically rather
than averaging out. It is easy to mistake for a degradation signal: in a long
ageing run with LAM and SEI submodels active it shows up as LLI [%] that does
not reconcile with the model's own loss terms.
Both ingredients are required. Uniform mesh: exact. Constant diffusivity: exact.
Cause
FickianDiffusion sets the surface boundary condition as a gradient, dividing
the interfacial current by the boundary value of the diffusivity
(fickian_diffusion.py:209):
f"{Domain} {phase_name}particle bc [mol.m-4]": -j * R_nondim / self.param.F / pybamm.surf(D_eff)The flux that condition actually produces is N = -D_eff * grad(c), where D_eff
is interpolated from nodes onto the surface edge. That interpolation is
FiniteVolume.shift's arithmetic_mean, whose boundary rows are hard-coded
[1.5, -0.5] and [-0.5, 1.5] (finite_volume.py:1551-1560) — a linear
extrapolation that assumes uniformly spaced nodes. pybamm.surf extrapolates
using the real node positions.
On a uniform mesh the two agree identically, so the imposed flux is exactly j/F.
On a graded mesh they do not, and the imposed flux is off by
D_node->edge / surf(D_eff).
Note that harmonic_mean in the same method already weights by the true
control-volume widths (beta = dx1 / (dx1 + dx2)); only arithmetic_mean
hard-codes uniform spacing.
Steps to Reproduce
1. Full cell, stock parameter sets. Leak is the percentage drop in
Total lithium in particles [mol] over three 1C cycles (negative = lithium
created):
import pybamm
def leak(pset, mesh, n_r=16):
model = pybamm.lithium_ion.SPMe()
st = model.default_submesh_types
if mesh == "exponential":
for dom in ("negative particle", "positive particle"):
st[dom] = pybamm.MeshGenerator(
pybamm.Exponential1DSubMesh, submesh_params={"side": "right"}
)
sol = pybamm.Simulation(
model,
parameter_values=pybamm.ParameterValues(pset),
experiment=pybamm.Experiment(
["Discharge at 1C until 3.0V", "Charge at 1C until 4.2V"] * 3
),
var_pts={"x_n": 8, "x_s": 8, "x_p": 8, "r_n": n_r, "r_p": n_r},
submesh_types=st,
).solve(initial_soc=1.0)
n = sol["Total lithium in particles [mol]"](sol.t)
return 100 * (1 - n[-1] / n[0])
for pset in ("ORegan2022", "Ecker2015", "Chen2020"):
for n_r in (16, 64):
print(pset, n_r, leak(pset, "uniform", n_r), leak(pset, "exponential", n_r))| parameter set | r_n | uniform % | exponential % |
|---|---|---|---|
| ORegan2022 | 16 | -0.0000 | -25.9351 |
| ORegan2022 | 64 | 0.0000 | -0.3540 |
| Ecker2015 | 16 | 0.0000 | 0.4274 |
| Ecker2015 | 64 | -0.0000 | 0.0244 |
| Chen2020 (constant D) | 16 | 0.0000 | 0.0000 |
| Chen2020 (constant D) | 64 | -0.0000 | 0.0000 |
2. Minimal reproduction, no battery model. Pure diffusion in a sphere with a
constant surface flux and D linear in c, using the same BC form the Fickian
particle uses. The exact answer is d<c>/dt = -3*N_surf/R:
import numpy as np
import pybamm
R, j_over_F, T = 1e-5, 1e-7, 3600.0
for label, gen in [
("uniform", pybamm.Uniform1DSubMesh),
("exponential(right)", pybamm.MeshGenerator(
pybamm.Exponential1DSubMesh, submesh_params={"side": "right"})),
]:
for npts in (16, 48, 128):
model = pybamm.BaseModel()
c = pybamm.Variable("c", domain="negative particle")
D = 1e-14 * (1 + 9 * c / 1000.0)
model.rhs = {c: -pybamm.div(-D * pybamm.grad(c))}
model.boundary_conditions = {
c: {"left": (pybamm.Scalar(0), "Neumann"),
"right": (-j_over_F / pybamm.surf(D), "Neumann")}
}
model.initial_conditions = {c: pybamm.Scalar(500.0)}
model.variables = {"mean": pybamm.r_average(c)}
r = pybamm.SpatialVariable(
"r_n", domain=["negative particle"], coord_sys="spherical polar")
geom = {"negative particle": {r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(R)}}}
mesh = pybamm.Mesh(geom, {"negative particle": gen}, {r: npts})
pybamm.Discretisation(
mesh, {"negative particle": pybamm.FiniteVolume()}).process_model(model)
sol = pybamm.IDAKLUSolver(rtol=1e-10, atol=1e-10).solve(model, np.linspace(0, T, 50))
got = float(sol["mean"](T) - sol["mean"](0))
exact = -3 * j_over_F * T / R
print(label, npts, 100 * (got - exact) / exact)| mesh | npts=16 | npts=48 | npts=128 |
|---|---|---|---|
| uniform | 0.000 % | 0.000 % | 0.000 % |
| exponential(right) | -0.002 % | -0.000 % | -0.000 % |
Uniform is exact at every resolution; graded is not, and converges under
refinement — i.e. it is discretisation error, not integration error. Tightening
the solver from rtol=atol=1e-6 to 1e-9 changes the full-cell leak by less
than 0.001 points.
3. The two extrapolations, compared directly. This is the whole bug in one
number:
import numpy as np
import pybamm
R = 1e-5
r = pybamm.SpatialVariable("r_n", domain=["negative particle"], coord_sys="spherical polar")
geom = {"negative particle": {r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(R)}}}
for label, gen in [
("uniform", pybamm.Uniform1DSubMesh),
("exponential(right)", pybamm.MeshGenerator(
pybamm.Exponential1DSubMesh, submesh_params={"side": "right"})),
]:
for npts in (16, 64):
mesh = pybamm.Mesh(geom, {"negative particle": gen}, {r: npts})
fv = pybamm.FiniteVolume()
disc = pybamm.Discretisation(mesh, {"negative particle": fv})
c = pybamm.Variable("c", domain="negative particle")
disc.set_variable_slices([c])
D = 1e-14 * (1 + 99 * c / 1000.0)
y = 1000.0 * mesh["negative particle"].nodes / R
surf_val = float(np.ravel(disc.process_symbol(pybamm.surf(D)).evaluate(y=y))[-1])
edge_val = float(np.ravel(
fv.node_to_edge(disc.process_symbol(D)).evaluate(y=y))[-1])
print(label, npts, 100 * (edge_val / surf_val - 1))| mesh | npts | surf(D) | node->edge D | imposed flux error |
|---|---|---|---|---|
| uniform | 16 | 1.000000e-12 | 1.000000e-12 | 0.000 % |
| uniform | 64 | 1.000000e-12 | 1.000000e-12 | 0.000 % |
| exponential(right) | 16 | 1.000000e-12 | 1.000659e-12 | 0.066 % |
| exponential(right) | 64 | 1.000000e-12 | 1.000037e-12 | 0.004 % |
Suggested fix
Make arithmetic_mean's node-to-edge weights use the real node spacing, as
harmonic_mean already does, so that the interior weights are distance-weighted
and the boundary rows extrapolate from the true node positions. That would make
the interpolated edge value agree with pybamm.surf and restore conservation on
any mesh.
A narrower alternative is to have FickianDiffusion impose the flux directly
rather than as a gradient divided by surf(D_eff), which removes the dependence
on the two extrapolations agreeing.
Happy to open a PR for either.
Activity
rtimms commentedon Sep 16, 2026
Correction to the cause described above: I named
arithmetic_meanas the culprit, but that is only half of it and not the half that matters here.process_binary_operatorssends the operand to the edges with the harmonic mean when the other child is aGradient, i.e. exactly theD_eff * grad(c)flux case.harmonic_meanweights its interior edges correctly (beta = left_dx / (left_dx + right_dx)), but its two exterior edge rows are the same hard-coded[1.5, -0.5]/[-0.5, 1.5]asarithmetic_mean's.So the conservation bug is the exterior edge rows, in both means. Fixing
arithmetic_meanalone movesORegan2022from -25.94 % to only -23.33 %; fixing the exterior rows in both takes it to -0.0000 %.The interior
[0.5, 0.5]weights inarithmetic_meanare a real but separate first-order accuracy error on non-uniform meshes.Fix in #5755.