Repository: pybamm-team/PyBaMM (python)

Issue:
Implement delta function source terms
Lead-acid model for recharge needs a source term that acts as a delta function between negative electrode and separator.

Pull request:
#618 implement basic delta function
# Description

Implement delta function on one side of a domain
Fixes #618 

## Type of change

- [x] New feature (non-breaking change which adds functionality)

# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks: 

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works
- [x] Any dependent changes have been merged and published in downstream modules


Patch:
diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst
index 9052d05e6..ed0bc19af 100644
--- a/docs/source/expression_tree/unary_operator.rst
+++ b/docs/source/expression_tree/unary_operator.rst
@@ -43,6 +43,9 @@ Unary Operators
 .. autoclass:: pybamm.BoundaryIntegral
   :members:
 
+.. autoclass:: pybamm.DeltaFunction
+  :members:
+
 .. autoclass:: pybamm.BoundaryOperator
   :members:
 
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index ca8d2e761..5348f770f 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -99,6 +99,7 @@ from .expression_tree.unary_operators import (
     IndefiniteIntegral,
     DefiniteIntegralVector,
     BoundaryIntegral,
+    DeltaFunction,
     grad,
     div,
     laplacian,
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index fc460c2d4..bf712ba68 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -627,6 +627,9 @@ class Discretisation(object):
                     )
                 return symbol
 
+            elif isinstance(symbol, pybamm.DeltaFunction):
+                return spatial_method.delta_function(symbol, disc_child)
+
             elif isinstance(symbol, pybamm.BoundaryOperator):
                 # if boundary operator applied on "negative tab" or
                 # "positive tab" *and* the mesh is 1D then change side to
diff --git a/pybamm/expression_tree/broadcasts.py b/pybamm/expression_tree/broadcasts.py
index a022511b0..5395b880c 100644
--- a/pybamm/expression_tree/broadcasts.py
+++ b/pybamm/expression_tree/broadcasts.py
@@ -124,9 +124,6 @@ class PrimaryBroadcast(Broadcast):
         child node
     broadcast_domain : iterable of str
         Primary domain for broadcast. This will become the domain of the symbol
-    broadcast_type : str, optional
-        Whether to broadcast to the full domain (primary and secondary) or only in the
-        primary direction. Default is "full".
     name : str
         name of the node
 
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index e8b24dc8f..a24ef0c4d 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -595,9 +595,7 @@ class BoundaryIntegral(SpatialOperator):
     def set_id(self):
         """ See :meth:`pybamm.Symbol.set_id()` """
         self._id = hash(
-            (self.__class__, self.name)
-            + (self.children[0].id,)
-            + tuple(self.domain)
+            (self.__class__, self.name) + (self.children[0].id,) + tuple(self.domain)
         )
 
     def _unary_simplify(self, simplified_child):
@@ -619,6 +617,48 @@ class BoundaryIntegral(SpatialOperator):
         return False
 
 
+class DeltaFunction(SpatialOperator):
+    """Delta function. Currently can only be implemented at the edge of a domain
+
+    Parameters
+    ----------
+    child : :class:`pybamm.Symbol`
+        The variable that sets the strength of the delta function
+    side : str
+        Which side of the domain to implement the delta function on
+
+    **Extends:** :class:`SpatialOperator`
+    """
+
+    def __init__(self, child, side, domain):
+        self.side = side
+        if child.domain != []:
+            auxiliary_domains = {"secondary": child.domain}
+        else:
+            auxiliary_domains = {}
+        super().__init__("delta function", child, domain, auxiliary_domains)
+
+    def set_id(self):
+        """ See :meth:`pybamm.Symbol.set_id()` """
+        self._id = hash(
+            (self.__class__, self.name, self.side, self.children[0].id)
+            + tuple(self.domain)
+            + tuple([(k, tuple(v)) for k, v in self.auxiliary_domains.items()])
+        )
+
+    def evaluates_on_edges(self):
+        """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
+        return False
+
+    def _unary_simplify(self, simplified_child):
+        """ See :meth:`UnaryOperator._unary_simplify()`. """
+        return self.__class__(simplified_child, self.side, self.domain)
+
+    def _unary_new_copy(self, child):
+        """ See :meth:`UnaryOperator._unary_new_copy()`. """
+        return self.__class__(child, self.side, self.domain)
+
+
 class BoundaryOperator(SpatialOperator):
     """A node in the expression tree which gets the boundary value of a variable.
 
