Repository: tinosulzer/PyBaMM

Issue:
discritisation of boundary conditions
I was just looking at the discretisation of the boundary conditions. At the moment this seems to be mainly done in the `gradient` and `divergence` functions of `MatrixVectorDiscretisation`, which both do this:

```python
       # Add boundary conditions if defined
        if symbol.id in boundary_conditions:
            lbc, rbc = boundary_conditions[symbol.id]
            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
        gradient_matrix = self.gradient_matrix(domain)
        return gradient_matrix * discretised_symbol
```
So it looks like extra ghost nodes are added to the current state vector to implement the bc. Should this be done instead when the `pybamm.Variable` node is processed in `process_symbol`? Also, the way the boundary conditions are implemented would depend on the method, so this should probably be done in `FiniteVolumeDiscretisation`? @tinosulzer? 



Discussion:

Patch:
diff --git a/.flake8 b/.flake8
index 70990935..9acc35b7 100644
--- a/.flake8
+++ b/.flake8
@@ -5,6 +5,12 @@ exclude=
     problems,
     __init__.py,
     venv,
+    bin,
+    etc,
+    lib,
+    lib64,
+    share,
+    pyvenv.cfg
 ignore=
     # False positive for white space before ':' on list slice
     # black should format these correctly
diff --git a/.gitignore b/.gitignore
index 65bb1425..4db0408e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -46,3 +46,9 @@ pyproject.toml
 
 # virtual enviroment
 venv/
+bin/
+etc/
+lib/
+lib64
+share/
+pyvenv.cfg
diff --git a/docs/source/discretisations/base_discretisation.rst b/docs/source/discretisations/base_discretisation.rst
new file mode 100644
index 00000000..95531920
--- /dev/null
+++ b/docs/source/discretisations/base_discretisation.rst
@@ -0,0 +1,5 @@
+Base Discretisation
+===================
+
+.. autoclass:: pybamm.BaseDiscretisation
+  :members:
diff --git a/docs/source/discretisations/core.rst b/docs/source/discretisations/core.rst
deleted file mode 100644
index 1306f77e..00000000
--- a/docs/source/discretisations/core.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-Core Discretisation
-===================
-
-.. autoclass:: pybamm.BaseDiscretisation
-  :members:
-
-.. autoclass:: pybamm.MatrixVectorDiscretisation
-  :members:
diff --git a/docs/source/discretisations/index.rst b/docs/source/discretisations/index.rst
index 50c9e2bc..d25791de 100644
--- a/docs/source/discretisations/index.rst
+++ b/docs/source/discretisations/index.rst
@@ -3,6 +3,6 @@ Discretisations
 
 .. toctree::
 
-  core
+  base_discretisation
   finite_volume_discretisations
   finite_volume_meshes
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index c392e99c..59856547 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -102,10 +102,7 @@ from .parameters import standard_parameters
 #
 # Mesh and Discretisation classes
 #
-from .discretisations.base_discretisation import (
-    BaseDiscretisation,
-    MatrixVectorDiscretisation,
-)
+from .discretisations.base_discretisation import BaseDiscretisation
 from .discretisations.finite_volume_discretisations import FiniteVolumeDiscretisation
 from .discretisations.base_mesh import KNOWN_DOMAINS
 from .discretisations.base_mesh import BaseMesh, BaseSubmesh
diff --git a/pybamm/discretisations/base_discretisation.py b/pybamm/discretisations/base_discretisation.py
index 7ae9efde..a19743e4 100644
--- a/pybamm/discretisations/base_discretisation.py
+++ b/pybamm/discretisations/base_discretisation.py
@@ -294,78 +294,3 @@ class BaseDiscretisation(object):
 
     def concatenate(self, *symbols):
         return pybamm.NumpyConcatenation(*symbols)
