Fix DomainError for hysteresis/MSMR OCP with particle size distribution and thermal

Combining hysteresis / MSMR / single open-circuit-potential submodels with
particle-size distribution and a non-isothermal thermal submodel raised a
DomainError when the thermal submodel built the hysteresis-heating term.

The root cause was that the "equilibrium open-circuit potential [V]" and
"OCP hysteresis [V]" variables were published on the particle-size domain
under PSD, whereas the parallel "open-circuit potential [V]" was correctly
size-averaged onto the electrode domain. base_thermal then subtracted them
and concatenated across electrodes, hitting a domain mismatch.

- base_hysteresis_ocp: size-average U_eq and H under PSD, publish
  distribution-named variables parallel to the existing "open-circuit
  potential distribution [V]", and properly x-average U_eq_x_av / H_x_av
  (previously they were the same expression as U_eq / H, never averaged).
- single_ocp and msmr_ocp: after _get_standard_ocp_variables has shaped
  the main OCP variables, reuse them for the equilibrium variables so
  both sides line up regardless of PSD.

Adds unit and integration regression tests for SPM, SPMe, DFN, MPM, and
Newman-Tobias covering:
- hysteresis OCP + PSD + lumped thermal
- single OCP + PSD + lumped thermal (pins equilibrium OCP domain)
- MSMR + PSD + lumped thermal

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: rtimms

src/pybamm/models/submodels/interface/open_circuit_potential/base_hysteresis_ocp.py

@@ -82,21 +82,44 @@ def _get_coupled_variables(self, variables):
                 variables
             )
             U_eq = self.phase_param.U(sto_surf, T)
-            U_eq_x_av = self.phase_param.U(sto_surf, T)
             U_lith = self.phase_param.U(sto_surf, T, "lithiation")
             U_lith_bulk = self.phase_param.U(sto_bulk, T_bulk, "lithiation")
             U_delith = self.phase_param.U(sto_surf, T, "delithiation")
             U_delith_bulk = self.phase_param.U(sto_bulk, T_bulk, "delithiation")
 
             H = U_lith - U_delith
-            H_x_av = pybamm.x_average(H)
+
+            # Under particle-size distribution the OCPs above live on the
+            # "<domain> particle size" domain. Publish the distribution-named
+            # variables at that shape, and size-average before publishing under
+            # the electrode-level names so they line up with other
+            # electrode-domain variables (e.g. in the thermal submodel).
+            if domain_options["particle size"] == "distribution":
+                variables.update(
+                    {
+                        f"{Domain} electrode {phase_name}OCP hysteresis distribution [V]": H,
+                        f"X-averaged {domain} electrode {phase_name}OCP hysteresis distribution [V]": pybamm.x_average(
+                            H
+                        ),
+                        f"{Domain} electrode {phase_name}equilibrium open-circuit potential distribution [V]": U_eq,
+                        f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential distribution [V]": pybamm.x_average(
+                            U_eq
+                        ),
+                    }
+                )
+                U_eq = pybamm.size_average(U_eq)
+                H = pybamm.size_average(H)
 
             variables.update(
                 {
                     f"{Domain} electrode {phase_name}OCP hysteresis [V]": H,
-                    f"X-averaged {domain} electrode {phase_name}OCP hysteresis [V]": H_x_av,
+                    f"X-averaged {domain} electrode {phase_name}OCP hysteresis [V]": pybamm.x_average(
+                        H
+                    ),
                     f"{Domain} electrode {phase_name}equilibrium open-circuit potential [V]": U_eq,
-                    f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential [V]": U_eq_x_av,
+                    f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential [V]": pybamm.x_average(
+                        U_eq
+                    ),
                 }
             )
 

src/pybamm/models/submodels/interface/open_circuit_potential/msmr_ocp.py

@@ -69,12 +69,31 @@ def get_coupled_variables(self, variables):
             dUdT = self.phase_param.dUdT(sto_surf)
 
         variables.update(self._get_standard_ocp_variables(ocp_surf, ocp_bulk, dUdT))
