Repository: tinosulzer/PyBaMM

Issue:
Add mass matrix
Allow construction of a non-trivial mass matrix when discretising the tree to allow for different discretisations e.g. 1D piece-wise finite element, or control volumes. This then needs to be passed to the solver.


Discussion:

Patch:
diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst
index 6884637f..45b50e6b 100644
--- a/docs/source/expression_tree/unary_operator.rst
+++ b/docs/source/expression_tree/unary_operator.rst
@@ -19,7 +19,7 @@ Unary Operators
 .. autoclass:: pybamm.Divergence
   :members:
 
-.. autoclass:: pybamm.BoundaryValue
+.. autoclass:: pybamm.SurfaceValue
   :members:
 
 .. autoclass:: pybamm.Integral
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index e294ff64..464d9cb5 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -75,7 +75,7 @@ from .expression_tree.unary_operators import (
     SpatialOperator,
     Gradient,
     Divergence,
-    BoundaryValue,
+    SurfaceValue,
     Integral,
     grad,
     div,
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index bf7991a1..ac38d53d 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -7,6 +7,7 @@ import pybamm
 
 import copy
 import numpy as np
+from scipy.sparse import block_diag, csr_matrix
 
 
 class Discretisation(object):
@@ -36,6 +37,8 @@ class Discretisation(object):
         self._bcs = {}
         self._y_slices = {}
 
+        mesh.add_ghost_meshes()
+
     @property
     def mesh(self):
         return self._mesh
@@ -75,6 +78,9 @@ class Discretisation(object):
             model.events[idx] = self.process_symbol(event)
         model.concatenated_events = self.concatenate(*model.events)
 
+        # Create mass matrix
+        self.create_mass_matrix(model)
+
         # Check that resulting model makes sense
         self.check_model(model)
 
@@ -152,6 +158,43 @@ class Discretisation(object):
         model.algebraic = self.process_dict(model.algebraic)
         model.concatenated_algebraic = self.concatenate(*model.algebraic.values())
 
+    def create_mass_matrix(self, model):
+        """Creates mass matrix of the discretised model.
+
+        Parameters
+        ----------
+        model : :class:`pybamm.BaseModel` (or subclass)
+            Model to dicretise. Must have attributes rhs, initial_conditions and
+            boundary_conditions (all dicts of {variable: equation})
+        """
+        # Create list of mass matrices for each equation to be put into block
+        # diagonal mass matrix for the model
+        mass_list = []
+
+        # Process mass matrices for the differential equations
+        for var in model.rhs.keys():
+            if var.domain == []:
+                # If variable domain empty then mass matrix is just 1
+                mass_list.append(1.0)
+            else:
+                mass_list.append(
+                    self._spatial_methods[var.domain[0]]
+                    .mass_matrix(var, self._bcs)
+                    .entries
+                )
+
+        # Create lumped mass matrix (of zeros) of the correct shape for the
+        # discretised algebraic equations
+        if model.algebraic.keys():
+            y0 = model.concatenated_initial_conditions
+            mass_algebraic_size = model.concatenated_algebraic.evaluate(0, y0).shape[0]
+            mass_algebraic = csr_matrix((mass_algebraic_size, mass_algebraic_size))
+            mass_list.append(mass_algebraic)
+
+        # Create block diagonal (sparse) mass matrix
+        mass_matrix = block_diag(mass_list)
+        model.mass_matrix = pybamm.Matrix(mass_matrix)
+
     def process_dict(self, var_eqn_dict):
         """Discretise a dictionary of {variable: equation}, broadcasting if necessary
         (can be model.rhs, model.initial_conditions or model.variables).
@@ -235,11 +278,11 @@ class Discretisation(object):
                 )
             return symbol
 
-        elif isinstance(symbol, pybamm.BoundaryValue):
+        elif isinstance(symbol, pybamm.SurfaceValue):
             child = symbol.children[0]
             discretised_child = self.process_symbol(child)
-            return self._spatial_methods[child.domain[0]].boundary_value(
-                discretised_child, symbol.side
+            return self._spatial_methods[child.domain[0]].surface_value(
+                discretised_child
             )
 
         elif isinstance(symbol, pybamm.BinaryOperator):
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index a6ccb55e..4d12e5b2 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -207,24 +207,17 @@ class Integral(SpatialOperator):
         return self._integration_variable
 
 
-class BoundaryValue(SpatialOperator):
-    """A node in the expression tree which gets the boundary value of a variable.
-
-    Parameters
-    ----------
-    child : `pybamm.Symbol`
-        The variable whose boundary value to take
-    side : string
-        Which side to take the boundary value on ("left" or "right")
+class SurfaceValue(SpatialOperator):
+    """A node in the expression tree which gets the surface value of a variable.
 
     **Extends:** :class:`SpatialOperator`
     """
 
-    def __init__(self, child, side):
-        super().__init__("boundary", child)
-        self.side = side
-        # Domain of BoundaryValue must be ([]) so that expressions can be formed
-        # of boundary values of variables in different domains
+    def __init__(self, child):
+        super().__init__("surf", child)
+
+        # Domain of SurfaceValue must be ([]) so that expressions can be formed
+        # of surface values of variables in different domains
         self.domain = []
 
 
@@ -292,4 +285,4 @@ def surf(variable):
         the surface value of ``variable``
     """
 
-    return BoundaryValue(variable, "right")
+    return SurfaceValue(variable)
diff --git a/pybamm/models/base_models.py b/pybamm/models/base_models.py
index 3a78d92a..e6ee3e6a 100644
--- a/pybamm/models/base_models.py
+++ b/pybamm/models/base_models.py
@@ -50,6 +50,7 @@ class BaseModel(object):
         self._events = []
         self._concatenated_rhs = None
         self._concatenated_initial_conditions = None
+        self._mass_matrix = None
 
         # Default parameter values, geometry, submesh, spatial methods and solver
         input_path = os.path.join(
@@ -190,6 +191,14 @@ class BaseModel(object):
     def concatenated_initial_conditions(self, concatenated_initial_conditions):
         self._concatenated_initial_conditions = concatenated_initial_conditions
 
+    @property
+    def mass_matrix(self):
+        return self._mass_matrix
+
+    @mass_matrix.setter
+    def mass_matrix(self, mass_matrix):
+        self._mass_matrix = mass_matrix
+
     def __getitem__(self, key):
         return self.rhs[key]
 
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index b3d66c27..e00c226e 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -179,10 +179,6 @@ class ParameterValues(dict):
             new_child = self.process_symbol(symbol.children[0])
             return pybamm.Integral(new_child, symbol.integration_variable)
 
-        elif isinstance(symbol, pybamm.BoundaryValue):
-            new_child = self.process_symbol(symbol.children[0])
-            return pybamm.BoundaryValue(new_child, symbol.side)
-
         elif isinstance(symbol, pybamm.UnaryOperator):
             new_child = self.process_symbol(symbol.children[0])
             return symbol.__class__(new_child)
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index a4624d61..057a5efe 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -49,20 +49,11 @@ class DaeSolver(pybamm.BaseSolver):
         def algebraic(t, y):
             return model.concatenated_algebraic.evaluate(t, y)
 
-        # Create event-dependent function to evaluate events
-        def event_fun(event):
-            def eval_event(t, y):
-                return event.evaluate(t, y)
-
-            return eval_event
-
-        events = [event_fun(event) for event in model.events]
-
         y0 = self.calculate_consistent_initial_conditions(
             rhs, algebraic, model.concatenated_initial_conditions
         )
 
-        self.t, self.y = self.integrate(residuals, y0, t_eval, events=events)
+        self.t, self.y = self.integrate(residuals, y0, t_eval)
 
     def calculate_consistent_initial_conditions(self, rhs, algebraic, y0_guess):
         """
@@ -103,7 +94,7 @@ class DaeSolver(pybamm.BaseSolver):
 
         return y0_consistent
 
-    def integrate(self, residuals, y0, t_eval, events=None):
+    def integrate(self, residuals, y0, t_eval, events=None, mass_matrix=None):
         """
         Solve a DAE model defined by residuals with initial conditions y0.
 
@@ -119,6 +110,7 @@ class DaeSolver(pybamm.BaseSolver):
         events : method, optional
             A function that takes in t and y and returns conditions for the solver to
             stop
-
+        mass_matrix : array_like
+            The (sparse) mass matrix for the chosen spatial method.
         """
         raise NotImplementedError
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 46dc8609..19c3ea6b 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -35,20 +35,12 @@ class OdeSolver(pybamm.BaseSolver):
         def dydt(t, y):
             return model.concatenated_rhs.evaluate(t, y)
 
-        # Create event-dependent function to evaluate events
-        def event_fun(event):
-            def eval_event(t, y):
-                return event.evaluate(t, y)
-
-            return eval_event
-
-        events = [event_fun(event) for event in model.events]
+        events = [lambda t, y: event.evaluate(t, y) for event in model.events]
 
         y0 = model.concatenated_initial_conditions
-
         self.t, self.y = self.integrate(dydt, y0, t_eval, events=events)
 
-    def integrate(self, derivs, y0, t_eval, events=None):
+    def integrate(self, derivs, y0, t_eval, events=None, mass_matrix=None):
         """
         Solve a model defined by dydt with initial conditions y0.
 
@@ -63,6 +55,7 @@ class OdeSolver(pybamm.BaseSolver):
         events : method, optional
             A function that takes in t and y and returns conditions for the solver to
             stop
-
+        mass_matrix : :class:`pybamm.Matrix`
+            The (sparse) mass matrix for the chosen spatial method.
         """
         raise NotImplementedError
diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py
index 83592860..164fb98d 100644
--- a/pybamm/solvers/scikits_dae_solver.py
+++ b/pybamm/solvers/scikits_dae_solver.py
@@ -43,7 +43,7 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
     def method(self, value):
         self._method = value
 
-    def integrate(self, residuals, y0, t_eval, events=None):
+    def integrate(self, residuals, y0, t_eval, events=None, mass_matrix=None):
         """
         Solve a DAE model defined by residuals with initial conditions y0.
 
@@ -59,7 +59,8 @@ class ScikitsDaeSolver(pybamm.DaeSolver):
         events : method, optional
             A function that takes in t and y and returns conditions for the solver to
             stop
-
+        mass_matrix : array_like
+            The (sparse) mass matrix for the chosen spatial method.
         """
 
         def eqsres(t, y, ydot, return_residuals):
diff --git a/pybamm/solvers/scikits_ode_solver.py b/pybamm/solvers/scikits_ode_solver.py
index 74b96e13..853cb8a7 100644
--- a/pybamm/solvers/scikits_ode_solver.py
+++ b/pybamm/solvers/scikits_ode_solver.py
@@ -43,7 +43,7 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
     def method(self, value):
         self._method = value
 
-    def integrate(self, derivs, y0, t_eval, events=None):
+    def integrate(self, derivs, y0, t_eval, events=None, mass_matrix=None):
         """
         Solve a model defined by dydt with initial conditions y0.
 
@@ -58,6 +58,8 @@ class ScikitsOdeSolver(pybamm.OdeSolver):
         events : method, optional
             A function that takes in t and y and returns conditions for the solver to
             stop
+        mass_matrix : array_like
+            The (sparse) mass matrix for the chosen spatial method.
 
         """
 
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index e47de6d7..504170d4 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -32,7 +32,7 @@ class ScipySolver(pybamm.OdeSolver):
     def method(self, value):
         self._method = value
 
-    def integrate(self, derivs, y0, t_eval, events=None):
+    def integrate(self, derivs, y0, t_eval, events=None, mass_matrix=None):
         """
         Solve a model defined by dydt with initial conditions y0.
 
@@ -48,6 +48,8 @@ class ScipySolver(pybamm.OdeSolver):
         events : method, optional
             A function that takes in t and y and returns conditions for the solver to
             stop
+        mass_matrix : array_like
+            The (sparse) mass matrix for the chosen spatial method.
 
         Returns
         -------
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 2352b69e..8e70f2ab 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -387,38 +387,75 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return new_discretised_symbol
 
-    def boundary_value(self, discretised_symbol, side):
+    def surface_value(self, discretised_symbol):
         """
-        Uses linear extrapolation to get the boundary value of a variable in the
+        Uses linear extrapolation to get the surface value of a variable in the
         Finite Volume Method.
 
         Parameters
         -----------
         discretised_symbol : :class:`pybamm.StateVector`
-            The discretised variable from which to calculate the boundary value
-        side : string
-            Which side to take the boundary value on ("left" or "right")
+            The discretised variable (a state vector) from which to calculate
+            the surface value.
 
         Returns
         -------
-        :class:`pybamm.Symbol`
-            The variable representing the boundary value.
+        :class:`pybamm.Variable`
+            The variable representing the surface value.
         """
+        # Better to make class similar NodeToEdge and pass function?
+        # def surface_value(array):
+        #     "Linear extrapolation for surface value"
+        #     array[-1] + (array[-1] - array[-2]) / 2
+        # ... or make StateVector and add?
+        y_slice_stop = discretised_symbol.y_slice.stop
+        last_node = pybamm.StateVector(slice(y_slice_stop - 1, y_slice_stop))
+        penultimate_node = pybamm.StateVector(slice(y_slice_stop - 2, y_slice_stop - 1))
+        surface_value = last_node + (last_node - penultimate_node) / 2
+        return surface_value
+
+    def mass_matrix(self, symbol, boundary_conditions):
+        """
+        Calculates the mass matrix for a spatial method.
 
-        def linear_extrapolation(array):
-            """Linearly extrapolates an array"""
-            if side == "left":
-                return array[0] + (array[0] - array[1]) / 2
-            elif side == "right":
-                return array[-1] + (array[-1] - array[-2]) / 2
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Variable`
+            The variable corresponding to the equation for which we are
+            calculating the mass matrix.
+        boundary_conditions : dict
+            The boundary conditions of the model
+            ({symbol.id: {"left": left bc, "right": right bc}})
 
-        return BoundaryValueEvaluated(discretised_symbol, linear_extrapolation)
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) mass matrix for the spatial method.
+        """
+        # NOTE: for different spatial methods the matrix may need to be adjusted
+        # to account for Dirichlet boundary conditions. Here, we just have that
+        # the mass matrix is the identity.
+
+        # Create appropriate submesh by combining submeshes in domain
+        submesh = self.mesh.combine_submeshes(*symbol.domain)
+
+        # Get number of points in primary dimension
+        n = submesh[0].npts
+
+        # Create mass matrix for primary dimension
+        prim_mass = eye(n)
+
+        # Get number of points in secondary dimension
+        sec_pts = len(submesh)
+
+        mass = kron(eye(sec_pts), prim_mass)
+        return pybamm.Matrix(mass)
 
     #######################################################
     # Can probably be moved outside of the spatial method
     ######################################################
 
-    def compute_diffusivity(self, discretised_symbol):
+    def compute_diffusivity(self, symbol):
         """
         Compute the diffusivity at cell edges, based on the diffusivity at cell nodes.
         For now we just take the arithemtic mean, though it may be better to take the
@@ -429,7 +466,7 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         Parameters
         ----------
-        discretised_symbol : :class:`pybamm.Symbol`
+        symbol : :class:`pybamm.Symbol`
             Symbol to be averaged. When evaluated, this symbol returns either a scalar
             or an array of shape (n,), where n is the number of points in the mesh for
             the symbol's domain (n = self.mesh[symbol.domain].npts)
@@ -445,37 +482,7 @@ class FiniteVolume(pybamm.SpatialMethod):
             """Calculate the arithemetic mean of an array"""
             return (array[1:] + array[:-1]) / 2
 
-        return pybamm.NodeToEdge(discretised_symbol, arithmetic_mean)
-
-
-class BoundaryValueEvaluated(pybamm.SpatialOperator):
-    """A node in the expression tree representing a unary operator that evaluates the
-    value of its child at a boundary.
-
-    Parameters
-    ----------
-    child : :class:`Symbol`
-        child node
-    boundary_function : method
-        the function used to calculate the boundary value
-
-    **Extends:** :class:`pybamm.SpatialOperator`
-    """
-
-    def __init__(self, child, boundary_function):
-        """ See :meth:`pybamm.UnaryOperator.__init__()`. """
-        super().__init__(
-            "boundary value ({})".format(boundary_function.__name__), child
-        )
-        self._boundary_function = boundary_function
-        # Domain of BoundaryValue must be ([]) so that expressions can be formed
-        # of boundary values of variables in different domains
-        self.domain = []
-
-    def evaluate(self, t=None, y=None):
-        """ See :meth:`pybamm.Symbol.evaluate()`. """
-        evaluated_child = self.children[0].evaluate(t, y)
-        return self._boundary_function(evaluated_child)
+        return pybamm.NodeToEdge(symbol, arithmetic_mean)
 
 
 class NodeToEdge(pybamm.SpatialOperator):
@@ -485,6 +492,8 @@ class NodeToEdge(pybamm.SpatialOperator):
     Parameters
     ----------
 
+    name : str
+        name of the node
     child : :class:`Symbol`
         child node
     node_to_edge_function : method
@@ -506,6 +515,7 @@ class NodeToEdge(pybamm.SpatialOperator):
         evaluated_child = self.children[0].evaluate(t, y)
         # If the evaluated child is a numpy array of shape (n,), do the averaging
         # NOTE: Doing this check every time might be slow?
+        # NOTE: Will need to deal with 2D arrays at some point
         if isinstance(evaluated_child, np.ndarray) and len(evaluated_child.shape) == 1:
             return self._node_to_edge_function(evaluated_child)
         # If not, no need to average
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 7e095388..72f05ba4 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -141,6 +141,26 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
+    def mass_matrix(self, symbol, boundary_conditions):
+        """
+        Calculates the mass matrix for a spatial method.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Variable`
+            The variable corresponding to the equation for which we are
+            calculating the mass matrix.
+        boundary_conditions : dict
+            The boundary conditions of the model
+            ({symbol.id: {"left": left bc, "right": right bc}})
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) mass matrix for the spatial method.
+        """
+        raise NotImplementedError
+
     # We could possibly move the following outside of SpatialMethod
     # depending on the requirements of the FiniteVolume
 


Test patch:
diff --git a/tests/shared.py b/tests/shared.py
index ac4206e7..6ed712be 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -47,6 +47,13 @@ class SpatialMethodForTesting(pybamm.SpatialMethod):
         divergence_matrix = pybamm.Matrix(np.eye(n))
         return divergence_matrix @ discretised_symbol
 
+    def mass_matrix(self, symbol, boundary_conditions):
+        n = 0
+        for domain in symbol.domain:
+            n += self.mesh[domain][0].npts
+        mass_matrix = pybamm.Matrix(np.eye(n))
+        return mass_matrix
+
     def compute_diffusivity(self, symbol):
         return symbol
 
diff --git a/tests/test_discretisations/test_discretisation.py b/tests/test_discretisations/test_discretisation.py
index 418a4af1..b09a0243 100644
--- a/tests/test_discretisations/test_discretisation.py
+++ b/tests/test_discretisations/test_discretisation.py
@@ -6,6 +6,7 @@ import pybamm
 import numpy as np
 import unittest
 from tests import get_mesh_for_testing, get_discretisation_for_testing
+from scipy.sparse import block_diag
 
 
 class TestDiscretise(unittest.TestCase):
@@ -343,10 +344,17 @@ class TestDiscretise(unittest.TestCase):
         y0 = model.concatenated_initial_conditions
         np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh[0].nodes))
         np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
+
         # grad and div are identity operators here
         np.testing.assert_array_equal(y0, model.variables["c"].evaluate(None, y0))
         np.testing.assert_array_equal(y0, model.variables["N"].evaluate(None, y0))
 
+        # mass matrix is identity
+        np.testing.assert_array_equal(
+            np.eye(combined_submesh[0].nodes.shape[0]),
+            model.mass_matrix.entries.toarray(),
+        )
+
         # several equations
         T = pybamm.Variable("T", domain=["negative electrode"])
         q = pybamm.grad(T)
@@ -386,6 +394,11 @@ class TestDiscretise(unittest.TestCase):
         )
         np.testing.assert_array_equal(S0 * T0, model.variables["ST"].evaluate(None, y0))
 
+        # mass matrix is identity
+        np.testing.assert_array_equal(
+            np.eye(np.size(y0)), model.mass_matrix.entries.toarray()
+        )
+
         # test that not enough initial conditions raises an error
         model = pybamm.BaseModel()
         model.rhs = {c: pybamm.div(N), T: pybamm.div(q), S: pybamm.div(p)}
@@ -439,6 +452,22 @@ class TestDiscretise(unittest.TestCase):
             np.zeros_like(combined_submesh[0].nodes),
         )
 
+        # mass matrix is identity upper left, zeros elsewhere
+        mass = block_diag(
+            (
+                np.eye(np.size(combined_submesh[0].nodes)),
+                np.zeros(
+                    (
+                        np.size(combined_submesh[0].nodes),
+                        np.size(combined_submesh[0].nodes),
+                    )
+                ),
+            )
+        )
+        np.testing.assert_array_equal(
+            mass.toarray(), model.mass_matrix.entries.toarray()
+        )
+
     def test_process_model_concatenation(self):
         # concatenation of variables as the key
         cn = pybamm.Variable("c", domain=["negative electrode"])
diff --git a/tests/test_models/test_li_ion.py b/tests/test_models/test_li_ion.py
index 904ab6ae..c289c5f8 100644
--- a/tests/test_models/test_li_ion.py
+++ b/tests/test_models/test_li_ion.py
@@ -17,7 +17,7 @@ class TestLiIonSPM(unittest.TestCase):
 
         modeltest.test_all()
 
-    def test_surface_concentration(self):
+    def test_surface_concentrartion(self):
         model = pybamm.li_ion.SPM()
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
@@ -26,16 +26,13 @@ class TestLiIonSPM(unittest.TestCase):
         # check surface concentration decreases in negative particle and
         # increases in positive particle for discharge
         np.testing.assert_array_less(
-            model.variables["cn_surf"].evaluate(T, Y)[1:],
-            model.variables["cn_surf"].evaluate(T, Y)[:-1],
+            model.variables["cn_surf"].evaluate(T, Y)[:, 1:],
+            model.variables["cn_surf"].evaluate(T, Y)[:, :-1],
         )
         np.testing.assert_array_less(
-            model.variables["cp_surf"].evaluate(T, Y)[:-1],
-            model.variables["cp_surf"].evaluate(T, Y)[1:],
+            model.variables["cp_surf"].evaluate(T, Y)[:, :-1],
+            model.variables["cp_surf"].evaluate(T, Y)[:, 1:],
         )
-        # test that surface concentrations are all positive
-        np.testing.assert_array_less(0, model.variables["cn_surf"].evaluate(T, Y))
-        np.testing.assert_array_less(0, model.variables["cp_surf"].evaluate(T, Y))
 
 
 if __name__ == "__main__":
diff --git a/tests/test_solvers/test_scikits_solvers.py b/tests/test_solvers/test_scikits_solvers.py
index e4cdc799..cbccac74 100644
--- a/tests/test_solvers/test_scikits_solvers.py
+++ b/tests/test_solvers/test_scikits_solvers.py
@@ -70,7 +70,7 @@ class TestScikitsSolver(unittest.TestCase):
         y0 = np.array([1])
         t_eval = np.linspace(0, 3, 100)
         t_sol, y_sol = solver.integrate(
-            exponential_growth, y0, t_eval, events=[ysq_eq_7, y_eq_9]
+            exponential_growth, y0, t_eval, events=[y_eq_9, ysq_eq_7]
         )
         self.assertLess(len(t_sol), len(t_eval))
         np.testing.assert_allclose(np.exp(t_sol), y_sol[0], rtol=1e-4)
@@ -213,10 +213,7 @@ class TestScikitsSolver(unittest.TestCase):
         var = pybamm.Variable("var", domain=whole_cell)
         model.rhs = {var: 0.1 * var}
         model.initial_conditions = {var: 1}
-        model.events = [
-            pybamm.Function(np.min, 2 * var - 2.5),
-            pybamm.Function(np.min, var - 1.5),
-        ]
+        model.events = [pybamm.Function(np.min, var - 1.5)]
         disc = StandardModelTest(model).disc
         disc.process_model(model)
 
@@ -226,7 +223,6 @@ class TestScikitsSolver(unittest.TestCase):
         solver.solve(model, t_eval)
         np.testing.assert_allclose(solver.y[0], np.exp(0.1 * solver.t))
         np.testing.assert_array_less(solver.y[0], 1.5)
-        np.testing.assert_array_less(solver.y[0], 1.25)
 
     def test_model_solver_dae(self):
         # Create model
@@ -286,7 +282,7 @@ class TestScikitsSolver(unittest.TestCase):
 
         # Solve
         solver = pybamm.ScikitsDaeSolver(tol=1e-8)
-        t_eval = np.linspace(0, 5, 100)
+        t_eval = np.linspace(0, 1, 100)
         solver.solve(model, t_eval)
         np.testing.assert_array_less(solver.y[0], 1.5)
         np.testing.assert_array_less(solver.y[-1], 2.5)
diff --git a/tests/test_spatial_methods/test_finite_volume.py b/tests/test_spatial_methods/test_finite_volume.py
index 17300738..8efa6ac5 100644
--- a/tests/test_spatial_methods/test_finite_volume.py
+++ b/tests/test_spatial_methods/test_finite_volume.py
@@ -5,6 +5,7 @@ import pybamm
 from tests import get_mesh_for_testing, get_p2d_mesh_for_testing
 
 import numpy as np
+from scipy.sparse import kron, eye
 import unittest
 
 
@@ -32,40 +33,14 @@ class TestFiniteVolume(unittest.TestCase):
         avd = pybamm.NodeToEdge(d, arithmetic_mean)
         np.testing.assert_array_equal(avd.evaluate(None, y_test), np.ones(9))
 
-    def test_extrapolate_left_right(self):
+    def test_surface_value(self):
         # create discretisation
         mesh = get_mesh_for_testing()
-        spatial_methods = {
-            "macroscale": pybamm.FiniteVolume,
-            "negative particle": pybamm.FiniteVolume,
-        }
+        spatial_methods = {"negative particle": pybamm.FiniteVolume}
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        macro_submesh = mesh.combine_submeshes(*whole_cell)
-        micro_submesh = mesh["negative particle"]
-
-        # Macroscale
-        # create variable
-        var = pybamm.Variable("var", domain=whole_cell)
-        # boundary value should work with something more complicated than a variable
-        extrap_left = pybamm.BoundaryValue(2 * var, "left")
-        extrap_right = pybamm.BoundaryValue(4 - var, "right")
-        disc.set_variable_slices([var])
-        extrap_left_disc = disc.process_symbol(extrap_left)
-        extrap_right_disc = disc.process_symbol(extrap_right)
-
-        # check constant extrapolates to constant
-        constant_y = np.ones_like(macro_submesh[0].nodes)
-        self.assertEqual(extrap_left_disc.evaluate(None, constant_y), 2)
-        self.assertEqual(extrap_right_disc.evaluate(None, constant_y), 3)
-
-        # check linear variable extrapolates correctly
-        linear_y = macro_submesh[0].nodes
-        self.assertEqual(extrap_left_disc.evaluate(None, linear_y), 0)
-        self.assertEqual(extrap_right_disc.evaluate(None, linear_y), 3)
+        combined_submesh = mesh.combine_submeshes("negative particle")
 
-        # Microscale
         # create variable
         var = pybamm.Variable("var", domain="negative particle")
         surf_eqn = pybamm.surf(var)
@@ -73,12 +48,12 @@ class TestFiniteVolume(unittest.TestCase):
         surf_eqn_disc = disc.process_symbol(surf_eqn)
 
         # check constant extrapolates to constant
-        constant_y = np.ones_like(micro_submesh[0].nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         self.assertEqual(surf_eqn_disc.evaluate(None, constant_y), 1)
 
         # check linear variable extrapolates correctly
-        linear_y = micro_submesh[0].nodes
-        y_surf = micro_submesh[0].nodes[-1] + micro_submesh[0].d_nodes[-1] / 2
+        linear_y = combined_submesh[0].nodes
+        y_surf = combined_submesh[0].nodes[-1] + combined_submesh[0].d_nodes[-1] / 2
         self.assertEqual(surf_eqn_disc.evaluate(None, linear_y), y_surf)
 
     def test_discretise_diffusivity_times_spatial_operator(self):
@@ -693,10 +668,10 @@ class TestFiniteVolume(unittest.TestCase):
         integral_eqn_disc = disc.process_symbol(integral_eqn)
 
         combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
-        constant_y = np.ones_like(combined_submesh.nodes)
-        constant_y_edges = np.ones_like(combined_submesh.edges)
-        linear_y = combined_submesh.nodes
-        linear_y_edges = combined_submesh.edges
+        constant_y = np.ones_like(combined_submesh[0].nodes)
+        constant_y_edges = np.ones_like(combined_submesh[0].edges)
+        linear_y = combined_submesh[0].nodes
+        linear_y_edges = combined_submesh[0].edges
         np.testing.assert_array_equal(
             integral_eqn_disc.evaluate(None, constant_y_edges), linear_y
         )
@@ -723,7 +698,7 @@ class TestFiniteVolume(unittest.TestCase):
         np.testing.assert_array_almost_equal(
             integral_eqn_disc.evaluate(None, linear_y), (linear_y ** 2 - (ln) ** 2) / 2
         )
-        cos_y = np.cos(combined_submesh.nodes)
+        cos_y = np.cos(combined_submesh[0].nodes)
         np.testing.assert_array_almost_equal(
             integral_eqn_disc.evaluate(None, cos_y),
             np.sin(linear_y) - np.sin(ln),
@@ -737,15 +712,15 @@ class TestFiniteVolume(unittest.TestCase):
         disc.set_variable_slices([var])
         integral_eqn_disc = disc.process_symbol(integral_eqn)
 
-        constant_y = np.ones_like(mesh["negative particle"].nodes)
+        constant_y = np.ones_like(mesh["negative particle"][0].nodes)
         np.testing.assert_array_equal(
             integral_eqn_disc.evaluate(None, constant_y), np.pi
         )
-        linear_y = mesh["negative particle"].nodes
+        linear_y = mesh["negative particle"][0].nodes
         np.testing.assert_array_almost_equal(
             integral_eqn_disc.evaluate(None, linear_y), 2 * np.pi / 3, places=4
         )
-        one_over_y = 1 / mesh["negative particle"].nodes
+        one_over_y = 1 / mesh["negative particle"][0].nodes
         np.testing.assert_array_equal(
             integral_eqn_disc.evaluate(None, one_over_y), 2 * np.pi
         )
@@ -1054,6 +1029,60 @@ class TestFiniteVolume(unittest.TestCase):
             r_disc.evaluate(), 3 * disc.mesh["negative particle"][0].nodes
         )
 
+    def test_mass_matrix_shape(self):
+        """
+        Test mass matrix shape
+        """
+        # one equation
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        c = pybamm.Variable("c", domain=whole_cell)
+        N = pybamm.grad(c)
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N)}
+        model.initial_conditions = {c: pybamm.Scalar(0)}
+        model.boundary_conditions = {
+            c: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
+        }
+        model.variables = {"c": c, "N": N}
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        disc.process_model(model)
+
+        # mass matrix
+        mass = np.eye(combined_submesh[0].npts)
+        np.testing.assert_array_equal(mass, model.mass_matrix.entries.toarray())
+
+    def test_p2d_mass_matrix_shape(self):
+        """
+        Test mass matrix shape in the pseudo 2-dimensional case
+        """
+        c = pybamm.Variable("c", domain=["negative particle"])
+        N = pybamm.grad(c)
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N)}
+        model.initial_conditions = {c: pybamm.Scalar(0)}
+        model.boundary_conditions = {
+            c: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
+        }
+        model.variables = {"c": c, "N": N}
+        mesh = get_p2d_mesh_for_testing()
+        spatial_methods = {"negative particle": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        prim_pts = mesh["negative particle"][0].npts
+        sec_pts = len(mesh["negative particle"])
+        mass_local = eye(prim_pts)
+        mass = kron(eye(sec_pts), mass_local)
+        np.testing.assert_array_equal(
+            mass.toarray(), model.mass_matrix.entries.toarray()
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_spatial_methods/test_spatial_methods.py b/tests/test_spatial_methods/test_spatial_methods.py
index 51f59f0a..6f4686b3 100644
--- a/tests/test_spatial_methods/test_spatial_methods.py
+++ b/tests/test_spatial_methods/test_spatial_methods.py
@@ -24,6 +24,8 @@ class TestSpatialMethod(unittest.TestCase):
             spatial_method.surface_value(None)
         with self.assertRaises(NotImplementedError):
             spatial_method.compute_diffusivity()
+        with self.assertRaises(NotImplementedError):
+            spatial_method.mass_matrix(None, None)
 
 
 if __name__ == "__main__":