-
-
-class MatrixVectorDiscretisation(BaseDiscretisation):
-    """
-    A base class for discretisations where the gradient and divergence operators are
-    always matrix-vector multiplications.
-
-    Parameters
-    ----------
-     mesh : :class:`pybamm.BaseMesh` (or subclass)
-            The underlying mesh for discretisation
-
-    **Extends:** :class:`pybamm.BaseDiscretisation`
-    """
-
-    def __init__(self, mesh):
-        super().__init__(mesh)
-
-    def gradient(self, symbol, y_slices, boundary_conditions):
-        """Matrix-vector multiplication to implement the gradient operator.
-        See :meth:`pybamm.BaseDiscretisation.gradient`
-        """
-        # Check that boundary condition keys are hashes (ids)
-        for key in boundary_conditions.keys():
-            assert isinstance(key, int), TypeError(
-                "boundary condition keys should be hashes, not {}".format(type(key))
-            )
-        # Discretise symbol
-        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
-        # Add boundary conditions if defined
-        if symbol.id in boundary_conditions:
-            lbc = boundary_conditions[symbol.id]["left"]
-            rbc = boundary_conditions[symbol.id]["right"]
-            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
-        gradient_matrix = self.gradient_matrix(symbol.domain)
-        return gradient_matrix * discretised_symbol
-
-    def gradient_matrix(self, domain):
-        """The gradient matrix
-
-        Parameters
-        ----------
-        domain : list
-            The domain(s) in which to compute the gradient matrix
-        """
-        raise NotImplementedError
-
-    def divergence(self, symbol, y_slices, boundary_conditions):
-        """Matrix-vector multiplication to implement the divergence operator.
-        See :meth:`pybamm.BaseDiscretisation.gradient`
-        """
-        # Check that boundary condition keys are hashes (ids)
-        for key in boundary_conditions.keys():
-            assert isinstance(key, int), TypeError(
-                "boundary condition keys should be hashes, not {}".format(type(key))
-            )
-        # Discretise symbol
-        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
-        # Add boundary conditions if defined
-        if symbol.id in boundary_conditions:
-            lbc = boundary_conditions[symbol.id]["left"]
-            rbc = boundary_conditions[symbol.id]["right"]
-            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
-        divergence_matrix = self.divergence_matrix(symbol.domain)
-        return divergence_matrix * discretised_symbol
-
-    def divergence_matrix(self, domain):
-        """The divergence matrix
-
-        Parameters
-        ----------
-        domain : list
-            The domain(s) in which to compute the divergence matrix
-        """
-        raise NotImplementedError
diff --git a/pybamm/discretisations/finite_volume_discretisations.py b/pybamm/discretisations/finite_volume_discretisations.py
index bfa015d6..b97eae23 100644
--- a/pybamm/discretisations/finite_volume_discretisations.py
+++ b/pybamm/discretisations/finite_volume_discretisations.py
@@ -9,33 +9,112 @@ import numpy as np
 from scipy.sparse import spdiags
 
 
-class FiniteVolumeDiscretisation(pybamm.MatrixVectorDiscretisation):
+class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
     """Discretisation using Finite Volumes.
-    Inherits from :class:`pybamm.MatrixVectorDiscretisation`, so we only need to
-    implement the gradient and divergence matrices
 
     Parameters
     ----------
-     mesh : :class:`BaseMesh` (or subclass)
-            The underlying mesh for discretisation
+    mesh : :class:`pybamm.BaseMesh` (or subclass)
+        The underlying mesh for discretisation
 
-    **Extends:** :class:`pybamm.MatrixVectorDiscretisation`
+    **Extends:** :class:`pybamm.BaseDiscretisation`
     """
 
     def __init__(self, mesh):
         super().__init__(mesh)
 
