Repository: pybamm-team/PyBaMM (python)

Issue:
Add handling of derivatives of derived variables
In SPMe we have explicit expressions for `phi_s`. We want to take gradients of `phi_s` to obtain the ohmic heating. However `phi_s` is not in the boundary conditions (and shouldn't be) so when we take the gradient the resulting symbol is of the wrong size.


Pull request description:
Issue 481 add handling derivatives of derived variables
# Description
Added derivatives of derived variables. Gradients can now be taken of phi_e_n etc for SPMe. Additionally, taking gradients of subdomain variables has been extended to the 2+1D models. Ohmic heating terms in all the thermal models have also been added. 

Still failing 1 debug mode test for 2+1 models. 

Fixes #481 

## Type of change

Please delete options that are not relevant.

- [x] New feature (non-breaking change which adds functionality)


# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [ ] Debug mode tests work. 
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks: 

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works
- [x] Any dependent changes have been merged and published in downstream modules


Patch:
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 67a3a1b4a..f90098f29 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -86,8 +86,8 @@ class Discretisation(object):
         self.set_variable_slices(variables)
 
         # set boundary conditions (only need key ids for boundary_conditions)
-        model = self.set_internal_boundary_conditions(model)
         self.bcs = self.process_boundary_conditions(model)
+        self.set_internal_boundary_conditions(model)
 
         # set up inplace vs not inplace
         if inplace:
@@ -177,23 +177,23 @@ class Discretisation(object):
         Note: this method modifies the state of self.boundary_conditions.
         """
 
-        def boundary_gradient(left_orphan, right_orphan):
+        def boundary_gradient(left_symbol, right_symbol):
 
-            left_domain = left_orphan.domain[0]
-            right_domain = right_orphan.domain[0]
+            left_domain = left_symbol.domain[0]
+            right_domain = right_symbol.domain[0]
 
-            left_submesh = self._spatial_methods[left_domain].mesh[left_domain][0]
-            right_submesh = self._spatial_methods[right_domain].mesh[right_domain][0]
+            left_mesh = self._spatial_methods[left_domain].mesh[left_domain]
+            right_mesh = self._spatial_methods[right_domain].mesh[right_domain]
 
-            x_left = left_submesh.nodes[-1]
-            x_right = right_submesh.nodes[0]
+            left_symbol_disc = self.process_symbol(left_symbol)
+            right_symbol_disc = self.process_symbol(right_symbol)
 
-            dy = pybamm.Index(right_orphan, 0) - pybamm.Index(left_orphan, -1)
-            dx = x_right - x_left
-
-            return dy / dx
+            return self._spatial_methods[left_domain].internal_neumann_condition(
+                left_symbol_disc, right_symbol_disc, left_mesh, right_mesh
+            )
 
-        bc_key_ids = [key.id for key in list(model.boundary_conditions.keys())]
+        # bc_key_ids = [key.id for key in list(model.boundary_conditions.keys())]
+        bc_key_ids = list(self.bcs.keys())
 
         internal_bcs = {}
         for var in model.boundary_conditions.keys():
@@ -205,11 +205,11 @@ class Discretisation(object):
                 next_child = children[1]
                 next_orphan = next_child.new_copy()
 
-                lbc = model.boundary_conditions[var]["left"]
+                lbc = self.bcs[var.id]["left"]
                 rbc = (boundary_gradient(first_orphan, next_orphan), "Neumann")
 
                 if first_child.id not in bc_key_ids:
-                    internal_bcs.update({first_child: {"left": lbc, "right": rbc}})
+                    internal_bcs.update({first_child.id: {"left": lbc, "right": rbc}})
 
                 for i, _ in enumerate(children[1:-1]):
                     current_child = next_child
@@ -221,16 +221,15 @@ class Discretisation(object):
                     rbc = (boundary_gradient(current_orphan, next_orphan), "Neumann")
                     if current_child.id not in bc_key_ids:
                         internal_bcs.update(
-                            {current_child: {"left": lbc, "right": rbc}}
+                            {current_child.id: {"left": lbc, "right": rbc}}
                         )
 
                 lbc = rbc
-                rbc = model.boundary_conditions[var]["right"]
+                rbc = self.bcs[var.id]["right"]
                 if children[-1].id not in bc_key_ids:
-                    internal_bcs.update({children[-1]: {"left": lbc, "right": rbc}})
+                    internal_bcs.update({children[-1].id: {"left": lbc, "right": rbc}})
 
-        model.boundary_conditions.update(internal_bcs)
-        return model
+        self.bcs.update(internal_bcs)
 
     def process_initial_conditions(self, model):
         """Discretise model initial_conditions.
@@ -275,7 +274,11 @@ class Discretisation(object):
             Dictionary of processed boundary conditions
 
         """
+
         processed_bcs = {}
+
+        # process and set pybamm.variables first incase required
+        # in discrisation of other boundary conditions
         for key, bcs in model.boundary_conditions.items():
             processed_bcs[key.id] = {}
             for side, bc in bcs.items():
@@ -286,6 +289,16 @@ class Discretisation(object):
 
         return processed_bcs
 
+    def _process_bc_entry(self, key, bcs):
+        processed_entry = {key.id: {}}
+        for side, bc in bcs.items():
+            eqn, typ = bc
+            pybamm.logger.debug("Discretise {} ({} bc)".format(key, side))
+            processed_eqn = self.process_symbol(eqn)
+            processed_entry[key.id][side] = (processed_eqn, typ)
+
+        return processed_entry
+
     def process_rhs_and_algebraic(self, model):
         """Discretise model equations - differential ('rhs') and algebraic.
 
@@ -415,6 +428,7 @@ class Discretisation(object):
             # note we are sending in the key.id here so we don't have to
             # keep calling .id
             pybamm.logger.debug("Discretise {!r}".format(eqn_key))
+
             new_var_eqn_dict[eqn_key] = self.process_symbol(eqn)
 
             new_var_eqn_dict[eqn_key].test_shape()
diff --git a/pybamm/expression_tree/simplify.py b/pybamm/expression_tree/simplify.py
index 45f5ec5ea..c9d785e64 100644
--- a/pybamm/expression_tree/simplify.py
+++ b/pybamm/expression_tree/simplify.py
@@ -215,14 +215,14 @@ def simplify_multiplication_division(myclass, left, right):
     groups of constant children (that produce a value) and simplify them
 
     The purpose of this function is to simplify expressions of the type (1 * c / 2),
-    which should simplify to (0.5 * c). The former expression consists of a Divsion,
+    which should simplify to (0.5 * c). The former expression consists of a Division,
     with a left child of a Multiplication containing a Scalar and a Parameter, and a
     right child consisting of a Scalar. For this case, this function will first flatten
     the expression to a list of the bottom level children on the numerator (i.e.
     [Scalar(1), Parameter(c)]) and their operators (i.e. [None, Multiplication]), as
     well as those children on the denominator (i.e. [Scalar(2)]. After this, all the
     constant children on the numerator and denominator (i.e. Scalar(1) and Scalar(2))
-    will be combined appropriatly, in this case to Scalar(0.5), and combined with the
+    will be combined appropriately, in this case to Scalar(0.5), and combined with the
     nonconstant children (i.e. Parameter(c))
 
     Note that this function will flatten the expression tree until a symbol is found
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index f93a57c9b..9ceb58315 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -122,7 +122,6 @@ class Mesh(dict):
                 raise pybamm.DomainError(
                     "trying to combine two meshes in different coordinate systems"
                 )
-
         submeshes = [None] * len(self[submeshnames[0]])
         for i in range(len(self[submeshnames[0]])):
             combined_submesh_edges = np.concatenate(
diff --git a/pybamm/models/full_battery_models/lead_acid/loqs.py b/pybamm/models/full_battery_models/lead_acid/loqs.py
index c488e6ca9..7ee5971c7 100644
--- a/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -91,6 +91,7 @@ class LOQS(BaseModel):
             self.submodels[
                 "electrolyte conductivity"
             ] = electrolyte.conductivity.LeadingOrder(self.param)
+
         elif self.options["surface form"] == "differential":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
@@ -98,6 +99,7 @@ class LOQS(BaseModel):
                 ] = surf_form.LeadingOrderDifferential(
                     self.param, domain, self.reactions
                 )
+
         elif self.options["surface form"] == "algebraic":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
diff --git a/pybamm/models/submodels/electrode/base_electrode.py b/pybamm/models/submodels/electrode/base_electrode.py
index 1fb3c1b67..3ad95c0a3 100644
--- a/pybamm/models/submodels/electrode/base_electrode.py
+++ b/pybamm/models/submodels/electrode/base_electrode.py
@@ -49,6 +49,7 @@ class BaseElectrode(pybamm.BaseSubModel):
             phi_s_av_dim = (
                 param.U_p_ref - param.U_n_ref + param.potential_scale * phi_s_av
             )
+
             v = pybamm.BoundaryValue(phi_s, "right")
             delta_phi_s = phi_s - pybamm.Broadcast(
                 v, ["positive electrode"], broadcast_type="primary"
@@ -71,6 +72,9 @@ class BaseElectrode(pybamm.BaseSubModel):
             "Average "
             + self.domain.lower()
             + " electrode ohmic losses [V]": delta_phi_s_av_dim,
+            "Gradient of "
+            + self.domain.lower()
+            + " electrode potential": pybamm.grad(phi_s),
         }
 
         return variables
diff --git a/pybamm/models/submodels/electrode/ohm/base_ohm.py b/pybamm/models/submodels/electrode/ohm/base_ohm.py
index f7755d550..bebee1ca4 100644
--- a/pybamm/models/submodels/electrode/ohm/base_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/base_ohm.py
@@ -23,6 +23,26 @@ class BaseModel(BaseElectrode):
     def __init__(self, param, domain, reactions=None):
         super().__init__(param, domain, reactions)
 
+    def set_boundary_conditions(self, variables):
+
+        phi_s = variables[self.domain + " electrode potential"]
+        eps = variables[self.domain + " electrode porosity"]
+        i_boundary_cc = variables["Current collector current density"]
+
+        if self.domain == "Negative":
+            lbc = (pybamm.Scalar(0), "Dirichlet")
+            rbc = (pybamm.Scalar(0), "Neumann")
+
+        elif self.domain == "Positive":
+            lbc = (pybamm.Scalar(0), "Neumann")
+            sigma_eff = self.param.sigma_p * (1 - eps) ** self.param.b
+            rbc = (
+                i_boundary_cc / pybamm.boundary_value(-sigma_eff, "right"),
+                "Neumann",
+            )
+
+        self.boundary_conditions[phi_s] = {"left": lbc, "right": rbc}
+
     @property
     def default_solver(self):
         """
diff --git a/pybamm/models/submodels/electrode/ohm/leading_ohm.py b/pybamm/models/submodels/electrode/ohm/leading_ohm.py
index ed7a61377..bb9b956c7 100644
--- a/pybamm/models/submodels/electrode/ohm/leading_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/leading_ohm.py
@@ -60,6 +60,15 @@ class LeadingOrder(BaseModel):
 
         return variables
 
+    def set_boundary_conditions(self, variables):
+
+        phi_s = variables[self.domain + " electrode potential"]
+
+        lbc = (pybamm.Scalar(0), "Neumann")
+        rbc = (pybamm.Scalar(0), "Neumann")
+
+        self.boundary_conditions[phi_s] = {"left": lbc, "right": rbc}
+
     @property
     def default_solver(self):
         """
diff --git a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py
index d8ce9927f..0bf11ffcd 100644
--- a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py
@@ -69,6 +69,10 @@ class BaseElectrolyteConductivity(pybamm.BaseSubModel):
             + pot_scale * phi_e_p_av,
             "Average electrolyte overpotential": eta_e_av,
             "Average electrolyte overpotential [V]": pot_scale * eta_e_av,
+            "Gradient of negative electrolyte potential": pybamm.grad(phi_e_n),
+            "Gradient of separator electrolyte potential": pybamm.grad(phi_e_s),
+            "Gradient of positive electrolyte potential": pybamm.grad(phi_e_p),
+            "Gradient of electrolyte potential": pybamm.grad(phi_e),
         }
 
         return variables
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py
index 198d6bf1c..1da015c25 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_stefan_maxwell_conductivity.py
@@ -2,6 +2,7 @@
 # Base class for electrolyte conductivity employing stefan-maxwell
 #
 from ...base_electrolyte_conductivity import BaseElectrolyteConductivity
+import pybamm
 
 
 class BaseModel(BaseElectrolyteConductivity):
@@ -19,3 +20,12 @@ class BaseModel(BaseElectrolyteConductivity):
 
     def __init__(self, param, domain=None, reactions=None):
         super().__init__(param, domain, reactions)
+
+    def set_boundary_conditions(self, variables):
+        phi_e = variables["Electrolyte potential"]
+        self.boundary_conditions = {
+            phi_e: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
index d4d8c34eb..3dcd749da 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
@@ -58,16 +58,6 @@ class Full(BaseModel):
 
         self.algebraic = {phi_e: pybamm.div(i_e) - sum_j}
 
-    def set_boundary_conditions(self, variables):
-        phi_e = variables["Electrolyte potential"]
-
-        self.boundary_conditions = {
-            phi_e: {
-                "left": (pybamm.Scalar(0), "Neumann"),
-                "right": (pybamm.Scalar(0), "Neumann"),
-            }
-        }
-
     def set_initial_conditions(self, variables):
         phi_e = variables["Electrolyte potential"]
         self.initial_conditions = {phi_e: -self.param.U_n(self.param.c_n_init)}
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
index a19613c07..a6b927fe2 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
@@ -106,6 +106,17 @@ class BaseModel(BaseStefanMaxwellConductivity):
             c_e: {"left": lbc_c_e, "right": rbc_c_e},
         }
 
+        if self.domain == "Negative":
+            phi_e = variables["Electrolyte potential"]
+            self.boundary_conditions.update(
+                {
+                    phi_e: {
+                        "left": (pybamm.Scalar(0), "Neumann"),
+                        "right": (pybamm.Scalar(0), "Neumann"),
+                    }
+                }
+            )
+
     def _get_conductivities(self, variables):
         param = self.param
         eps = variables[self.domain + " electrode porosity"]
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 c1260cddd..b257cc5d7 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
@@ -2,6 +2,7 @@
 # Class for leading-order surface form electrolyte conductivity employing stefan-maxwell
 #
 import pybamm
+
 from ..leading_stefan_maxwell_conductivity import LeadingOrder
 
 
@@ -62,6 +63,16 @@ class BaseLeadingOrderSurfaceForm(LeadingOrder):
         else:
             return variables
 
+    def set_boundary_conditions(self, variables):
+        if self.domain == "Negative":
+            phi_e = variables["Electrolyte potential"]
+            self.boundary_conditions = {
+                phi_e: {
+                    "left": (pybamm.Scalar(0), "Neumann"),
+                    "right": (pybamm.Scalar(0), "Neumann"),
+                }
+            }
+
 
 class LeadingOrderDifferential(BaseLeadingOrderSurfaceForm):
     """Leading-order model for conservation of charge in the electrolyte employing the
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py
index 0bb8ec46e..2fd5146bc 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/base_stefan_maxwell_diffusion.py
@@ -2,6 +2,7 @@
 # Base class for electrolyte diffusion employing stefan-maxwell
 #
 from ...base_electrolyte_diffusion import BaseElectrolyteDiffusion
+import pybamm
 
 
 class BaseModel(BaseElectrolyteDiffusion):
@@ -19,3 +20,14 @@ class BaseModel(BaseElectrolyteDiffusion):
 
     def __init__(self, param, reactions=None):
         super().__init__(param, reactions)
+
+    def set_boundary_conditions(self, variables):
+
+        c_e = variables["Electrolyte concentration"]
+
+        self.boundary_conditions = {
+            c_e: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
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 31eddfef8..70d62c418 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
@@ -77,17 +77,6 @@ class Full(BaseModel):
             * (-pybamm.div(N_e) / param.C_e + source_terms - c_e * deps_dt)
         }
 
-    def set_boundary_conditions(self, variables):
-
-        c_e = variables["Electrolyte concentration"]
-
-        self.boundary_conditions = {
-            c_e: {
-                "left": (pybamm.Scalar(0), "Neumann"),
-                "right": (pybamm.Scalar(0), "Neumann"),
-            }
-        }
-
     def set_initial_conditions(self, variables):
 
         c_e = variables["Electrolyte concentration"]
diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py
index 909bf77b5..bc04f42a5 100644
--- a/pybamm/models/submodels/thermal/base_thermal.py
+++ b/pybamm/models/submodels/thermal/base_thermal.py
@@ -59,18 +59,22 @@ class BaseModel(pybamm.BaseSubModel):
         dUdT_n = variables["Negative electrode entropic change"]
         dUdT_p = variables["Positive electrode entropic change"]
 
-        # i_e = variables["Electrolyte current density"]
-        # phi_e = variables["Electrolyte potential"]
-
-        # phi_s_n = variables["Negative electrode potential"]
-        # phi_s_s = pybamm.Broadcast(0, ["separator"])
-        # phi_s_p = variables["Positive electrode potential"]
-        # phi_s = pybamm.Concatenation(phi_s_n, phi_s_s, phi_s_p)
-
-        # TODO: add ohmic heating from solid and electrolyte
-        # Q_ohm = -i_s * pybamm.grad(phi_s) - i_e * pybamm.grad(phi_e)
-        # Q_ohm = -pybamm.inner(i_e, pybamm.grad(phi_e))
-        Q_ohm = pybamm.Scalar(0)
+        i_e = variables["Electrolyte current density"]
+        phi_e = variables["Electrolyte potential"]
+
+        i_s_n = variables["Negative electrode current density"]
+        i_s_p = variables["Positive electrode current density"]
+        phi_s_n = variables["Negative electrode potential"]
+        phi_s_p = variables["Positive electrode potential"]
+
+        Q_ohm_s_n = -pybamm.inner(i_s_n, pybamm.grad(phi_s_n))
+        Q_ohm_s_s = pybamm.Broadcast(0, ["separator"])
+        Q_ohm_s_p = -pybamm.inner(i_s_p, pybamm.grad(phi_s_p))
+        Q_ohm_s = pybamm.Concatenation(Q_ohm_s_n, Q_ohm_s_s, Q_ohm_s_p)
+
+        Q_ohm_e = -pybamm.inner(i_e, pybamm.grad(phi_e))
+
+        Q_ohm = Q_ohm_s + Q_ohm_e
 
         Q_rxn_n = j_n * eta_r_n
         Q_rxn_p = j_p * eta_r_p
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index fa6f0cda3..801427b6c 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -133,6 +133,7 @@ class OdeSolver(pybamm.BaseSolver):
         if model.use_jacobian:
             # Create Jacobian from simplified rhs
             y = pybamm.StateVector(slice(0, np.size(y0)))
+
             pybamm.logger.info("Calculating jacobian")
             jac_rhs = concatenated_rhs.jac(y)
             if model.use_simplify:
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index c9c4518bd..b9e3c71da 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -69,7 +69,9 @@ class FiniteVolume(pybamm.SpatialMethod):
         if symbol.id in boundary_conditions:
             bcs = boundary_conditions[symbol.id]
             # add ghost nodes
-            discretised_symbol = self.add_ghost_nodes(symbol, discretised_symbol, bcs)
+            discretised_symbol = self.add_ghost_nodes(
+                symbol, discretised_symbol, bcs
+            )
             # edit domain
             domain = (
                 [domain[0] + "_left ghost cell"]
@@ -321,6 +323,52 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return pybamm.Matrix(matrix)
 
+    def internal_neumann_condition(
+        self, left_symbol_disc, right_symbol_disc, left_mesh, right_mesh
+    ):
+        """
+        A method to find the internal neumann conditions between two symbols
+        on adjacent subdomains.
+
+        Parameters
+        ----------
+        left_symbol_disc : :class:`pybamm.Symbol`
+            The discretised symbol on the left subdomain
+        right_symbol_disc : :class:`pybamm.Symbol`
+            The discretised symbol on the right subdomain
+        left_mesh : list
+            The mesh on the left subdomain
+        right_mesh : list
+            The mesh on the right subdomain
+        """
+
+        left_npts = left_mesh[0].npts
+        right_npts = right_mesh[0].npts
+
+        sec_pts = len(left_mesh)
+
+        if sec_pts != len(right_mesh):
+            raise pybamm.DomainError(
+                """Number of secondary points in subdomains do not match"""
+            )
+
+        left_sub_matrix = np.zeros((1, left_npts))
+        left_sub_matrix[0][left_npts - 1] = 1
+        left_matrix = pybamm.Matrix(csr_matrix(kron(eye(sec_pts), left_sub_matrix)))
+
+        right_sub_matrix = np.zeros((1, right_npts))
+        right_sub_matrix[0][0] = 1
+        right_matrix = pybamm.Matrix(csr_matrix(kron(eye(sec_pts), right_sub_matrix)))
+
+        right_copy = right_symbol_disc.new_copy()
+        left_copy = left_symbol_disc.new_copy()
+        right_copy.domain = []
+        left_copy.domain = []
+        dy = right_matrix @ right_copy - left_matrix @ left_copy
+        dx = right_mesh[0].nodes[0] - left_mesh[0].nodes[-1]
+
+        return dy / dx
+
     def indefinite_integral_matrix_nodes(self, domain):
         """
         Matrix for finite-volume implementation of the indefinite integral where the
@@ -541,8 +589,11 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         # inner product takes fluxes from edges to nodes
         if isinstance(bin_op, pybamm.Inner):
-            disc_left = self.edge_to_node(disc_left)
-            disc_right = self.edge_to_node(disc_right)
+            if left_evaluates_on_edges:
+                disc_left = self.edge_to_node(disc_left)
+            if right_evaluates_on_edges:
+                disc_right = self.edge_to_node(disc_right)
+
         # If neither child evaluates on edges, or both children have gradients,
         # no need to do any averaging
         elif left_evaluates_on_edges == right_evaluates_on_edges:
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 7ddf64e5b..4df976253 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -210,6 +210,27 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
+    def internal_neumann_condition(
+        self, left_symbol_disc, right_symbol_disc, left_mesh, right_mesh
+    ):
+        """
+        A method to find the internal neumann conditions between two symbols
+        on adjacent subdomains.
+
+        Parameters
+        ----------
+        left_symbol_disc : :class:`pybamm.Symbol`
+            The discretised symbol on the left subdomain
+        right_symbol_disc : :class:`pybamm.Symbol`
+            The discretised symbol on the right subdomain
+        left_mesh : list
+            The mesh on the left subdomain
+        right_mesh : list
+            The mesh on the right subdomain
+        """
+
+        raise NotImplementedError
+
     def boundary_value_or_flux(self, symbol, discretised_child):
         """
         Returns the boundary value or flux using the approriate expression for the


Test patch:
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 62b329f54..bf6303a97 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -448,6 +448,11 @@ class PotentialTests(BaseOutputTest):
             "Positive electrode surface potential difference [V]"
         ]
 
+        self.grad_phi_e = variables["Gradient of electrolyte potential"]
+        self.grad_phi_e_n = variables["Gradient of negative electrolyte potential"]
+        self.grad_phi_e_s = variables["Gradient of separator electrolyte potential"]
+        self.grad_phi_e_p = variables["Gradient of positive electrolyte potential"]
+
     def test_negative_electrode_potential_profile(self):
         """Test that negative electrode potential is zero on left boundary. Test
         average negative electrode potential is less than or equal to zero."""
@@ -482,6 +487,20 @@ class PotentialTests(BaseOutputTest):
 
         np.testing.assert_array_less(-self.phi_s_p(self.t, self.x_p), 0)
 
+    def test_gradient_splitting(self):
+
+        t, x_n, x_s, x_p, x = self.t, self.x_n, self.x_s, self.x_p, self.x
+        grad_phi_e_combined = np.concatenate(
+            (
+                self.grad_phi_e_n(t, x_n),
+                self.grad_phi_e_s(t, x_s),
+                self.grad_phi_e_p(t, x_p),
+            ),
+            axis=0,
+        )
+
+        np.testing.assert_array_equal(self.grad_phi_e(t, x), grad_phi_e_combined)
+
     def test_all(self):
         self.test_negative_electrode_potential_profile()
         self.test_positive_electrode_potential_profile()
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index fcd9218ac..736b991c1 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -50,75 +50,6 @@ class TestDFN(unittest.TestCase):
         modeltest = tests.StandardModelTest(model, var_pts=var_pts)
         modeltest.test_all()
 
-    def test_gradient_of_subdomain_var(self):
-        # a hacky test that will be replaced by a more general one in
-        # standard output tests when #481 is completed
-
-        options = {"thermal": None}
-        model = pybamm.lithium_ion.DFN(options)
-        c_e = model.variables["Electrolyte concentration"]
-        c_e_n, c_e_s, c_e_p = c_e.orphans
-        phi_e = model.variables["Electrolyte potential"]
-        phi_e_n, phi_e_s, phi_e_p = phi_e.orphans
-        model.variables.update(
-            {
-                "Test grad c_e": pybamm.grad(c_e),
-                "Test grad c_e_n": pybamm.grad(c_e_n),
-                "Test grad c_e_s": pybamm.grad(c_e_s),
-                "Test grad c_e_p": pybamm.grad(c_e_p),
-                "Test grad phi_e": pybamm.grad(phi_e),
-                "Test grad phi_e_n": pybamm.grad(phi_e_n),
-                "Test grad phi_e_s": pybamm.grad(phi_e_s),
-                "Test grad phi_e_p": pybamm.grad(phi_e_p),
-            }
-        )
-        geometry = model.default_geometry
-        param = model.default_parameter_values
-        param.process_model(model)
-        param.process_geometry(geometry)
-        mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
-        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-        disc.process_model(model)
-
-        t_eval = np.linspace(0, 2, 100)
-        solution = model.default_solver.solve(model, t_eval)
-        t = solution.t
-        x_n = disc.mesh["negative electrode"][0].nodes
-        x_s = disc.mesh["separator"][0].nodes
-        x_p = disc.mesh["positive electrode"][0].nodes
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        x = disc.mesh.combine_submeshes(*whole_cell)[0].nodes
-
-        variables = pybamm.post_process_variables(
-            model.variables, solution.t, solution.y, mesh=disc.mesh
-        )
-
-        grad_c_e_n = variables["Test grad c_e_n"]
-        grad_c_e_s = variables["Test grad c_e_s"]
-        grad_c_e_p = variables["Test grad c_e_p"]
-        grad_c_e = variables["Test grad c_e"]
-
-        grad_c_e_combined = np.concatenate(
-            (grad_c_e_n(t, x_n), grad_c_e_s(t, x_s), grad_c_e_p(t, x_p)), axis=0
-        )
-
-        np.testing.assert_array_almost_equal(
-            grad_c_e(t, x), grad_c_e_combined, decimal=6
-        )
-
-        grad_phi_e_n = variables["Test grad phi_e_n"]
-        grad_phi_e_s = variables["Test grad phi_e_s"]
-        grad_phi_e_p = variables["Test grad phi_e_p"]
-        grad_phi_e = variables["Test grad phi_e"]
-
-        grad_phi_e_combined = np.concatenate(
-            (grad_phi_e_n(t, x_n), grad_phi_e_s(t, x_s), grad_phi_e_p(t, x_p)), axis=0
-        )
-
-        np.testing.assert_array_almost_equal(
-            grad_phi_e(t, x), grad_phi_e_combined, decimal=6
-        )
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/shared.py b/tests/shared.py
index 7d06e002c..acb499db7 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -25,6 +25,11 @@ class SpatialMethodForTesting(pybamm.SpatialMethod):
         divergence_matrix = pybamm.Matrix(eye(n))
         return divergence_matrix @ discretised_symbol
 
+    def internal_neumann_condition(
+        self, left_symbol_disc, right_symbol_disc, left_mesh, right_mesh
+    ):
+        return pybamm.Scalar(0)
+
     def mass_matrix(self, symbol, boundary_conditions):
         n = 0
         for domain in symbol.domain:
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index aa90603d4..ca9a397f5 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -10,6 +10,8 @@ from tests import (
     get_discretisation_for_testing,
     get_1p1d_discretisation_for_testing,
 )
+from tests.shared import SpatialMethodForTesting
+
 from scipy.sparse import block_diag
 
 
@@ -46,18 +48,23 @@ class TestDiscretise(unittest.TestCase):
 
     def test_add_internal_boundary_conditions(self):
         model = pybamm.BaseModel()
-        c_e = pybamm.standard_variables.c_e
+        c_e_n = pybamm.Broadcast(0, ["negative electrode"])
+        c_e_s = pybamm.Broadcast(0, ["separator"])
+        c_e_p = pybamm.Broadcast(0, ["positive electrode"])
+        c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
         lbc = (pybamm.Scalar(0), "Neumann")
         rbc = (pybamm.Scalar(0), "Neumann")
         model.boundary_conditions = {c_e: {"left": lbc, "right": rbc}}
 
         mesh = get_mesh_for_testing()
-        spatial_methods = {"macroscale": pybamm.SpatialMethod}
+        spatial_methods = {"macroscale": SpatialMethodForTesting}
+
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        model = disc.set_internal_boundary_conditions(model)
+        disc.bcs = disc.process_boundary_conditions(model)
+        disc.set_internal_boundary_conditions(model)
 
         for child in c_e.children:
-            self.assertTrue(child in model.boundary_conditions.keys())
+            self.assertTrue(child.id in disc.bcs.keys())
 
     def test_discretise_slicing(self):
         # create discretisation
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_base_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_base_stefan_maxwell_conductivity.py
index b6110f048..fea3f3acf 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_base_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_base_stefan_maxwell_conductivity.py
@@ -11,7 +11,9 @@ class TestBaseModel(unittest.TestCase):
     def test_public_functions(self):
         submodel = pybamm.electrolyte.stefan_maxwell.conductivity.BaseModel(None)
         std_tests = tests.StandardSubModelTests(submodel)
-        std_tests.test_all()
+
+        with self.assertRaises(KeyError):
+            std_tests.test_all()
 
 
 if __name__ == "__main__":
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
index 1d6f0b613..ab63f0662 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
@@ -19,6 +19,7 @@ class TestFull(unittest.TestCase):
             "Negative electrode porosity": a_n,
             "Negative electrolyte concentration": a_n,
             "Negative electrode interfacial current density": a_n,
+            "Electrolyte potential": pybamm.Concatenation(a_n, a_s, a_p),
         }
         icd = " interfacial current density"
         reactions = {
diff --git a/tests/unit/test_models/test_submodels/test_thermal/test_full_thermal.py b/tests/unit/test_models/test_submodels/test_thermal/test_full_thermal.py
index dc1cc06c2..9aafd527f 100644
--- a/tests/unit/test_models/test_submodels/test_thermal/test_full_thermal.py
+++ b/tests/unit/test_models/test_submodels/test_thermal/test_full_thermal.py
@@ -11,6 +11,7 @@ class TestFull(unittest.TestCase):
     def test_public_functions(self):
         param = pybamm.standard_parameters_lithium_ion
         a_n = pybamm.Broadcast(pybamm.Scalar(0), ["negative electrode"])
+        a_s = pybamm.Broadcast(pybamm.Scalar(0), ["separator"])
         a_p = pybamm.Broadcast(pybamm.Scalar(0), ["positive electrode"])
         variables = {
             "Negative electrode interfacial current density": a_n,
@@ -19,6 +20,12 @@ class TestFull(unittest.TestCase):
             "Positive electrode reaction overpotential": a_p,
             "Negative electrode entropic change": a_n,
             "Positive electrode entropic change": a_p,
+            "Electrolyte potential": pybamm.Concatenation(a_n, a_s, a_p),
+            "Electrolyte current density": pybamm.Concatenation(a_n, a_s, a_p),
+            "Negative electrode potential": a_n,
+            "Negative electrode current density": a_n,
+            "Positive electrode potential": a_p,
+            "Positive electrode current density": a_p,
         }
         submodel = pybamm.thermal.Full(param)
         std_tests = tests.StandardSubModelTests(submodel, variables)
diff --git a/tests/unit/test_models/test_submodels/test_thermal/test_lumped_thermal.py b/tests/unit/test_models/test_submodels/test_thermal/test_lumped_thermal.py
index 7c75b3ebb..d01e86ea7 100644
--- a/tests/unit/test_models/test_submodels/test_thermal/test_lumped_thermal.py
+++ b/tests/unit/test_models/test_submodels/test_thermal/test_lumped_thermal.py
@@ -11,6 +11,7 @@ class TestLumped(unittest.TestCase):
     def test_public_functions(self):
         param = pybamm.standard_parameters_lithium_ion
         a_n = pybamm.Broadcast(pybamm.Scalar(0), ["negative electrode"])
+        a_s = pybamm.Broadcast(pybamm.Scalar(0), ["separator"])
         a_p = pybamm.Broadcast(pybamm.Scalar(0), ["positive electrode"])
         variables = {
             "Negative electrode interfacial current density": a_n,
@@ -19,6 +20,12 @@ class TestLumped(unittest.TestCase):
             "Positive electrode reaction overpotential": a_p,
             "Negative electrode entropic change": a_n,
             "Positive electrode entropic change": a_p,
+            "Electrolyte potential": pybamm.Concatenation(a_n, a_s, a_p),
+            "Electrolyte current density": pybamm.Concatenation(a_n, a_s, a_p),
+            "Negative electrode potential": a_n,
+            "Negative electrode current density": a_n,
+            "Positive electrode potential": a_p,
+            "Positive electrode current density": a_p,
         }
         submodel = pybamm.thermal.Lumped(param)
         std_tests = tests.StandardSubModelTests(submodel, variables)