fix: sum the electrochemical heat sources over every particle phase

CHANGELOG.md

@@ -7,6 +7,7 @@
 
 ## Bug fixes
 
+- 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))
 - `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))

packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py

@@ -136,7 +136,8 @@ def _get_standard_coupled_variables(self, variables):
             phase_names = ["primary ", "secondary "]
         ocp_options = self.options.positive["open-circuit potential"]
         if isinstance(ocp_options, str):
-            ocp_options = [ocp_options]
+            # a single option applies to every phase
+            ocp_options = (ocp_options,) * len(phase_names)
 
         Q_rxn_p, Q_rev_p, Q_hys_p = 0, 0, 0
         T_p = variables["Positive electrode temperature [K]"]
@@ -171,7 +172,8 @@ def _get_standard_coupled_variables(self, variables):
             phase_names = ["primary ", "secondary "]
         ocp_options = self.options.negative["open-circuit potential"]
         if isinstance(ocp_options, str):
-            ocp_options = [ocp_options]
+            # a single option applies to every phase
+            ocp_options = (ocp_options,) * len(phase_names)
 
         if self.options.electrode_types["negative"] == "planar":
             i_n = variables["Lithium metal total interfacial current density [A.m-2]"]

packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py

@@ -176,3 +176,64 @@ def test_compare_SPMe_silicon_graphite(self):
     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):
+        """
+        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",
+        }
+        t_eval = [0, 3600]
+        t_interp = np.linspace(0, 3600)
+
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        sol = pybamm.Simulation(
+            model_class(options), parameter_values=parameter_values
+        ).solve(t_eval=t_eval, t_interp=t_interp)
+
+        parameter_values_two_phase = pybamm.ParameterValues("Chen2020")
+        for parameter in [
+            "Negative electrode OCP [V]",
+            "Negative electrode OCP entropic change [V.K-1]",
+            "Maximum concentration in negative electrode [mol.m-3]",
+            "Initial concentration in negative electrode [mol.m-3]",
+            "Negative particle radius [m]",
+            "Negative particle diffusivity [m2.s-1]",
+            "Negative electrode exchange-current density [A.m-2]",
+            "Negative electrode active material volume fraction",
+        ]:
+            value = parameter_values_two_phase[parameter]
+            if parameter.endswith("active material volume fraction"):
+                value = value / 2
+            parameter_values_two_phase.update(
+                {f"Primary: {parameter}": value, f"Secondary: {parameter}": value}
+            )
+            del parameter_values_two_phase[parameter]
+
+        sol_two_phase = pybamm.Simulation(
+            model_class({"particle phases": ("2", "1"), **options}),
+            parameter_values=parameter_values_two_phase,
+        ).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 total heating [W.m-3]",
+            "Voltage [V]",
+        ]:
+            np.testing.assert_allclose(
+                sol[variable](t_interp),
+                sol_two_phase[variable](t_interp),
+                rtol=1e-2,
+                atol=1e-8,
+            )
+
+    def test_compare_heat_sources_SPM_graphite_graphite(self):
+        self.compare_heat_sources_two_phase_graphite_graphite(pybamm.lithium_ion.SPM)
+
+    def test_compare_heat_sources_DFN_graphite_graphite(self):
+        self.compare_heat_sources_two_phase_graphite_graphite(pybamm.lithium_ion.DFN)