+    def gradient(self, symbol, y_slices, boundary_conditions):
+        """Matrix-vector multiplication to implement the gradient operator.
+        See :meth:`pybamm.BaseDiscretisation.gradient`
+        """
+        # Check that boundary condition keys are hashes (ids)
+        for key in boundary_conditions.keys():
+            assert isinstance(key, int), TypeError(
+                "boundary condition keys should be hashes, not {}".format(type(key))
+            )
+        # Discretise symbol
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        domain = symbol.domain
+        # Add Dirichlet boundary conditions, if defined
+        if symbol.id in boundary_conditions:
+            lbc = boundary_conditions[symbol.id]["left"]
+            rbc = boundary_conditions[symbol.id]["right"]
+            discretised_symbol = self.add_ghost_nodes(discretised_symbol, lbc, rbc)
+            domain = (
+                [domain[0] + "_left ghost cell"]
+                + domain
+                + [domain[-1] + "_right ghost cell"]
+            )
+        gradient_matrix = self.gradient_matrix(domain)
+        return gradient_matrix * discretised_symbol
+
+    def add_ghost_nodes(self, discretised_symbol, lbc, rbc):
+        """
+        Add Dirichlet boundary conditions via ghost nodes.
+
+        For a boundary condition "y = a at the left-hand boundary",
+        we concatenate a ghost node to the start of the vector y with value "2*a - y1"
+        where y1 is the value of the first node.
+        Similarly for the right-hand boundary condition.
+
+        Currently, Dirichlet boundary conditions can only be applied on state
+        variables (e.g. concentration, temperature), and not on expressions.
+        To access the value of the first node (y1), we create a "first_node" object
+        which is a StateVector whose y_slice is the start of the y_slice of
+        discretised_symbol.
+        Similarly, the last node is a StateVector whose y_slice is the end of the
+        y_slice of discretised_symbol
+
+        Parameters
+        ----------
+        discretised_symbol : :class:`pybamm.StateVector` (size n)
+            The discretised variable (a state vector) to which to add ghost nodes
+        lbc : :class:`pybamm.Scalar`
+            Dirichlet bouncary condition on the left-hand side
+        rbc : :class:`pybamm.Scalar`
+            Dirichlet bouncary condition on the right-hand side
+
+        Returns
+        -------
+        :class:`pybamm.Concatenation` (size n+2)
+            Concatenation of the variable (a state vector) and ghost nodes
+
+        """
+        assert isinstance(discretised_symbol, pybamm.StateVector), NotImplementedError(
+            """discretised_symbol must be a StateVector, not {}""".format(
+                type(discretised_symbol)
+            )
+        )
+        # left ghost cell
+        y_slice_start = discretised_symbol.y_slice.start
+        first_node = pybamm.StateVector(slice(y_slice_start, y_slice_start + 1))
+        left_ghost_cell = 2 * lbc - first_node
+        # right ghost cell
+        y_slice_stop = discretised_symbol.y_slice.stop
+        last_node = pybamm.StateVector(slice(y_slice_stop - 1, y_slice_stop))
+        right_ghost_cell = 2 * rbc - last_node
+        # concatenate
+        return self.concatenate(left_ghost_cell, discretised_symbol, right_ghost_cell)
+
     def gradient_matrix(self, domain):
         """
         Gradient matrix for finite volumes in the appropriate domain.
         Equivalent to grad(y) = (y[1:] - y[:-1])/dx
-        See :meth:`pybamm.MatrixVectorDiscretisation.gradient_matrix()`
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) in which to compute the gradient matrix
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) finite volume gradient matrix for the domain
         """
-        assert len(domain) == 1
-        # TODO: implement for when there are several domains
+        # Create appropriate submesh by combining submeshes in domain
+        submesh = self.mesh.combine_submeshes(*domain)
 
