#746 remove outers and some primary broadcasts

pybamm/expression_tree/binary_operators.py

@@ -116,12 +116,16 @@ def format(self, left, right):
         # Do some broadcasting in special cases
         if (
             left.domain != right.domain
+            and left.domain != []
+            and 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 left.domain != []
+            and right.domain != []
             and "secondary" in left.auxiliary_domains
             and right.domain == left.auxiliary_domains["secondary"]
         ):

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/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))

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py

@@ -162,8 +162,7 @@ def _get_neg_pos_coupled_variables(self, variables):
         i_e = conductivity * (
             ((1 + param.Theta * T) * param.chi(c_e) / c_e) * pybamm.grad(c_e)
             + pybamm.grad(delta_phi)
-            + pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast)
-            / sigma_eff
+            + i_boundary_cc / sigma_eff
         )
         variables.update(self._get_domain_current_variables(i_e))
 
@@ -193,13 +192,9 @@ def _get_sep_coupled_variables(self, variables):
         chi_e_s = param.chi(c_e_s)
         kappa_s_eff = param.kappa_e(c_e_s, T) * tor_s
 
-        phi_e_s = pybamm.PrimaryBroadcast(
-            pybamm.boundary_value(phi_e_n, "right"), "separator"
-        ) + pybamm.IndefiniteIntegral(
+        phi_e_s = pybamm.boundary_value(phi_e_n, "right") + pybamm.IndefiniteIntegral(
             (1 + param.Theta * T) * chi_e_s / c_e_s * pybamm.grad(c_e_s)
-            - param.C_e
-            * pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast)
-            / kappa_s_eff,
+            - param.C_e * i_boundary_cc / kappa_s_eff,
             x_s,
         )
 
@@ -227,7 +222,7 @@ def nasty_hack_to_get_phi_s(self, variables):
         i_boundary_cc = variables["Current collector current density"]
         i_e = variables[self.domain + " electrolyte current density"]
 
-        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
@@ -240,12 +235,9 @@ def nasty_hack_to_get_phi_s(self, variables):
 
             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")
             )
 
         return phi_s

pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/composite_stefan_maxwell_diffusion.py

@@ -38,12 +38,7 @@ def get_coupled_variables(self, variables):
 
         param = self.param
 
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        N_e_diffusion = (
-            -tor_0
-            * pybamm.PrimaryBroadcast(param.D_e(c_e_0_av, T_0), whole_cell)
-            * pybamm.grad(c_e)
-        )
+        N_e_diffusion = -tor_0 * param.D_e(c_e_0_av, T_0) * pybamm.grad(c_e)
         # N_e_migration = (param.C_e * param.t_plus) / param.gamma_e * i_e
         # N_e_convection = c_e * v_box_0
 
@@ -52,7 +47,7 @@ def get_coupled_variables(self, variables):
         if v_box_0.id == pybamm.Scalar(0).id:
             N_e = N_e_diffusion
         else:
-            N_e = N_e_diffusion + pybamm.outer(v_box_0, c_e)
+            N_e = N_e_diffusion + v_box_0 * c_e
 
         variables.update(self._get_standard_flux_variables(N_e))
 

pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/first_order_stefan_maxwell_diffusion.py

@@ -78,23 +78,15 @@ def get_coupled_variables(self, variables):
         D_e_p = tor_p_0 * param.D_e(c_e_0, T_0)
 
         # Fluxes
-        N_e_n_1 = -pybamm.outer(rhs_n, x_n)
-        N_e_s_1 = -(
-            pybamm.outer(rhs_s, (x_s - l_n))
-            + pybamm.PrimaryBroadcast(rhs_n * l_n, "separator")
-        )
-        N_e_p_1 = -pybamm.outer(rhs_p, (x_p - 1))
+        N_e_n_1 = -rhs_n * x_n
+        N_e_s_1 = -(rhs_s * (x_s - l_n) + rhs_n * l_n)
+        N_e_p_1 = -rhs_p * (x_p - 1)
 
         # Concentrations