@@ -666,6 +706,7 @@ class BoundaryOperator(SpatialOperator):
         self._id = hash(
             (self.__class__, self.name, self.side, self.children[0].id)
             + tuple(self.domain)
+            + tuple([(k, tuple(v)) for k, v in self.auxiliary_domains.items()])
         )
 
     def _unary_simplify(self, simplified_child):
@@ -798,6 +839,7 @@ def grad_squared(expression):
 
     return Gradient_Squared(expression)
 
+
 #
 # Method to call SurfaceValue
 #
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 99a55c731..9698b4ffd 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -338,6 +338,47 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return pybamm.Matrix(matrix)
 
+    def delta_function(self, symbol, discretised_symbol):
+        """
+        Delta function. Implemented as a vector whose only non-zero element is the
+        first (if symbol.side = "left") or last (if symbol.side = "right"), with
+        appropriate value so that the integral of the delta function across the whole
+        domain is the same as the integral of the discretised symbol across the whole
+        domain.
+
+        See :meth:`pybamm.SpatialMethod.delta_function`
+        """
+        # Find the number of submeshes
+        submesh_list = self.mesh.combine_submeshes(*symbol.domain)
+
+        prim_pts = submesh_list[0].npts
+        sec_pts = len(submesh_list)
+
+        # Create submatrix to compute delta function as a flux
+        if symbol.side == "left":
+            dx = submesh_list[0].d_nodes[0]
+            sub_matrix = csr_matrix(([1], ([0], [0])), shape=(prim_pts, 1))
+        elif symbol.side == "right":
+            dx = submesh_list[0].d_nodes[-1]
+            sub_matrix = csr_matrix(([1], ([prim_pts - 1], [0])), shape=(prim_pts, 1))
+
+        # Calculate domain width, to make sure that the integral of the delta function
+        # is the same as the integral of the child
+        domain_width = submesh_list[0].edges[-1] - submesh_list[0].edges[0]
+        # Generate full matrix from the submatrix
+        # Convert to csr_matrix so that we can take the index (row-slicing), which is
+        # not supported by the default kron format
+        # Note that this makes column-slicing inefficient, but this should not be an
+        # issue
+        matrix = kron(eye(sec_pts), sub_matrix).toarray()
+
+        # Return delta function, keep domains
+        delta_fn = pybamm.Matrix(domain_width / dx * matrix) * discretised_symbol
+        delta_fn.domain = symbol.domain
+        delta_fn.auxiliary_domains = symbol.auxiliary_domains
+
+        return delta_fn
+
     def internal_neumann_condition(
         self, left_symbol_disc, right_symbol_disc, left_mesh, right_mesh
     ):
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 878b6e23c..7919a8204 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -244,6 +244,19 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
+    def delta_function(self, symbol, discretised_symbol):
+        """
+        Implements the delta function on the approriate side for a spatial method.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol to which is being integrated
+        discretised_symbol: :class:`pybamm.Symbol`
+            The discretised symbol of the correct size
+        """
+        raise NotImplementedError
+
     def internal_neumann_condition(
         self, left_symbol_disc, right_symbol_disc, left_mesh, right_mesh
     ):
@@ -281,7 +294,7 @@ class SpatialMethod:
 
         Returns
         -------