-        # implementation for a single domain
-        submesh = self.mesh[domain[0]]
+        # Create matrix using submesh
         n = submesh.npts
         e = 1 / submesh.d_nodes
         data = np.vstack(
@@ -45,17 +124,44 @@ class FiniteVolumeDiscretisation(pybamm.MatrixVectorDiscretisation):
         matrix = spdiags(data, diags, n - 1, n)
         return pybamm.Matrix(matrix)
 
+    def divergence(self, symbol, y_slices, boundary_conditions):
+        """Matrix-vector multiplication to implement the divergence operator.
+        See :meth:`pybamm.BaseDiscretisation.gradient`
+        """
+        # Check that boundary condition keys are hashes (ids)
+        for key in boundary_conditions.keys():
+            assert isinstance(key, int), TypeError(
+                "boundary condition keys should be hashes, not {}".format(type(key))
+            )
+        # Discretise symbol
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        # Add Neumann boundary conditions if defined
+        if symbol.id in boundary_conditions:
+            lbc = boundary_conditions[symbol.id]["left"]
+            rbc = boundary_conditions[symbol.id]["right"]
+            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
+        divergence_matrix = self.divergence_matrix(symbol.domain)
+        return divergence_matrix * discretised_symbol
+
     def divergence_matrix(self, domain):
         """
         Divergence matrix for finite volumes in the appropriate domain.
         Equivalent to div(N) = (N[1:] - N[:-1])/dx
-        See :meth:`pybamm.MatrixVectorDiscretisation.divergence_matrix()`
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) in which to compute the divergence matrix
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) finite volume divergence matrix for the domain
         """
-        assert len(domain) == 1
-        # TODO: implement for when there are several domains
+        # Create appropriate submesh by combining submeshes in domain
+        submesh = self.mesh.combine_submeshes(*domain)
 
-        # implementation for a single domain
-        submesh = self.mesh[domain[0]]
+        # Create matrix using submesh
         n = submesh.npts + 1
         e = 1 / submesh.d_edges
         data = np.vstack(
diff --git a/pybamm/discretisations/finite_volume_meshes.py b/pybamm/discretisations/finite_volume_meshes.py
index 5a36edf2..08ca3091 100644
--- a/pybamm/discretisations/finite_volume_meshes.py
+++ b/pybamm/discretisations/finite_volume_meshes.py
@@ -51,15 +51,65 @@ class FiniteVolumeMacroMesh(pybamm.BaseMesh):
         self.total_mesh_points = (
             self.neg_mesh_points + (self.sep_mesh_points - 2) + self.pos_mesh_points
         )
-        self["whole cell"] = self.submeshclass(
-            np.concatenate(
-                [
-                    self["negative electrode"].edges,
-                    self["separator"].edges[1:-1],
-                    self["positive electrode"].edges,
-                ]
-            )
+        self["whole cell"] = self.combine_submeshes(
+            "negative electrode", "separator", "positive electrode"
+        )
+
+        # Add ghost meshes for ghost cells for Dirichlet boundary conditions
+        self.add_ghost_meshes()
+
+    def combine_submeshes(self, *submeshnames):
+        """Combine submeshes into a new submesh, using self.submeshclass
+        Raises pybamm.DomainError if submeshes to be combined do not match up (edges are
+        not aligned).
+
+        Parameters
+        ----------
+        submeshnames: list of str
+            The names of the submeshes to be combined
+
+        Returns
+        -------
+        submesh: :class:`self.submeshclass`
+            A new submesh with the class defined by self.submeshclass
+        """
+        # Check that the final edge of each submesh is the same as the first edge of the
+        # next submesh
+        edges_aligned = all(
+            [
+                self[submeshnames[i]].edges[-1] == self[submeshnames[i + 1]].edges[0]
+                for i in range(len(submeshnames) - 1)
+            ]
         )
+        if edges_aligned:
+            # Combine submeshes, being careful not to double-count repeated edges at the
+            # intersection of submeshes
+            combined_submesh_edges = np.concatenate(
+                [self[submeshnames[0]].edges]
+                + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]]
+            )
+            return self.submeshclass(combined_submesh_edges)
+        else:
+            raise pybamm.DomainError("submesh edges are not aligned")
+
+    def add_ghost_meshes(self):
+        """
+        Create meshes for potential ghost nodes on either side of each submesh, using
+        self.submeshclass
+        This will be useful for calculating the gradient with Dirichlet BCs.
+        """
+        # Get all submeshes relating to space (i.e. exclude time)
+        submeshes = [
+            (name, submesh) for name, submesh in self.items() if name != "time"
+        ]
+        for submeshname, submesh in submeshes:
+            edges = submesh.edges
+            # left ghost cell: two edges, one node, to the left of existing submesh
+            lgs_edges = np.array([2 * edges[0] - edges[1], edges[0]])
+            self[submeshname + "_left ghost cell"] = self.submeshclass(lgs_edges)
+            # right ghost cell: two edges, one node, to the right of existing submesh
+            rgs_edges = np.array([edges[-1], 2 * edges[-1] - edges[-2]])
+            self[submeshname + "_right ghost cell"] = self.submeshclass(rgs_edges)
 
 
 class FiniteVolumeSubmesh:


Test patch:
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index adfe5b30..7ff51bea 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -15,19 +15,23 @@ class MeshForTesting(pybamm.BaseMesh):
         self["positive electrode"] = self.submeshclass(self["whole cell"].nodes[40:])
 
 
-class DiscretisationForTesting(pybamm.MatrixVectorDiscretisation):
-    """Interpolating operators."""
+class DiscretisationForTesting(pybamm.BaseDiscretisation):
+    """Identity operators, no boundary conditions."""
 
     def __init__(self, mesh):
         super().__init__(mesh)
 
-    def gradient_matrix(self, domain):
-        n = self.mesh[domain[0]].npts
-        return pybamm.Matrix(np.eye(n))
+    def gradient(self, symbol, y_slices, boundary_conditions):
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        n = self.mesh[symbol.domain[0]].npts
+        gradient_matrix = pybamm.Matrix(np.eye(n))
+        return gradient_matrix * discretised_symbol
 
-    def divergence_matrix(self, domain):
-        n = self.mesh[domain[0]].npts
-        return pybamm.Matrix(np.eye(n))
+    def divergence(self, symbol, y_slices, boundary_conditions):
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        n = self.mesh[symbol.domain[0]].npts
+        divergence_matrix = pybamm.Matrix(np.eye(n))
+        return divergence_matrix * discretised_symbol
 
 
 class TestDiscretise(unittest.TestCase):
@@ -165,18 +169,13 @@ class TestDiscretise(unittest.TestCase):
             disc.gradient(None, None, {})
         with self.assertRaises(NotImplementedError):
             disc.divergence(None, None, {})
-        disc = pybamm.MatrixVectorDiscretisation(None)
-        with self.assertRaises(NotImplementedError):
-            disc.gradient_matrix(None)
-        with self.assertRaises(NotImplementedError):
-            disc.divergence_matrix(None)
 
     def test_process_initial_conditions(self):
         # one equation
         c = pybamm.Variable("c", domain=["whole cell"])
         initial_conditions = {c: pybamm.Scalar(3)}
         mesh = MeshForTesting()
-        disc = DiscretisationForTesting(mesh)
+        disc = pybamm.BaseDiscretisation(mesh)
         y0 = disc.process_initial_conditions(initial_conditions)
         np.testing.assert_array_equal(
             y0[c].evaluate(0, None), 3 * np.ones_like(mesh["whole cell"].nodes)
diff --git a/tests/test_discretisations/test_finite_volume_discretisations.py b/tests/test_discretisations/test_finite_volume_discretisations.py
index 1d00b816..f0c444ed 100644
--- a/tests/test_discretisations/test_finite_volume_discretisations.py
+++ b/tests/test_discretisations/test_finite_volume_discretisations.py
@@ -8,7 +8,92 @@ import unittest
 
 
 class TestFiniteVolumeDiscretisation(unittest.TestCase):
-    def test_grad_div_shapes(self):
+    def test_add_ghost_nodes(self):
+        # Set up
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+        # Add ghost nodes
+        var = pybamm.Variable("var", domain=["whole cell"])
+        y_slices = disc.get_variable_slices([var])
+        discretised_symbol = pybamm.StateVector(y_slices[var.id])
+        lbc = pybamm.Scalar(0)
+        rbc = pybamm.Scalar(3)
+        symbol_plus_ghost = disc.add_ghost_nodes(discretised_symbol, lbc, rbc)
+
+        # Test
+        y_test = np.ones_like(mesh["whole cell"].nodes)
+        np.testing.assert_array_equal(
+            symbol_plus_ghost.evaluate(None, y_test)[1:-1],
+            discretised_symbol.evaluate(None, y_test),
+        )
+        self.assertEqual(
+            (
+                symbol_plus_ghost.evaluate(None, y_test)[0]
+                + symbol_plus_ghost.evaluate(None, y_test)[1]
+            )
+            / 2,
+            0,
+        )
+        self.assertEqual(
+            (
+                symbol_plus_ghost.evaluate(None, y_test)[-2]
+                + symbol_plus_ghost.evaluate(None, y_test)[-1]
+            )
+            / 2,
+            3,
+        )
+
+    def test_grad_div_shapes_Dirichlet_bcs(self):
+        """
+        Test grad and div with Dirichlet boundary conditions (applied by grad on var)
+        """
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+        # grad
+        var = pybamm.Variable("var", domain=["whole cell"])
+        grad_eqn = pybamm.grad(var)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
+        }
+        y_slices = disc.get_variable_slices([var])
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+
+        constant_y = np.ones_like(mesh["whole cell"].nodes)
+        np.testing.assert_array_equal(
+            grad_eqn_disc.evaluate(None, constant_y),
+            np.zeros_like(mesh["whole cell"].edges),
+        )
+
+        # div: test on linear y (should have laplacian zero) so change bcs
+        linear_y = mesh["whole cell"].nodes
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}
+        }
+
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            grad_eqn_disc.evaluate(None, linear_y),
+            np.ones_like(mesh["whole cell"].edges),
+        )
+
+        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(None, linear_y),
+            np.zeros_like(mesh["whole cell"].nodes),
+        )
+
+    def test_grad_div_shapes_Neumann_bcs(self):
+        """Test grad and div with Neumann boundary conditions (applied by div on N)"""
         param = pybamm.ParameterValues(
             base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
         )
