Repository: pybamm-team/PyBaMM

Issue:
[Bug]: Binary operators on concatenations assume that each concatenation has the same number of orphans
### PyBaMM Version

develop

### Python Version

3.12

### Describe the bug

If a binary operator, such as `+`, is invoked with 2 concatenations (i.e. `isinstance(left, pybamm.Concatenation) and isinstance(right, pybamm.Concatenation) == True`), and the concatenations have a different number of orphans, the "extra orphans" are dropped. This can be easily seen [here](https://github.com/aabills/PyBaMM/blob/8ce4a20090fee6bae2251d54fd28b4b0d49abdfb/src/pybamm/expression_tree/binary_operators.py#L859).

This is causing reaction heating to not be present in DFN half-cell models (see comments in #4557)






### Steps to Reproduce

no MWE yet, but should be self evident in the code.

### Relevant log output

_No response_

Discussion:
aabills: Tests are going to fail on this due to an unrelated issue that I am working on concurrently, so opening this for now and once we merge the other one we can merge this one.
aabills: Hi all, this PR now disables thermal half cell models. Previously, thermal half cell models were giving incorrect results. I am working on a long-term fix for this problem, but for now, I think it's better to throw an error than give incorrect results. Specifically, reaction heating from both electrodes has been ignored for at least the last few PyBaMM releases (and likely all time). 
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4562?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.26%. Comparing base [(`f05cae2`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/f05cae2ecda531cde049ccfc69d470cdb845bf98?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`30a9547`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/30a9547ec5932ab650265831dede6440a673b4fa?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #4562      +/-   ##
===========================================
- Coverage    99.27%   99.26%   -0.02%     
===========================================
  Files          300      300              
  Lines        22795    22801       +6     
===========================================
+ Hits         22630    22633       +3     
- Misses         165      168       +3     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4562?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

kratman: @aabills Can you update the bug fix section of the change log
aabills: > @aabills Can you update the bug fix section of the change log

Done, I realized I haven't done that for the other PR's either so I did all 3.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index bf376c1a1..eee7de05d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+- Added `CoupledVariable` which provides a placeholder variable whose equation can be elsewhere in the model. ([#4556](https://github.com/pybamm-team/PyBaMM/pull/4556))
 - Adds support to `pybamm.Experiment` for the `output_variables` option in the `IDAKLUSolver`. ([#4534](https://github.com/pybamm-team/PyBaMM/pull/4534))
 - Adds an option "voltage as a state" that can be "false" (default) or "true". If "true" adds an explicit algebraic equation for the voltage. ([#4507](https://github.com/pybamm-team/PyBaMM/pull/4507))
 - Improved `QuickPlot` accuracy for simulations with Hermite interpolation. ([#4483](https://github.com/pybamm-team/PyBaMM/pull/4483))
@@ -21,6 +22,8 @@
 - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416))
 
 ## Bug Fixes
+- Added error for binary operators on two concatenations with different numbers of children. Previously, the extra children were dropped. Also fixed bug where Q_rxn was dropped from the total heating term in half-cell models. ([#4562](https://github.com/pybamm-team/PyBaMM/pull/4562))
+- Fixed bug where Q_rxn was set to 0 for the negative electrode in half-cell models. ([#4557](https://github.com/pybamm-team/PyBaMM/pull/4557))
 - Fixed bug in post-processing solutions with infeasible experiments using the (`IDAKLUSolver`). ([#4541](https://github.com/pybamm-team/PyBaMM/pull/4541))
 - Disabled IREE on MacOS due to compatibility issues and added the CasADI
   path to the environment to resolve issues on MacOS and Linux. Windows
diff --git a/src/pybamm/expression_tree/binary_operators.py b/src/pybamm/expression_tree/binary_operators.py
index 42f2f74e2..be3df653a 100644
--- a/src/pybamm/expression_tree/binary_operators.py
+++ b/src/pybamm/expression_tree/binary_operators.py
@@ -856,12 +856,17 @@ def _simplified_binary_broadcast_concatenation(
         elif isinstance(right, pybamm.Concatenation) and not isinstance(
             right, pybamm.ConcatenationVariable
         ):
-            return left.create_copy(
-                [
-                    operator(left_child, right_child)
-                    for left_child, right_child in zip(left.orphans, right.orphans)
-                ]
-            )
+            if len(left.orphans) == len(right.orphans):
+                return left.create_copy(
+                    [
+                        operator(left_child, right_child)
+                        for left_child, right_child in zip(left.orphans, right.orphans)
+                    ]
+                )
+            else:
+                raise AssertionError(
+                    "Concatenations must have the same number of children"
+                )
     if isinstance(right, pybamm.Concatenation) and not isinstance(
         right, pybamm.ConcatenationVariable
     ):
diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py
index a445f47d1..153c447ed 100644
--- a/src/pybamm/models/full_battery_models/base_battery_model.py
+++ b/src/pybamm/models/full_battery_models/base_battery_model.py
@@ -989,7 +989,6 @@ class BaseBatteryModel(pybamm.BaseModel):
             raise pybamm.OptionError(
                 f"must use surface formulation to solve {self!s} with hydrolysis"
             )
-
         self._options = options
 
     def set_standard_output_variables(self):
diff --git a/src/pybamm/models/submodels/thermal/base_thermal.py b/src/pybamm/models/submodels/thermal/base_thermal.py
index 384cb1723..d5a39435f 100644
--- a/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -88,6 +88,12 @@ class BaseThermal(pybamm.BaseSubModel):
         # TODO: change full stefan-maxwell conductivity so that i_e is always
         # a Concatenation
         i_e = variables["Electrolyte current density [A.m-2]"]
+        # Special case for half cell -- i_e has to be a concatenation for this to work due to a mismatch with Q_ohm, so we make a new i_e which is a concatenation.
+        if (not isinstance(i_e, pybamm.Concatenation)) and (
+            self.options.electrode_types["negative"] == "planar"
+        ):
+            i_e = self._get_i_e_for_half_cell_thermal(variables)
+
         phi_e = variables["Electrolyte potential [V]"]
         if isinstance(i_e, pybamm.Concatenation):
             # compute by domain if possible
@@ -411,3 +417,23 @@ class BaseThermal(pybamm.BaseSubModel):
             return pybamm.z_average(var)
         elif self.options["dimensionality"] == 2:
             return pybamm.yz_average(var)
+
+    def _get_i_e_for_half_cell_thermal(self, variables):
+        c_e_s = variables["Separator electrolyte concentration [mol.m-3]"]
+        c_e_p = variables["Positive electrolyte concentration [mol.m-3]"]
+        grad_phi_e_s = variables["Gradient of separator electrolyte potential [V.m-1]"]
+        grad_phi_e_p = variables["Gradient of positive electrolyte potential [V.m-1]"]
+        grad_c_e_s = pybamm.grad(c_e_s)
+        grad_c_e_p = pybamm.grad(c_e_p)
+        T_s = variables["Separator temperature [K]"]
+        T_p = variables["Positive electrode temperature [K]"]
+        tor_s = variables["Separator electrolyte transport efficiency"]
+        tor_p = variables["Positive electrolyte transport efficiency"]
+        i_e_s = (self.param.kappa_e(c_e_s, T_s) * tor_s) * (
+            self.param.chiRT_over_Fc(c_e_s, T_s) * grad_c_e_s - grad_phi_e_s
+        )
+        i_e_p = (self.param.kappa_e(c_e_p, T_p) * tor_p) * (
+            self.param.chiRT_over_Fc(c_e_p, T_p) * grad_c_e_p - grad_phi_e_p
+        )
+        i_e = pybamm.concatenation(i_e_s, i_e_p)
+        return i_e


Test patch:
diff --git a/tests/unit/test_expression_tree/test_concatenations.py b/tests/unit/test_expression_tree/test_concatenations.py
index 9d653a4d5..e8fad71ce 100644
--- a/tests/unit/test_expression_tree/test_concatenations.py
+++ b/tests/unit/test_expression_tree/test_concatenations.py
@@ -445,3 +445,18 @@ class TestConcatenations:
 
         # test _from_json
         assert pybamm.NumpyConcatenation._from_json(np_json) == conc_np
+
+    def test_same_number_of_children(self):
+        a = pybamm.Variable("y", domain=["1", "2"])
+        b = pybamm.Variable("z", domain=["3"])
+
+        d1 = pybamm.Variable("d1", domain=["1"])
+        d2 = pybamm.Variable("d2", domain=["2"])
+        d3 = pybamm.Variable("d3", domain=["3"])
+
+        d_concat = pybamm.concatenation(pybamm.sin(d1), pybamm.sin(d2), pybamm.sin(d3))
+        a_concat = pybamm.concatenation(pybamm.sin(a), pybamm.sin(b))
+        with pytest.raises(
+            AssertionError, match="Concatenations must have the same number of children"
+        ):
+            a_concat + d_concat