Repository: pybamm-team/PyBaMM

Issue:
[Bug]: SPMe with particle size distribution
### PyBaMM Version

23.5

### Python Version

3.9

### Describe the bug

Trying to use a particle size distribution with SPMe fails due to a `DomainError`. DFN works fine. Looks like a broadcasting or averaging issue.

### Steps to Reproduce

import pybamm

options = {
    "particle size": "distribution",
    "surface form": "algebraic",
}
model = pybamm.lithium_ion.SPMe(options)

param = pybamm.ParameterValues("Marquis2019")
param = pybamm.get_size_distribution_parameters(param)

sim = pybamm.Simulation(model, parameter_values=param)
sim.solve([0, 3600])
sim.plot()


### Relevant log output

```shell
File "/Users/robertwtimms/Documents/PyBaMM/pybamm/models/submodels/interface/kinetics/base_kinetics.py", line 98, in get_coupled_variables
    eta_r = delta_phi - ocp
  File "/Users/robertwtimms/Documents/PyBaMM/pybamm/expression_tree/symbol.py", line 551, in __sub__
    return pybamm.simplified_subtraction(self, other)
  File "/Users/robertwtimms/Documents/PyBaMM/pybamm/expression_tree/binary_operators.py", line 928, in simplified_subtraction
    return pybamm.simplify_if_constant(Subtraction(left, right))
  File "/Users/robertwtimms/Documents/PyBaMM/pybamm/expression_tree/binary_operators.py", line 222, in __init__
    super().__init__("-", left, right)
  File "/Users/robertwtimms/Documents/PyBaMM/pybamm/expression_tree/binary_operators.py", line 58, in __init__
    domains = self.get_children_domains([left, right])
  File "/Users/robertwtimms/Documents/PyBaMM/pybamm/expression_tree/symbol.py", line 371, in get_children_domains
    raise pybamm.DomainError(
pybamm.expression_tree.exceptions.DomainError: children must have same or empty domains, not ['negative electrode'] and ['current collector']
```


Discussion:

Patch:
diff --git a/examples/scripts/cycling_ageing.py b/examples/scripts/cycling_ageing.py
index 29e2c0d16..6ac0da8b6 100644
--- a/examples/scripts/cycling_ageing.py
+++ b/examples/scripts/cycling_ageing.py
@@ -1,7 +1,7 @@
 import pybamm as pb
 
 pb.set_logging_level("NOTICE")