@@ -46,7 +131,54 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             np.zeros_like(mesh["whole cell"].nodes),
         )
 
-    def test_grad_convergence(self):
+    def test_grad_div_shapes_mixed_domain(self):
+        """
+        Test grad and div with Dirichlet boundary conditions (applied by grad on var)
+        """
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+        # grad
+        var = pybamm.Variable("var", domain=["negative electrode", "separator"])
+        grad_eqn = pybamm.grad(var)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
+        }
+        y_slices = disc.get_variable_slices([var])
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+
+        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        constant_y = np.ones_like(combined_submesh.nodes)
+        np.testing.assert_array_equal(
+            grad_eqn_disc.evaluate(None, constant_y),
+            np.zeros_like(combined_submesh.edges),
+        )
+
+        # div: test on linear y (should have laplacian zero) so change bcs
+        linear_y = combined_submesh.nodes
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            var.id: {
+                "left": pybamm.Scalar(0),
+                "right": pybamm.Scalar(combined_submesh.edges[-1]),
+            }
+        }
+
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            grad_eqn_disc.evaluate(None, linear_y), np.ones_like(combined_submesh.edges)
+        )
+
+        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(None, linear_y), np.zeros_like(combined_submesh.nodes)
+        )
+
+    def test_grad_convergence_without_bcs(self):
         # Convergence
         param = pybamm.ParameterValues(
             base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
@@ -54,10 +186,8 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         var = pybamm.Variable("var", domain=["whole cell"])
         grad_eqn = pybamm.grad(var)
 
-        # Prepare convergence testing
-        ns = [50, 100, 200]
-        errs = [0] * len(ns)
-        for i, n in enumerate(ns):
+        # Function for convergence testing
+        def get_l2_error(n):
             # Set up discretisation
             mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
             disc = pybamm.FiniteVolumeDiscretisation(mesh)
@@ -72,14 +202,96 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             grad_approx = grad_eqn_disc.evaluate(None, y)
 
             # Calculate errors
-            errs[i] = np.linalg.norm(grad_approx - grad_exact) / np.linalg.norm(
-                grad_exact
-            )
+            return np.linalg.norm(grad_approx - grad_exact) / np.linalg.norm(grad_exact)
 
-        # Expect h**2 convergence
-        [self.assertLess(errs[i + 1] / errs[i], 0.26) for i in range(len(errs) - 1)]
+        # Get errors
+        ns = 100 * (2 ** np.arange(1, 7))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
+
+        # Get rates: expect h**2 convergence
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
+
+    def test_grad_convergence_with_bcs(self):
+        # Convergence
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        var = pybamm.Variable("var", domain=["whole cell"])
+        grad_eqn = pybamm.grad(var)
+        boundary_conditions = {
+            var.id: {
+                "left": pybamm.Scalar(np.sin(0)),
+                "right": pybamm.Scalar(np.sin(1)),
+            }
+        }
+
+        # Function for convergence testing
+        def get_l2_error(n):
+            # Set up discretisation
+            mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
+            disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+            # Define exact solutions
+            y = np.sin(mesh["whole cell"].nodes)
+            grad_exact = np.cos(mesh["whole cell"].edges)
+
+            # Discretise and evaluate
+            y_slices = disc.get_variable_slices([var])
+            grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+            grad_approx = grad_eqn_disc.evaluate(None, y)
+
+            # Calculate errors
+            return np.linalg.norm(grad_approx - grad_exact) / np.linalg.norm(grad_exact)
+
+        # Get errors
+        ns = 100 * (2 ** np.arange(1, 7))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
+
+        # Get rates: expect h**1.5 convergence
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.49 * np.ones_like(rates), rates)
+
+    def test_div_convergence_internal(self):
+        # Convergence
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        var = pybamm.Variable("var", domain=["whole cell"])
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            N.id: {"left": pybamm.Scalar(np.cos(0)), "right": pybamm.Scalar(np.cos(1))}
+        }
+
+        # Function for convergence testing
+        def get_l2_error(n):
+            # Set up discretisation
+            mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
+            disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+            # Define exact solutions
+            y = np.sin(mesh["whole cell"].nodes)
+            div_exact_internal = -np.sin(mesh["whole cell"].nodes[1:-1])
+
+            # Discretise and evaluate
+            y_slices = disc.get_variable_slices([var])
+            div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+            div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
+
+            # Calculate errors
+            return np.linalg.norm(
+                div_approx_internal - div_exact_internal
+            ) / np.linalg.norm(div_exact_internal)
+
+        # Get errors
+        ns = 10 * (2 ** np.arange(1, 6))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
+
+        # Get rates: expect h**2 convergence
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
 
