#746 reformat lithium-ion models

pybamm/expression_tree/binary_operators.py

@@ -77,19 +77,8 @@ class BinaryOperator(pybamm.Symbol):
     """
 
     def __init__(self, name, left, right):
-        # Turn numbers into scalars
-        if isinstance(left, numbers.Number):
-            left = pybamm.Scalar(left)
-        if isinstance(right, numbers.Number):
-            right = pybamm.Scalar(right)
+        left, right = self.format(left, right)
 
-        # Check both left and right are pybamm Symbols
-        if not (isinstance(left, pybamm.Symbol) and isinstance(right, pybamm.Symbol)):
-            raise NotImplementedError(
-                """'{}' not implemented for symbols of type {} and {}""".format(
-                    self.__class__.__name__, type(left), type(right)
-                )
-            )
         # Check and process domains, except for Outer symbol which takes the outer
         # product of two smbols in different domains, and gives it the domain of the
         # right child.
@@ -108,6 +97,38 @@ def __init__(self, name, left, right):
         self.left = self.children[0]
         self.right = self.children[1]
 
+    def format(self, left, right):
+        "Format children left and right into compatible form"
+        # Turn numbers into scalars
+        if isinstance(left, numbers.Number):
+            left = pybamm.Scalar(left)
+        if isinstance(right, numbers.Number):
+            right = pybamm.Scalar(right)
+
+        # Check both left and right are pybamm Symbols
+        if not (isinstance(left, pybamm.Symbol) and isinstance(right, pybamm.Symbol)):
+            raise NotImplementedError(
+                """'{}' not implemented for symbols of type {} and {}""".format(
+                    self.__class__.__name__, type(left), type(right)
+                )
+            )
+
+        # Do some broadcasting in special cases
+        if (
+            left.domain != right.domain
+            and "secondary" in right.auxiliary_domains
+            and left.domain == right.auxiliary_domains["secondary"]
+        ):
+            left = pybamm.PrimaryBroadcast(left, right.domain)
+        if (
+            right.domain != left.domain
+            and "secondary" in left.auxiliary_domains
+            and right.domain == left.auxiliary_domains["secondary"]
+        ):
+            right = pybamm.PrimaryBroadcast(right, left.domain)
+
+        return left, right
+
     def __str__(self):
         """ See :meth:`pybamm.Symbol.__str__()`. """
         return "{!s} {} {!s}".format(self.left, self.name, self.right)

pybamm/expression_tree/independent_variable.py

@@ -24,12 +24,14 @@ class IndependentVariable(pybamm.Symbol):
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=[]):
-        super().__init__(name, domain=domain)
+    def __init__(self, name, domain=None, auxiliary_domains=None):
+        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
 
     def evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
-        return pybamm.evaluate_for_shape_using_domain(self.domain)
+        return pybamm.evaluate_for_shape_using_domain(
+            self.domain, self.auxiliary_domains
+        )
 
     def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
@@ -77,9 +79,9 @@ class SpatialVariable(IndependentVariable):
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, coord_sys=None):
+    def __init__(self, name, domain=None, auxiliary_domains=None, coord_sys=None):
         self.coord_sys = coord_sys
-        super().__init__(name, domain=domain)
+        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
         domain = self.domain
 
         if name not in KNOWN_SPATIAL_VARS:
@@ -109,7 +111,9 @@ def __init__(self, name, domain=None, coord_sys=None):
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return SpatialVariable(self.name, self.domain, self.coord_sys)
+        return SpatialVariable(
+            self.name, self.domain, self.auxiliary_domains, self.coord_sys
+        )
 
 
 # the independent variable time

pybamm/geometry/standard_spatial_vars.py

@@ -4,35 +4,77 @@
 
 # Domains at cell centres
 x_n = pybamm.SpatialVariable(
-    "x_n", domain=["negative electrode"], coord_sys="cartesian"
+    "x_n",
+    domain=["negative electrode"],
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
+)
+x_s = pybamm.SpatialVariable(
+    "x_s",
+    domain=["separator"],
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
 )
-x_s = pybamm.SpatialVariable("x_s", domain=["separator"], coord_sys="cartesian")
 x_p = pybamm.SpatialVariable(
-    "x_p", domain=["positive electrode"], coord_sys="cartesian"
+    "x_p",
+    domain=["positive electrode"],
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
+)
+x = pybamm.SpatialVariable(
+    "x",
+    domain=whole_cell,
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
 )
-x = pybamm.SpatialVariable("x", domain=whole_cell, coord_sys="cartesian")
 
 y = pybamm.SpatialVariable("y", domain="current collector", coord_sys="cartesian")
 z = pybamm.SpatialVariable("z", domain="current collector", coord_sys="cartesian")
 
 r_n = pybamm.SpatialVariable(
-    "r_n", domain=["negative particle"], coord_sys="spherical polar"
+    "r_n",
+    domain=["negative particle"],
+    auxiliary_domains={
+        "secondary": "negative electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="spherical polar",
 )
 r_p = pybamm.SpatialVariable(
-    "r_p", domain=["positive particle"], coord_sys="spherical polar"
+    "r_p",
+    domain=["positive particle"],
+    auxiliary_domains={
+        "secondary": "positive electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="spherical polar",
 )
 
 # Domains at cell edges
 x_n_edge = pybamm.SpatialVariable(
-    "x_n_edge", domain=["negative electrode"], coord_sys="cartesian"
+    "x_n_edge",
+    domain=["negative electrode"],
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
 )
 x_s_edge = pybamm.SpatialVariable(
-    "x_s_edge", domain=["separator"], coord_sys="cartesian"
+    "x_s_edge",
+    domain=["separator"],
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
 )
 x_p_edge = pybamm.SpatialVariable(
-    "x_p_edge", domain=["positive electrode"], coord_sys="cartesian"
+    "x_p_edge",
+    domain=["positive electrode"],
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
+)
+x_edge = pybamm.SpatialVariable(
+    "x_edge",
+    domain=whole_cell,
+    auxiliary_domains={"secondary": "current collector"},
+    coord_sys="cartesian",
 )
-x_edge = pybamm.SpatialVariable("x_edge", domain=whole_cell, coord_sys="cartesian")
 
 y_edge = pybamm.SpatialVariable(
     "y_edge", domain="current collector", coord_sys="cartesian"
@@ -42,8 +84,20 @@
 )
 
 r_n_edge = pybamm.SpatialVariable(
-    "r_n_edge", domain=["negative particle"], coord_sys="spherical polar"
+    "r_n_edge",
+    domain=["negative particle"],
+    auxiliary_domains={
+        "secondary": "negative electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="spherical polar",
 )
 r_p_edge = pybamm.SpatialVariable(
-    "r_p_edge", domain=["positive particle"], coord_sys="spherical polar"
+    "r_p_edge",
+    domain=["positive particle"],
+    auxiliary_domains={
+        "secondary": "positive electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="spherical polar",
 )

pybamm/models/submodels/electrode/ohm/composite_ohm.py

@@ -41,12 +41,10 @@ def get_coupled_variables(self, variables):
 
         if self._domain == "Negative":
             sigma_eff_0 = self.param.sigma_n * tor_0
-            phi_s = pybamm.PrimaryBroadcast(
-                phi_s_cn, "negative electrode"
-            ) + pybamm.outer(
-                i_boundary_cc_0 / sigma_eff_0, x_n * (x_n - 2 * l_n) / (2 * l_n)
+            phi_s = phi_s_cn + (i_boundary_cc_0 / sigma_eff_0) * (
+                x_n * (x_n - 2 * l_n) / (2 * l_n)
             )
-            i_s = pybamm.outer(i_boundary_cc_0, 1 - x_n / l_n)
+            i_s = i_boundary_cc_0 * (1 - x_n / l_n)
 
         elif self.domain == "Positive":
             delta_phi_p_av = variables[
@@ -62,12 +60,10 @@ def get_coupled_variables(self, variables):
                 + (i_boundary_cc_0 / sigma_eff_0) * (1 - l_p / 3)
             )
 
-            phi_s = pybamm.PrimaryBroadcast(
-                const, ["positive electrode"]
-            ) - pybamm.outer(
-                i_boundary_cc_0 / sigma_eff_0, x_p + (x_p - 1) ** 2 / (2 * l_p)
+            phi_s = const - (i_boundary_cc_0 / sigma_eff_0) * (
+                x_p + (x_p - 1) ** 2 / (2 * l_p)
             )
-            i_s = pybamm.outer(i_boundary_cc_0, 1 - (1 - x_p) / l_p)
+            i_s = i_boundary_cc_0 * (1 - (1 - x_p) / l_p)
 
         variables.update(self._get_standard_potential_variables(phi_s))
         variables.update(self._get_standard_current_variables(i_s))

pybamm/models/submodels/electrode/ohm/leading_ohm.py

@@ -41,7 +41,7 @@ def get_coupled_variables(self, variables):
 
         if self.domain == "Negative":
             phi_s = pybamm.PrimaryBroadcast(phi_s_cn, "negative electrode")
-            i_s = pybamm.outer(i_boundary_cc, 1 - x_n / l_n)
+            i_s = i_boundary_cc * (1 - x_n / l_n)
 
         elif self.domain == "Positive":
             # recall delta_phi = phi_s - phi_e
@@ -53,7 +53,7 @@ def get_coupled_variables(self, variables):
             v = delta_phi_p_av + phi_e_p_av
 
             phi_s = pybamm.PrimaryBroadcast(v, ["positive electrode"])
-            i_s = pybamm.outer(i_boundary_cc, 1 - (1 - x_p) / l_p)
+            i_s = i_boundary_cc * (1 - (1 - x_p) / l_p)
 
         variables.update(self._get_standard_potential_variables(phi_s))
         variables.update(self._get_standard_current_variables(i_s))

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py

@@ -30,12 +30,6 @@ def unpack(self, variables):
         raise NotImplementedError
 
     def get_coupled_variables(self, variables):
-        # NOTE: the heavy use of Broadcast and outer in this method is mainly so
-        # that products are handled correctly when using 1 or 2D current collector
-        # models. In standard 1D battery models outer behaves as a normal multiply.
-        # In the future, multiply will automatically handle switching between
-        # normal multiply and outer products as appropriate.
-
         c_e_av = self.unpack(variables)
 
         i_boundary_cc_0 = variables["Leading-order current collector current density"]
@@ -78,9 +72,9 @@ def get_coupled_variables(self, variables):
             chi_av_p = chi_av
 
         # electrolyte current
-        i_e_n = pybamm.outer(i_boundary_cc_0, x_n / l_n)
+        i_e_n = i_boundary_cc_0 * x_n / l_n
         i_e_s = pybamm.PrimaryBroadcast(i_boundary_cc_0, "separator")
-        i_e_p = pybamm.outer(i_boundary_cc_0, (1 - x_p) / l_p)
+        i_e_p = i_boundary_cc_0 * (1 - x_p) / l_p
         i_e = pybamm.Concatenation(i_e_n, i_e_s, i_e_p)
 
         # electrolyte potential
@@ -90,11 +84,7 @@ def get_coupled_variables(self, variables):
             - (
                 chi_av
                 * (1 + param.Theta * T_av)
-                * pybamm.x_average(
-                    self._higher_order_macinnes_function(
-                        c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode")
-                    )
-                )
+                * pybamm.x_average(self._higher_order_macinnes_function(c_e_n / c_e_av))
             )
             - (
                 (i_boundary_cc_0 * param.C_e * l_n / param.gamma_e)
@@ -103,55 +93,39 @@ def get_coupled_variables(self, variables):
         )
 
         phi_e_n = (
-            pybamm.PrimaryBroadcast(phi_e_const, "negative electrode")
+            phi_e_const
             + (
                 chi_av_n
                 * (1 + param.Theta * T_av_n)
-                * self._higher_order_macinnes_function(
-                    c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode")
-                )
-            )
-            - pybamm.outer(
-                i_boundary_cc_0 * (param.C_e / param.gamma_e) / kappa_n_av,
-                (x_n ** 2 - l_n ** 2) / (2 * l_n),
-            )
-            - pybamm.PrimaryBroadcast(
-                i_boundary_cc_0 * l_n * (param.C_e / param.gamma_e) / kappa_s_av,
-                "negative electrode",
+                * self._higher_order_macinnes_function(c_e_n / c_e_av)
             )
+            - (i_boundary_cc_0 * (param.C_e / param.gamma_e) / kappa_n_av)
+            * (x_n ** 2 - l_n ** 2)
+            / (2 * l_n)
+            - i_boundary_cc_0 * l_n * (param.C_e / param.gamma_e) / kappa_s_av
         )
 
         phi_e_s = (
-            pybamm.PrimaryBroadcast(phi_e_const, "separator")
+            phi_e_const
             + (
                 chi_av_s
                 * (1 + param.Theta * T_av_s)
-                * self._higher_order_macinnes_function(
-                    c_e_s / pybamm.PrimaryBroadcast(c_e_av, "separator")
-                )
-            )
-            - pybamm.outer(
-                i_boundary_cc_0 * param.C_e / param.gamma_e / kappa_s_av, x_s
+                * self._higher_order_macinnes_function(c_e_s / c_e_av)
             )
+            - (i_boundary_cc_0 * param.C_e / param.gamma_e / kappa_s_av) * x_s
         )
 
         phi_e_p = (
-            pybamm.PrimaryBroadcast(phi_e_const, "positive electrode")
+            phi_e_const
             + (
                 chi_av_p
                 * (1 + param.Theta * T_av_p)
-                * self._higher_order_macinnes_function(
-                    c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode")
-                )
-            )
-            - pybamm.outer(
-                i_boundary_cc_0 * (param.C_e / param.gamma_e) / kappa_p_av,
-                (x_p * (2 - x_p) + l_p ** 2 - 1) / (2 * l_p),
-            )
-            - pybamm.PrimaryBroadcast(
-                i_boundary_cc_0 * (1 - l_p) * (param.C_e / param.gamma_e) / kappa_s_av,
-                "positive electrode",
+                * self._higher_order_macinnes_function(c_e_p / c_e_av)
             )
+            - (i_boundary_cc_0 * (param.C_e / param.gamma_e) / kappa_p_av)
+            * (x_p * (2 - x_p) + l_p ** 2 - 1)
+            / (2 * l_p)
+            - i_boundary_cc_0 * (1 - l_p) * (param.C_e / param.gamma_e) / kappa_s_av
         )
 
         phi_e = pybamm.Concatenation(phi_e_n, phi_e_s, phi_e_p)
@@ -162,16 +136,8 @@ def get_coupled_variables(self, variables):
             chi_av
             * (1 + param.Theta * T_av)
             * (
-                pybamm.x_average(
-                    self._higher_order_macinnes_function(
-                        c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode")
-                    )
-                )
-                - pybamm.x_average(
-                    self._higher_order_macinnes_function(
-                        c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode")
-                    )
-                )
+                pybamm.x_average(self._higher_order_macinnes_function(c_e_p / c_e_av))
+                - pybamm.x_average(self._higher_order_macinnes_function(c_e_n / c_e_av))
             )
         )