Repository: pybamm-team/PyBaMM

Issue:
Lead-acid convection
Add models for convection in lead-acid batteries

Discussion:

Patch:
diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py
index ea26a0b9b..e15a3b255 100644
--- a/examples/scripts/compare_lead_acid.py
+++ b/examples/scripts/compare_lead_acid.py
@@ -5,9 +5,9 @@ pybamm.set_logging_level("INFO")
 
 # load models
 models = [
-    pybamm.lead_acid.LOQS(),
-    pybamm.lead_acid.Composite(),
-    pybamm.lead_acid.NewmanTiedemann(),
+    pybamm.lead_acid.LOQS({"convection": True}),
+    pybamm.lead_acid.Composite({"convection": True}),
+    pybamm.lead_acid.NewmanTiedemann({"convection": True}),
 ]
 
 # create geometry
@@ -44,7 +44,7 @@ output_variables = [
     "Electrolyte concentration [mol.m-3]",
     "Current [A]",
     "Porosity",
-    "Electrolyte potential [V]",
+    "Volume-averaged velocity [m.s-1]",
     "Terminal voltage [V]",
 ]
 plot = pybamm.QuickPlot(models, mesh, solvers, output_variables)
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index a6131e9c3..1b24ff382 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -145,6 +145,7 @@ from .models.submodels import (
     particle,
     porosity,
     potential,
+    velocity,
     vertical,
 )
 
@@ -170,9 +171,8 @@ from .geometry.geometry import (
     Geometry1p1DMicro,
 )
 
-from .expression_tree.independent_variable import KNOWN_SPATIAL_VARS
+from .expression_tree.independent_variable import KNOWN_SPATIAL_VARS, KNOWN_COORD_SYS
 from .geometry import standard_spatial_vars
-from .geometry.standard_spatial_vars import KNOWN_COORD_SYS
 
 #
 # Mesh and Discretisation classes
diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index cec22982a..348370e8a 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -3,7 +3,10 @@
 #
 import pybamm
 
+KNOWN_COORD_SYS = ["cartesian", "spherical polar"]
 KNOWN_SPATIAL_VARS = ["x", "y", "z", "r", "x_n", "x_s", "x_p", "r_n", "r_p"]
+KNOWN_SPATIAL_VARS_EXTENDED = [v + "_edge" for v in KNOWN_SPATIAL_VARS]
+KNOWN_SPATIAL_VARS.extend(KNOWN_SPATIAL_VARS_EXTENDED)
 
 
 class IndependentVariable(pybamm.Symbol):
diff --git a/pybamm/geometry/standard_spatial_vars.py b/pybamm/geometry/standard_spatial_vars.py
index d23dd757a..231d515dd 100644
--- a/pybamm/geometry/standard_spatial_vars.py
+++ b/pybamm/geometry/standard_spatial_vars.py
@@ -1,10 +1,8 @@
 import pybamm
 
-KNOWN_COORD_SYS = ["cartesian", "spherical polar"]
-KNOWN_SPATIAL_VARS = ["x", "y", "z", "r", "x_n", "x_s", "x_p", "r_n", "r_p"]
-
 whole_cell = ["negative electrode", "separator", "positive electrode"]
 
+# Domains at cell centres
 x_n = pybamm.SpatialVariable(
     "x_n", domain=["negative electrode"], coord_sys="cartesian"
 )