-    @unittest.skip("div errors do not converge as expected")
     def test_div_convergence(self):
         # Convergence
         param = pybamm.ParameterValues(
@@ -92,28 +304,33 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             N.id: {"left": pybamm.Scalar(np.cos(0)), "right": pybamm.Scalar(np.cos(1))}
         }
 
-        # Prepare convergence testing
-        ns = [50, 100, 200]
-        errs = [0] * len(ns)
-        for i, n in enumerate(ns):
+        # Function for convergence testing
+        def get_l2_error(n):
             # Set up discretisation
             mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
             disc = pybamm.FiniteVolumeDiscretisation(mesh)
 
             # Define exact solutions
             y = np.sin(mesh["whole cell"].nodes)
-            div_exact = -np.sin(mesh["whole cell"].nodes)
+            div_exact_internal = -np.sin(mesh["whole cell"].nodes[1:-1])
 
             # Discretise and evaluate
             y_slices = disc.get_variable_slices([var])
             div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
-            div_approx = div_eqn_disc.evaluate(None, y)
+            div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
 
             # Calculate errors
-            errs[i] = np.linalg.norm(div_approx - div_exact) / np.linalg.norm(div_exact)
+            return np.linalg.norm(
+                div_approx_internal - div_exact_internal
+            ) / np.linalg.norm(div_exact_internal)
+
+        # Get errors
+        ns = 10 * (2 ** np.arange(1, 6))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
 
