Repository: pybamm-team/PyBaMM (python)

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_

Pull request description:
Throw an error when concatenations have different number of children
# Description

Please include a summary of the change and which issue is fixed. Please also include relevant motivation and context. List any dependencies that are required for this change.

Fixes #4561 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [ ] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [x] Bug fix (non-breaking change which fixes an issue)

# Key checklist:

- [ ] No style issues: `$ pre-commit run` (or `$ nox -s pre-commit`) (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) for how to set this up to run automatically when committing locally, in just two lines of code)
- [ ] All tests pass: `$ python run-tests.py --all` (or `$ nox -s tests`)
- [ ] The documentation builds: `$ python run-tests.py --doctest` (or `$ nox -s doctests`)

You can run integration tests, unit tests, and doctests together at once, using `$ python run-tests.py --quick` (or `$ nox -s quick`).

## Further checks:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


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