+
+        # MSMR has no hysteresis, so the equilibrium OCP equals the OCP.
+        # Reuse the already-shaped variables published by
+        # _get_standard_ocp_variables so the equilibrium variable has the
+        # same domain regardless of particle-size distribution (it is
+        # consumed by the thermal hysteresis-heating term).
         variables.update(
             {
-                f"{Domain} electrode {phase_name}equilibrium open-circuit potential [V]": ocp_surf,
-                f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential [V]": pybamm.x_average(
-                    ocp_surf
-                ),
+                f"{Domain} electrode {phase_name}equilibrium open-circuit potential [V]": variables[
+                    f"{Domain} electrode {phase_name}open-circuit potential [V]"
+                ],
+                f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential [V]": variables[
+                    f"X-averaged {domain} electrode {phase_name}open-circuit potential [V]"
+                ],
             }
         )
+        if domain_options["particle size"] == "distribution":
+            variables.update(
+                {
+                    f"{Domain} electrode {phase_name}equilibrium open-circuit potential distribution [V]": variables[
+                        f"{Domain} electrode {phase_name}open-circuit potential distribution [V]"
+                    ],
+                    f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential distribution [V]": variables[
+                        f"X-averaged {domain} electrode {phase_name}open-circuit potential distribution [V]"
+                    ],
+                }
+            )
         return variables

src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py

@@ -9,6 +9,7 @@
 class SingleOpenCircuitPotential(BaseOpenCircuitPotential):
     def get_coupled_variables(self, variables):
         domain, Domain = self.domain_Domain
+        domain_options = getattr(self.options, domain)
         phase_name = self.phase_name
 
         if self.reaction == "lithium-ion main":
@@ -21,15 +22,6 @@ def get_coupled_variables(self, variables):
 
             dUdT = self.phase_param.dUdT(sto_surf)
 
-            variables.update(
-                {
-                    f"{Domain} electrode {phase_name}equilibrium open-circuit potential [V]": ocp_surf,
-                    f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential [V]": pybamm.x_average(
-                        ocp_surf
-                    ),
-                }
-            )
-
         elif self.reaction == "lithium metal plating":
             T = variables[f"{Domain} electrode temperature [K]"]
             ocp_surf = 0 * T
@@ -52,4 +44,32 @@ def get_coupled_variables(self, variables):
             dUdT = pybamm.Scalar(0)
 
         variables.update(self._get_standard_ocp_variables(ocp_surf, ocp_bulk, dUdT))
+
+        if self.reaction == "lithium-ion main":
+            # For "single" OCP the equilibrium OCP equals the OCP. Reuse the
+            # already-shaped variables published by _get_standard_ocp_variables
+            # so the equilibrium variable has the same domain regardless of
+            # particle-size distribution (it is consumed by the thermal
+            # hysteresis-heating term).
+            variables.update(
+                {
+                    f"{Domain} electrode {phase_name}equilibrium open-circuit potential [V]": variables[
+                        f"{Domain} electrode {phase_name}open-circuit potential [V]"
+                    ],
+                    f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential [V]": variables[
+                        f"X-averaged {domain} electrode {phase_name}open-circuit potential [V]"
+                    ],
+                }
+            )
+            if domain_options["particle size"] == "distribution":
+                variables.update(
+                    {
+                        f"{Domain} electrode {phase_name}equilibrium open-circuit potential distribution [V]": variables[
+                            f"{Domain} electrode {phase_name}open-circuit potential distribution [V]"
+                        ],
+                        f"X-averaged {domain} electrode {phase_name}equilibrium open-circuit potential distribution [V]": variables[
+                            f"X-averaged {domain} electrode {phase_name}open-circuit potential distribution [V]"
+                        ],
+                    }
+                )
         return variables

tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py

@@ -316,6 +316,35 @@ def test_both_swelling_only(self):
         parameter_values = pybamm.ParameterValues("Ai2020")
         self.run_basic_processing_test(options, parameter_values=parameter_values)
 