-        # Expect h**2 convergence
-        [self.assertLess(errs[i + 1] / errs[i], 0.26) for i in range(len(errs) - 1)]
+        # Get rates: expect h**1.5 convergence because of boundary conditions
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.49 * np.ones_like(rates), rates)
 
 
 if __name__ == "__main__":
diff --git a/tests/test_discretisations/test_finite_volume_meshes.py b/tests/test_discretisations/test_finite_volume_meshes.py
index 983a83b4..42284f68 100644
--- a/tests/test_discretisations/test_finite_volume_meshes.py
+++ b/tests/test_discretisations/test_finite_volume_meshes.py
@@ -50,6 +50,52 @@ class TestFiniteVolumeMesh(unittest.TestCase):
         self.assertEqual(len(mesh["positive electrode"].edges), mesh.pos_mesh_points)
         self.assertEqual(len(mesh["whole cell"].edges), mesh.total_mesh_points)
 
+    def test_combine_submeshes(self):
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.01, "Ls": 0.5, "Lp": 0.12}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 50)
+        submesh = mesh.combine_submeshes("negative electrode", "separator")
+        self.assertEqual(submesh.edges[0], 0)
+        self.assertEqual(submesh.edges[-1], mesh["separator"].edges[-1])
+        self.assertAlmostEqual(
+            np.linalg.norm(
+                submesh.nodes
+                - np.concatenate(
+                    [mesh["negative electrode"].nodes, mesh["separator"].nodes]
+                )
+            ),
+            0,
+        )
+        with self.assertRaises(pybamm.DomainError):
+            submesh = mesh.combine_submeshes("negative electrode", "positive electrode")
+
+    def test_ghost_cells(self):
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.01, "Ls": 0.5, "Lp": 0.12}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 50)
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges,
+            mesh["whole cell_left ghost cell"].edges,
+        )
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges[1],
+            mesh["negative electrode"].edges[0],
+        )
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges[0],
+            -mesh["negative electrode"].edges[1],
+        )
+        np.testing.assert_array_equal(
+            mesh["positive electrode_right ghost cell"].edges,
+            mesh["whole cell_right ghost cell"].edges,
+        )
+        np.testing.assert_array_equal(
+            mesh["positive electrode_right ghost cell"].edges[0],
+            mesh["positive electrode"].edges[-1],
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")