Issue 746 broadcasting

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,43 @@ 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, to avoid having to do this manually
+        if (
+            not isinstance(self, (Outer, Kron))
+            and left.domain != []
+            and right.domain != []
+        ):
+            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)
@@ -641,14 +667,12 @@ class Outer(BinaryOperator):
 
     def __init__(self, left, right):
         """ See :meth:`pybamm.BinaryOperator.__init__()`. """
-        # cannot have Variable, StateVector or Matrix in the right symbol, as these
+        # cannot have certain types of objects in the right symbol, as these
         # can already be 2D objects (so we can't take an outer product with them)
         if right.has_symbol_of_classes(
-            (pybamm.Variable, pybamm.StateVector, pybamm.Matrix)
+            (pybamm.Variable, pybamm.StateVector, pybamm.Matrix, pybamm.SpatialVariable)
         ):
-            raise TypeError(
-                "right child must only contain SpatialVariable and scalars" ""
-            )
+            raise TypeError("right child must only contain Vectors and Scalars" "")
 
         super().__init__("outer product", left, 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/convection/base_convection.py

@@ -135,8 +135,6 @@ def _separator_velocity(self, variables):
                 auxiliary_domains={"secondary": "current collector"},
             ),
         )
-        v_box_s = pybamm.outer(d_vbox_s__dx, (x_s - l_n)) + pybamm.PrimaryBroadcast(
-            v_box_n_right, "separator"
-        )
+        v_box_s = d_vbox_s__dx * (x_s - l_n) + v_box_n_right
 
         return v_box_s, dVbox_dz

pybamm/models/submodels/convection/composite_convection.py

@@ -36,8 +36,8 @@ def get_coupled_variables(self, variables):
         j_p_av = variables["X-averaged positive electrode interfacial current density"]
 
         # Volume-averaged velocity
-        v_box_n = param.beta_n * pybamm.outer(j_n_av, x_n)
-        v_box_p = param.beta_p * pybamm.outer(j_p_av, x_p - 1)
+        v_box_n = param.beta_n * j_n_av * x_n
+        v_box_p = param.beta_p * j_p_av * (x_p - 1)
 
         v_box_s, dVbox_dz = self._separator_velocity(variables)
         v_box = pybamm.Concatenation(v_box_n, v_box_s, v_box_p)

pybamm/models/submodels/convection/leading_convection.py

@@ -30,8 +30,8 @@ def get_coupled_variables(self, variables):
         j_p_av = variables["X-averaged positive electrode interfacial current density"]
 
         # Volume-averaged velocity
-        v_box_n = param.beta_n * pybamm.outer(j_n_av, x_n)
-        v_box_p = param.beta_p * pybamm.outer(j_p_av, x_p - 1)
+        v_box_n = param.beta_n * j_n_av * x_n
+        v_box_p = param.beta_p * j_p_av * (x_p - 1)
 
         v_box_s, dVbox_dz = self._separator_velocity(variables)
         v_box = pybamm.Concatenation(v_box_n, v_box_s, v_box_p)

pybamm/models/submodels/electrode/base_electrode.py

@@ -54,7 +54,7 @@ def _get_standard_potential_variables(self, phi_s):
             )
 
             v = pybamm.boundary_value(phi_s, "right")
-            delta_phi_s = phi_s - pybamm.PrimaryBroadcast(v, ["positive electrode"])
+            delta_phi_s = phi_s - v
         delta_phi_s_av = pybamm.x_average(delta_phi_s)
         delta_phi_s_dim = delta_phi_s * param.potential_scale
         delta_phi_s_av_dim = delta_phi_s_av * param.potential_scale

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/electrode/ohm/surface_form_ohm.py

@@ -34,26 +34,21 @@ def get_coupled_variables(self, variables):
         tor = variables[self.domain + " electrode tortuosity"]
         phi_s_cn = variables["Negative current collector potential"]
 
-        i_s = pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast) - i_e
+        i_s = i_boundary_cc - i_e
 
         if self.domain == "Negative":
             conductivity = param.sigma_n * tor
-            phi_s = pybamm.PrimaryBroadcast(
-                phi_s_cn, "negative electrode"
-            ) - pybamm.IndefiniteIntegral(i_s / conductivity, x_n)
+            phi_s = phi_s_cn - pybamm.IndefiniteIntegral(i_s / conductivity, x_n)
 
         elif self.domain == "Positive":
 
             phi_e_s = variables["Separator electrolyte potential"]
             delta_phi_p = variables["Positive electrode surface potential difference"]
 
             conductivity = param.sigma_p * tor
-            phi_s = -pybamm.IndefiniteIntegral(
-                i_s / conductivity, x_p
-            ) + pybamm.PrimaryBroadcast(
+            phi_s = -pybamm.IndefiniteIntegral(i_s / conductivity, x_p) + (
                 pybamm.boundary_value(phi_e_s, "right")
-                + pybamm.boundary_value(delta_phi_p, "left"),
-                "positive electrode",
+                + pybamm.boundary_value(delta_phi_p, "left")
             )
 
         variables.update(self._get_standard_potential_variables(phi_s))