[mleot]

Description

Fixes #5792

"heat of mixing": "true" could not be used with a composite (multi-phase) electrode.
BaseThermal._heat_of_mixing read the single-phase variable names and took its
material properties from param.n.prim / param.p.prim only, so any model built with
"particle phases": "2" failed at build time:

import pybamm

pybamm.lithium_ion.DFN(options={
    "particle phases": ("2", "1"),
    "thermal": "lumped",
    "heat of mixing": "true",
})

ModelError: Missing variable for submodel 'thermal': 'Negative particle radius [m]'.

The same failure occurred for SPM, SPMe and DFN, and for a composite positive as well as
a composite negative.

Changes

• Sum over phases. _phase_heat_of_mixing evaluates Richardson's term for one phase;
  _electrode_heat_of_mixing sums it over the electrode's phases, looked up through
  param.domain_params and options.phases.
• Effective diffusivity. The heat of mixing is the dissipation of the particle flux,
  so it now reads the particle's "effective diffusivity" variable instead of
  re-evaluating the base phase_param.D. They differ under "diffusivity": "current sigmoid" and "stress-induced diffusion": "true"; with lithiation and delithiation
  diffusivities of D/3 and 3D the heat sources were off the energy balance by 33% of the
  heat.
• Uniform profile. A uniform particle profile has no gradient, so its heat of mixing
  is zero outright. A DFN with "particle": "uniform profile" and heat of mixing now
  builds; it raised a shape error before.
• Which OCP. Unchanged: the equilibrium OCP, which the docstring now explains. The
  reaction and hysteresis heats sum to a j (phi_s - phi_e - U_eq), so the heat of
  mixing must use the same U_eq for energy to balance; the hysteresis-weighted OCP
  that sets the voltage would not.

Default results are unchanged. SPM, SPMe and DFN with "thermal": "lumped" and
"heat of mixing": "true" give bit-identical voltage, temperature and heat of mixing on
Chen2020 and Marquis2019 (1800 s, 1C).

Tests

• Energy balance (test_thermal_models.py, new): over a 1C SPM discharge, the fall
  in stored chemical energy -F * sum int U_eq dc must equal electrical work plus heat
  to 0.5% of the heat, for a one-phase cell and for a composite positive whose second
  phase has current-sigmoid hysteresis (branches tilted away from its OCP) and
  current-sigmoid diffusivity. The SPM keeps the electrolyte uniform and an isothermal
  model with zero entropic change keeps temperature out of the stored energy. Residuals
  are about -0.1% while the heat of mixing is ~45% of the heat. Summing only the primary
  phase fails (+28%), as do the base diffusivity (-18%) and the hysteresis-weighted OCP
  (-2%).
• Composite vs one phase (test_compare_outputs_two_phase.py, extending fix: sum the electrochemical heat sources over every particle phase #5770's
  test): a graphite electrode split into two identical halves must reproduce the
  one-phase heat sources, including the heat of mixing, for SPM, SPMe and DFN under
  every "thermal" option (isothermal heat source, lumped, x-full, x-lumped with 1+1D
  current collectors). It uses non-zero entropic change (Ai2020 for Chen2020;
  Marquis2019's own for the pouch cases), pins "surface form": "algebraic" on both
  sides (composite models default to it, one-phase models to "false") and solves at
  rtol 1e-6 so the ~300 K temperature state does not drift by more than the differences
  tested.
• Si/graphite smoke test: Chen2020_composite with a lumped model heats the cell
  through a positive heat of mixing.
• Well-posedness for "particle phases" "2", ("2", "1"), ("1", "2").

Not addressed here

• The heat of mixing differentiates the OCP symbolically, which fails for tabulated
  OCPs: a cubic Interpolant crashes in CasADi at solve time and a linear one at build
  time (issues to follow). The tests use analytic OCPs; Chen2020_composite's graphite
  OCP is swapped for the analytic Chen2020 fit where needed.
• Chen2020_composite has no "Negative electrode density [kg.m-3]", so no thermal
  model runs with it as shipped; the tests supply one (issue to follow).
• "thermal": "lumped" with 1+1D current collectors does not discretise, with or
  without the heat of mixing (issue to follow), so the 1+1D case uses "x-lumped".

Type of change

Bug fix: two entries under ## Bug fixes in CHANGELOG.md (phase summation; effective
diffusivity and uniform profile).

Important checks:

• No style issues: nox -s pre-commit (every hook passes on the branch diff,
  including comment-slop and repo-review)