-        :class:`pybamm.Variable`
+        :class:`pybamm.MatrixMultiplication`
             The variable representing the surface value.
         """
         if any(len(self.mesh[dom]) > 1 for dom in discretised_child.domain):
diff --git a/setup.py b/setup.py
index ae509d5e6..23ebed662 100644
--- a/setup.py
+++ b/setup.py
@@ -38,6 +38,7 @@ setup(
         "docs": ["sphinx>=1.5", "guzzle-sphinx-theme"],  # For doc generation
         "dev": [
             "flake8>=3",  # For code style checking
+            "black",  # For code style auto-formatting
             "jupyter",  # For documentation and testing
         ],
     },


Test patch:
diff --git a/tests/unit/test_expression_tree/test_simplify.py b/tests/unit/test_expression_tree/test_simplify.py
index 372b44b42..f53beff41 100644
--- a/tests/unit/test_expression_tree/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_simplify.py
@@ -80,6 +80,12 @@ class TestSimplify(unittest.TestCase):
             (pybamm.boundary_value(v_neg, "right")).simplify(), pybamm.BoundaryValue
         )
 
+        # Delta function
+        self.assertIsInstance(
+            (pybamm.DeltaFunction(v_neg, "right", None)).simplify(),
+            pybamm.DeltaFunction,
+        )
+
         # addition
         self.assertIsInstance((a + b).simplify(), pybamm.Scalar)
         self.assertEqual((a + b).simplify().evaluate(), 1)
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index 387cd8489..d7d7e73b0 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -355,9 +355,10 @@ class TestSymbol(unittest.TestCase):
         state = 2 * pybamm.StateVector(slice(100000))
         self.assertEqual(state.shape, (100000, 1))
 
-    def test_shape_for_testing(self):
+    def test_shape_and_size_for_testing(self):
         scal = pybamm.Scalar(1)
         self.assertEqual(scal.shape_for_testing, scal.shape)
+        self.assertEqual(scal.size_for_testing, scal.size)
 
         state = pybamm.StateVector(slice(10, 25))
         self.assertEqual(state.shape_for_testing, state.shape)
@@ -372,6 +373,7 @@ class TestSymbol(unittest.TestCase):
         self.assertEqual(concat.shape_for_testing, (0,))
         concat = pybamm.Concatenation(state, state)
         self.assertEqual(concat.shape_for_testing, (30, 1))
+        self.assertEqual(concat.size_for_testing, 30)
 
         var = pybamm.Variable("var", domain="negative electrode")
         broadcast = pybamm.Broadcast(0, "negative electrode")
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 156812525..2e4f08cbf 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -181,6 +181,13 @@ class TestUnaryOperators(unittest.TestCase):
         un5 = pybamm.UnaryOperator("test", d)
         self.assertNotEqual(un1.id, un5.id)
 
+    def test_delta_function(self):
+        a = pybamm.Symbol("a")
+        delta_a = pybamm.DeltaFunction(a, "right", "some domain")
+        self.assertEqual(delta_a.side, "right")
+        self.assertEqual(delta_a.child.id, a.id)
+        self.assertFalse(delta_a.evaluates_on_edges())
+
     def test_boundary_operators(self):
         a = pybamm.Symbol("a", domain="some domain")
         boundary_a = pybamm.BoundaryOperator("boundary", a, "right")
diff --git a/tests/unit/test_geometry/test_geometry.py b/tests/unit/test_geometry/test_geometry.py
index b4030bc0e..0cea2d1ff 100644
--- a/tests/unit/test_geometry/test_geometry.py
+++ b/tests/unit/test_geometry/test_geometry.py
@@ -168,60 +168,50 @@ class TestGeometryxp1p1DMicro(unittest.TestCase):
 
 class TestGeometry(unittest.TestCase):
     def test_add_domain(self):
-        L = pybamm.Parameter('L')
+        L = pybamm.Parameter("L")
         zero = pybamm.Scalar(0)
-        x = pybamm.SpatialVariable('x', domain=["negative electrode"])
+        x = pybamm.SpatialVariable("x", domain=["negative electrode"])
         geometry = pybamm.Geometry()
 
-        geometry.add_domain("name_of_domain",
-                            {"primary": {x: {"min": zero, "max": L}}}
-                            )
+        geometry.add_domain("name_of_domain", {"primary": {x: {"min": zero, "max": L}}})
 
-        geometry.add_domain("name_of_domain",
-                            {"tabs": {"negative": {"z_centre": L}}}
-                            )
+        geometry.add_domain("name_of_domain", {"tabs": {"negative": {"z_centre": L}}})
 
         # name must be a string
         with self.assertRaisesRegex(ValueError, "name must be a string"):
-            geometry.add_domain(123,
-                                {"primary": {x: {"min": zero, "max": L}}}
-                                )
+            geometry.add_domain(123, {"primary": {x: {"min": zero, "max": L}}})
 
         # keys of geometry must be either \"primary\" or \"secondary\"
         with self.assertRaisesRegex(ValueError, "primary.*secondary"):
-            geometry.add_domain("name_of_domain",
-                                {"primaryy": {x: {"min": zero, "max": L}}}
-                                )
+            geometry.add_domain(
+                "name_of_domain", {"primaryy": {x: {"min": zero, "max": L}}}
+            )
 
         # inner dict of geometry must have pybamm.SpatialVariable as keys
         with self.assertRaisesRegex(ValueError, "pybamm\.SpatialVariable as keys"):
-            geometry.add_domain("name_of_domain",
-                                {"primary": {L: {"min": zero, "max": L}}}
-                                )
+            geometry.add_domain(
+                "name_of_domain", {"primary": {L: {"min": zero, "max": L}}}
+            )
 
         # no minimum extents for variable
         with self.assertRaisesRegex(ValueError, "minimum"):
-            geometry.add_domain("name_of_domain",
-                                {"primary": {x: {"max": L}}}
-                                )
+            geometry.add_domain("name_of_domain", {"primary": {x: {"max": L}}})
 
         # no maximum extents for variable
         with self.assertRaisesRegex(ValueError, "maximum"):
-            geometry.add_domain("name_of_domain",
-                                {"primary": {x: {"min": zero}}}
-                                )
+            geometry.add_domain("name_of_domain", {"primary": {x: {"min": zero}}})
 
         # tabs region must be \"negative\" or \"positive\"
         with self.assertRaisesRegex(ValueError, "negative.*positive"):
-            geometry.add_domain("name_of_domain",
-                                {"tabs": {"negativee": {"z_centre": L}}}
-                                )
+            geometry.add_domain(
+                "name_of_domain", {"tabs": {"negativee": {"z_centre": L}}}
+            )
 
         # tabs region params must be \"y_centre\", "\"z_centre\" or \"width\"
         with self.assertRaisesRegex(ValueError, "y_centre.*z_centre.*width"):
-            geometry.add_domain("name_of_domain",
-                                {"tabs": {"negative": {"z_centree": L}}}
-                                )
+            geometry.add_domain(
+                "name_of_domain", {"tabs": {"negative": {"z_centree": L}}}
+            )
 
     def test_combine_geometries(self):
         geometry1Dmacro = pybamm.Geometry1DMacro()
diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py
index 26366ef4b..266cc1dde 100644
--- a/tests/unit/test_meshes/test_meshes.py
+++ b/tests/unit/test_meshes/test_meshes.py
@@ -33,7 +33,7 @@ class TestMesh(unittest.TestCase):
         self.assertEqual(len(mesh["negative particle"][0].nodes), var_pts[r])
         self.assertEqual(
             len(mesh["negative particle"][0].edges),
-            len(mesh["negative particle"][0].nodes) + 1
+            len(mesh["negative particle"][0].nodes) + 1,
         )
 
     def test_mesh_creation(self):
diff --git a/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_many_particles.py b/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_many_particles.py
index c63052193..cb842fc9a 100644
--- a/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_many_particles.py
+++ b/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_many_particles.py
@@ -18,17 +18,13 @@ class TestManyParticles(unittest.TestCase):
             pybamm.Scalar(0), "positive electrode", {"secondary": "current collector"}
         )
 
-        variables = {
-            "Negative electrode interfacial current density": a_n,
-        }
+        variables = {"Negative electrode interfacial current density": a_n}
 
         submodel = pybamm.particle.fast.ManyParticles(param, "Negative")
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
-        variables = {
-            "Positive electrode interfacial current density": a_p,
-        }
+        variables = {"Positive electrode interfacial current density": a_p}
         submodel = pybamm.particle.fast.ManyParticles(param, "Positive")
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
diff --git a/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_single_particle.py b/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_single_particle.py
index 2a2d1cf78..6d1a3a21a 100644
--- a/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_single_particle.py
+++ b/tests/unit/test_models/test_submodels/test_particle/test_fast/test_fast_single_particle.py
@@ -12,17 +12,13 @@ class TestSingleParticle(unittest.TestCase):
         param = pybamm.standard_parameters_lithium_ion
 
         a = pybamm.PrimaryBroadcast(pybamm.Scalar(0), "current collector")
-        variables = {
-            "X-averaged negative electrode interfacial current density": a,
-        }
+        variables = {"X-averaged negative electrode interfacial current density": a}
 
         submodel = pybamm.particle.fast.SingleParticle(param, "Negative")
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
-        variables = {
-            "X-averaged positive electrode interfacial current density": a,
-        }
+        variables = {"X-averaged positive electrode interfacial current density": a}
         submodel = pybamm.particle.fast.SingleParticle(param, "Positive")
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
diff --git a/tests/unit/test_models/test_submodels/test_thermal/test_isothermal/test_base_isothermal.py b/tests/unit/test_models/test_submodels/test_thermal/test_isothermal/test_base_isothermal.py
index 2bfc3bd41..254feca74 100644
--- a/tests/unit/test_models/test_submodels/test_thermal/test_isothermal/test_base_isothermal.py
+++ b/tests/unit/test_models/test_submodels/test_thermal/test_isothermal/test_base_isothermal.py
@@ -23,4 +23,3 @@ if __name__ == "__main__":
         debug = True
     pybamm.settings.debug_mode = True
     unittest.main()
-
diff --git a/tests/unit/test_models/test_submodels/test_thermal/test_x_lumped/test_x_lumped_0D_current_collector.py b/tests/unit/test_models/test_submodels/test_thermal/test_x_lumped/test_x_lumped_0D_current_collector.py
index 6d82c71aa..36efd6d74 100644
--- a/tests/unit/test_models/test_submodels/test_thermal/test_x_lumped/test_x_lumped_0D_current_collector.py
+++ b/tests/unit/test_models/test_submodels/test_thermal/test_x_lumped/test_x_lumped_0D_current_collector.py
@@ -31,4 +31,3 @@ if __name__ == "__main__":
         debug = True
     pybamm.settings.debug_mode = True
     unittest.main()
-
diff --git a/tests/unit/test_models/test_submodels/test_thermal/test_xyz_lumped/test_base_xyz_lumped.py b/tests/unit/test_models/test_submodels/test_thermal/test_xyz_lumped/test_base_xyz_lumped.py
index 33596f5c5..db94559a7 100644
--- a/tests/unit/test_models/test_submodels/test_thermal/test_xyz_lumped/test_base_xyz_lumped.py
+++ b/tests/unit/test_models/test_submodels/test_thermal/test_xyz_lumped/test_base_xyz_lumped.py
@@ -23,4 +23,3 @@ if __name__ == "__main__":
         debug = True
     pybamm.settings.debug_mode = True
     unittest.main()
-
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 20509eb90..10d09886d 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -91,6 +91,16 @@ class TestParameterValues(unittest.TestCase):
         self.assertIsInstance(processed_grad.children[0], pybamm.Scalar)
         self.assertEqual(processed_grad.children[0].value, 1)
 
+        # process delta function
+        aa = pybamm.Parameter("a")
+        delta_aa = pybamm.DeltaFunction(aa, "left", "some domain")
+        processed_delta_aa = parameter_values.process_symbol(delta_aa)
+        self.assertIsInstance(processed_delta_aa, pybamm.DeltaFunction)
+        self.assertEqual(processed_delta_aa.side, "left")
+        processed_a = processed_delta_aa.children[0]
+        self.assertIsInstance(processed_a, pybamm.Scalar)
+        self.assertEqual(processed_a.value, 1)
+
         # process boundary operator (test for BoundaryValue)
         aa = pybamm.Parameter("a", domain=["negative electrode"])
         x = pybamm.SpatialVariable("x", domain=["negative electrode"])
diff --git a/tests/unit/test_processed_variable.py b/tests/unit/test_processed_variable.py
index 08c4ae69b..2f360180b 100644
--- a/tests/unit/test_processed_variable.py
+++ b/tests/unit/test_processed_variable.py
@@ -58,10 +58,10 @@ class TestProcessedVariable(unittest.TestCase):
     def test_processed_variable_2D_unknown_domain(self):
         x = pybamm.SpatialVariable("x", domain="SEI layer", coord_sys="cartesian")
         geometry = pybamm.Geometry()
-        geometry.add_domain("SEI layer",
-                            {"primary":
-                                {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
-                            )
+        geometry.add_domain(
+            "SEI layer",
+            {"primary": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}},
+        )
 
         submesh_types = {"SEI layer": pybamm.Uniform1DSubMesh}
         var_pts = {x: 100}
@@ -74,13 +74,11 @@ class TestProcessedVariable(unittest.TestCase):
             np.zeros((var_pts[x], nt)),
             np.linspace(0, 1, 1),
             np.zeros((var_pts[x])),
-            "test"
+            "test",
         )
 
         c = pybamm.StateVector(slice(0, var_pts[x]), domain=["SEI layer"])
-        pybamm.ProcessedVariable(
-            c, solution.t, solution.y, mesh
-        )
+        pybamm.ProcessedVariable(c, solution.t, solution.y, mesh)
 
     def test_processed_variable_3D(self):
         var = pybamm.Variable("var", domain=["negative particle"])
@@ -135,8 +133,7 @@ class TestProcessedVariable(unittest.TestCase):
 
         processed_var = pybamm.ProcessedVariable(var_sol, t_sol, u_sol, mesh=disc.mesh)
         np.testing.assert_array_equal(
-            processed_var.entries,
-            np.reshape(u_sol, [len(y_sol), len(z_sol)]),
+            processed_var.entries, np.reshape(u_sol, [len(y_sol), len(z_sol)])
         )
 
     def test_processed_var_1D_interpolation(self):
diff --git a/tests/unit/test_spatial_methods/test_base_spatial_method.py b/tests/unit/test_spatial_methods/test_base_spatial_method.py
index af0521c5b..1bd9101e5 100644
--- a/tests/unit/test_spatial_methods/test_base_spatial_method.py
+++ b/tests/unit/test_spatial_methods/test_base_spatial_method.py
@@ -26,6 +26,8 @@ class TestSpatialMethod(unittest.TestCase):
             spatial_method.indefinite_integral(None, None)
         with self.assertRaises(NotImplementedError):
             spatial_method.boundary_integral(None, None, None)
+        with self.assertRaises(NotImplementedError):
+            spatial_method.delta_function(None, None)
         with self.assertRaises(NotImplementedError):
             spatial_method.internal_neumann_condition(None, None, None, None)
 
diff --git a/tests/unit/test_spatial_methods/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume.py
index fbe5ba1a7..cdee8c22e 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume.py
@@ -1335,7 +1335,6 @@ class TestFiniteVolume(unittest.TestCase):
         }
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
-        # add ghost nodes
         c_s_n = pybamm.Variable("c_s_n", domain=["negative particle"])
         c_s_p = pybamm.Variable("c_s_p", domain=["positive particle"])
 
@@ -1357,6 +1356,49 @@ class TestFiniteVolume(unittest.TestCase):
         self.assertEqual(c_s_n_surf_disc.domain, ["negative electrode"])
         self.assertEqual(c_s_p_surf_disc.domain, ["positive electrode"])
 
+    def test_delta_function(self):
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var")
+        delta_fn_left = pybamm.DeltaFunction(var, "left", "negative electrode")
+        delta_fn_right = pybamm.DeltaFunction(var, "right", "negative electrode")
+        disc.set_variable_slices([var])
+        delta_fn_left_disc = disc.process_symbol(delta_fn_left)
+        delta_fn_right_disc = disc.process_symbol(delta_fn_right)
+
+        # Basic shape and type tests
+        y = np.ones_like(mesh["negative electrode"][0].nodes[:, np.newaxis])
+        # Left
+        self.assertEqual(delta_fn_left_disc.domain, delta_fn_left.domain)
+        self.assertEqual(
+            delta_fn_left_disc.auxiliary_domains, delta_fn_left.auxiliary_domains
+        )
+        self.assertIsInstance(delta_fn_left_disc, pybamm.Multiplication)
+        self.assertIsInstance(delta_fn_left_disc.left, pybamm.Matrix)
+        np.testing.assert_array_equal(delta_fn_left_disc.left.evaluate()[:, 1:], 0)
+        self.assertEqual(delta_fn_left_disc.shape, y.shape)
+        # Right
+        self.assertEqual(delta_fn_right_disc.domain, delta_fn_right.domain)
+        self.assertEqual(
+            delta_fn_right_disc.auxiliary_domains, delta_fn_right.auxiliary_domains
+        )
+        self.assertIsInstance(delta_fn_right_disc, pybamm.Multiplication)
+        self.assertIsInstance(delta_fn_right_disc.left, pybamm.Matrix)
+        np.testing.assert_array_equal(delta_fn_right_disc.left.evaluate()[:, :-1], 0)
+        self.assertEqual(delta_fn_right_disc.shape, y.shape)
+
+        # Value tests
+        # Delta function should integrate to the same thing as variable
+        var_disc = disc.process_symbol(var)
+        x = pybamm.standard_spatial_vars.x_n
+        delta_fn_int_disc = disc.process_symbol(pybamm.Integral(delta_fn_left, x))
+        np.testing.assert_array_equal(
+            var_disc.evaluate(y=y) * mesh["negative electrode"][0].edges[-1],
+            np.sum(delta_fn_int_disc.evaluate(y=y)),
+        )
+
     def test_grad_div_with_bcs_on_tab(self):
         # 2d macroscale
         mesh = get_1p1d_mesh_for_testing()
diff --git a/tests/unit/test_timer.py b/tests/unit/test_timer.py
index acaa0efed..74d8eabcd 100644
--- a/tests/unit/test_timer.py
+++ b/tests/unit/test_timer.py
@@ -20,7 +20,7 @@ class TestTimer(unittest.TestCase):
         t = pybamm.Timer()
         a = t.time()
         self.assertGreaterEqual(a, 0)
-        for i in range(10):
+        for _ in range(100):
             self.assertGreater(t.time(), a)
         a = t.time()
         t.reset()