From bb24fb6e4daf4d52d2739f14bae56a3a50945d04 Mon Sep 17 00:00:00 2001
From: mleot <tleo@silanano.com>
Date: Mon, 21 Sep 2026 12:34:48 -0700
Subject: [PATCH 1/9] fix: sum the heat of mixing over every particle phase

`"heat of mixing": "true"` could not be used with a composite electrode.
`BaseThermal._heat_of_mixing` read the single-phase variable names
("Negative particle radius [m]", "X-averaged negative particle
concentration [mol.m-3]") and took its material properties from
`param.n.prim` / `param.p.prim` alone, so any model built with
`"particle phases": "2"` failed at build time with

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

Richardson's term is per-particle-population, so the fix is to evaluate it
per phase and add the contributions, which is what the reaction, reversible
and hysteresis heat sources in the same class already do. The phase loop
reuses the `phase_names = ["primary ", "secondary "]` idiom used there.

The single-phase path is untouched: it takes the `num_phases == 1` branch,
which builds the same expression from the same unprefixed variable names as
before. A DFN with Chen2020 and "heat of mixing": "true" returns voltage,
temperature and "Volume-averaged heat of mixing [W.m-3]" bit-identical to
the previous implementation.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 CHANGELOG.md                                  |  1 +
 .../models/submodels/thermal/base_thermal.py  | 94 +++++++++----------
 .../base_lithium_ion_tests.py                 | 22 +++++
 .../test_lithium_ion/test_newman_tobias.py    |  4 +
 4 files changed, 74 insertions(+), 47 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e6e78c74ee..83bed6ff1e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -17,6 +17,7 @@
 - `BasicDFN2D` now reports `"Total lithium [mol]"` and the `"Negative/Positive/Total solid lithium [mol]"` variables in mol for the whole cell: the integrals over `x` and `z` were in mol per metre of width, and are now scaled by the electrode width `L_y`. The volumetric interfacial current density `a * j`, previously mislabelled `"Current density [A.m-2]"`, is renamed to `"Sum of volumetric interfacial current densities [A.m-3]"` to match the full DFN. ([#5775](https://github.com/pybamm-team/PyBaMM/pull/5775))
 - The irreversible, reversible and hysteresis heat sources now sum over every particle phase. `"open-circuit potential"` stays a plain string when it is left at its default, so on a composite electrode it was zipped against the two phase names and silently truncated the loop after the primary phase: a graphite electrode split into two identical halves reported 56% of the irreversible heating of the equivalent single-phase electrode. Passing the option once per phase was a workaround. ([#5770](https://github.com/pybamm-team/PyBaMM/pull/5770))
 - Discretisation no longer mutates symbols it has already cached, so it no longer depends on the order outputs are processed in. An isothermal DFN with a `pybamm.Interpolant` of `pybamm.t` as `"Ambient temperature [K]"` and concentration-dependent electrolyte parameters raised a `ShapeError`. Copying a `StateVectorDot` also keeps it a `StateVectorDot`; `create_copy()` used to return a `StateVector`, so a copied time derivative read the state instead. ([#3630](https://github.com/pybamm-team/PyBaMM/issues/3630), [#4670](https://github.com/pybamm-team/PyBaMM/issues/4670), [#5771](https://github.com/pybamm-team/PyBaMM/pull/5771))
+- The heat of mixing now sums over every particle phase, so `"heat of mixing": "true"` works on a composite electrode. It previously read the single-phase variable names (`"Negative particle radius [m]"`) and `param.n.prim` only, so any model with `"particle phases": "2"` raised `ModelError: Missing variable for submodel 'thermal': 'Negative particle radius [m]'` at build time. Single-phase results are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
 - `pybamm.Function(scipy.special.erf, ...)` converts to CasADi again; CasADi 3.8 stopped dispatching scipy's `erf` ufunc onto symbolic values, so the generic `Function` path now maps it to `casadi.erf` explicitly. `pybamm.erf` was unaffected. ([#5761](https://github.com/pybamm-team/PyBaMM/pull/5761))
 - The dimensional minimum particle surface concentration is registered as `"Minimum <domain> particle surface concentration [mol.m-3]"`. An implicit string concatenation repeated its first line, so the variable was only reachable under the doubled name `"Minimum negative particle Minimum negative particle surface concentration [mol.m-3]"`, and the intended name resolved to nothing. The `Maximum` sibling was unaffected. ([#5759](https://github.com/pybamm-team/PyBaMM/pull/5759))
 - `FiniteVolume` node-to-edge shifts now use the true node spacing instead of weights that assume a uniform mesh, so the exterior edge values agree with `boundary_value`. On a non-uniform mesh they did not, which broke lithium conservation wherever a flux is imposed as a gradient divided by a surface value, as the Fickian particle does: with a concentration-dependent particle diffusivity on an `Exponential1DSubMesh`, `ORegan2022` gained 26% of its particle lithium over three 1C cycles. Uniform-mesh results are unchanged. ([#5755](https://github.com/pybamm-team/PyBaMM/pull/5755))
diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
index fa5e0b590f..4d601339b7 100644
--- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -383,65 +383,65 @@ def _current_collector_heating(self, variables):
                 Q_s_cp = self.param.p.sigma_cc * pybamm.grad_squared(phi_s_cp)
         return Q_s_cn, Q_s_cp
 
+    def _phase_heat_of_mixing(self, variables, domain, phase_name, phase_param):
+        """Compute the heat of mixing source term for a single particle phase."""
+        Domain = domain.capitalize()
+        F = pybamm.constants.F.value
+        a = variables[
+            f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]"
+        ]
+        R = variables[f"{Domain} {phase_name}particle radius [m]"]
+        N = a / (4 * np.pi * R**2)
+        if self.x_average:
+            c = variables[
+                f"X-averaged {domain} {phase_name}particle concentration [mol.m-3]"
+            ]
+            T = variables[f"X-averaged {domain} electrode temperature [K]"]
+        else:
+            c = variables[f"{Domain} {phase_name}particle concentration [mol.m-3]"]
+            T = variables[f"{Domain} electrode temperature [K]"]
+        T_part = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle"])
+        dc_dr2 = pybamm.inner(pybamm.grad(c), pybamm.grad(c))
+        D = phase_param.D(c, T_part)
+        dUeq = phase_param.U(c / phase_param.c_max, T_part).diff(c)
+        integrand_r = D * dc_dr2 * dUeq
+        integration_variable_r = [
+            pybamm.SpatialVariable("r", domain=integrand_r.domain)
+        ]
+        integral_r = pybamm.Integral(integrand_r, integration_variable_r)
+        return -F * N * integral_r
+
+    def _electrode_heat_of_mixing(self, variables, domain):
+        """Sum the heat of mixing over the particle phases of one electrode."""
+        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,
+            )
+        )
+
     def _heat_of_mixing(self, variables):
         """Compute heat of mixing source terms."""
 
         if self.options["heat of mixing"] == "true":
-            F = pybamm.constants.F.value
-            pi = np.pi
-
             # Compute heat of mixing in negative electrode
             if self.options.electrode_types["negative"] == "planar":
                 Q_mix_s_n = pybamm.FullBroadcast(
                     0, ["negative electrode"], "current collector"
                 )
             else:
-                a_n = variables["Negative electrode surface area to volume ratio [m-1]"]
-                R_n = variables["Negative particle radius [m]"]
-                N_n = a_n / (4 * pi * R_n**2)
-                if self.x_average:
-                    c_n = variables[
-                        "X-averaged negative particle concentration [mol.m-3]"
-                    ]
-                    T_n = variables["X-averaged negative electrode temperature [K]"]
-                else:
-                    c_n = variables["Negative particle concentration [mol.m-3]"]
-                    T_n = variables["Negative electrode temperature [K]"]
-                T_n_part = pybamm.PrimaryBroadcast(T_n, ["negative particle"])
-                dc_n_dr2 = pybamm.inner(pybamm.grad(c_n), pybamm.grad(c_n))
-                D_n = self.param.n.prim.D(c_n, T_n_part)
-                dUeq_n = self.param.n.prim.U(
-                    c_n / self.param.n.prim.c_max, T_n_part
-                ).diff(c_n)
-                integrand_r_n = D_n * dc_n_dr2 * dUeq_n
-                integration_variable_r_n = [
-                    pybamm.SpatialVariable("r", domain=integrand_r_n.domain)
-                ]
-                integral_r_n = pybamm.Integral(integrand_r_n, integration_variable_r_n)
-                Q_mix_s_n = -F * N_n * integral_r_n
+                Q_mix_s_n = self._electrode_heat_of_mixing(variables, "negative")
 
             # Compute heat of mixing in positive electrode
-            a_p = variables["Positive electrode surface area to volume ratio [m-1]"]
-            R_p = variables["Positive particle radius [m]"]
-            N_p = a_p / (4 * pi * R_p**2)
-            if self.x_average:
-                c_p = variables["X-averaged positive particle concentration [mol.m-3]"]
-                T_p = variables["X-averaged positive electrode temperature [K]"]
-            else:
-                c_p = variables["Positive particle concentration [mol.m-3]"]
-                T_p = variables["Positive electrode temperature [K]"]
-            T_p_part = pybamm.PrimaryBroadcast(T_p, ["positive particle"])
-            dc_p_dr2 = pybamm.inner(pybamm.grad(c_p), pybamm.grad(c_p))
-            D_p = self.param.p.prim.D(c_p, T_p_part)
-            dUeq_p = self.param.p.prim.U(c_p / self.param.p.prim.c_max, T_p_part).diff(
-                c_p
-            )
-            integrand_r_p = D_p * dc_p_dr2 * dUeq_p
-            integration_variable_r_p = [
-                pybamm.SpatialVariable("r", domain=integrand_r_p.domain)
-            ]
-            integral_r_p = pybamm.Integral(integrand_r_p, integration_variable_r_p)
-            Q_mix_s_p = -F * N_p * integral_r_p
+            Q_mix_s_p = self._electrode_heat_of_mixing(variables, "positive")
             Q_mix_s_s = pybamm.FullBroadcast(0, ["separator"], "current collector")
         else:
             Q_mix_s_n = pybamm.FullBroadcast(
diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index 965f5b8949..ca5f9124c3 100644
--- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -533,6 +533,28 @@ def test_well_posed_particle_phases_thermal(self):
         options = {"particle phases": "2", "thermal": "lumped"}
         self.check_well_posedness(options)
 
+    def test_well_posed_particle_phases_heat_of_mixing(self):
+        options = {
+            "particle phases": "2",
+            "thermal": "lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+        options = {
+            "particle phases": ("2", "1"),
+            "thermal": "lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+        options = {
+            "particle phases": ("1", "2"),
+            "thermal": "lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_particle_phases_sei(self):
         options = {"particle phases": "2", "SEI": "ec reaction limited"}
         self.check_well_posedness(options)
diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py
index 306be41d8d..013ed24f6d 100644
--- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py
+++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py
@@ -24,6 +24,10 @@ def test_well_posed_particle_phases(self):
     def test_well_posed_particle_phases_thermal(self):
         pass  # Skip this test
 
+    @pytest.mark.skip(reason="Test currently not implemented")
+    def test_well_posed_particle_phases_heat_of_mixing(self):
+        pass  # skip this test
+
     @pytest.mark.skip(reason="Test currently not implemented")
     def test_well_posed_particle_phases_sei(self):
         pass  # skip this test

From 4ba0a9e6f3e126c3859d21d124e5cb7ece0ef532 Mon Sep 17 00:00:00 2001
From: tleo <tleo@silanano.com>
Date: Thu, 24 Sep 2026 18:32:02 -0700
Subject: [PATCH 2/9] test: check composite heat of mixing values, not just
 well-posedness

The composite heat-of-mixing tests only built the model. Extend the
graphite/graphite heat-source comparison from #5770 so that splitting one
graphite phase into two identical halves must also reproduce the one-phase
"Volume-averaged heat of mixing [W.m-3]":

- run it with "heat of mixing": "true", for isothermal (heat source only)
  and lumped thermal models, comparing the temperature rise for the latter
- give Chen2020 the Ai2020 entropic changes; with zero dU/dT the reversible
  heat and the temperature dependence of U inside the heat of mixing were
  compared as zero against zero
- solve at rtol 1e-6: at the default rtol the ~300 K temperature state
  drifts by ~0.03 K, which is 1% of the rise and failed SPM + lumped

A heat of mixing that sums only the primary phase fails all four cases.

Add a Si/graphite smoke test (Chen2020_composite, lumped) checking the heat
of mixing stays positive and heats the cell. It supplies the per-layer
"Negative electrode density [kg.m-3]" the set lacks and swaps the cubic
Interpolant graphite OCP for the analytic Chen2020 fit, whose derivative
the solver can convert.

Also look up the electrode and phase parameters through
param.domain_params and options.phases instead of getattr.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
---
 .../models/submodels/thermal/base_thermal.py  |  14 +--
 .../test_compare_outputs_two_phase.py         | 112 ++++++++++++++++--
 2 files changed, 107 insertions(+), 19 deletions(-)

diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
index 4d601339b7..90fbfb27c8 100644
--- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -413,19 +413,17 @@ def _phase_heat_of_mixing(self, variables, domain, phase_name, phase_param):
 
     def _electrode_heat_of_mixing(self, variables, domain):
         """Sum the heat of mixing over the particle phases of one electrode."""
-        electrode_param = getattr(self.param, domain[0])
-        num_phases = int(getattr(self.options, domain)["particle phases"])
-        if num_phases == 1:
+        electrode_param = self.param.domain_params[domain]
+        phases = self.options.phases[domain]
+        if len(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,
+            self._phase_heat_of_mixing(
+                variables, domain, f"{phase} ", electrode_param.phase_params[phase]
             )
+            for phase in phases
         )
 
     def _heat_of_mixing(self, variables):
diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
index dea71aca7a..ab67b8cf0f 100644
--- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
+++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
@@ -177,24 +177,38 @@ def test_compare_DFN_silicon_graphite(self):
         model_class = pybamm.lithium_ion.DFN
         self.compare_outputs_two_phase_silicon_graphite(model_class)
 
-    def compare_heat_sources_two_phase_graphite_graphite(self, model_class):
+    def compare_heat_sources_two_phase_graphite_graphite(self, model_class, thermal):
         """
         Every heat source is a sum over the particle phases, so splitting one graphite
         phase into two identical halves must reproduce the one-phase heat sources
         """
-        options = {
-            "thermal": "isothermal",
-            "calculate heat source for isothermal models": "true",
-        }
+        if thermal == "isothermal":
+            options = {
+                "thermal": "isothermal",
+                "calculate heat source for isothermal models": "true",
+            }
+        else:
+            options = {"thermal": thermal}
+        options["heat of mixing"] = "true"
         t_eval = [0, 3600]
         t_interp = np.linspace(0, 3600)
 
+        # Chen2020 has zero entropic change, which would leave the reversible heat
+        # and the temperature dependence of the heat of mixing untested
+        entropic_change = pybamm.ParameterValues("Ai2020")
         parameter_values = pybamm.ParameterValues("Chen2020")
+        for domain in ["Negative", "Positive"]:
+            name = f"{domain} electrode OCP entropic change [V.K-1]"
+            parameter_values[name] = entropic_change[name]
+
+        # the default rtol lets a temperature state near 300 K drift by ~0.03 K,
+        # which is comparable to the one- and two-phase differences being tested
+        solver = pybamm.IDAKLUSolver(rtol=1e-6, atol=1e-8)
         sol = pybamm.Simulation(
-            model_class(options), parameter_values=parameter_values
+            model_class(options), parameter_values=parameter_values, solver=solver
         ).solve(t_eval=t_eval, t_interp=t_interp)
 
-        parameter_values_two_phase = pybamm.ParameterValues("Chen2020")
+        parameter_values_two_phase = parameter_values.copy()
         for parameter in [
             "Negative electrode OCP [V]",
             "Negative electrode OCP entropic change [V.K-1]",
@@ -216,24 +230,100 @@ def compare_heat_sources_two_phase_graphite_graphite(self, model_class):
         sol_two_phase = pybamm.Simulation(
             model_class({"particle phases": ("2", "1"), **options}),
             parameter_values=parameter_values_two_phase,
+            solver=solver,
         ).solve(t_eval=t_eval, t_interp=t_interp)
 
         for variable in [
             "Volume-averaged irreversible electrochemical heating [W.m-3]",
             "Volume-averaged reversible heating [W.m-3]",
             "Volume-averaged hysteresis electrochemical heating [W.m-3]",
+            "Volume-averaged heat of mixing [W.m-3]",
             "Volume-averaged total heating [W.m-3]",
             "Voltage [V]",
         ]:
+            one_phase = sol[variable](t_interp)
+            # scale atol to the variable, since the reversible heat changes sign
             np.testing.assert_allclose(
-                sol[variable](t_interp),
                 sol_two_phase[variable](t_interp),
+                one_phase,
                 rtol=1e-2,
-                atol=1e-8,
+                atol=1e-3 * np.nanmax(np.abs(one_phase)) + 1e-8,
+            )
+
+        if thermal != "isothermal":
+            temperature = "Volume-averaged cell temperature [K]"
+            temperature_rise = sol[temperature](t_interp) - sol[temperature](0)
+            temperature_rise_two_phase = sol_two_phase[temperature](
+                t_interp
+            ) - sol_two_phase[temperature](0)
+            np.testing.assert_allclose(
+                temperature_rise_two_phase,
+                temperature_rise,
+                rtol=1e-2,
+                atol=1e-3 * np.nanmax(temperature_rise),
             )
 
     def test_compare_heat_sources_SPM_graphite_graphite(self):
-        self.compare_heat_sources_two_phase_graphite_graphite(pybamm.lithium_ion.SPM)
+        self.compare_heat_sources_two_phase_graphite_graphite(
+            pybamm.lithium_ion.SPM, "isothermal"
+        )
 
     def test_compare_heat_sources_DFN_graphite_graphite(self):
-        self.compare_heat_sources_two_phase_graphite_graphite(pybamm.lithium_ion.DFN)
+        self.compare_heat_sources_two_phase_graphite_graphite(
+            pybamm.lithium_ion.DFN, "isothermal"
+        )
+
+    def test_compare_heat_sources_SPM_graphite_graphite_lumped(self):
+        self.compare_heat_sources_two_phase_graphite_graphite(
+            pybamm.lithium_ion.SPM, "lumped"
+        )
+
+    def test_compare_heat_sources_DFN_graphite_graphite_lumped(self):
+        self.compare_heat_sources_two_phase_graphite_graphite(
+            pybamm.lithium_ion.DFN, "lumped"
+        )
+
+    def heat_of_mixing_silicon_graphite(self, model_class):
+        # Check the heat of mixing of two genuinely different phases heats the cell
+        options = {
+            "particle phases": ("2", "1"),
+            "open-circuit potential": (("single", "current sigmoid"), "single"),
+            "thermal": "lumped",
+        }
+        parameter_values = pybamm.ParameterValues("Chen2020_composite")
+        # thermal properties are per layer, but the set only has per-phase densities
+        parameter_values.update(
+            {"Negative electrode density [kg.m-3]": 1657.0},
+            check_already_exists=False,
+        )
+        # the shipped graphite OCP is a cubic Interpolant, whose derivative cannot
+        # yet be converted for the solver, so use the analytic Chen2020 fit
+        parameter_values["Primary: Negative electrode OCP [V]"] = (
+            pybamm.ParameterValues("Chen2020")["Negative electrode OCP [V]"]
+        )
+        t_eval = [0, 3600]
+        t_interp = np.linspace(0, 3600)
+
+        solutions = {
+            heat_of_mixing: pybamm.Simulation(
+                model_class({**options, "heat of mixing": heat_of_mixing}),
+                parameter_values=parameter_values,
+            ).solve(t_eval=t_eval, t_interp=t_interp)
+            for heat_of_mixing in ["false", "true"]
+        }
+
+        heat_of_mixing = solutions["true"]["Volume-averaged heat of mixing [W.m-3]"](
+            t_interp
+        )
+        heat_of_mixing = heat_of_mixing[~np.isnan(heat_of_mixing)]
+        assert np.all(heat_of_mixing > 0)
+        temperature = "Volume-averaged cell temperature [K]"
+        assert solutions["true"][temperature](t_interp[-2]) > solutions["false"][
+            temperature
+        ](t_interp[-2])
+
+    def test_heat_of_mixing_SPM_silicon_graphite(self):
+        self.heat_of_mixing_silicon_graphite(pybamm.lithium_ion.SPM)
+
+    def test_heat_of_mixing_DFN_silicon_graphite(self):
+        self.heat_of_mixing_silicon_graphite(pybamm.lithium_ion.DFN)

From 016fc432c40f6bf6b334fbccb0ff0008f1b1d6d8 Mon Sep 17 00:00:00 2001
From: tleo <tleo@silanano.com>
Date: Thu, 24 Sep 2026 18:38:49 -0700
Subject: [PATCH 3/9] test: compare composite heat sources under every thermal
 option

Parametrise the graphite/graphite heat-source comparison over each value
of the "thermal" option, for SPM and DFN:

- isothermal (heat source only) and lumped, on Chen2020 with the Ai2020
  entropic changes, as before
- x-full, which resolves the temperature in x, so the heat of mixing reads
  the local electrode temperature instead of the x-averaged one
- x-lumped with 1+1D current collectors

The pouch cases use Marquis2019, which ships the tab and heat transfer
parameters, a non-zero entropic change and an analytic, temperature
dependent particle diffusivity. A heat of mixing that sums only the primary
phase fails all eight cases.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
---
 .../test_compare_outputs_two_phase.py         | 90 ++++++++++++-------
 1 file changed, 56 insertions(+), 34 deletions(-)

diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
index ab67b8cf0f..67a929ccf4 100644
--- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
+++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
@@ -2,6 +2,7 @@
 # Tests for the surface formulation
 #
 import numpy as np
+import pytest
 
 import pybamm
 
@@ -177,29 +178,25 @@ def test_compare_DFN_silicon_graphite(self):
         model_class = pybamm.lithium_ion.DFN
         self.compare_outputs_two_phase_silicon_graphite(model_class)
 
-    def compare_heat_sources_two_phase_graphite_graphite(self, model_class, thermal):
+    def compare_heat_sources_two_phase_graphite_graphite(
+        self, model_class, thermal_options, parameter_set="Chen2020"
+    ):
         """
         Every heat source is a sum over the particle phases, so splitting one graphite
         phase into two identical halves must reproduce the one-phase heat sources
         """
-        if thermal == "isothermal":
-            options = {
-                "thermal": "isothermal",
-                "calculate heat source for isothermal models": "true",
-            }
-        else:
-            options = {"thermal": thermal}
-        options["heat of mixing"] = "true"
+        options = {**thermal_options, "heat of mixing": "true"}
         t_eval = [0, 3600]
         t_interp = np.linspace(0, 3600)
 
-        # Chen2020 has zero entropic change, which would leave the reversible heat
-        # and the temperature dependence of the heat of mixing untested
-        entropic_change = pybamm.ParameterValues("Ai2020")
-        parameter_values = pybamm.ParameterValues("Chen2020")
-        for domain in ["Negative", "Positive"]:
-            name = f"{domain} electrode OCP entropic change [V.K-1]"
-            parameter_values[name] = entropic_change[name]
+        parameter_values = pybamm.ParameterValues(parameter_set)
+        if parameter_set == "Chen2020":
+            # Chen2020 has zero entropic change, which would leave the reversible
+            # heat and the temperature dependence of the heat of mixing untested
+            entropic_change = pybamm.ParameterValues("Ai2020")
+            for domain in ["Negative", "Positive"]:
+                name = f"{domain} electrode OCP entropic change [V.K-1]"
+                parameter_values[name] = entropic_change[name]
 
         # the default rtol lets a temperature state near 300 K drift by ~0.03 K,
         # which is comparable to the one- and two-phase differences being tested
@@ -250,7 +247,7 @@ def compare_heat_sources_two_phase_graphite_graphite(self, model_class, thermal)
                 atol=1e-3 * np.nanmax(np.abs(one_phase)) + 1e-8,
             )
 
-        if thermal != "isothermal":
+        if options["thermal"] != "isothermal":
             temperature = "Volume-averaged cell temperature [K]"
             temperature_rise = sol[temperature](t_interp) - sol[temperature](0)
             temperature_rise_two_phase = sol_two_phase[temperature](
@@ -263,24 +260,49 @@ def compare_heat_sources_two_phase_graphite_graphite(self, model_class, thermal)
                 atol=1e-3 * np.nanmax(temperature_rise),
             )
 
-    def test_compare_heat_sources_SPM_graphite_graphite(self):
-        self.compare_heat_sources_two_phase_graphite_graphite(
-            pybamm.lithium_ion.SPM, "isothermal"
-        )
-
-    def test_compare_heat_sources_DFN_graphite_graphite(self):
-        self.compare_heat_sources_two_phase_graphite_graphite(
-            pybamm.lithium_ion.DFN, "isothermal"
-        )
-
-    def test_compare_heat_sources_SPM_graphite_graphite_lumped(self):
-        self.compare_heat_sources_two_phase_graphite_graphite(
-            pybamm.lithium_ion.SPM, "lumped"
-        )
-
-    def test_compare_heat_sources_DFN_graphite_graphite_lumped(self):
+    @pytest.mark.parametrize(
+        ("model_class", "thermal_options", "parameter_set"),
+        [
+            pytest.param(model_class, options, parameter_set, id=f"{name}-{label}")
+            for name, model_class in [
+                ("SPM", pybamm.lithium_ion.SPM),
+                ("DFN", pybamm.lithium_ion.DFN),
+            ]
+            for label, options, parameter_set in [
+                (
+                    "isothermal",
+                    {
+                        "thermal": "isothermal",
+                        "calculate heat source for isothermal models": "true",
+                    },
+                    "Chen2020",
+                ),
+                ("lumped", {"thermal": "lumped"}, "Chen2020"),
+                # x-full resolves the temperature in x, so the heat of mixing uses
+                # the local rather than the x-averaged electrode temperature
+                (
+                    "x-full",
+                    {"thermal": "x-full", "cell geometry": "pouch"},
+                    "Marquis2019",
+                ),
+                (
+                    "x-lumped-1plus1D",
+                    {
+                        "thermal": "x-lumped",
+                        "cell geometry": "pouch",
+                        "current collector": "potential pair",
+                        "dimensionality": 1,
+                    },
+                    "Marquis2019",
+                ),
+            ]
+        ],
+    )
+    def test_compare_heat_sources_graphite_graphite(
+        self, model_class, thermal_options, parameter_set
+    ):
         self.compare_heat_sources_two_phase_graphite_graphite(
-            pybamm.lithium_ion.DFN, "lumped"
+            model_class, thermal_options, parameter_set
         )
 
     def heat_of_mixing_silicon_graphite(self, model_class):

From 08d401b82aaf2cef2a8b003db0b9842df4f4cea3 Mon Sep 17 00:00:00 2001
From: tleo <tleo@silanano.com>
Date: Thu, 24 Sep 2026 18:41:06 -0700
Subject: [PATCH 4/9] test: pin the surface form in the composite heat-source
 comparison

Composite models default to "surface form": "algebraic" while one-phase
models default to "false", so the comparison was between two electrolyte
formulations. Pin "algebraic" for both, which makes it like with like and
lets SPMe join SPM and DFN: with the defaults SPMe differs by 1.4e-3 in
voltage from the formulation change alone.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
---
 .../test_lithium_ion/test_compare_outputs_two_phase.py   | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
index 67a929ccf4..57dbe40c36 100644
--- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
+++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py
@@ -185,7 +185,13 @@ def compare_heat_sources_two_phase_graphite_graphite(
         Every heat source is a sum over the particle phases, so splitting one graphite
         phase into two identical halves must reproduce the one-phase heat sources
         """
-        options = {**thermal_options, "heat of mixing": "true"}
+        # composite models default to the algebraic surface form, so pin it for the
+        # one-phase model too and compare like with like
+        options = {
+            **thermal_options,
+            "heat of mixing": "true",
+            "surface form": "algebraic",
+        }
         t_eval = [0, 3600]
         t_interp = np.linspace(0, 3600)
 
@@ -266,6 +272,7 @@ def compare_heat_sources_two_phase_graphite_graphite(
             pytest.param(model_class, options, parameter_set, id=f"{name}-{label}")
             for name, model_class in [
                 ("SPM", pybamm.lithium_ion.SPM),
+                ("SPMe", pybamm.lithium_ion.SPMe),
                 ("DFN", pybamm.lithium_ion.DFN),
             ]
             for label, options, parameter_set in [

From dfd3e57a1c232f89afd2b47791b677488a56ddf2 Mon Sep 17 00:00:00 2001
From: tleo <tleo@silanano.com>
Date: Thu, 24 Sep 2026 18:56:10 -0700
Subject: [PATCH 5/9] test: check the heat sources close the energy balance
 with heat of mixing

pybamm had no check that its heat sources conserve energy. The reaction and
hysteresis heats sum to a j (phi_s - phi_e - U_eq) at the particle surface,
so the stored chemical energy is -F * int U_eq dc over the particles, and
integrating by parts leaves exactly the heat of mixing taken with U_eq.
The new test integrates a 1C SPM discharge and checks that the fall in that
stored energy equals electrical work plus heat to 0.5% of the heat, for a
one-phase cell and for a composite positive electrode whose second phase
has current-sigmoid hysteresis with branches tilted away from its OCP. The
SPM keeps the electrolyte uniform and an isothermal model with zero
entropic change keeps temperature out of the stored energy, so the
particles are the only store.

Residuals are about -0.1%. The heat of mixing is ~45% of the heat, and two
plausible mistakes fail: summing only the primary phase (+19%) and taking
the heat of mixing from the hysteresis-weighted branch instead of U_eq
(-1.4%). The heat-of-mixing docstring now says it must share U_eq with the
hysteresis heat.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
---
 .../models/submodels/thermal/base_thermal.py  |   8 +-
 .../test_lithium_ion/test_thermal_models.py   | 119 ++++++++++++++++++
 2 files changed, 126 insertions(+), 1 deletion(-)

diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
index 90fbfb27c8..0c615883b0 100644
--- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -384,7 +384,13 @@ def _current_collector_heating(self, variables):
         return Q_s_cn, Q_s_cp
 
     def _phase_heat_of_mixing(self, variables, domain, phase_name, phase_param):
-        """Compute the heat of mixing source term for a single particle phase."""
+        """
+        Compute the heat of mixing source term for a single particle phase.
+
+        The OCP is the phase's equilibrium OCP, the same one the hysteresis heat is
+        referenced to, not the hysteresis-weighted OCP that sets the voltage. Energy
+        only balances if both heat sources share that reference.
+        """
         Domain = domain.capitalize()
         F = pybamm.constants.F.value
         a = variables[
diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
index 7f2609c4e8..5be0565c3f 100644
--- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
+++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
@@ -222,3 +222,122 @@ def R_contact(T):
             # cell heats up above the reference, so R (and hence the contact
             # overpotential magnitude) is larger than the constant case
             np.testing.assert_array_less(np.abs(dphi_const[1:]), np.abs(dphi_fn[1:]))
+
+    @pytest.mark.parametrize("composite", [False, True], ids=["one-phase", "composite"])
+    def test_heat_of_mixing_energy_balance(self, composite):
+        """
+        Electrical work plus heat must equal the fall in stored chemical energy,
+        -F * sum over the particles of int_0^c U_eq dc, where U_eq is the OCP the
+        hysteresis heat is referenced to. The SPM keeps the electrolyte uniform and
+        the isothermal model with zero entropic change keeps the stored energy
+        independent of temperature, so the particles are the only energy store
+        """
+        options = {
+            "thermal": "isothermal",
+            "calculate heat source for isothermal models": "true",
+            "heat of mixing": "true",
+        }
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        if composite:
+            # Split the positive electrode, which carries most of the heat of mixing,
+            # into two identical halves, and give the second hysteresis about its OCP.
+            # The branches are tilted so that their slopes differ from the OCP's
+            options["particle phases"] = ("1", "2")
+            options["open-circuit potential"] = (
+                "single",
+                ("single", "current sigmoid"),
+            )
+            for name in [
+                "Positive electrode OCP [V]",
+                "Positive electrode OCP entropic change [V.K-1]",
+                "Maximum concentration in positive electrode [mol.m-3]",
+                "Initial concentration in positive electrode [mol.m-3]",
+                "Positive particle radius [m]",
+                "Positive particle diffusivity [m2.s-1]",
+                "Positive electrode exchange-current density [A.m-2]",
+                "Positive electrode active material volume fraction",
+            ]:
+                value = parameter_values[name]
+                if name.endswith("volume fraction"):
+                    value = value / 2
+                parameter_values.update(
+                    {f"Primary: {name}": value, f"Secondary: {name}": value}
+                )
+                del parameter_values[name]
+            mean_ocp = parameter_values["Secondary: Positive electrode OCP [V]"]
+            parameter_values.update(
+                {
+                    "Secondary: Positive electrode lithiation OCP [V]": lambda sto: (
+                        mean_ocp(sto) - 0.005 - 0.1 * sto
+                    ),
+                    "Secondary: Positive electrode delithiation OCP [V]": lambda sto: (
+                        mean_ocp(sto) + 0.005 + 0.1 * sto
+                    ),
+                },
+                check_already_exists=False,
+            )
+            particles = [
+                ("negative", ""),
+                ("positive", "primary "),
+                ("positive", "secondary "),
+            ]
+        else:
+            particles = [("negative", ""), ("positive", "")]
+
+        sim = pybamm.Simulation(
+            pybamm.lithium_ion.SPM(options),
+            parameter_values=parameter_values,
+            solver=pybamm.IDAKLUSolver(rtol=1e-9, atol=1e-11),
+        )
+        solution = sim.solve([0, 3000])
+        t = np.linspace(0, solution.t[-1], 2001)
+
+        F = pybamm.constants.F.value
+        area = (
+            parameter_values["Electrode width [m]"]
+            * parameter_values["Electrode height [m]"]
+        )
+        stored_energy = 0
+        for domain, phase in particles:
+            Domain = domain.capitalize()
+            prefix = phase.strip().capitalize() + ": " if phase else ""
+            c_max = parameter_values[
+                f"{prefix}Maximum concentration in {domain} electrode [mol.m-3]"
+            ]
+            # int_0^sto U_eq, avoiding the OCP asymptotes at the ends
+            sto_table = np.linspace(1e-4, 1 - 1e-4, 4001)
+            ocp = parameter_values[f"{prefix}{Domain} electrode OCP [V]"]
+            ocp_table = np.ravel(ocp(pybamm.Vector(sto_table)).evaluate())
+            integral_table = np.concatenate(
+                [
+                    [0],
+                    np.cumsum(
+                        (ocp_table[1:] + ocp_table[:-1]) / 2 * np.diff(sto_table)
+                    ),
+                ]
+            )
+            c = solution[
+                f"X-averaged {domain} {phase}particle concentration [mol.m-3]"
+            ](t=t)
+            energy_density = (
+                -F * c_max * np.interp(c / c_max, sto_table, integral_table)
+            )
+            # finite volume cells are spherical shells
+            edges = sim.mesh[f"{domain} {phase}particle"].edges
+            shell_fraction = (edges[1:] ** 3 - edges[:-1] ** 3) / edges[-1] ** 3
+            solid_volume = (
+                parameter_values[
+                    f"{prefix}{Domain} electrode active material volume fraction"
+                ]
+                * parameter_values[f"{Domain} electrode thickness [m]"]
+                * area
+            )
+            stored_energy += solid_volume * (shell_fraction @ energy_density)
+
+        released = stored_energy[0] - stored_energy[-1]
+        work = np.trapezoid(solution["Current [A]"](t) * solution["Voltage [V]"](t), t)
+        heat = np.trapezoid(solution["Total heating [W]"](t), t)
+        heat_of_mixing = np.trapezoid(solution["Heat of mixing [W]"](t), t)
+        np.testing.assert_allclose(released, work + heat, rtol=0, atol=5e-3 * heat)
+        # without the heat of mixing the balance fails, so the check has teeth
+        assert abs(released - (work + heat - heat_of_mixing)) > 0.2 * heat

From 20fdd0151745dd565fa2c6e348c62d6456c75f82 Mon Sep 17 00:00:00 2001
From: tleo <tleo@silanano.com>
Date: Thu, 24 Sep 2026 19:16:13 -0700
Subject: [PATCH 6/9] fix: take the heat of mixing's diffusivity from the
 particle flux

The heat of mixing is the dissipation of the particle flux,
-N . grad(mu) = F D |grad c|^2 dU_eq/dc, so it needs the D the flux uses.
It re-evaluated the base `phase_param.D(c, T)` instead, which ignores the
"current sigmoid" diffusivity branches and the stress-induced diffusion
factor. Read the particle's "effective diffusivity" variable, x-averaged
for the x-averaged models, as the particle submodels publish it.

With lithiation and delithiation diffusivities of D/3 and 3D on one phase
the heat sources missed the energy balance by 33% of the heat. The
energy-balance test now gives the composite positive's hysteretic phase a
current-sigmoid diffusivity (D/2, 2D): the old code fails it at -18%, the
new code passes at -0.2%.

A uniform particle profile publishes no diffusivity and has no gradient,
so its heat of mixing is now zero outright. That also makes a DFN with
"particle": "uniform profile" and heat of mixing build; it raised a shape
error before.

Default options give bit-identical voltage, temperature and heat of mixing
for SPM, SPMe and DFN on Chen2020 and Marquis2019.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md                                      |  1 +
 .../models/submodels/thermal/base_thermal.py      | 14 +++++++++++++-
 .../test_lithium_ion/test_thermal_models.py       | 15 +++++++++++++--
 3 files changed, 27 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 83bed6ff1e..8e38d6780e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -18,6 +18,7 @@
 - The irreversible, reversible and hysteresis heat sources now sum over every particle phase. `"open-circuit potential"` stays a plain string when it is left at its default, so on a composite electrode it was zipped against the two phase names and silently truncated the loop after the primary phase: a graphite electrode split into two identical halves reported 56% of the irreversible heating of the equivalent single-phase electrode. Passing the option once per phase was a workaround. ([#5770](https://github.com/pybamm-team/PyBaMM/pull/5770))
 - Discretisation no longer mutates symbols it has already cached, so it no longer depends on the order outputs are processed in. An isothermal DFN with a `pybamm.Interpolant` of `pybamm.t` as `"Ambient temperature [K]"` and concentration-dependent electrolyte parameters raised a `ShapeError`. Copying a `StateVectorDot` also keeps it a `StateVectorDot`; `create_copy()` used to return a `StateVector`, so a copied time derivative read the state instead. ([#3630](https://github.com/pybamm-team/PyBaMM/issues/3630), [#4670](https://github.com/pybamm-team/PyBaMM/issues/4670), [#5771](https://github.com/pybamm-team/PyBaMM/pull/5771))
 - The heat of mixing now sums over every particle phase, so `"heat of mixing": "true"` works on a composite electrode. It previously read the single-phase variable names (`"Negative particle radius [m]"`) and `param.n.prim` only, so any model with `"particle phases": "2"` raised `ModelError: Missing variable for submodel 'thermal': 'Negative particle radius [m]'` at build time. Single-phase results are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
+- The heat of mixing now uses the particle's effective diffusivity, the one its flux uses, instead of re-evaluating the base diffusivity. With `"diffusivity": "current sigmoid"` or `"stress-induced diffusion": "true"` the two differ, and the heat sources stopped conserving energy: a phase with lithiation and delithiation diffusivities of a third and three times the base value was off by 33% of the heat. A uniform particle profile now has zero heat of mixing, so a DFN with `"particle": "uniform profile"` and `"heat of mixing": "true"` builds again. Results with the default options are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
 - `pybamm.Function(scipy.special.erf, ...)` converts to CasADi again; CasADi 3.8 stopped dispatching scipy's `erf` ufunc onto symbolic values, so the generic `Function` path now maps it to `casadi.erf` explicitly. `pybamm.erf` was unaffected. ([#5761](https://github.com/pybamm-team/PyBaMM/pull/5761))
 - The dimensional minimum particle surface concentration is registered as `"Minimum <domain> particle surface concentration [mol.m-3]"`. An implicit string concatenation repeated its first line, so the variable was only reachable under the doubled name `"Minimum negative particle Minimum negative particle surface concentration [mol.m-3]"`, and the intended name resolved to nothing. The `Maximum` sibling was unaffected. ([#5759](https://github.com/pybamm-team/PyBaMM/pull/5759))
 - `FiniteVolume` node-to-edge shifts now use the true node spacing instead of weights that assume a uniform mesh, so the exterior edge values agree with `boundary_value`. On a non-uniform mesh they did not, which broke lithium conservation wherever a flux is imposed as a gradient divided by a surface value, as the Fickian particle does: with a concentration-dependent particle diffusivity on an `Exponential1DSubMesh`, `ORegan2022` gained 26% of its particle lithium over three 1C cycles. Uniform-mesh results are unchanged. ([#5755](https://github.com/pybamm-team/PyBaMM/pull/5755))
diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
index 0c615883b0..c74e2d1718 100644
--- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -391,6 +391,10 @@ def _phase_heat_of_mixing(self, variables, domain, phase_name, phase_param):
         referenced to, not the hysteresis-weighted OCP that sets the voltage. Energy
         only balances if both heat sources share that reference.
         """
+        phase_options = getattr(getattr(self.options, domain), phase_param.phase)
+        if phase_options["particle"] == "uniform profile":
+            # no concentration gradient, so nothing to mix
+            return pybamm.FullBroadcast(0, [f"{domain} electrode"], "current collector")
         Domain = domain.capitalize()
         F = pybamm.constants.F.value
         a = variables[
@@ -398,17 +402,25 @@ def _phase_heat_of_mixing(self, variables, domain, phase_name, phase_param):
         ]
         R = variables[f"{Domain} {phase_name}particle radius [m]"]
         N = a / (4 * np.pi * R**2)
+        # the diffusivity the particle flux actually uses, including any current
+        # sigmoid or stress-induced diffusion, so the dissipation matches the flux
         if self.x_average:
             c = variables[
                 f"X-averaged {domain} {phase_name}particle concentration [mol.m-3]"
             ]
             T = variables[f"X-averaged {domain} electrode temperature [K]"]
+            D = variables[
+                f"X-averaged {domain} {phase_name}particle effective "
+                "diffusivity [m2.s-1]"
+            ]
         else:
             c = variables[f"{Domain} {phase_name}particle concentration [mol.m-3]"]
             T = variables[f"{Domain} electrode temperature [K]"]
+            D = variables[
+                f"{Domain} {phase_name}particle effective diffusivity [m2.s-1]"
+            ]
         T_part = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle"])
         dc_dr2 = pybamm.inner(pybamm.grad(c), pybamm.grad(c))
-        D = phase_param.D(c, T_part)
         dUeq = phase_param.U(c / phase_param.c_max, T_part).diff(c)
         integrand_r = D * dc_dr2 * dUeq
         integration_variable_r = [
diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
index 5be0565c3f..36e497bf39 100644
--- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
+++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
@@ -240,13 +240,15 @@ def test_heat_of_mixing_energy_balance(self, composite):
         parameter_values = pybamm.ParameterValues("Chen2020")
         if composite:
             # Split the positive electrode, which carries most of the heat of mixing,
-            # into two identical halves, and give the second hysteresis about its OCP.
-            # The branches are tilted so that their slopes differ from the OCP's
+            # into two identical halves, and give the second hysteresis about its OCP
+            # and a direction-dependent diffusivity. The OCP branches are tilted so
+            # that their slopes differ from the OCP's
             options["particle phases"] = ("1", "2")
             options["open-circuit potential"] = (
                 "single",
                 ("single", "current sigmoid"),
             )
+            options["diffusivity"] = ("single", ("single", "current sigmoid"))
             for name in [
                 "Positive electrode OCP [V]",
                 "Positive electrode OCP entropic change [V.K-1]",
@@ -265,6 +267,9 @@ def test_heat_of_mixing_energy_balance(self, composite):
                 )
                 del parameter_values[name]
             mean_ocp = parameter_values["Secondary: Positive electrode OCP [V]"]
+            diffusivity = parameter_values[
+                "Secondary: Positive particle diffusivity [m2.s-1]"
+            ]
             parameter_values.update(
                 {
                     "Secondary: Positive electrode lithiation OCP [V]": lambda sto: (
@@ -273,6 +278,12 @@ def test_heat_of_mixing_energy_balance(self, composite):
                     "Secondary: Positive electrode delithiation OCP [V]": lambda sto: (
                         mean_ocp(sto) + 0.005 + 0.1 * sto
                     ),
+                    "Secondary: Positive particle lithiation diffusivity [m2.s-1]": (
+                        diffusivity / 2
+                    ),
+                    "Secondary: Positive particle delithiation diffusivity [m2.s-1]": (
+                        diffusivity * 2
+                    ),
                 },
                 check_already_exists=False,
             )

From d469838097e55e31f31c52836299dd7b703dadcf Mon Sep 17 00:00:00 2001
From: tleo <tleo@silanano.com>
Date: Thu, 24 Sep 2026 19:44:23 -0700
Subject: [PATCH 7/9] refactor: sum the heat of mixing without a one-phase
 branch

Apply the review suggestion: build the phase-name prefixes from
options.phases, with "" for a one-phase electrode, and sum over them, so
the one-phase case needs no branch of its own. One-phase results are
bit-identical.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
---
 .../src/pybamm/models/submodels/thermal/base_thermal.py  | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
index c74e2d1718..2cf7fd0f0c 100644
--- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -433,15 +433,12 @@ def _electrode_heat_of_mixing(self, variables, domain):
         """Sum the heat of mixing over the particle phases of one electrode."""
         electrode_param = self.param.domain_params[domain]
         phases = self.options.phases[domain]
-        if len(phases) == 1:
-            return self._phase_heat_of_mixing(
-                variables, domain, "", electrode_param.prim
-            )
+        names = [""] if len(phases) == 1 else [f"{phase} " for phase in phases]
         return sum(
             self._phase_heat_of_mixing(
-                variables, domain, f"{phase} ", electrode_param.phase_params[phase]
+                variables, domain, name, electrode_param.phase_params[phase]
             )
-            for phase in phases
+            for name, phase in zip(names, phases, strict=True)
         )
 
     def _heat_of_mixing(self, variables):

From 49e02543f6cad0417a5c5d5bebac5343390a7a5a Mon Sep 17 00:00:00 2001
From: rtimms <rob@ionworks.com>
Date: Thu, 8 Oct 2026 19:06:53 +0100
Subject: [PATCH 8/9] docs: move the heat of mixing CHANGELOG entries to
 Unreleased

Merging main put them under the released v26.9.0.0 section.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7f0c23e068..a11c7538e0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,8 @@
 
 ## Bug fixes
 
+- The heat of mixing now sums over every particle phase, so `"heat of mixing": "true"` works on a composite electrode. It previously read the single-phase variable names (`"Negative particle radius [m]"`) and `param.n.prim` only, so any model with `"particle phases": "2"` raised `ModelError: Missing variable for submodel 'thermal': 'Negative particle radius [m]'` at build time. Single-phase results are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
+- The heat of mixing now uses the particle's effective diffusivity, the one its flux uses, instead of re-evaluating the base diffusivity. With `"diffusivity": "current sigmoid"` or `"stress-induced diffusion": "true"` the two differ, and the heat sources stopped conserving energy: a phase with lithiation and delithiation diffusivities of a third and three times the base value was off by 33% of the heat. A uniform particle profile now has zero heat of mixing, so a DFN with `"particle": "uniform profile"` and `"heat of mixing": "true"` builds again. Results with the default options are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
 - A single `IDAKLUSolver` solve no longer runs inside a one-thread OpenMP parallel region. OpenMP code the model calls during the solve, such as a CasADi external function, now runs at the top level instead of nested, so it takes the top-level thread count rather than the nested-level one (for example 8, not 1, under `OMP_NUM_THREADS=8,1`). ([#5837](https://github.com/pybamm-team/PyBaMM/pull/5837))
 - `NonlinearSolver` now logs its Newton iterations when `pybamm.logger` is at `DEBUG` level; previously they were silently dropped. ([#5849](https://github.com/pybamm-team/PyBaMM/pull/5849))
 - Fixed `FiniteVolume2D` boundary values, gradients and integrals raising an `IndexError` on 2D meshes with fewer than 3 nodes in either direction, even when the extrapolation only needed the other direction. Meshes too coarse for the requested extrapolation order now raise a `DiscretisationError`. ([#5839](https://github.com/pybamm-team/PyBaMM/pull/5839))
@@ -88,8 +90,6 @@
 - `BasicDFN2D` now reports `"Total lithium [mol]"` and the `"Negative/Positive/Total solid lithium [mol]"` variables in mol for the whole cell: the integrals over `x` and `z` were in mol per metre of width, and are now scaled by the electrode width `L_y`. The volumetric interfacial current density `a * j`, previously mislabelled `"Current density [A.m-2]"`, is renamed to `"Sum of volumetric interfacial current densities [A.m-3]"` to match the full DFN. ([#5775](https://github.com/pybamm-team/PyBaMM/pull/5775))
 - The irreversible, reversible and hysteresis heat sources now sum over every particle phase. `"open-circuit potential"` stays a plain string when it is left at its default, so on a composite electrode it was zipped against the two phase names and silently truncated the loop after the primary phase: a graphite electrode split into two identical halves reported 56% of the irreversible heating of the equivalent single-phase electrode. Passing the option once per phase was a workaround. ([#5770](https://github.com/pybamm-team/PyBaMM/pull/5770))
 - Discretisation no longer mutates symbols it has already cached, so it no longer depends on the order outputs are processed in. An isothermal DFN with a `pybamm.Interpolant` of `pybamm.t` as `"Ambient temperature [K]"` and concentration-dependent electrolyte parameters raised a `ShapeError`. Copying a `StateVectorDot` also keeps it a `StateVectorDot`; `create_copy()` used to return a `StateVector`, so a copied time derivative read the state instead. ([#3630](https://github.com/pybamm-team/PyBaMM/issues/3630), [#4670](https://github.com/pybamm-team/PyBaMM/issues/4670), [#5771](https://github.com/pybamm-team/PyBaMM/pull/5771))
-- The heat of mixing now sums over every particle phase, so `"heat of mixing": "true"` works on a composite electrode. It previously read the single-phase variable names (`"Negative particle radius [m]"`) and `param.n.prim` only, so any model with `"particle phases": "2"` raised `ModelError: Missing variable for submodel 'thermal': 'Negative particle radius [m]'` at build time. Single-phase results are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
-- The heat of mixing now uses the particle's effective diffusivity, the one its flux uses, instead of re-evaluating the base diffusivity. With `"diffusivity": "current sigmoid"` or `"stress-induced diffusion": "true"` the two differ, and the heat sources stopped conserving energy: a phase with lithiation and delithiation diffusivities of a third and three times the base value was off by 33% of the heat. A uniform particle profile now has zero heat of mixing, so a DFN with `"particle": "uniform profile"` and `"heat of mixing": "true"` builds again. Results with the default options are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
 - `pybamm.Function(scipy.special.erf, ...)` converts to CasADi again; CasADi 3.8 stopped dispatching scipy's `erf` ufunc onto symbolic values, so the generic `Function` path now maps it to `casadi.erf` explicitly. `pybamm.erf` was unaffected. ([#5761](https://github.com/pybamm-team/PyBaMM/pull/5761))
 - The dimensional minimum particle surface concentration is registered as `"Minimum <domain> particle surface concentration [mol.m-3]"`. An implicit string concatenation repeated its first line, so the variable was only reachable under the doubled name `"Minimum negative particle Minimum negative particle surface concentration [mol.m-3]"`, and the intended name resolved to nothing. The `Maximum` sibling was unaffected. ([#5759](https://github.com/pybamm-team/PyBaMM/pull/5759))
 - `FiniteVolume` node-to-edge shifts now use the true node spacing instead of weights that assume a uniform mesh, so the exterior edge values agree with `boundary_value`. On a non-uniform mesh they did not, which broke lithium conservation wherever a flux is imposed as a gradient divided by a surface value, as the Fickian particle does: with a concentration-dependent particle diffusivity on an `Exponential1DSubMesh`, `ORegan2022` gained 26% of its particle lithium over three 1C cycles. Uniform-mesh results are unchanged. ([#5755](https://github.com/pybamm-team/PyBaMM/pull/5755))

From add7b2a23e62ed026cdbd2222d90c903f724a154 Mon Sep 17 00:00:00 2001
From: Robert Timms <43040151+rtimms@users.noreply.github.com>
Date: Thu, 22 Oct 2026 13:05:07 +0100
Subject: [PATCH 9/9] test: shorten heat of mixing energy-balance comment

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../test_lithium_ion/test_thermal_models.py                 | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
index 36e497bf39..0ada0967fb 100644
--- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
+++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
@@ -239,10 +239,8 @@ def test_heat_of_mixing_energy_balance(self, composite):
         }
         parameter_values = pybamm.ParameterValues("Chen2020")
         if composite:
-            # Split the positive electrode, which carries most of the heat of mixing,
-            # into two identical halves, and give the second hysteresis about its OCP
-            # and a direction-dependent diffusivity. The OCP branches are tilted so
-            # that their slopes differ from the OCP's
+            # Tilted hysteresis branches and a direction-dependent diffusivity make
+            # the second phase's equilibrium OCP and flux differ from the first's
             options["particle phases"] = ("1", "2")
             options["open-circuit potential"] = (
                 "single",