Throw an error when concatenations have different number of children (#…

…4562)

* add error and test

* disable thermal half cells for now

* easier way to fix it

* clean up

* move half cell code to function

* style: pre-commit fixes

* changelog

* Update CHANGELOG.md

Co-authored-by: Eric G. Kratz <[email protected]>

---------

Co-authored-by: Eric G. Kratz <[email protected]>
Co-authored-by: Alexander Bills <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

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

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
     ):

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -989,7 +989,6 @@ def options(self, extra_options):
             raise pybamm.OptionError(
                 f"must use surface formulation to solve {self!s} with hydrolysis"
             )
-
         self._options = options
 
     def set_standard_output_variables(self):

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

@@ -88,6 +88,12 @@ def _get_standard_coupled_variables(self, variables):
         # 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 @@ def _yz_average(self, var):
             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

tests/unit/test_expression_tree/test_concatenations.py

@@ -445,3 +445,18 @@ def test_to_from_json(self, mocker):
 
         # 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