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/CHANGELOG.md b/CHANGELOG.md
index 17170ab05..678ada791 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,9 +9,11 @@
 
 ## Bug fixes
 
+- Fixed a bug where the "basic" lithium-ion models gave incorrect results when using nonlinear particle diffusivity ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
+- Particle size distributions now work with SPMe and NewmanTobias models ([#3207](https://github.com/pybamm-team/PyBaMM/pull/3207))
 - Fix to simulate c_rate steps with drive cycles ([#3186](https://github.com/pybamm-team/PyBaMM/pull/3186))
 - Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019.
-- Error generated when invalid parameter values are passed.
+- Error generated when invalid parameter values are passed. ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))
 - Thevenin() model is now constructed with standard variables: `Time [s], Time [min], Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
 - Fix SEI Example Notebook ([#3166](https://github.com/pybamm-team/PyBaMM/pull/3166))
 
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
index 9275436ef..5b3892669 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
@@ -170,14 +170,14 @@ class BasicDFN(BaseModel):
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_n / (param.F * param.n.prim.D(c_s_surf_n, T)),
+                -j_n / (param.F * pybamm.surf(param.n.prim.D(c_s_n, T))),
                 "Neumann",
             ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_p / (param.F * param.p.prim.D(c_s_surf_p, T)),
+                -j_p / (param.F * pybamm.surf(param.p.prim.D(c_s_p, T))),
                 "Neumann",
             ),
         }
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_composite.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn_composite.py
index a7aaf27d5..76889d28e 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_composite.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn_composite.py
@@ -219,15 +219,24 @@ class BasicDFNComposite(BaseModel):
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n_p1] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_n_p1 / param.F / param.n.prim.D(c_s_surf_n_p1, T), "Neumann"),
+            "right": (
+                -j_n_p1 / param.F / pybamm.surf(param.n.prim.D(c_s_n_p1, T)),
+                "Neumann",
+            ),
         }
         self.boundary_conditions[c_s_n_p2] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_n_p2 / param.F / param.n.sec.D(c_s_surf_n_p2, T), "Neumann"),
+            "right": (
+                -j_n_p2 / param.F / pybamm.surf(param.n.sec.D(c_s_n_p2, T)),
+                "Neumann",
+            ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_p / param.F / param.p.prim.D(c_s_surf_p, T), "Neumann"),
+            "right": (
+                -j_p / param.F / pybamm.surf(param.p.prim.D(c_s_p, T)),
+                "Neumann",
+            ),
         }
         self.initial_conditions[c_s_n_p1] = param.n.prim.c_init
         self.initial_conditions[c_s_n_p2] = param.n.sec.c_init
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py
index 786ae559e..6586a4518 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py
@@ -160,7 +160,7 @@ class BasicDFNHalfCell(BaseModel):
         # derivatives
         self.boundary_conditions[c_s_w] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-j_w / D_w(c_s_surf_w, T) / param.F, "Neumann"),
+            "right": (-j_w / pybamm.surf(D_w(c_s_w, T)) / param.F, "Neumann"),
         }
         self.initial_conditions[c_s_w] = c_w_init
 
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
index aa10cc61a..802f8037e 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
@@ -86,14 +86,14 @@ class BasicSPM(BaseModel):
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_n / param.F / param.n.prim.D(c_s_surf_n, T),
+                -j_n / (param.F * pybamm.surf(param.n.prim.D(c_s_n, T))),
                 "Neumann",
             ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -j_p / param.F / param.p.prim.D(c_s_surf_p, T),
+                -j_p / (param.F * pybamm.surf(param.p.prim.D(c_s_p, T))),
                 "Neumann",
             ),
         }
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index ccfefe6bf..6705998d4 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -57,6 +57,7 @@ class BaseKinetics(BaseInterface):
         reaction_name = self.reaction_name
         phase_name = self.phase_name
 
+        # Get surface potential difference
         if self.reaction == "lithium metal plating":  # li metal electrode (half-cell)
             delta_phi = variables[
                 "Lithium metal interface surface potential difference [V]"
@@ -79,7 +80,8 @@ class BaseKinetics(BaseInterface):
 
         # Get exchange-current density
         j0 = self._get_exchange_current_density(variables)
-        # Get open-circuit potential variables and reaction overpotential
+
+        # Get open-circuit potential
         if (
             domain_options["particle size"] == "distribution"
             and self.options.electrode_types[domain] == "porous"
@@ -92,9 +94,14 @@ class BaseKinetics(BaseInterface):
             ocp = variables[
                 f"{Domain} electrode {reaction_name}open-circuit potential [V]"
             ]
-        # If ocp was broadcast, take only the orphan.
-        if isinstance(ocp, pybamm.Broadcast):
+        # 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]
+
+        # Get reaction overpotential
         eta_r = delta_phi - ocp
 
         # Get average interfacial current density
@@ -129,8 +136,15 @@ class BaseKinetics(BaseInterface):
                 eta_sei = pybamm.Scalar(0)
             eta_r += eta_sei
 
+        # Broadcast j0 to match eta_r's domain, if necessary
+        if j0.secondary_domain == ["current collector"] and eta_r.secondary_domain == [
+            f"{domain} electrode"
+        ]:
+            j0 = pybamm.SecondaryBroadcast(j0, [f"{domain} electrode"])
+
         # Get number of electrons in reaction
         ne = self._get_number_of_electrons_in_reaction()
+
         # Get kinetics. Note: T and u must have the same domain as j0 and eta_r
         if self.options.electrode_types[domain] == "planar":
             T = variables["X-averaged cell temperature [K]"]


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py
index 51e96342e..58d6f1c20 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py
@@ -10,17 +10,18 @@ import unittest
 
 class TestCompareBasicModels(TestCase):
     def test_compare_dfns(self):
+        parameter_values = pybamm.ParameterValues("Ecker2015")
         basic_dfn = pybamm.lithium_ion.BasicDFN()
         dfn = pybamm.lithium_ion.DFN()
 
         # Solve basic DFN
-        basic_sim = pybamm.Simulation(basic_dfn)
+        basic_sim = pybamm.Simulation(basic_dfn, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         basic_sim.solve(t_eval)
         basic_sol = basic_sim.solution
 
         # Solve main DFN
-        sim = pybamm.Simulation(dfn)
+        sim = pybamm.Simulation(dfn, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         sim.solve(t_eval)
         sol = sim.solution
@@ -64,17 +65,18 @@ class TestCompareBasicModels(TestCase):
             )
 
     def test_compare_spms(self):
+        parameter_values = pybamm.ParameterValues("Ecker2015")
         basic_spm = pybamm.lithium_ion.BasicSPM()
         spm = pybamm.lithium_ion.SPM()
 
         # Solve basic SPM
-        basic_sim = pybamm.Simulation(basic_spm)
+        basic_sim = pybamm.Simulation(basic_spm, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         basic_sim.solve(t_eval)
         basic_sol = basic_sim.solution
 
         # Solve main SPM
-        sim = pybamm.Simulation(spm)
+        sim = pybamm.Simulation(spm, parameter_values=parameter_values)
         t_eval = np.linspace(0, 3600)
         sim.solve(t_eval)
         sol = sim.solution
@@ -84,7 +86,7 @@ class TestCompareBasicModels(TestCase):
         # Compare variables
         for name in basic_spm.variables:
             np.testing.assert_allclose(
-                basic_sol[name].entries, sol[name].entries, rtol=1e-5
+                basic_sol[name].entries, sol[name].entries, rtol=1e-4
             )
 
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index 1911f2989..d60dcda8b 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -367,3 +367,7 @@ class BaseUnitTestLithiumIon:
     def test_well_posed_current_sigmoid_diffusivity(self):
         options = {"diffusivity": "current sigmoid"}
         self.check_well_posedness(options)
+
+    def test_well_posed_psd(self):
+        options = {"particle size": "distribution", "surface form": "algebraic"}
+        self.check_well_posedness(options)