Fix SPMe Surface Form Discrepancy

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -435,16 +435,15 @@ def __init__(self, extra_options):
         # The "stress-induced diffusion" option will still be overridden by
         # extra_options if provided
 
-        # Change the default for surface form based on which particle
-        # phases option is provided.
-        # return "1" if option not given
         phases_option = extra_options.get("particle phases", "1")
-        if phases_option == "1":
+        if phases_option == "1" or phases_option == ("1", "1"):
             default_options["surface form"] = "false"
         else:
             default_options["surface form"] = "algebraic"
-        # The "surface form" option will still be overridden by
-        # extra_options if provided
+            pybamm.logger.warning(
+                "Multi-phase configuration detected. Setting 'surface form' to "
+                "'algebraic' by default. Specify 'surface form' explicitly to override."
+            )
 
         # Change default SEI model based on which lithium plating option is provided
         # return "none" if option not given
@@ -638,10 +637,16 @@ def __init__(self, extra_options):
                 and options["particle"] == "Fickian diffusion"
             ):
                 raise pybamm.OptionError(
-                    "If there are multiple particle phases: 'surface form' cannot be "
-                    "'false', 'particle size' must be 'single', 'particle' must be "
-                    "'Fickian diffusion'."
+                    "If there are multiple particle phases: 'surface form' must be "
+                    "'differential' or 'algebraic', 'particle size' must be 'single', "
+                    "'particle' must be 'Fickian diffusion'."
                 )
+        elif options["surface form"] == "algebraic":
+            pybamm.logger.warning(
+                "Surface form 'algebraic' is not necessary for single-phase "
+                "configurations ('particle phases' = '1' or ('1', '1')). "
+                "Consider using 'false' for equivalent results with simpler computation."
+            )
 
         if options["surface temperature"] == "lumped":
             if options["thermal"] not in ["isothermal", "lumped"]:

src/pybamm/models/full_battery_models/lithium_ion/spme.py

@@ -80,3 +80,44 @@
                 self.submodels[f"{domain} surface potential difference [V]"] = (
                     surf_model(self.param, domain, self.options)
                 )
+
+        if self.options["surface form"] in ["false", "algebraic"]:
+            for domain in ["negative", "positive"]:
+                if self.options.electrode_types[domain] == "porous":
+                    phi_s = self.variables.get(
+                        f"{domain.capitalize()} electrode potential [V]"
+                    )
+                    phi_e = self.variables.get(
+                        f"{domain.capitalize()} electrolyte potential [V]"
+                    )
+                    delta_phi = self.variables.get(
+                        f"{domain} surface potential difference [V]"
+                    )
+                    delta_phi_av = self.variables.get(
+                        f"X-averaged {domain} surface potential difference [V]"
+                    )
+                    if all([phi_s, phi_e, delta_phi, delta_phi_av]):
+                        pybamm.logger.debug(
+                            f"{domain.capitalize()} electrode for surface form "
+                            f"'{self.options['surface form']}': "
+                            f"phi_s={phi_s}, phi_e={phi_e}, delta_phi={delta_phi}, "
+                            f"delta_phi_av={delta_phi_av}"
+                        )
+            voltage = self.variables.get("Terminal voltage [V]")
+            if voltage:
+                pybamm.logger.debug(
+                    f"Terminal voltage for surface form '{self.options['surface form']}': {voltage}"
+                )
+            phases = self.options["particle phases"]
+            if phases == "1" or phases == ("1", "1"):
+                pybamm.logger.info(
+                    "Surface forms 'false' and 'algebraic' should produce identical "
+                    "results for single-phase configurations. If differences are observed, "
+                    "check logs for phi_s, phi_e, delta_phi, and current densities."
+                )
+            else:
+                pybamm.logger.warning(
+                    "Surface form 'algebraic' is designed for multi-phase configurations. "
+                    "Differences with 'false' may occur due to algebraic constraints in "
+                    "CompositeAlgebraic. Verify with logs if unexpected."
+                )

src/pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/composite_surface_form_conductivity.py

@@ -157,4 +157,24 @@ def set_algebraic(self, variables):
             f"X-averaged {domain} electrode surface potential difference [V]"
         ]
 
-        self.algebraic[delta_phi] = (sum_a_j_av - sum_a_j) / self.param.a_j_scale
+        pybamm.logger.debug(
+            f"CompositeAlgebraic ({domain}): "
+            f"delta_phi={delta_phi}, sum_a_j={sum_a_j}, sum_a_j_av={sum_a_j_av}"
+        )
+
+        phases = self.options.get("particle phases", "1")
+        if phases == "1" or phases == ("1", "1"):
+            phi_s = variables[f"{domain.capitalize()} electrode potential [V]"]
+            phi_e = variables[f"{domain.capitalize()} electrolyte potential [V]"]
+            delta_phi_calc = phi_s - phi_e
+            self.algebraic[delta_phi] = delta_phi - delta_phi_calc
+            pybamm.logger.debug(
+                f"CompositeAlgebraic ({domain}): Single-phase detected, using "
+                f"delta_phi = phi_s - phi_e = {delta_phi_calc}"
+            )
+        else:
+            self.algebraic[delta_phi] = (sum_a_j_av - sum_a_j) / self.param.a_j_scale
+            pybamm.logger.debug(
+                f"CompositeAlgebraic ({domain}): Multi-phase detected, using "
+                f"algebraic constraint (sum_a_j_av - sum_a_j) / a_j_scale"
+            )

src/pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/explicit_surface_form_conductivity.py

@@ -33,6 +33,12 @@ def get_coupled_variables(self, variables):
         phi_s = variables[f"{Domain} electrode potential [V]"]
         phi_e = variables[f"{Domain} electrolyte potential [V]"]
         delta_phi = phi_s - phi_e
+
+        pybamm.logger.debug(
+            f"Explicit surface form ({self.domain}): "
+            f"phi_s={phi_s}, phi_e={phi_e}, delta_phi={delta_phi}"
+        )
+
         variables.update(
             self._get_standard_surface_potential_difference_variables(delta_phi)
         )