+    def test_psd_hysteresis_thermal(self):
+        # Regression test: hysteresis OCP + particle size distribution + a
+        # non-isothermal thermal submodel previously raised a DomainError when
+        # the thermal submodel built the hysteresis heating term.
+        options = {
+            "open-circuit potential": ("one-state hysteresis", "single"),
+            "particle size": "distribution",
+            "surface form": "algebraic",
+            "thermal": "lumped",
+        }
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values = pybamm.get_size_distribution_parameters(parameter_values)
+        parameter_values.update(
+            {
+                "Negative electrode lithiation OCP [V]": lambda sto: (
+                    parameter_values["Negative electrode OCP [V]"](sto) - 0.1
+                ),
+                "Negative electrode delithiation OCP [V]": lambda sto: (
+                    parameter_values["Negative electrode OCP [V]"](sto) + 0.1
+                ),
+                "Negative particle lithiation hysteresis decay rate": 10,
+                "Negative particle delithiation hysteresis decay rate": 10,
+                "Initial hysteresis state in negative electrode": 0.0,
+            }
+        )
+        model = self.model(options)
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all(skip_output_tests=True)
+
     def test_composite_graphite_silicon(self):
         options = {
             "particle phases": ("2", "1"),
@@ -356,6 +385,7 @@ def test_basic_processing_msmr(self):
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all(skip_output_tests=True)
 
+
     def test_basic_processing_temperature_interpolant(self):
         times = np.arange(0, 4000, 10)
         tmax = max(times)

tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py

@@ -589,6 +589,54 @@ def test_well_posed_psd_swelling_only(self):
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_psd_hysteresis_thermal(self):
+        # Regression test: hysteresis OCP + particle size distribution + a
+        # non-isothermal thermal submodel previously raised a DomainError in
+        # the thermal hysteresis-heating term because the "equilibrium
+        # open-circuit potential [V]" variable was left on the particle-size
+        # domain while "open-circuit potential [V]" was size-averaged.
+        options = {
+            "open-circuit potential": "one-state hysteresis",
+            "particle size": "distribution",
+            "surface form": "algebraic",
+            "thermal": "lumped",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_psd_single_ocp_thermal(self):
+        # Regression test: for "single" OCP + particle size distribution, the
+        # "equilibrium open-circuit potential [V]" variable used to be
+        # published on the particle-size domain. The thermal submodel skips
+        # the hysteresis branch for "single" so this did not raise today, but
+        # the variable itself was on the wrong domain; this test pins the
+        # corrected electrode-domain shape.
+        options = {
+            "particle size": "distribution",
+            "surface form": "algebraic",
+            "thermal": "lumped",
+        }
+        model = self.model(options)
+        model.check_well_posedness()
+        v = model.variables[
+            "Positive electrode equilibrium open-circuit potential [V]"
+        ]
+        assert v.domains["primary"] == ["positive electrode"]
+
+    def test_well_posed_psd_msmr_thermal(self):
+        # Regression test: MSMR + particle size distribution + thermal
+        # previously failed building Q_hys because "equilibrium open-circuit
+        # potential [V]" was on the particle-size domain.
+        options = {
+            "open-circuit potential": "MSMR",
+            "particle": "MSMR",
+            "intercalation kinetics": "MSMR",
+            "number of MSMR reactions": ("6", "4"),
+            "particle size": "distribution",
+            "surface form": "differential",
+            "thermal": "lumped",
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_transport_efficiency_Bruggeman(self):
         options = {"transport efficiency": "Bruggeman"}
         self.check_well_posedness(options)

tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py

@@ -111,6 +111,17 @@ def test_one_state_hysteresis_ocp(self):
         model = pybamm.lithium_ion.MPM(options)
         model.check_well_posedness()
 
+    def test_one_state_hysteresis_ocp_lumped_thermal(self):
+        # Regression test: MPM always uses particle size distribution, and
+        # combining hysteresis OCP with a non-isothermal thermal submodel
+        # previously raised a DomainError in the hysteresis heating term.
+        options = {
+            "open-circuit potential": "one-state hysteresis",
+            "thermal": "lumped",
+        }
+        model = pybamm.lithium_ion.MPM(options)
+        model.check_well_posedness()
+
     def test_mpm_with_lithium_plating(self):
         options = {
             "lithium plating": "irreversible",

tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py

@@ -59,3 +59,10 @@ def test_well_posed_composite_differential_surface_form(self):
     @pytest.mark.skip(reason="Test currently not implemented")
     def test_well_posed_composite_algebraic_surface_form(self):
         pass  # skip this test
+
+    @pytest.mark.skip(
+        reason="Newman-Tobias + MSMR + particle size distribution is not "
+        "supported (pre-existing DomainError in MSMR Butler-Volmer kinetics)"
+    )
+    def test_well_posed_psd_msmr_thermal(self):
+        pass