-        c_e_n_1 = pybamm.outer(rhs_n / (2 * D_e_n), x_n ** 2 - l_n ** 2)
-        c_e_s_1 = pybamm.outer(rhs_s / 2, (x_s - l_n) ** 2) + pybamm.outer(
-            rhs_n * l_n / D_e_s, x_s - l_n
-        )
-        c_e_p_1 = pybamm.outer(
-            rhs_p / (2 * D_e_p), (x_p - 1) ** 2 - l_p ** 2
-        ) + pybamm.PrimaryBroadcast(
-            (rhs_s * l_s ** 2 / (2 * D_e_s)) + (rhs_n * l_n * l_s / D_e_s),
-            "positive electrode",
+        c_e_n_1 = (rhs_n / (2 * D_e_n)) * (x_n ** 2 - l_n ** 2)
+        c_e_s_1 = (rhs_s / 2) * ((x_s - l_n) ** 2) + (rhs_n * l_n / D_e_s) * (x_s - l_n)
+        c_e_p_1 = (rhs_p / (2 * D_e_p)) * ((x_p - 1) ** 2 - l_p ** 2) + (
+            (rhs_s * l_s ** 2 / (2 * D_e_s)) + (rhs_n * l_n * l_s / D_e_s)
         )
 
         # Correct for integral
@@ -108,18 +100,18 @@ def get_coupled_variables(self, variables):
         A_e = -(eps_n_0 * c_e_n_1_av + eps_s_0 * c_e_s_1_av + eps_p_0 * c_e_p_1_av) / (
             l_n * eps_n_0 + l_s * eps_s_0 + l_p * eps_p_0
         )
-        c_e_n_1 += pybamm.PrimaryBroadcast(A_e, "negative electrode")
-        c_e_s_1 += pybamm.PrimaryBroadcast(A_e, "separator")
-        c_e_p_1 += pybamm.PrimaryBroadcast(A_e, "positive electrode")
+        c_e_n_1 += A_e
+        c_e_s_1 += A_e
+        c_e_p_1 += A_e
         c_e_n_1_av += A_e
         c_e_s_1_av += A_e
         c_e_p_1_av += A_e
 
         # Update variables
         c_e = pybamm.Concatenation(
-            pybamm.PrimaryBroadcast(c_e_0, "negative electrode") + param.C_e * c_e_n_1,
-            pybamm.PrimaryBroadcast(c_e_0, "separator") + param.C_e * c_e_s_1,
-            pybamm.PrimaryBroadcast(c_e_0, "positive electrode") + param.C_e * c_e_p_1,
+            c_e_0 + param.C_e * c_e_n_1,
+            c_e_0 + param.C_e * c_e_s_1,
+            c_e_0 + param.C_e * c_e_p_1,
         )
         variables.update(self._get_standard_concentration_variables(c_e))
         # Update with analytical expressions for first-order x-averages

pybamm/models/submodels/interface/kinetics/base_first_order_kinetics.py

@@ -1,7 +1,6 @@
 #
 # First-order Butler-Volmer kinetics
 #
-import pybamm
 from .base_kinetics import BaseModel
 
 
@@ -52,11 +51,7 @@ def get_coupled_variables(self, variables):
             + self.reaction_name
             + " interfacial current density"
         ]
-        j_1 = (
-            pybamm.PrimaryBroadcast(dj_dc_0, self.domain_for_broadcast) * c_e_1
-            + pybamm.PrimaryBroadcast(dj_ddeltaphi_0, self.domain_for_broadcast)
-            * delta_phi_1
-        )
+        j_1 = dj_dc_0 * c_e_1 + dj_ddeltaphi_0 * delta_phi_1
         j = j_0 + self.param.C_e * j_1
         # Get exchange-current density
         j0 = self._get_exchange_current_density(variables)

pybamm/models/submodels/oxygen_diffusion/first_order_oxygen_diffusion.py

@@ -58,14 +58,15 @@ def get_coupled_variables(self, variables):
         # Fluxes
         N_ox_n_1 = pybamm.FullBroadcast(0, "negative electrode", "current collector")
         N_ox_s_1 = -pybamm.PrimaryBroadcast(sj_ox_p * l_p, "separator")
-        N_ox_p_1 = pybamm.outer(sj_ox_p, x_p - 1)
+        N_ox_p_1 = sj_ox_p * (x_p - 1)
 
         # Concentrations
         c_ox_n_1 = pybamm.FullBroadcast(0, "negative electrode", "current collector")
-        c_ox_s_1 = pybamm.outer(sj_ox_p * l_p / D_ox_s, x_s - l_n)
-        c_ox_p_1 = pybamm.outer(
-            -sj_ox_p / (2 * D_ox_p), (x_p - 1) ** 2 - l_p ** 2
-        ) + pybamm.PrimaryBroadcast(sj_ox_p * l_p * l_s / D_ox_s, "positive electrode")
+        c_ox_s_1 = sj_ox_p * l_p / D_ox_s * (x_s - l_n)
+        c_ox_p_1 = (
+            -sj_ox_p / (2 * D_ox_p) * ((x_p - 1) ** 2 - l_p ** 2)
+            + sj_ox_p * l_p * l_s / D_ox_s
+        )
 
         # Update variables
         c_ox = pybamm.Concatenation(

pybamm/models/submodels/particle/fickian/fickian_single_particle.py

@@ -45,7 +45,6 @@ def get_coupled_variables(self, variables):
             variables["X-averaged " + self.domain.lower() + " electrode temperature"],
             [self.domain.lower() + " particle"],
         )
-
         N_s_xav = self._flux_law(c_s_xav, T_k_xav)
         N_s = pybamm.PrimaryBroadcast(N_s_xav, [self._domain.lower() + " electrode"])
 

pybamm/models/submodels/thermal/base_thermal.py

@@ -173,8 +173,8 @@ def _get_standard_coupled_variables(self, variables):
                 "X-averaged irreversible electrochemical heating [W.m-3]": Q_rxn_av
                 * Q_scale,
                 "Volume-averaged irreversible electrochemical heating": Q_rxn_vol_av,
-                "Volume-averaged irreversible electrochemical heating [W.m-3]": Q_rxn_vol_av
-                * Q_scale,
+                "Volume-averaged irreversible electrochemical heating "
+                + "[W.m-3]": Q_rxn_vol_av * Q_scale,
                 "Reversible heating": Q_rev,
                 "Reversible heating [W.m-3]": Q_rev * Q_scale,
                 "X-averaged reversible heating": Q_rev_av,

pybamm/spatial_methods/finite_volume.py

@@ -639,9 +639,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         # https://github.com/Scottmar93/extrapolation-coefficents/tree/master
         if isinstance(symbol, pybamm.BoundaryValue):
 
-            if use_bcs and pybamm.has_bc_of_form(
-                child, symbol.side, bcs, "Dirichlet"
-            ):
+            if use_bcs and pybamm.has_bc_of_form(child, symbol.side, bcs, "Dirichlet"):
                 # just use the value from the bc: f(x*)
                 sub_matrix = csr_matrix((1, prim_pts))
                 additive = bcs[child.id][symbol.side][0]
@@ -655,7 +653,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                     if use_bcs and pybamm.has_bc_of_form(
                         child, symbol.side, bcs, "Neumann"
                     ):
-                        sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, prim_pts),)
+                        sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, prim_pts))
 
                         additive = -dx0 * bcs[child.id][symbol.side][0]
 
@@ -676,7 +674,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         alpha = -(dx0 * (dx0 + dx1)) / (2 * dx0 + dx1)
 
                         sub_matrix = csr_matrix(
-                            ([a, b], ([0, 0], [0, 1])), shape=(1, prim_pts),
+                            ([a, b], ([0, 0], [0, 1])), shape=(1, prim_pts)
                         )
                         additive = alpha * bcs[child.id][symbol.side][0]
 
@@ -686,7 +684,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         c = dx0 * (dx0 + dx1) / (dx2 * (dx1 + dx2))
 
                         sub_matrix = csr_matrix(
-                            ([a, b, c], ([0, 0, 0], [0, 1, 2])), shape=(1, prim_pts),
+                            ([a, b, c], ([0, 0, 0], [0, 1, 2])), shape=(1, prim_pts)
                         )
 
                         additive = pybamm.Scalar(0)
@@ -703,7 +701,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         # use formula:
                         # f(x*) = fN + dxN * f'(x*)
                         sub_matrix = csr_matrix(
-                            ([1], ([0], [prim_pts - 1]),), shape=(1, prim_pts),
+                            ([1], ([0], [prim_pts - 1])), shape=(1, prim_pts)
                         )
                         additive = dxN * bcs[child.id][symbol.side][0]
 
@@ -727,7 +725,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         b = -(dxN ** 2) / (2 * dxN * dxNm1 + dxNm1 ** 2)
                         alpha = dxN * (dxN + dxNm1) / (2 * dxN + dxNm1)
                         sub_matrix = csr_matrix(
-                            ([b, a], ([0, 0], [prim_pts - 2, prim_pts - 1]),),
+                            ([b, a], ([0, 0], [prim_pts - 2, prim_pts - 1])),
                             shape=(1, prim_pts),
                         )
 
@@ -755,9 +753,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
 
         elif isinstance(symbol, pybamm.BoundaryGradient):
 
-            if use_bcs and pybamm.has_bc_of_form(
-                child, symbol.side, bcs, "Neumann"
-            ):
+            if use_bcs and pybamm.has_bc_of_form(child, symbol.side, bcs, "Neumann"):
                 # just use the value from the bc: f'(x*)
                 sub_matrix = csr_matrix((1, prim_pts))
                 additive = bcs[child.id][symbol.side][0]