• All tests pass: nox -s tests. Locally the lithium-ion model, thermal and
  two-phase tests pass (353 integration tests in test_lithium_ion/ and
  test_submodels/test_thermal/). The only local failures need optional dependencies
  my environment lacks (skfem, matplotlib, pybtex, bpx) and fail the same way
  on main; CI runs them.
• The documentation builds: nox -s doctests (not run locally; no docs changed,
  and the readthedocs check passes)
• Code is commented for hard-to-understand areas
• Tests added that prove fix is effective or that feature works

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.10%. Comparing base (3df83e8) to head (49e0254).
⚠️ Report is 2 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #5769      +/-   ##
==========================================
- Coverage   98.10%   98.10%   -0.01%     
==========================================
  Files         346      346              
  Lines       35482    35478       -4     
==========================================
- Hits        34809    34805       -4     
  Misses        673      673              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[BradyPlanden]

Thanks @mleot - this looks good. One suggestion to tighten the phase summing.

[BradyPlanden]

suggestion: the below would avoid hard-coding the ['primary', 'secondary'] and [prim, sec] while removing the multi-line branching.

Suggested change

	electrode_param = getattr(self.param, domain[0])
	num_phases = int(getattr(self.options, domain)["particle phases"])
	if num_phases == 1:
	return self._phase_heat_of_mixing(
	variables, domain, "", electrode_param.prim
	)
	return sum(
	self._phase_heat_of_mixing(variables, domain, phase_name, phase_param)
	for phase_name, phase_param in zip(
	["primary ", "secondary "],
	[electrode_param.prim, electrode_param.sec],
	strict=True,
	)
	)
	phases = self.options.phases[domain]
	names = [""] if len(phases) == 1 else [f"{p} " for p in phases]
	return sum(
	self._phase_heat_of_mixing(variables, domain, n, electrode_param.phase_params[p])
	for n, p in zip(names, phases, strict=True)
	)

[mleot]

Thanks @BradyPlanden, I think I addressed this in d469838. Now the phase names now come from options.phases, with "" for a one-phase electrode, so there's no branching for the phases now.

[mleot]

@BradyPlanden thanks for the review! Since the review I rebased onto main (after #5770); the heat of mixing now uses the particle's effective diffusivity, which it needs to conserve energy under current-sigmoid or stress-induced diffusion.

There are some addition tests including an energy-balance check and composite-vs-one-phase comparisons under every thermal option.

I think it's ready for another review when you get the chance.

[github-actions]

Bencher Report

Project	PyBaMM
Branch	heat-of-mixing-composite-electrodes
Testbed	bare-metal

⚠️ WARNING: Truncated view!

The full continuous benchmarking report exceeds the maximum length allowed on this platform.

🐰 View full continuous benchmarking report in Bencher

[rtimms]

Thanks @mleot, this looks good. Summing per phase through options.phases and param.domain_params is clean, and switching to the particle's effective diffusivity is the right call. The energy-balance test is a great addition: it fails if only the primary phase is summed, so it really tests the fix.

I've pushed two commits to get this ready for the release:

• 8057d1a84 merges main.
• 49e02543f moves both CHANGELOG bullets to # [Unreleased]. Since v26.9.0.0 was cut, git silently put them under the released section, with no conflict to flag it.

Checks I ran (model build and discretisation, main vs this branch):

• SPM, SPMe and DFN with lumped and x-full thermal, and stress-induced diffusion (Ai2020): unchanged.
• Composite "2" and ("2", "1"): now build (ModelError before).
• DFN with "particle": "uniform profile": now builds (ValueError before).
• MPM and size distributions: still raise the clear NotImplementedError.

I've filed the follow-ups mentioned in the description, plus one I found:

• #5876: heat of mixing with a quadratic or quartic particle profile gives a ShapeError. This is already broken on main, so it doesn't block this PR.
• #5877: Chen2020_composite gives the electrode density per phase, but the thermal model reads one value per electrode.
• #5878: "lumped" thermal doesn't discretise with 1+1D current collectors.

[rtimms]

@mleot I've merged main into your branch (14803e8), since this had picked up conflicts in the CHANGELOG and test_compare_outputs_two_phase.py. I kept your parametrised test_compare_heat_sources_graphite_graphite, kept main's new test_compare_graphite_graphite_differential_surface_form, and dropped main's old test_compare_heat_sources_DFN_graphite_graphite, which your parametrisation replaces. Both CHANGELOG bullets are back under # [Unreleased]. #5935 builds on your _phase_heat_of_mixing to fix #5876 (heat of mixing with polynomial profiles), so it'll follow this one.

[mleot]

Thanks for your review, support, and fixes for this @rtimms ! Let me know if I need to attend to anything else before it's good to go.