[Bug]: Half-cell and custom models do not conserve lithium

PyBaMM Version

26.4.1; 26.7.1.0

Python Version

3.12.3

Describe the bug

The finite_volume spatial method has been designed with Neumann boundary conditions in mind (i.e. fluxes that are directly proportional to gradients). By default, internal Neumann conditions are applied to concatenated gradients, instead of flux conditions. This means that lithium is typically not conserved in models with fluxes that contain functional parameters.

Currently, flux boundary conditions can only be applied (inexactly) via extrapolation. PR #5524 adds the ability to apply flux boundary conditions during the discretisation of the divergence, but unfortunately highlights some limitations of the finite volume implementation. To overcome these limitations, PyBaMM's DFN applies extra boundary conditions, but the BasicDFN demonstrates the issue: it does not conserve lithium when the transference number varies with electrolyte concentration. Half-cell ("working electrode": "positive") simulations also do not preserve lithium (after taking lithium exchange with the metal into account).

These issues have been reported before and should be solved at source. I have created a branch building on PR #5524, in which I have disabled the internal Neumann conditions and ensured that all gradient vectors discretise to the same length (as they do before the addition of flux BCs). Users should be advised to concatenate before applying the gradient operator, but for a custom model (with a constant current through the separator), I also needed a FluxConcatenation class to concatenate fluxes (without the circular application of edge_to_node/node_to_node that results in artificial smoothing of the intended flux). I have validated this approach for custom models in which the (constant) current in the separator determines the flux on the interface.

The remaining issue is that any source term, such as the (1 - t_plus(c_e, T)) * a_j / F in the electrolyte continuity equation, requires a systematic application of the finite volume method, rather than naive evaluation on nodes, to ensure conservation of lithium.

Thanks for reading! If you would find it useful, I'd be happy to discuss these issues on a video call.

Steps to Reproduce

import pybamm
from pybamm import Parameter

for model in [
    pybamm.lithium_ion.BasicDFN(),
    pybamm.lithium_ion.DFN({"working electrode": "positive"}),
]:
    # Parameter set with functional electrolyte parameters
    if model.options.get("working electrode") == "positive":
        parameter_values = pybamm.ParameterValues("Xu2019")
    else:
        parameter_values = pybamm.ParameterValues("ORegan2022")
        parameter_values["Negative particle diffusivity [m2.s-1]"] = 3.3e-14

        """The BasicDFN does not support a functional transference number"""
        # ce0 = parameter_values["Initial concentration in electrolyte [mol.m-3]"]
        # T = parameter_values["Ambient temperature [K]"]
        # parameter_values["Cation transference number"] = parameter_values[
        #     "Cation transference number"
        # ](ce0, T)

    # Define the total lithium, if necessary
    if "Total lithium [mol]" not in model.variables.keys():
        total_lithium_per_area = (
            pybamm.x_average(
                model.variables["Separator electrolyte concentration [mol.m-3]"]
                * Parameter("Separator porosity")
            )
            * pybamm.Parameter("Separator thickness [m]")
            + (
                pybamm.x_average(
                    model.variables["Positive electrolyte concentration [mol.m-3]"]
                    * Parameter("Positive electrode porosity")
                )
                + pybamm.x_average(
                    Parameter("Positive electrode active material volume fraction")
                    * pybamm.r_average(
                        model.variables["Positive particle concentration [mol.m-3]"]
                    )
                )
            )
            * pybamm.Parameter("Positive electrode thickness [m]")
        ) / pybamm.constants.F.value
        if model.options.get("working electrode") != "positive":
            total_lithium_per_area += (
                (
                    pybamm.x_average(
                        Parameter("Negative electrode active material volume fraction")
                        * pybamm.r_average(
                            model.variables["Negative particle concentration [mol.m-3]"]
                        )
                    )
                    + pybamm.x_average(
                        model.variables["Negative electrolyte concentration [mol.m-3]"]
                        * Parameter("Negative electrode porosity")
                    )
                )
                * Parameter("Negative electrode thickness [m]")
            ) / pybamm.constants.F.value
        model.variables["Total lithium [mol]"] = (
            total_lithium_per_area
            * pybamm.Parameter("Electrode width [m]")
            * pybamm.Parameter("Electrode height [m]")
            * pybamm.Parameter(
                "Number of electrodes connected in parallel to make a cell"
            )
        )

    # Define the quantity that should be conserved
    plot_variables = model.default_quick_plot_variables
    if model.options.get("working electrode") == "positive":
        model.variables["Lithium conservation [mol]"] = (
            model.variables["Total lithium [mol]"]
            - model.variables["Discharge capacity [A.h]"]
            * 3600
            / pybamm.constants.F.value
        )
        plot_variables.append("Total lithium [mol]")
        plot_variables.append("Lithium conservation [mol]")
    else:
        model.variables["Lithium conservation [mol]"] = model.variables[
            "Total lithium [mol]"
        ]
        plot_variables.append("Lithium conservation [mol]")

    # Define a test protocol
    init_soc = 0.45
    experiment = pybamm.Experiment(
        [
            "Rest for 10 minutes (5 seconds period)",
            "Discharge at 1C for 20 minutes (5 seconds period)",
            "Rest for 10 minutes (5 seconds period)",
            "Charge at 1C for 20 minutes (5 seconds period)",
            "Rest for 20 minutes (5 seconds period)",
        ],
    )

    # Run example simulation
    solution = pybamm.Simulation(
        model, parameter_values=parameter_values, experiment=experiment
    ).solve(initial_soc=init_soc)
    solution.plot(plot_variables)

