#546 some progress with Newman-Tiedemann

pybamm/expression_tree/unary_operators.py

@@ -449,7 +449,7 @@ class IndefiniteIntegral(Integral):
     **Extends:** :class:`Integral`
     """
 
-    def __init__(self, child, integration_variable):
+    def __init__(self, child, integration_variable, domain=None):
         if isinstance(integration_variable, list):
             if len(integration_variable) > 1:
                 raise NotImplementedError(

pybamm/models/standard_variables.py

@@ -122,9 +122,18 @@
 )
 
 # Electrolyte pressure
-pressure_n = pybamm.Variable("Negative electrolyte pressure", ["negative electrode"])
-pressure_s = pybamm.Variable("Separator electrolyte pressure", ["separator"])
-pressure_p = pybamm.Variable("Positive electrolyte pressure", ["positive electrode"])
+pressure_n = pybamm.SecondaryBroadcast(
+    pybamm.Variable("Negative electrolyte pressure", ["negative electrode"]),
+    "current collector",
+)
+pressure_s = pybamm.SecondaryBroadcast(
+    pybamm.Variable("Separator electrolyte pressure", ["separator"]),
+    "current collector",
+)
+pressure_p = pybamm.SecondaryBroadcast(
+    pybamm.Variable("Positive electrolyte pressure", ["positive electrode"]),
+    "current collector",
+)
 pressure = pybamm.Concatenation(pressure_n, pressure_s, pressure_p)
 
 # Temperature

pybamm/models/submodels/convection/base_convection.py

@@ -123,9 +123,15 @@ def _separator_velocity(self, variables):
 
         # Simple formula for velocity in the separator
         dVbox_dz = pybamm.Concatenation(
-            pybamm.Broadcast(0, "negative electrode"),
-            pybamm.PrimaryBroadcast(-d_vbox_s__dx, "separator"),
-            pybamm.Broadcast(0, "positive electrode"),
+            pybamm.SecondaryBroadcast(
+                pybamm.Broadcast(0, "negative electrode"), "current collector"
+            ),
+            pybamm.SecondaryBroadcast(
+                pybamm.PrimaryBroadcast(-d_vbox_s__dx, "separator"), "current collector"
+            ),
+            pybamm.SecondaryBroadcast(
+                pybamm.Broadcast(0, "positive electrode"), "current collector"
+            ),
         )
         v_box_s = pybamm.outer(d_vbox_s__dx, (x_s - l_n)) + pybamm.PrimaryBroadcast(
             v_box_n_right, "separator"

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

@@ -33,6 +33,9 @@ def get_coupled_variables(self, variables):
         i_e = variables[self.domain + " electrolyte current density"]
         eps = variables[self.domain + " electrode porosity"]
 
+        if isinstance(i_boundary_cc, pybamm.Broadcast):
+            i_boundary_cc = i_boundary_cc.orphans[0]
+
         i_s = i_boundary_cc - i_e
 
         if self.domain == "Negative":
@@ -45,10 +48,12 @@ def get_coupled_variables(self, variables):
             delta_phi_p = variables["Positive electrode surface potential difference"]
 
             conductivity = param.sigma_p * (1 - eps) ** param.b
-            phi_s = (
-                -pybamm.IndefiniteIntegral(i_s / conductivity, x_p)
-                + pybamm.boundary_value(phi_e_s, "right")
-                + pybamm.boundary_value(delta_phi_p, "left")
+            phi_s = -pybamm.IndefiniteIntegral(
+                i_s / conductivity, x_p
+            ) + pybamm.PrimaryBroadcast(
+                pybamm.boundary_value(phi_e_s, "right")
+                + pybamm.boundary_value(delta_phi_p, "left"),
+                "positive electrode",
             )
 
         variables.update(self._get_standard_potential_variables(phi_s))

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

@@ -51,7 +51,9 @@ def set_algebraic(self, variables):
         sum_j = sum(
             pybamm.Concatenation(
                 variables[reaction["Negative"]["aj"]],
-                pybamm.Broadcast(0, "separator"),
+                pybamm.SecondaryBroadcast(
+                    pybamm.Broadcast(0, "separator"), "current collector"
+                ),
                 variables[reaction["Positive"]["aj"]],
             )
             for reaction in self.reactions.values()

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

@@ -74,6 +74,9 @@ def set_boundary_conditions(self, variables):
         c_e = variables[self.domain + " electrolyte concentration"]
         delta_phi = variables[self.domain + " electrode surface potential difference"]
 
+        if isinstance(i_boundary_cc, pybamm.Broadcast):
+            i_boundary_cc = i_boundary_cc.orphans[0]
+
         if self.domain == "Negative":
             c_e_flux = pybamm.BoundaryFlux(c_e, "right")
             flux_left = -pybamm.BoundaryValue(i_boundary_cc / sigma_eff, "left")
@@ -151,6 +154,9 @@ def _get_neg_pos_coupled_variables(self, variables):
         c_e = variables[self.domain + " electrolyte concentration"]
         delta_phi = variables[self.domain + " electrode surface potential difference"]
 
+        if isinstance(i_boundary_cc, pybamm.Broadcast):
+            i_boundary_cc = i_boundary_cc.orphans[0]
+
         i_e = conductivity * (
             (param.chi(c_e) / c_e) * pybamm.grad(c_e)
             + pybamm.grad(delta_phi)
@@ -180,10 +186,15 @@ def _get_sep_coupled_variables(self, variables):
         phi_e_n = variables["Negative electrolyte potential"]
         eps_s = variables["Separator porosity"]
 
+        if isinstance(i_boundary_cc, pybamm.Broadcast):
+            i_boundary_cc = i_boundary_cc.orphans[0]
+
         chi_e_s = param.chi(c_e_s)
         kappa_s_eff = param.kappa_e(c_e_s) * (eps_s ** param.b)
 
-        phi_e_s = pybamm.boundary_value(phi_e_n, "right") + pybamm.IndefiniteIntegral(
+        phi_e_s = pybamm.PrimaryBroadcast(
+            pybamm.boundary_value(phi_e_n, "right"), "separator"
+        ) + pybamm.IndefiniteIntegral(
             chi_e_s / c_e_s * pybamm.grad(c_e_s)
             - param.C_e * i_boundary_cc / kappa_s_eff,
             x_s,
@@ -213,6 +224,9 @@ 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"]
 
+        if isinstance(i_boundary_cc, pybamm.Broadcast):
+            i_boundary_cc = i_boundary_cc.orphans[0]
+
         i_s = i_boundary_cc - i_e
 
         if self.domain == "Negative":
@@ -225,10 +239,13 @@ def nasty_hack_to_get_phi_s(self, variables):
             delta_phi_p = variables["Positive electrode surface potential difference"]
 
             conductivity = param.sigma_p * (1 - eps) ** param.b
-            phi_s = (
-                -pybamm.IndefiniteIntegral(i_s / conductivity, x_p)
-                + pybamm.boundary_value(phi_e_s, "right")
-                + pybamm.boundary_value(delta_phi_p, "left")
+
+            phi_s = -pybamm.IndefiniteIntegral(
+                i_s / conductivity, x_p
+            ) + pybamm.PrimaryBroadcast(
+                pybamm.boundary_value(phi_e_s, "right")
+                + pybamm.boundary_value(delta_phi_p, "left"),
+                "positive electrode",
             )
 
         return phi_s

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

@@ -60,7 +60,6 @@ def set_rhs(self, variables):
         N_e = variables["Electrolyte flux"]
         # i_e = variables["Electrolyte current density"]
 
-        # TODO: check lead acid version in new form
         # source_term = ((param.s - param.t_plus) / param.gamma_e) * pybamm.div(i_e)
         # source_term = pybamm.div(i_e) / param.gamma_e  # lithium-ion
         source_terms = sum(

pybamm/models/submodels/interface/lead_acid_oxygen.py

@@ -42,7 +42,7 @@ def _get_exchange_current_density(self, variables):
         """
         c_e = variables[self.domain + " electrolyte concentration"]
         # If c_e was broadcast, take only the orphan
-        if isinstance(c_e, pybamm.Broadcast):
+        if isinstance(c_e, pybamm.Broadcast) and c_e.broadcast_type != "secondary":
             c_e = c_e.orphans[0]
 
         if self.domain == "Negative":

pybamm/models/submodels/porosity/full_reaction_driven_porosity.py

@@ -32,7 +32,9 @@ def get_coupled_variables(self, variables):
         j_p = variables["Positive electrode interfacial current density"]
 
         deps_dt_n = -self.param.beta_surf_n * j_n
-        deps_dt_s = pybamm.Broadcast(0, ["separator"])
+        deps_dt_s = pybamm.SecondaryBroadcast(
+            pybamm.Broadcast(0, ["separator"]), "current collector"
+        )
         deps_dt_p = -self.param.beta_surf_p * j_p
 
         deps_dt = pybamm.Concatenation(deps_dt_n, deps_dt_s, deps_dt_p)