@@ -23,3 +21,29 @@ r_n = pybamm.SpatialVariable(
 r_p = pybamm.SpatialVariable(
     "r_p", domain=["positive particle"], coord_sys="spherical polar"
 )
+
+# Domains at cell edges
+x_n_edge = pybamm.SpatialVariable(
+    "x_n_edge", domain=["negative electrode"], coord_sys="cartesian"
+)
+x_s_edge = pybamm.SpatialVariable(
+    "x_s_edge", domain=["separator"], coord_sys="cartesian"
+)
+x_p_edge = pybamm.SpatialVariable(
+    "x_p_edge", domain=["positive electrode"], 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"
+)
+z_edge = pybamm.SpatialVariable(
+    "z_edge", domain="current collector", coord_sys="cartesian"
+)
+
+r_n_edge = pybamm.SpatialVariable(
+    "r_n_edge", domain=["negative particle"], coord_sys="spherical polar"
+)
+r_p_edge = pybamm.SpatialVariable(
+    "r_p_edge", domain=["positive particle"], coord_sys="spherical polar"
+)
diff --git a/pybamm/models/base_models.py b/pybamm/models/base_models.py
index 5556b9a02..67e7030e4 100644
--- a/pybamm/models/base_models.py
+++ b/pybamm/models/base_models.py
@@ -463,7 +463,11 @@ class StandardBatteryBaseModel(BaseModel):
 
     @property
     def options(self):
-        default_options = {"capacitance": False, "bc_options": {"dimensionality": 0}}
+        default_options = {
+            "bc_options": {"dimensionality": 0},
+            "capacitance": False,
+            "convection": False,
+        }
         if self._extra_options is None:
             options = default_options
         else:
diff --git a/pybamm/models/lead_acid/composite.py b/pybamm/models/lead_acid/composite.py
index fcddbaf1c..c821ca00e 100644
--- a/pybamm/models/lead_acid/composite.py
+++ b/pybamm/models/lead_acid/composite.py
@@ -36,6 +36,7 @@ class Composite(pybamm.LeadAcidBaseModel):
         self.set_diffusion_submodel()
         self.set_electrolyte_current_model(int_curr_model)
         self.set_current_variables()
+        self.set_convection_variables()
 
     def set_boundary_conditions(self, bc_variables):
         """Set boundary conditions, dependent on self.options"""
@@ -163,6 +164,17 @@ class Composite(pybamm.LeadAcidBaseModel):
         voltage = self.variables["Terminal voltage"]
         self.events.append(voltage - param.voltage_low_cut)
 
+    def set_convection_variables(self):
+        velocity_model = pybamm.velocity.Velocity(self.set_of_parameters)
+        if self.options["convection"] is not False:
+            velocity_vars = velocity_model.get_explicit_composite(self.variables)
+            self.variables.update(velocity_vars)
+        else:
+            whole_cell = ["negative electrode", "separator", "positive electrode"]
+            v_box = pybamm.Broadcast(0, whole_cell)
+            dVbox_dz = pybamm.Broadcast(0, whole_cell)
+            self.variables.update(velocity_model.get_variables(v_box, dVbox_dz))
+
     @property
     def default_solver(self):
         """
diff --git a/pybamm/models/lead_acid/loqs.py b/pybamm/models/lead_acid/loqs.py
index a6f72f8da..5b64a968f 100644
--- a/pybamm/models/lead_acid/loqs.py
+++ b/pybamm/models/lead_acid/loqs.py
@@ -25,6 +25,7 @@ class LOQS(pybamm.LeadAcidBaseModel):
         self.set_porosity_submodel()
         self.set_diffusion_submodel()
         self.set_current_variables()
+        self.set_convection_variables()
 
         # ODEs only (don't use jacobian, use base spatial method)
         self.use_jacobian = False
@@ -156,7 +157,6 @@ class LOQS(pybamm.LeadAcidBaseModel):
         self.variables.update({**ocp_vars, **eta_r_vars})
 
     def set_porosity_submodel(self):
-        " Set model for porosity"
         param = self.set_of_parameters
         porosity_model = pybamm.porosity.Standard(param)
         porosity_model.set_leading_order_system(self.variables)
@@ -197,6 +197,17 @@ class LOQS(pybamm.LeadAcidBaseModel):
             voltage = pybamm.boundary_value(voltage, "right")
         self.events.append(voltage - param.voltage_low_cut)
 
+    def set_convection_variables(self):
+        velocity_model = pybamm.velocity.Velocity(self.set_of_parameters)
+        if self.options["convection"] is not False:
+            velocity_vars = velocity_model.get_explicit_leading_order(self.variables)
+            self.variables.update(velocity_vars)
+        else:
+            whole_cell = ["negative electrode", "separator", "positive electrode"]
+            v_box = pybamm.Broadcast(0, whole_cell)
+            dVbox_dz = pybamm.Broadcast(0, whole_cell)
+            self.variables.update(velocity_model.get_variables(v_box, dVbox_dz))
+
     @property
     def default_spatial_methods(self):
         # ODEs only in the macroscale, so use base spatial method
diff --git a/pybamm/models/lead_acid/newman_tiedemann.py b/pybamm/models/lead_acid/newman_tiedemann.py
index 65eb6e5cf..9b606404d 100644
--- a/pybamm/models/lead_acid/newman_tiedemann.py
+++ b/pybamm/models/lead_acid/newman_tiedemann.py
@@ -24,6 +24,7 @@ class NewmanTiedemann(pybamm.LeadAcidBaseModel):
         self.set_boundary_conditions(None)
         self.set_interfacial_variables()
         self.set_porosity_submodel()
+        self.set_convection_submodel()
         self.set_diffusion_submodel()
         self.set_current_submodels()
 
@@ -117,6 +118,18 @@ class NewmanTiedemann(pybamm.LeadAcidBaseModel):
             "Porosity change"
         ]
 
+    def set_convection_submodel(self):
+        velocity_model = pybamm.velocity.Velocity(self.set_of_parameters)
+        if self.options["convection"] is not False:
+            self.variables["Electrolyte pressure"] = pybamm.standard_variables.pressure
+            velocity_model.set_algebraic_system(self.variables)
+            self.update(velocity_model)
+        else:
+            whole_cell = ["negative electrode", "separator", "positive electrode"]
+            v_box = pybamm.Broadcast(0, whole_cell)
+            dVbox_dz = pybamm.Broadcast(0, whole_cell)
+            self.variables.update(velocity_model.get_variables(v_box, dVbox_dz))
+
     def set_diffusion_submodel(self):
         param = self.set_of_parameters
         reactions = self.reactions
diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py
index f62e3abe9..7360d98a1 100644
--- a/pybamm/models/standard_variables.py
+++ b/pybamm/models/standard_variables.py
@@ -46,3 +46,9 @@ eps_piecewise_constant = pybamm.Concatenation(
     pybamm.Broadcast(eps_s_pc, ["separator"]),
     pybamm.Broadcast(eps_p_pc, ["positive electrode"]),
 )
+
+# 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 = pybamm.Concatenation(pressure_n, pressure_s, pressure_p)
diff --git a/pybamm/models/submodels/electrolyte_diffusion.py b/pybamm/models/submodels/electrolyte_diffusion.py
index 2c92dd1f3..21dfd059d 100644
--- a/pybamm/models/submodels/electrolyte_diffusion.py
+++ b/pybamm/models/submodels/electrolyte_diffusion.py
@@ -43,10 +43,9 @@ class StefanMaxwell(pybamm.SubModel):
         except KeyError:
             epsilon = param.epsilon
             deps_dt = pybamm.Scalar(0)
+
         # Use convection velocity if it exists, otherwise set to zero
         v_box = variables.get("Volume-averaged velocity", pybamm.Scalar(0))
-        # # Use variable for div_v_box, to avoid having an extra divergence in rhs
-        # div_v_box = variables.get("div(volume-averaged velocity)", pybamm.Scalar(0))
 
         # Flux
         N_e_diff = -(epsilon ** param.b) * param.D_e(c_e) * pybamm.grad(c_e)
diff --git a/pybamm/models/submodels/velocity.py b/pybamm/models/submodels/velocity.py
new file mode 100644
index 000000000..50a234d18
--- /dev/null
+++ b/pybamm/models/submodels/velocity.py
@@ -0,0 +1,168 @@
+#
+# Equation classes for the electrolyte velocity
+#
+import pybamm
+
+
+class Velocity(pybamm.SubModel):
+    """Electrolyte velocity
+
+    Parameters
+    ----------
+    set_of_parameters : parameter class
+        The parameters to use for this submodel
+
+    *Extends:* :class:`pybamm.SubModel`
+    """
+
+    def __init__(self, set_of_parameters):
+        super().__init__(set_of_parameters)
+
+    def set_algebraic_system(self, variables):
+        """
+        Algebraic system for pressure and volume-averaged velocity in the electrolyte.
+
+        Parameters
+        ----------
+        variables : dict
+            Dictionary of symbols to use in the model
+        """
+        # Unpack variables
+        param = self.set_of_parameters
+        p = variables["Electrolyte pressure"]
+        # c_e = variables["Electrolyte concentration"]
+
+        # Set up reactions
+        j = variables["Interfacial current density"]
+        v_mass = -pybamm.grad(p)
+        v_box = v_mass
+        _, dVbox_dz = self.get_separator_velocities(variables)
+
+        # Build model
+        self.algebraic = {p: pybamm.div(v_box) + dVbox_dz - param.beta * j}
+        self.initial_conditions = {p: 0}
+        self.boundary_conditions = {
+            p: {"left": (0, "Dirichlet"), "right": (0, "Neumann")}
+        }
+        self.variables = self.get_variables(v_box, dVbox_dz)
+
+    def get_explicit_leading_order(self, variables):
+        """
+        Provides explicit velocity for the leading-order models, as a post-processing
+        step.
+
+        Parameters
+        ----------
+        variables : dict
+            Dictionary of symbols to use in the model
+        """
+        # Set up
+        param = self.set_of_parameters
+        x_n = pybamm.standard_spatial_vars.x_n
+        x_p = pybamm.standard_spatial_vars.x_p
+
+        j_n = variables["Negative electrode interfacial current density"]
+        j_p = variables["Positive electrode interfacial current density"]
+
+        # Volume-averaged velocity
+        v_box_n = param.beta_n * pybamm.outer(j_n, x_n)
+        v_box_p = param.beta_p * pybamm.outer(j_p, x_p - 1)
+
+        v_box_s, dVbox_dz = self.get_separator_velocities(variables)
+        v_box = pybamm.Concatenation(v_box_n, v_box_s, v_box_p)
+
+        return self.get_variables(v_box, dVbox_dz)
+
+    def get_explicit_composite(self, variables):
+        """
+        Provides explicit velocity for the composite models, as a post-processing step.
+
+        Parameters
+        ----------
+        variables : dict
+            Dictionary of symbols to use in the model
+        """
+        # Set up
+        param = self.set_of_parameters
+        x_n = pybamm.standard_spatial_vars.x_n
+        x_p = pybamm.standard_spatial_vars.x_p
+        j_n = variables["Negative electrode interfacial current density"]
+        j_p = variables["Positive electrode interfacial current density"]
+
+        # Volume-averaged velocity
+        v_box_n = param.beta_n * pybamm.IndefiniteIntegral(j_n, x_n)
+        # Shift v_box_p to be equal to 0 at x_p = 1
+        v_box_p = param.beta_p * (
+            pybamm.IndefiniteIntegral(j_p, x_p) - pybamm.Integral(j_p, x_p)
+        )
+
+        v_box_s, dVbox_dz = self.get_separator_velocities(variables)
+        v_box = pybamm.Concatenation(v_box_n, v_box_s, v_box_p)
+
+        return self.get_variables(v_box, dVbox_dz)
+
+    def get_separator_velocities(self, variables):
+        """
+        Calculate x- and z-components of velocity in the separator
+
+        Parameters
+        ----------
+        variables : dict
+            Dictionary of symbols to use in the model
+
+        Returns
+        -------
+        v_box_s : :class:`pybamm.Symbol`
+            The x-component of velocity in the separator
+        dVbox_dz : :class:`pybamm.Symbol`
+            The z-component of velocity in the separator
+        """
+        # Set up
+        param = self.set_of_parameters
+        l_n = pybamm.geometric_parameters.l_n
+        l_s = pybamm.geometric_parameters.l_s
+        x_s = pybamm.standard_spatial_vars.x_s
+
+        # Difference in negative and positive electrode velocities determines the
+        # velocity in the separator
+        i_boundary_cc = variables["Current collector current density"]
+        v_box_n_right = param.beta_n * i_boundary_cc
+        v_box_p_left = param.beta_p * i_boundary_cc
+        d_vbox_s__dx = (v_box_p_left - v_box_n_right) / l_s
+
+        # Simple formula for velocity in the separator
+        dVbox_dz = pybamm.Concatenation(
+            pybamm.Broadcast(0, "negative electrode"),
+            pybamm.Broadcast(-d_vbox_s__dx, "separator"),
+            pybamm.Broadcast(0, "positive electrode"),
+        )
+        v_box_s = d_vbox_s__dx * (x_s - l_n) + v_box_n_right
+
+        return v_box_s, dVbox_dz
+
+    def get_variables(self, v_box, dVbox_dz):
+        """
+        Calculate dimensionless and dimensional variables for the electrolyte current
+        submodel
+
+        Parameters
+        ----------
+        v_box : :class:`pybamm.Symbol`
+            Volume-averaged velocity in the x-direction
+        dVbox_dz : :class:`pybamm.Symbol`
+            Volume-averaged acceleration in the z-direction (z-derivative of velocity)
+
+        Returns
+        -------
+        dict
+            Dictionary {string: :class:`pybamm.Symbol`} of relevant variables
+        """
+        vel_scale = self.set_of_parameters.velocity_scale
+        L_z = self.set_of_parameters.L_z
+
+        return {
+            "Volume-averaged velocity": v_box,
+            "Volume-averaged velocity [m.s-1]": vel_scale * v_box,
+            "Vertical volume-averaged acceleration": dVbox_dz,
+            "Vertical volume-averaged acceleration [m.s-2]": vel_scale / L_z * dVbox_dz,
+        }
diff --git a/pybamm/parameters/standard_parameters_lead_acid.py b/pybamm/parameters/standard_parameters_lead_acid.py
index 19b09aeb7..e964be578 100644
--- a/pybamm/parameters/standard_parameters_lead_acid.py
+++ b/pybamm/parameters/standard_parameters_lead_acid.py
@@ -319,6 +319,11 @@ beta_liq_n = -c_e_typ * DeltaVliq_n / ne_n_S  # Molar volume change (electrolyte
 beta_liq_p = -c_e_typ * DeltaVliq_p / ne_p_S  # Molar volume change (electrolyte, pos)
 beta_n = beta_surf_n + beta_liq_n  # Total molar volume change (neg)
 beta_p = beta_surf_p + beta_liq_p  # Total molar volume change (pos)
+beta = pybamm.Concatenation(
+    pybamm.Broadcast(beta_n, ["negative electrode"]),
+    pybamm.Broadcast(0, ["separator"]),
+    pybamm.Broadcast(beta_p, ["positive electrode"]),
+)
 beta_Ox = -c_e_typ * (s_plus_Ox * V_plus + s_w_Ox * V_w + s_ox_Ox * V_ox)
 beta_Hy = -c_e_typ * (s_plus_Hy * V_plus + s_hy_Hy * V_hy)
 
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 98667f564..1349108c5 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -34,7 +34,8 @@ class SpatialMethod:
     def spatial_variable(self, symbol):
         """
         Convert a :class:`pybamm.SpatialVariable` node to a linear algebra object that
-        can be evaluated (here, a :class:`pybamm.Vector` on the nodes).
+        can be evaluated (here, a :class:`pybamm.Vector` on either the nodes or the
+        edges).
 
         Parameters
         -----------
@@ -47,7 +48,10 @@ class SpatialMethod:
             Contains the discretised spatial variable
         """
         symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
-        return pybamm.Vector(symbol_mesh[0].nodes, domain=symbol.domain)
+        if symbol.name.endswith("_edge"):
+            return pybamm.Vector(symbol_mesh[0].edges, domain=symbol.domain)
+        else:
+            return pybamm.Vector(symbol_mesh[0].nodes, domain=symbol.domain)
 
     def broadcast(self, symbol, domain):
         """


Test patch:
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 5bc1b04f8..330b70efa 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -47,6 +47,9 @@ class StandardOutputTests(object):
         if self.chemistry == "Lithium-ion":
             self.run_test_class(ParticleConcentrationTests)
 
+        if self.model.options["convection"] is not False:
+            self.run_test_class(VelocityTests)
+
 
 class BaseOutputTest(object):
     def __init__(self, model, param, disc, solver, operating_condition):
@@ -371,8 +374,8 @@ class ElectrolyteConcentrationTests(BaseOutputTest):
         """Test that the internal boundary fluxes are continuous. Test current
         collector fluxes are zero."""
         t, x = self.t, self.x_edge
-        np.testing.assert_array_equal(self.N_e_hat(t, x[0]), 0)
-        np.testing.assert_array_equal(self.N_e_hat(t, x[-1]), 0)
+        np.testing.assert_array_almost_equal(self.N_e_hat(t, x[0]), 0)
+        np.testing.assert_array_almost_equal(self.N_e_hat(t, x[-1]), 0)
 
     def test_splitting(self):
         """Test that when splitting the concentrations and fluxes by negative electrode,
@@ -504,3 +507,44 @@ class CurrentTests(BaseOutputTest):
         # current density will be affected slightly by capacitance effects
         if self.model.options["capacitance"] != "differential":
             self.test_interfacial_current_average()
+
+
+class VelocityTests(BaseOutputTest):
+    def __init__(self, model, param, disc, solver, operating_condition):
+        super().__init__(model, param, disc, solver, operating_condition)
+
+        self.v_box = self.get_var("Volume-averaged velocity")
+        self.i_e = self.get_var("Electrolyte current density")
+        self.dVbox_dz = self.get_var("Vertical volume-averaged acceleration")
+
+    def test_velocity_boundaries(self):
+        """Test the boundary values of the current densities"""
+        np.testing.assert_array_almost_equal(self.v_box(self.t, 0), 0)
+        np.testing.assert_array_almost_equal(self.v_box(self.t, 1), 0)
+
+    def test_vertical_velocity(self):
+        """Test the boundary values of the current densities"""
+        np.testing.assert_array_equal(self.dVbox_dz(self.t, 0), 0)
+        np.testing.assert_array_less(self.dVbox_dz(self.t, 0.5), 0)
+        np.testing.assert_array_equal(self.dVbox_dz(self.t, 1), 0)
+
+    def test_velocity_vs_current(self):
+        """Test the boundary values of the current densities"""
+        t, x_n, x_p = self.t, self.x_n, self.x_p
+
+        beta_n = pybamm.standard_parameters_lead_acid.beta_n
+        beta_n = self.param.process_symbol(beta_n).evaluate()
+        beta_p = pybamm.standard_parameters_lead_acid.beta_p
+        beta_p = self.param.process_symbol(beta_p).evaluate()
+
+        np.testing.assert_array_almost_equal(
+            self.v_box(t, x_n), beta_n * self.i_e(t, x_n)
+        )
+        np.testing.assert_array_almost_equal(
+            self.v_box(t, x_p), beta_p * self.i_e(t, x_p)
+        )
+
+    def test_all(self):
+        self.test_velocity_boundaries()
+        self.test_vertical_velocity()
+        self.test_velocity_vs_current()
diff --git a/tests/integration/test_models/test_lead_acid/test_composite.py b/tests/integration/test_models/test_lead_acid/test_composite.py
index 2e7ba0b7e..e2fea15f4 100644
--- a/tests/integration/test_models/test_lead_acid/test_composite.py
+++ b/tests/integration/test_models/test_lead_acid/test_composite.py
@@ -15,6 +15,11 @@ class TestLeadAcidComposite(unittest.TestCase):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_basic_processing_with_convection(self):
+        model = pybamm.lead_acid.Composite({"convection": True})
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
     def test_optimisations(self):
         model = pybamm.lead_acid.Composite()
         optimtest = tests.OptimisationsTest(model)
diff --git a/tests/integration/test_models/test_lead_acid/test_loqs.py b/tests/integration/test_models/test_lead_acid/test_loqs.py
index 1666e2c0d..e5581a247 100644
--- a/tests/integration/test_models/test_lead_acid/test_loqs.py
+++ b/tests/integration/test_models/test_lead_acid/test_loqs.py
@@ -40,6 +40,11 @@ class TestLeadAcidLOQS(unittest.TestCase):
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all()
 
+    def test_basic_processing_with_convection(self):
+        model = pybamm.lead_acid.LOQS({"convection": True})
+        modeltest = tests.StandardModelTest(model)
+        modeltest.test_all()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/integration/test_models/test_lead_acid/test_newman_tiedemann.py b/tests/integration/test_models/test_lead_acid/test_newman_tiedemann.py
index 81d422659..20499d492 100644
--- a/tests/integration/test_models/test_lead_acid/test_newman_tiedemann.py
+++ b/tests/integration/test_models/test_lead_acid/test_newman_tiedemann.py
@@ -16,6 +16,13 @@ class TestLeadAcidNewmanTiedemann(unittest.TestCase):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all(t_eval=np.linspace(0, 0.6))
 
+    def test_basic_processing_with_convection(self):
+        model = pybamm.lead_acid.NewmanTiedemann({"convection": True})
+        var = pybamm.standard_spatial_vars
+        var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10}
+        modeltest = tests.StandardModelTest(model, var_pts=var_pts)
+        modeltest.test_all(t_eval=np.linspace(0, 0.6))
+
     def test_optimisations(self):
         model = pybamm.lead_acid.NewmanTiedemann()
         optimtest = tests.OptimisationsTest(model)
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 334ad0dc7..287c4a15e 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -705,36 +705,6 @@ class TestDiscretise(unittest.TestCase):
         )[:, np.newaxis]
         np.testing.assert_allclose(eqn_disc.evaluate(), expected_vector)
 
-    def test_discretise_spatial_variable(self):
-        # create discretisation
-        disc = get_discretisation_for_testing()
-
-        # space
-        x1 = pybamm.SpatialVariable("x", ["negative electrode"])
-        x1_disc = disc.process_symbol(x1)
-        self.assertIsInstance(x1_disc, pybamm.Vector)
-        np.testing.assert_array_equal(
-            x1_disc.evaluate(), disc.mesh["negative electrode"][0].nodes[:, np.newaxis]
-        )
-
-        x2 = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
-        x2_disc = disc.process_symbol(x2)
-        self.assertIsInstance(x2_disc, pybamm.Vector)
-        np.testing.assert_array_equal(
-            x2_disc.evaluate(),
-            disc.mesh.combine_submeshes("negative electrode", "separator")[0].nodes[
-                :, np.newaxis
-            ],
-        )
-
-        r = 3 * pybamm.SpatialVariable("r", ["negative particle"])
-        r_disc = disc.process_symbol(r)
-        self.assertIsInstance(r_disc.children[1], pybamm.Vector)
-        np.testing.assert_array_equal(
-            r_disc.evaluate(),
-            3 * disc.mesh["negative particle"][0].nodes[:, np.newaxis],
-        )
-
     def test_exceptions(self):
         c_n = pybamm.Variable("c", domain=["negative electrode"])
         N_n = pybamm.grad(c_n)
diff --git a/tests/unit/test_models/test_lead_acid/test_composite.py b/tests/unit/test_models/test_lead_acid/test_composite.py
index ffb2223f4..d7b43207d 100644
--- a/tests/unit/test_models/test_lead_acid/test_composite.py
+++ b/tests/unit/test_models/test_lead_acid/test_composite.py
@@ -10,6 +10,10 @@ class TestLeadAcidComposite(unittest.TestCase):
         model = pybamm.lead_acid.Composite()
         model.check_well_posedness()
 
+    def test_well_posed_with_convection(self):
+        model = pybamm.lead_acid.Composite({"convection": True})
+        model.check_well_posedness()
+
 
 class TestLeadAcidCompositeCapacitance(unittest.TestCase):
     def test_well_posed_differential(self):
diff --git a/tests/unit/test_models/test_lead_acid/test_loqs.py b/tests/unit/test_models/test_lead_acid/test_loqs.py
index 5066ebbe1..9a0858092 100644
--- a/tests/unit/test_models/test_lead_acid/test_loqs.py
+++ b/tests/unit/test_models/test_lead_acid/test_loqs.py
@@ -10,6 +10,10 @@ class TestLeadAcidLOQS(unittest.TestCase):
         model = pybamm.lead_acid.LOQS()
         model.check_well_posedness()
 
+    def test_well_posed_with_convection(self):
+        model = pybamm.lead_acid.LOQS({"convection": True})
+        model.check_well_posedness()
+
     def test_default_geometry(self):
         model = pybamm.lead_acid.LOQS()
         self.assertIsInstance(model.default_geometry, pybamm.Geometry)
diff --git a/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py b/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py
index 4416176f1..af1c36c04 100644
--- a/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py
+++ b/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py
@@ -10,6 +10,10 @@ class TestLeadAcidNewmanTiedemann(unittest.TestCase):
         model = pybamm.lead_acid.NewmanTiedemann()
         model.check_well_posedness()
 
+    def test_well_posed_with_convection(self):
+        model = pybamm.lead_acid.NewmanTiedemann({"convection": True})
+        model.check_well_posedness()
+
     def test_default_solver(self):
         model = pybamm.lead_acid.NewmanTiedemann()
         self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
diff --git a/tests/unit/test_spatial_methods/test_spatial_methods.py b/tests/unit/test_spatial_methods/test_base_spatial_method.py
similarity index 61%
rename from tests/unit/test_spatial_methods/test_spatial_methods.py
rename to tests/unit/test_spatial_methods/test_base_spatial_method.py
index 5145fc41c..61491e119 100644
--- a/tests/unit/test_spatial_methods/test_spatial_methods.py
+++ b/tests/unit/test_spatial_methods/test_base_spatial_method.py
@@ -1,10 +1,10 @@
 #
 # Test for the base Spatial Method class
 #
+import numpy as np
 import pybamm
-from tests import get_mesh_for_testing, get_1p1d_mesh_for_testing
-
 import unittest
+from tests import get_mesh_for_testing, get_1p1d_mesh_for_testing
 
 
 class TestSpatialMethod(unittest.TestCase):
@@ -21,6 +21,33 @@ class TestSpatialMethod(unittest.TestCase):
         with self.assertRaises(NotImplementedError):
             spatial_method.indefinite_integral(None, None, None)
 
+    def test_discretise_spatial_variable(self):
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_method = pybamm.SpatialMethod(mesh)
+
+        # centre
+        x1 = pybamm.SpatialVariable("x", ["negative electrode"])
+        x2 = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
+        r = pybamm.SpatialVariable("r", ["negative particle"])
+        for var in [x1, x2, r]:
+            var_disc = spatial_method.spatial_variable(var)
+            self.assertIsInstance(var_disc, pybamm.Vector)
+            np.testing.assert_array_equal(
+                var_disc.evaluate()[:, 0], mesh.combine_submeshes(*var.domain)[0].nodes
+            )
+
+        # edges
+        x1_edge = pybamm.SpatialVariable("x_edge", ["negative electrode"])
+        x2_edge = pybamm.SpatialVariable("x_edge", ["negative electrode", "separator"])
+        r_edge = pybamm.SpatialVariable("r_edge", ["negative particle"])
+        for var in [x1_edge, x2_edge, r_edge]:
+            var_disc = spatial_method.spatial_variable(var)
+            self.assertIsInstance(var_disc, pybamm.Vector)
+            np.testing.assert_array_equal(
+                var_disc.evaluate()[:, 0], mesh.combine_submeshes(*var.domain)[0].edges
+            )
+
     def test_broadcast_checks(self):
         child = pybamm.Symbol("sym", domain=["negative electrode"])
         symbol = pybamm.BoundaryFlux(child, "left")