Output

@aabills

@aabills not sure if this is in any way relevant to your work on higher dimension FVMs (the 1D version has not been edited for quite some time), otherwise @SarahRo can take it.

@SarahRo

Yeah, one thing I'd say (which I already told Ferran on slack) is that it might be worth considering deprecating all of the old spatial methods in favor of FiniteVolumeUnstructured (#5688), not sure how much of a lift it would be to add 1D, radial coordinates, and any other features we need and don't have right now.

Also if you have time @SarahRo and @NicolaCourtier could y'all look at the 2D and 3D spatial methods and se if this bug exists there too?

@aabills

@aabills I have used a modified version of the code @NicolaCourtier posted in the issue description to check the 2D version on the main branch (BasicDFN2D) and the version for the unstructured finite volume method (BasicDFN2DUnstructured - see PR #5690 ).

It looks like the bug is also present for those. For the code to run for the unstructured 2D DFN model I had to replace 'None' with something else in the following two lines of code:
https://github.com/pybamm-team/PyBaMM/blob/ufv-4-models/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_2d_unstructured.py#L215
https://github.com/pybamm-team/PyBaMM/blob/ufv-4-models/packages/pybamm/src/pybamm/models/full_battery_models/lithium_ion/basic_dfn_2d_unstructured.py#L216

Steps to reproduce

import pybamm
from pybamm import Parameter

for model in [
    pybamm.lithium_ion.BasicDFN2D(),
    # pybamm.lithium_ion.BasicDFN2DUnstructured()
]:
    # Parameter set with functional electrolyte parameters
    parameter_values = pybamm.ParameterValues("ORegan2022")
    parameter_values["Negative particle diffusivity [m2.s-1]"] = 3.3e-14


    # Define the total lithium, if necessary
    if "Total lithium [mol]" not in model.variables.keys():
        total_lithium_per_area = (
            pybamm.z_average(pybamm.x_average(
                model.variables["Separator electrolyte concentration [mol.m-3]"]
                * Parameter("Separator porosity")
            ))
            * pybamm.Parameter("Separator thickness [m]")
            + (
                pybamm.z_average(pybamm.x_average(
                    model.variables["Positive electrolyte concentration [mol.m-3]"]
                    * Parameter("Positive electrode porosity")
                )
                + pybamm.z_average(pybamm.x_average(
                    Parameter("Positive electrode active material volume fraction")
                    * pybamm.r_average(
                        model.variables["Positive particle concentration [mol.m-3]"]
                    )
                )))

            )
            * pybamm.Parameter("Positive electrode thickness [m]")
        ) / pybamm.constants.F.value
        total_lithium_per_area += (
            (
                pybamm.z_average(pybamm.x_average(
                    Parameter("Negative electrode active material volume fraction")
                    * pybamm.r_average(
                        model.variables["Negative particle concentration [mol.m-3]"]
                    )
                ))
                + pybamm.z_average(pybamm.x_average(
                    model.variables["Negative electrolyte concentration [mol.m-3]"]
                    * Parameter("Negative electrode porosity")
                ))
            )
            * Parameter("Negative electrode thickness [m]")
        ) / pybamm.constants.F.value
        model.variables["Total lithium [mol]"] = (
            total_lithium_per_area
            * pybamm.Parameter("Electrode width [m]")
            * pybamm.Parameter("Electrode height [m]")
            * pybamm.Parameter(
                "Number of electrodes connected in parallel to make a cell"
            )
        )

    # Define the quantity that should be conserved
    plot_variables = model.default_quick_plot_variables
    model.variables["Lithium conservation [mol]"] = model.variables[
        "Total lithium [mol]"
    ]
    plot_variables.append("Lithium conservation [mol]")

    # Define a test protocol
    init_soc = 0.45
    experiment = pybamm.Experiment(
        [
            "Rest for 10 minutes (5 seconds period)",
            "Discharge at 1C for 20 minutes (5 seconds period)",
            "Rest for 10 minutes (5 seconds period)",
            "Charge at 1C for 20 minutes (5 seconds period)",
            "Rest for 20 minutes (5 seconds period)",
        ],
    )

    # Run example simulation
    solution = pybamm.Simulation(
        model, parameter_values=parameter_values, experiment=experiment,
    ).solve(initial_soc=init_soc)

Ouput

I ran the two 1D cases from the description against v26.4.1 (d9e5afb), main (258fdc8) and the PR #5524 head (719763b), summing the lithium inventory directly from the finite-volume cell volumes: porosity × concentration × cell width for the electrolyte, active fraction × spherical-shell volume fraction for the particles, read at the mesh nodes. The sum uses PyBaMM's mesh, parameters and solution; it is independent only of the averaging and total-lithium diagnostics (x_average, r_average, Total lithium [mol]), which it does not call. Protocol as in the description, initial SOC 0.45, 20 cells per x-domain, 20 radial shells, IDAKLUSolver(rtol=atol=1e-8).

Observed (main and PR; v26.4.1 agrees to 1e-11 mol on BasicDFN)

Case	max abs ΔN, mol	relative to N(0)
BasicDFN, ORegan2022, D_n = 3.3e-14, variable t+	6.52e-3	2.44e-2
same, t+ fixed at its initial value	5.0e-15	1.9e-14
DFN full cell, same parameters	3.6e-15	1.4e-14
DFN `{"working electrode": "positive"}`, Xu2019, minus ∫I/F dt	1.64e-11	1.7e-7

BasicDFN refinement, max residual in mol: 6.51839e-3 at 20/20 for tol 1e-6, 1e-8 and 1e-10; 6.51612e-3 at 40/40; 6.51555e-3 at 80/80 (tol 1e-10). The drift does not go away with tolerance or mesh.
For all 24 configurations that completed on both main and the PR, the saved time series are byte-identical. Neither model declares a "Flux" boundary condition (BasicDFN: 11 Neumann, 1 Dirichlet; half-cell DFN: 9 Neumann, 1 Dirichlet), so as written they do not exercise the new code path.
Half-cell diagnostic: Total lithium in electrolyte [mol] is L_x * A * x_average(eps * c_e) with L_x = L_n + L_s + L_p (767 µm in Xu2019), but the electrolyte domain is separator + positive (67 µm). The built-in Total lithium [mol] equals N_particles + 11.448 × N_electrolyte (checked against the raw arrays to 4e-19 mol). With that diagnostic the half-cell residual is 1.80e-10 mol; with the direct sum it is 1.64e-11 mol.
Half-cell balance: the right-edge electrolyte flux is exactly zero, so dN/dt − I/F equals A × (left-edge electrolyte flux) − I/F. During the current steps this is a sustained bias of about 1.3e-14 mol/s (peak 7.1e-14) against I/F = 2.49e-8 mol/s at 20 cells, and about 3.0e-15 mol/s at 40 cells. Its sign follows the current.

Interpretation supported by the source

BasicDFN writes eps dc_e/dt = -div(-tor D grad c_e) + (1 - t+(c_e,T)) a_j/F. Integrating this with the particle equations over the cell gives dN/dt = −A ∫ t+(c_e) a_j/F dx exactly. The conservative form used by the modular Full electrolyte submodel, -div(-tor D grad c_e + t+ i_e/F) + a_j/F, differs by the product-rule term -i_e ∂x t+ / F, whose integral cancels that drift. The measured inventory change matches the time integral of −A ∫ t+ a_j/F dx to 1.2e-6 mol (5 s trapezoid), and fixing t+ removes it. So the 2.4% BasicDFN drift is a property of the reduced continuum equation for concentration-dependent t+, not of how the source is evaluated on nodes, and a boundary-condition change alone would not remove it.
The half-cell built-in diagnostic overweights the electrolyte by L_x/(L_s + L_p). Either Total lithium in electrolyte [mol] should use the length of whole_cell_domains when the negative electrode is absent, or the variable should be documented as full-cell only. Currently it is neither.

Not established here

Which extrapolation (D_e, tor or t+ at the boundary edge versus boundary_value in the Neumann expression) produces the half-cell left-edge flux bias. I did not decompose it.
Half-cell spatial convergence. 20/20 and 80/80 at tol 1e-10 fail in IDAKLU on main and PR (IDA_BAD_K, IDA_CONV_FAIL), so there is no fixed-tolerance mesh sequence. The 40/40 value is one point.
Whether the PR's flux boundary conditions fix these cases once the models are rewritten to use them. Nothing here tests that.
Any real-cell consequence. These are conservation checks on the default parameter sets only.

Minor: the reproducer in the description divides the concentration × length sum by F before multiplying by area, so its "mol" axis is scaled by 1/F. The relative drift is unaffected.

Inventory code, a wrapper that reproduces the table, the per-run outcomes for all 84 configurations, and the raw per-step arrays for the PR 20/20, 1e-8 runs: https://github.com/repowazdogz-droid/pybamm-5700-lithium-inventory (commit adaa877).

@repowazdogz-droid

Hi @repowazdogz-droid, thanks for sharing! Let me reply to a few of your points:

I ran the two 1D cases from the description against v26.4.1 (d9e5afb), main (258fdc8) and the PR

Did you see my branch here? It should work for constant t+.

Half-cell diagnostic: Total lithium in electrolyte [mol] is L_x * A * x_average(eps * c_e) with L_x = L_n + L_s + L_p (767 µm in Xu2019), but the electrolyte domain is separator + positive (67 µm).

This is true, and could be fixed as a separate issue.

So the 2.4% BasicDFN drift is a property of the reduced continuum equation for concentration-dependent t+, not of how the source is evaluated on nodes, and a boundary-condition change alone would not remove it.

Well done for identifying the missing -∂t+/∂x i_e / F in the BasicDFN or, put another way, the incorrect assumption that t+ i_e / F can be moved out of the flux (which only works for constant t+). Please also raise this as a separate issue.
I am concerned that there is still an underlying issue with the FV discretisation for a nonlinear source term (e.g. with a spatially varying active material volume fraction), but it should be at least first-order accurate.

Minor: the reproducer in the description divides the concentration × length sum by F before multiplying by area, so its "mol" axis is scaled by 1/F. The relative drift is unaffected.

True, good catch! Since the relative error is unaffected I'll leave the extraneous factor of F in the original post so the code matches the plots.

Thanks Nicola , that’s very helpful, particularly the distinction between the BasicDFN equation issue and the remaining FV discretisation question. I’ll raise the variable-t_+ BasicDFN term and the half-cell diagnostic separately as suggested. And agreed on the extra 1/F: it doesn’t affect the relative drift, so leaving the original reproducer unchanged makes sense.

Raised separately as suggested: #5745 (BasicDFN with a concentration-dependent t+) and #5746 (Total lithium in electrolyte [mol] using L_x in half-cell models).

Related: #5755

Now that the related conservation work has landed (including #5755 and #5765, plus the BasicDFN and unstructured-model changes), could we confirm whether these resolve the original scope of this issue? If so, I suggest closing #5700; otherwise, could the remaining reproducible defect be captured in a focused follow-up issue?

This issue is not yet resolved. If you're going to ask me to rewrite it, it might be better to do this after flux boundary conditions have been added.

Navigation Menu

Repository navigation

[Bug]: Half-cell and custom models do not conserve lithium #5700

PyBaMM Version

Python Version

Describe the bug

Steps to Reproduce

Output

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

[Bug]: Half-cell and custom models do not conserve lithium #5700

Description

PyBaMM Version

Python Version

Describe the bug

Steps to Reproduce

Output

Activity

brosaplanella commented on Aug 6, 2026

aabills commented on Aug 6, 2026

SarahRo commented on Aug 24, 2026

repowazdogz-droid commented on Sep 5, 2026

NicolaCourtier commented on Sep 6, 2026

repowazdogz-droid commented on Sep 6, 2026

repowazdogz-droid commented on Sep 6, 2026

rtimms commented on Sep 21, 2026

rtimms commented on Sep 25, 2026

NicolaCourtier commented on Sep 30, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions