Repository: pybamm-team/PyBaMM

Issue:
Average current densities
See #510 comments. 
Separate issue for how we deal with average interfacial current densities.

Discussion:

Patch:
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index c4d0e3fe9..5820edd3d 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -281,22 +281,30 @@ class BaseBatteryModel(pybamm.BaseModel):
     def build_model(self):
 
         # Get the fundamental variables
-        for submodel in self.submodels.values():
+        for name, submodel in self.submodels.items():
+            pybamm.logger.debug("Getting fundamental variables for {}".format(name))
             self.variables.update(submodel.get_fundamental_variables())
 
         # Get coupled variables
-        for submodel in self.submodels.values():
+        for name, submodel in self.submodels.items():
+            pybamm.logger.debug("Getting coupled variables for {}".format(name))
             self.variables.update(submodel.get_coupled_variables(self.variables))
 
-        # Set model equations
-        for submodel in self.submodels.values():
+            # Set model equations
+        for name, submodel in self.submodels.items():
+            pybamm.logger.debug("Setting rhs for {}".format(name))
             submodel.set_rhs(self.variables)
+            pybamm.logger.debug("Setting algebraic for {}".format(name))
             submodel.set_algebraic(self.variables)
+            pybamm.logger.debug("Setting boundary conditions for {}".format(name))
             submodel.set_boundary_conditions(self.variables)
+            pybamm.logger.debug("Setting initial conditions for {}".format(name))
             submodel.set_initial_conditions(self.variables)
             submodel.set_events(self.variables)
+            pybamm.logger.debug("Updating {}".format(name))
             self.update(submodel)
 
+        pybamm.logger.debug("Setting voltage variables")
         self.set_voltage_variables()
 
     def set_thermal_submodel(self):
diff --git a/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py b/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
index 49412f4f1..365bed078 100644
--- a/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/base_electrolyte_diffusion.py
@@ -19,7 +19,7 @@ class BaseElectrolyteDiffusion(pybamm.BaseSubModel):
     def __init__(self, param):
         super().__init__(param)
 
-    def _get_standard_concentration_variables(self, c_e, c_e_av):
+    def _get_standard_concentration_variables(self, c_e):
         """
         A private function to obtain the standard variables which
         can be derived from the concentration in the electrolyte.
@@ -40,6 +40,7 @@ class BaseElectrolyteDiffusion(pybamm.BaseSubModel):
 
         c_e_typ = self.param.c_e_typ
         c_e_n, c_e_s, c_e_p = c_e.orphans
+        c_e_av = pybamm.average(c_e)
 
         variables = {
             "Electrolyte concentration": c_e,
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py
index 56691b6a5..a07163fdc 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py
@@ -96,8 +96,10 @@ class LeadingOrderDifferential(BaseLeadingOrderSurfaceForm):
         param = self.param
 
         j = variables[self.domain + " electrode interfacial current density"].orphans[0]
-        j_av = variables[
-            "Average " + self.domain.lower() + " electrode interfacial current density"
+        j_tot_av = variables[
+            "Average "
+            + self.domain.lower()
+            + " electrode total interfacial current density"
         ]
         delta_phi = variables[
             "Average " + self.domain.lower() + " electrode surface potential difference"
@@ -108,7 +110,7 @@ class LeadingOrderDifferential(BaseLeadingOrderSurfaceForm):
         elif self.domain == "Positive":
             C_dl = param.C_dl_p
 
-        self.rhs[delta_phi] = 1 / C_dl * (j_av - j)
+        self.rhs[delta_phi] = 1 / C_dl * (j_tot_av - j)
 
 
 class LeadingOrderAlgebraic(BaseLeadingOrderSurfaceForm):
@@ -133,11 +135,13 @@ class LeadingOrderAlgebraic(BaseLeadingOrderSurfaceForm):
             return
 
         j = variables[self.domain + " electrode interfacial current density"].orphans[0]
-        j_av = variables[
-            "Average " + self.domain.lower() + " electrode interfacial current density"
+        j_tot_av = variables[
+            "Average "
+            + self.domain.lower()
+            + " electrode total interfacial current density"
         ]
         delta_phi = variables[
             "Average " + self.domain.lower() + " electrode surface potential difference"
         ]
 
-        self.algebraic[delta_phi] = j_av - j
+        self.algebraic[delta_phi] = j_tot_av - j
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py
index d8f3add84..f6fdb14af 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/constant_stefan_maxwell_diffusion.py
@@ -24,13 +24,12 @@ class ConstantConcentration(BaseModel):
         super().__init__(param)
 
     def get_fundamental_variables(self):
-        c_e_av = pybamm.Scalar(1)
-        c_e_n = pybamm.Broadcast(c_e_av, ["negative electrode"])
-        c_e_s = pybamm.Broadcast(c_e_av, ["separator"])
-        c_e_p = pybamm.Broadcast(c_e_av, ["positive electrode"])
+        c_e_n = pybamm.Broadcast(1, ["negative electrode"])
+        c_e_s = pybamm.Broadcast(1, ["separator"])
+        c_e_p = pybamm.Broadcast(1, ["positive electrode"])
         c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
 
-        variables = self._get_standard_concentration_variables(c_e, c_e_av)
+        variables = self._get_standard_concentration_variables(c_e)
 
         N_e = pybamm.Broadcast(
             0, ["negative electrode", "separator", "positive electrode"]
@@ -42,4 +41,3 @@ class ConstantConcentration(BaseModel):
 
     def set_boundary_conditions(self, variables):
         return None
-
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
index 99d1c4e78..61b93fd16 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
@@ -25,9 +25,8 @@ class Full(BaseModel):
 
     def get_fundamental_variables(self):
         c_e = pybamm.standard_variables.c_e
-        c_e_av = pybamm.average(c_e)
 
-        return self._get_standard_concentration_variables(c_e, c_e_av)
+        return self._get_standard_concentration_variables(c_e)
 
     def get_coupled_variables(self, variables):
 
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py
index e28fd6f9d..00b5edd11 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/leading_stefan_maxwell_diffusion.py
@@ -31,9 +31,7 @@ class LeadingOrder(BaseModel):
         c_e_p = pybamm.Broadcast(c_e_av, ["positive electrode"])
         c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
 
-        c_e_n.evaluate_for_shape()
-
-        return self._get_standard_concentration_variables(c_e, c_e_av)
+        return self._get_standard_concentration_variables(c_e)
 
     def get_coupled_variables(self, variables):
 
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index cfbe69d4f..64df69b9c 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -21,25 +21,27 @@ class BaseInterface(pybamm.BaseSubModel):
     def __init__(self, param, domain):
         super().__init__(param, domain)
 
-    def _get_average_interfacial_current_density(self, variables):
+    def _get_average_total_interfacial_current_density(self, variables):
         """