-model = pb.lithium_ion.SPM(
+model = pb.lithium_ion.DFN(
     {
         "SEI": "ec reaction limited",
         "SEI film resistance": "distributed",
diff --git a/noxfile.py b/noxfile.py
index 264a22a5f..2f1231214 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -114,13 +114,13 @@ def set_dev(session):
     session.install("cmake")
     if sys.platform == "linux" or sys.platform == "darwin":
         session.run(
-            "echo",
-            "export",
-            f"LD_LIBRARY_PATH={PYBAMM_ENV['LD_LIBRARY_PATH']}",
-            ">>",
-            f"{envbindir}/activate",
-            external=True,  # silence warning about echo being an external command
-        )
+        "echo",
+        "export",
+        f"LD_LIBRARY_PATH={PYBAMM_ENV['LD_LIBRARY_PATH']}",
+        ">>",
+        f"{envbindir}/activate",
+        external=True,  # silence warning about echo being an external command
+    )
 
 
 @nox.session(name="tests")
diff --git a/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py b/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
index e90b8d964..b06d6664a 100644
--- a/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
+++ b/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py
@@ -268,8 +268,12 @@ class ElectrodeSOHSolver:
         sol_dict = {key: sol[key].data[0] for key in sol.all_models[0].variables.keys()}
 
         # Calculate theoretical energy
+        x_0 = sol_dict["x_0"]
+        y_0 = sol_dict["y_0"]
+        x_100 = sol_dict["x_100"]
+        y_100 = sol_dict["y_100"]
         energy = pybamm.lithium_ion.electrode_soh.theoretical_energy_integral(
-            self.parameter_values, sol_dict
+            self.parameter_values, x_100, x_0, y_100, y_0
         )
         sol_dict.update({"Maximum theoretical energy [W.h]": energy})
         return sol_dict
@@ -594,7 +598,7 @@ def get_min_max_stoichiometries(
     return esoh_solver.get_min_max_stoichiometries()
 
 
-def theoretical_energy_integral(parameter_values, inputs, points=100):
+def theoretical_energy_integral(parameter_values, n_i, n_f, p_i, p_f, points=100):
     """
     Calculate maximum energy possible from a cell given OCV, initial soc, and final soc
     given voltage limits, open-circuit potentials, etc defined by parameter_values
@@ -614,25 +618,20 @@ def theoretical_energy_integral(parameter_values, inputs, points=100):
     E
         The total energy of the cell in Wh
     """
-    x_0 = inputs["x_0"]
-    y_0 = inputs["y_0"]
-    x_100 = inputs["x_100"]
-    y_100 = inputs["y_100"]
-    Q_p = inputs["Q_p"]
-    x_vals = np.linspace(x_100, x_0, num=points)
-    y_vals = np.linspace(y_100, y_0, num=points)
+    n_vals = np.linspace(n_i, n_f, num=points)
+    p_vals = np.linspace(p_i, p_f, num=points)
     # Calculate OCV at each stoichiometry
     param = pybamm.LithiumIonParameters()
     T = param.T_amb(0)
-    Vs = np.empty(x_vals.shape)
-    for i in range(x_vals.size):
+    Vs = np.empty(n_vals.shape)
+    for i in range(n_vals.size):
         Vs[i] = (
-            parameter_values.evaluate(param.p.prim.U(y_vals[i], T)).item()
-            - parameter_values.evaluate(param.n.prim.U(x_vals[i], T)).item()
+            parameter_values.evaluate(param.p.prim.U(p_vals[i], T)).item()
+            - parameter_values.evaluate(param.n.prim.U(n_vals[i], T)).item()
         )
     # Calculate dQ
-    Q = Q_p * (y_0 - y_100)
-    dQ = Q / (points - 1)
+    Q_p = parameter_values.evaluate(param.p.prim.Q_init) * (p_f - p_i)
+    dQ = Q_p / (points - 1)
     # Integrate and convert to W-h
     E = np.trapz(Vs, dx=dQ)
     return E
@@ -664,10 +663,7 @@ def calculate_theoretical_energy(
     # Get initial and final stoichiometric values.
     x_100, y_100 = get_initial_stoichiometries(initial_soc, parameter_values)
     x_0, y_0 = get_initial_stoichiometries(final_soc, parameter_values)
-    Q_p = parameter_values.evaluate(pybamm.LithiumIonParameters().p.prim.Q_init)
     E = theoretical_energy_integral(
-        parameter_values,
-        {"x_100": x_100, "x_0": x_0, "y_100": y_100, "y_0": y_0, "Q_p": Q_p},
-        points=points,
+        parameter_values, x_100, x_0, y_100, y_0, points=points
     )
     return E
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index 6705998d4..9b3746789 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -94,12 +94,14 @@ class BaseKinetics(BaseInterface):
             ocp = variables[
                 f"{Domain} electrode {reaction_name}open-circuit potential [V]"
             ]
-        # If ocp was broadcast, and delta_phi's secondary domain is "current collector",
-        # then take only the orphan.
-        if isinstance(ocp, pybamm.Broadcast) and delta_phi.domains["secondary"] == [
-            "current collector"
-        ]:
-            ocp = ocp.orphans[0]
+        # If ocp was broadcast, and the reaction is lithium metal plating OR
+        # delta_phi's secondary domain is "current collector", then take only the
+        # orphan.
+        if isinstance(ocp, pybamm.Broadcast):
+            if self.reaction == "lithium metal plating":
+                ocp = ocp.orphans[0]
+            elif delta_phi.domains["secondary"] == ["current collector"]:
+                ocp = ocp.orphans[0]
 
         # Get reaction overpotential
         eta_r = delta_phi - ocp


Test patch:
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py
index b1c320909..f0db6a3a6 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_electrode_soh.py
@@ -36,12 +36,12 @@ class TestElectrodeSOH(TestCase):
                 self.assertAlmostEqual(sol[key], sol_split[key].data[0], places=5)
             else:
                 # theoretical_energy is not present in sol_split
-                inputs = {
-                    k: sol_split[k].data[0]
-                    for k in ["x_0", "y_0", "x_100", "y_100", "Q_p"]
-                }
+                x_0 = sol_split["x_0"].data[0]
+                y_0 = sol_split["y_0"].data[0]
+                x_100 = sol_split["x_100"].data[0]
+                y_100 = sol_split["y_100"].data[0]
                 energy = pybamm.lithium_ion.electrode_soh.theoretical_energy_integral(
-                    parameter_values, inputs
+                    parameter_values, x_100, x_0, y_100, y_0
                 )
                 self.assertAlmostEqual(sol[key], energy, places=5)