-        Method to obtain the average interfacial current density.
+        Method to obtain the average total interfacial current density.
         """
 
         i_boundary_cc = variables["Current collector current density"]
 
         if self.domain == "Negative":
-            j_av = i_boundary_cc / pybamm.geometric_parameters.l_n
+            j_total_average = i_boundary_cc / pybamm.geometric_parameters.l_n
 
         elif self.domain == "Positive":
-            j_av = -i_boundary_cc / pybamm.geometric_parameters.l_p
+            j_total_average = -i_boundary_cc / pybamm.geometric_parameters.l_p
 
-        return j_av
+        return j_total_average
 
-    def _get_standard_interfacial_current_variables(self, j, j_av):
+    def _get_standard_interfacial_current_variables(self, j):
 
         i_typ = self.param.i_typ
-        # Broadcast if necessary
+
+        # Average, and broadcast if necessary
+        j_av = pybamm.average(j)
         if j.domain in [[], ["current collector"]]:
             j = pybamm.Broadcast(j, self.domain_for_broadcast)
 
@@ -56,6 +58,21 @@ class BaseInterface(pybamm.BaseSubModel):
 
         return variables
 
+    def _get_standard_total_interfacial_current_variables(self, j_tot_av):
+
+        i_typ = self.param.i_typ
+
+        variables = {
+            "Average "
+            + self.domain.lower()
+            + " electrode total interfacial current density": j_tot_av,
+            "Average "
+            + self.domain.lower()
+            + " electrode total interfacial current density [A.m-2]": i_typ * j_tot_av,
+        }
+
+        return variables
+
     def _get_standard_whole_cell_interfacial_current_variables(self, variables):
 
         i_typ = self.param.i_typ
diff --git a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
index c51941057..491f9c5be 100644
--- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
+++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py
@@ -29,8 +29,8 @@ class BaseInverseButlerVolmer(BaseInterface):
         ocp, dUdT = self._get_open_circuit_potential(variables)
 
         j0 = self._get_exchange_current_density(variables)
-        j_av = self._get_average_interfacial_current_density(variables)
-        j = pybamm.Broadcast(j_av, [self.domain.lower() + " electrode"])
+        j_tot_av = self._get_average_total_interfacial_current_density(variables)
+        j = pybamm.Broadcast(j_tot_av, [self.domain.lower() + " electrode"])
 
         if self.domain == "Negative":
             ne = self.param.ne_n
@@ -41,7 +41,10 @@ class BaseInverseButlerVolmer(BaseInterface):
 
         delta_phi = eta_r + ocp
 
-        variables.update(self._get_standard_interfacial_current_variables(j, j_av))
+        variables.update(self._get_standard_interfacial_current_variables(j))
+        variables.update(
+            self._get_standard_total_interfacial_current_variables(j_tot_av)
+        )
         variables.update(self._get_standard_exchange_current_variables(j0))
         variables.update(self._get_standard_overpotential_variables(eta_r))
         variables.update(
diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
index a20440e61..671cd42e1 100644
--- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py
+++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py
@@ -55,10 +55,13 @@ class BaseModel(BaseInterface):
             ne = self.param.ne_p
 
         j = self._get_kinetics(j0, ne, eta_r)
-        j_av = self._get_average_interfacial_current_density(variables)
-        # j = j_av + (j - pybamm.average(j))  # enforce true average
+        j_tot_av = self._get_average_total_interfacial_current_density(variables)
+        # j = j_tot_av + (j - pybamm.average(j))  # enforce true average
 
-        variables.update(self._get_standard_interfacial_current_variables(j, j_av))
+        variables.update(self._get_standard_interfacial_current_variables(j))
+        variables.update(
+            self._get_standard_total_interfacial_current_variables(j_tot_av)
+        )
         variables.update(self._get_standard_exchange_current_variables(j0))
         variables.update(self._get_standard_overpotential_variables(eta_r))
         variables.update(self._get_standard_ocp_variables(ocp, dUdT))
diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py
index cf16ae95c..909bf77b5 100644
--- a/pybamm/models/submodels/thermal/base_thermal.py
+++ b/pybamm/models/submodels/thermal/base_thermal.py
@@ -19,9 +19,10 @@ class BaseModel(pybamm.BaseSubModel):
     def __init__(self, param):
         super().__init__(param)
 
-    def _get_standard_fundamental_variables(self, T, T_av):
+    def _get_standard_fundamental_variables(self, T):
         param = self.param
         T_n, T_s, T_p = T.orphans
+        T_av = pybamm.average(T)
 
         q = self._flux_law(T)
 
@@ -131,4 +132,3 @@ class BaseModel(pybamm.BaseSubModel):
         T, _, _ = self._unpack(variables)
 
         self.initial_conditions = {T: self.param.T_init}
-
diff --git a/pybamm/models/submodels/thermal/full_thermal.py b/pybamm/models/submodels/thermal/full_thermal.py
index d7f4ca462..7563c9011 100644
--- a/pybamm/models/submodels/thermal/full_thermal.py
+++ b/pybamm/models/submodels/thermal/full_thermal.py
@@ -24,9 +24,8 @@ class Full(BaseModel):
     def get_fundamental_variables(self):
 
         T = pybamm.standard_variables.T
-        T_av = pybamm.average(T)
 
-        variables = self._get_standard_fundamental_variables(T, T_av)
+        variables = self._get_standard_fundamental_variables(T)
         return variables
 
     def get_coupled_variables(self, variables):
diff --git a/pybamm/models/submodels/thermal/isothermal.py b/pybamm/models/submodels/thermal/isothermal.py
index bc01fb96e..704567be0 100644
--- a/pybamm/models/submodels/thermal/isothermal.py
+++ b/pybamm/models/submodels/thermal/isothermal.py
@@ -23,15 +23,13 @@ class Isothermal(BaseModel):
 
     def get_fundamental_variables(self):
 
-        T_av = pybamm.Scalar(0)
-
-        T_n = pybamm.Broadcast(T_av, ["negative electrode"])
-        T_s = pybamm.Broadcast(T_av, ["separator"])
-        T_p = pybamm.Broadcast(T_av, ["positive electrode"])
+        T_n = pybamm.Broadcast(0, ["negative electrode"])
+        T_s = pybamm.Broadcast(0, ["separator"])
+        T_p = pybamm.Broadcast(0, ["positive electrode"])
 
         T = pybamm.Concatenation(T_n, T_s, T_p)
 
-        variables = self._get_standard_fundamental_variables(T, T_av)
+        variables = self._get_standard_fundamental_variables(T)
         return variables
 
     def _flux_law(self, T):
@@ -43,4 +41,3 @@ class Isothermal(BaseModel):
 
     def set_initial_conditions(self, variables):
         return {}
-
diff --git a/pybamm/models/submodels/thermal/lumped_thermal.py b/pybamm/models/submodels/thermal/lumped_thermal.py
index cfea2b99d..949b1e42a 100644
--- a/pybamm/models/submodels/thermal/lumped_thermal.py
+++ b/pybamm/models/submodels/thermal/lumped_thermal.py
@@ -29,7 +29,7 @@ class Lumped(BaseModel):
         T_p = pybamm.Broadcast(T_av, ["positive electrode"])
         T = pybamm.Concatenation(T_n, T_s, T_p)
 
-        variables = self._get_standard_fundamental_variables(T, T_av)
+        variables = self._get_standard_fundamental_variables(T)
         return variables
 
     def get_coupled_variables(self, variables):


Test patch:
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py
index 15d03f6e3..17011808d 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_leading_surface_form_stefan_maxwell_conductivity.py
@@ -19,7 +19,7 @@ class TestLeadingOrderModel(unittest.TestCase):
             "Negative electrode porosity": a_n,
             "Negative electrolyte concentration": a_n,
             "Negative electrode interfacial current density": a_n,
-            "Average negative electrode interfacial current density": a,
+            "Average negative electrode total interfacial current density": a,
         }
 
         spf = pybamm.electrolyte.stefan_maxwell.conductivity.surface_potential_form
@@ -39,7 +39,7 @@ class TestLeadingOrderModel(unittest.TestCase):
             "Positive electrode porosity": a_p,
             "Positive electrolyte concentration": a_p,
             "Positive electrode interfacial current density": a_p,
-            "Average positive electrode interfacial current density": a,
+            "Average positive electrode total interfacial current density": a,
         }
         submodel = spf.LeadingOrderAlgebraic(param, "Positive")
         std_tests = tests.StandardSubModelTests(submodel, variables)