Repository: pybamm-team/PyBaMM

Issue:
Move `Discretisation.process_binary_operators` to the spatial method
This is a little bit slow (if we're being picky) and entirely unneeded in cases where the spatial method is something other than Finite Volumes (e.g. ODE models)

Discussion:

Patch:
diff --git a/docs/tutorials/add-spatial-method.rst b/docs/tutorials/add-spatial-method.rst
index 652982691..09f14382f 100644
--- a/docs/tutorials/add-spatial-method.rst
+++ b/docs/tutorials/add-spatial-method.rst
@@ -51,18 +51,18 @@ and add the class to ``pybamm/__init__.py``:
 You can then start implementing the spatial method by adding functions to the class.
 In particular, any spatial method *must* have the following functions (from the base class :class:`pybamm.SpatialMethod`):
 
-- :meth:`pybamm.SpatialMethod.spatial_variable`
 - :meth:`pybamm.SpatialMethod.gradient`
 - :meth:`pybamm.SpatialMethod.divergence`
 - :meth:`pybamm.SpatialMethod.integral`
 - :meth:`pybamm.SpatialMethod.indefinite integral`
-- :meth:`pybamm.SpatialMethod.boundary_value`
 
 Optionally, a new spatial method can also overwrite the default behaviour for the following functions:
 
+- :meth:`pybamm.SpatialMethod.spatial_variable`
 - :meth:`pybamm.SpatialMethod.broadcast`
 - :meth:`pybamm.SpatialMethod.mass_matrix`
-- :meth:`pybamm.SpatialMethod.compute_diffusivity`
+- :meth:`pybamm.SpatialMethod.process_binary_operators`
+- :meth:`pybamm.SpatialMethod.boundary_value`
 
 For an example of an existing spatial method implementation, see the Finite Volume
 `API docs <https://pybamm.readthedocs.io/en/latest/source/spatial_methods/finite_volume.html>`_
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 09b90fa0a..337f1d89b 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -304,7 +304,19 @@ class Discretisation(object):
             Discretised symbol
 
         """
-        if isinstance(symbol, pybamm.Gradient):
+        if isinstance(symbol, pybamm.BinaryOperator):
+            # Pre-process children
+            left, right = symbol.children
+            disc_left = self.process_symbol(left)
+            disc_right = self.process_symbol(right)
+            if symbol.domain == []:
+                return symbol.__class__(disc_left, disc_right)
+            else:
+                return self._spatial_methods[symbol.domain[0]].process_binary_operators(
+                    symbol, left, right, disc_left, disc_right
+                )
+
+        elif isinstance(symbol, pybamm.Gradient):
             child = symbol.children[0]
             discretised_child = self.process_symbol(child)
             return self._spatial_methods[child.domain[0]].gradient(
@@ -352,9 +364,6 @@ class Discretisation(object):
                 child, discretised_child, symbol.side
             )
 
-        elif isinstance(symbol, pybamm.BinaryOperator):
-            return self.process_binary_operators(symbol)
-
         elif isinstance(symbol, pybamm.Function):
             new_child = self.process_symbol(symbol.children[0])
             return pybamm.Function(symbol.func, new_child)
@@ -391,69 +400,6 @@ class Discretisation(object):
                 "Cannot discretise symbol of type '{}'".format(type(symbol))
             )
 
-    def process_binary_operators(self, bin_op):
-        """Discretise binary operators in model equations.  Performs appropriate
-        averaging of diffusivities if one of the children is a gradient operator, so
-        that discretised sizes match up.  This is mainly an issue for the Finite Volume
-        Discretisation: see
-        :meth:`pybamm.FiniteVolumeDiscretisation.compute_diffusivity()`
-
-        Parameters
-        ----------
-        bin_op : :class:`pybamm.BinaryOperator`
-            Binary operator to discretise
-
-        Returns
-        -------
-        :class:`pybamm.BinaryOperator`
-            Discretised binary operator
-
-        """
-        # Pre-process children
-        left, right = bin_op.children
-        new_left = self.process_symbol(left)
-        new_right = self.process_symbol(right)
-        # Post-processing to make sure discretised dimensions match
-        # If neither child has gradients, or both children have gradients
-        # no need to do any averaging
-        if (
-            left.has_gradient_and_not_divergence()
-            == right.has_gradient_and_not_divergence()
-        ):
-            pass
-        # If only left child has gradient, compute diffusivity for right child
-        elif (
-            left.has_gradient_and_not_divergence()
-            and not right.has_gradient_and_not_divergence()
-        ):
-            # Extrapolate at either end depending on the ghost cells (from gradient)
-            extrapolate_left = any(
-                [x.has_left_ghost_cell for x in new_left.pre_order()]
-            )
-            extrapolate_right = any(
-                [x.has_right_ghost_cell for x in new_left.pre_order()]
-            )
-            new_right = self._spatial_methods[bin_op.domain[0]].compute_diffusivity(
-                new_right, extrapolate_left, extrapolate_right
-            )
-        # If only right child has gradient, compute diffusivity for left child
-        elif (
-            right.has_gradient_and_not_divergence()
-            and not left.has_gradient_and_not_divergence()
-        ):
-            # Extrapolate at either end depending on the ghost cells (from gradient)
-            extrapolate_left = any(
-                [x.has_left_ghost_cell for x in new_right.pre_order()]
-            )
-            extrapolate_right = any(
-                [x.has_right_ghost_cell for x in new_right.pre_order()]
-            )
-            new_left = self._spatial_methods[bin_op.domain[0]].compute_diffusivity(
-                new_left, extrapolate_left, extrapolate_right
-            )
-        # Return new binary operator with appropriate class
-        return bin_op.__class__(new_left, new_right)
-
     def concatenate(self, *symbols):
         return pybamm.NumpyConcatenation(*symbols)
 
diff --git a/pybamm/models/lead_acid/loqs.py b/pybamm/models/lead_acid/loqs.py
index 3ebe7ab71..40a03eb8a 100644
--- a/pybamm/models/lead_acid/loqs.py
+++ b/pybamm/models/lead_acid/loqs.py
@@ -88,5 +88,11 @@ class LOQS(pybamm.LeadAcidBaseModel):
 
         "-----------------------------------------------------------------------------"
         "Settings"
-        # ODEs only (don't use jacobian)
+        # ODEs only (don't use jacobian, use base spatial method)
         self.use_jacobian = False
+
+        self.default_spatial_methods = {
+            "macroscale": pybamm.SpatialMethod,
+            "negative particle": pybamm.SpatialMethod,
+            "positive particle": pybamm.SpatialMethod,
+        }
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 546396c8f..4573ed997 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -26,11 +26,11 @@ class FiniteVolume(pybamm.SpatialMethod):
     """
 
     def __init__(self, mesh):
+        super().__init__(mesh)
         # add npts_for_broadcast to mesh domains for this particular discretisation
         for dom in mesh.keys():
             for i in range(len(mesh[dom])):
                 mesh[dom][i].npts_for_broadcast = mesh[dom][i].npts
-        super().__init__(mesh)
 
     def spatial_variable(self, symbol):
         """
@@ -556,6 +556,70 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return boundary_value
 
+    def process_binary_operators(self, bin_op, left, right, disc_left, disc_right):
+        """Discretise binary operators in model equations.  Performs appropriate
+        averaging of diffusivities if one of the children is a gradient operator, so
+        that discretised sizes match up.
+
+        Parameters
+        ----------
+        bin_op : :class:`pybamm.BinaryOperator`
+            Binary operator to discretise
+        left : :class:`pybamm.Symbol`
+            The left child of `bin_op`
+        right : :class:`pybamm.Symbol`
+            The right child of `bin_op`
+        disc_left : :class:`pybamm.Symbol`
+            The discretised left child of `bin_op`
+        disc_right : :class:`pybamm.Symbol`
+            The discretised right child of `bin_op`
+        Returns
+        -------
+        :class:`pybamm.BinaryOperator`
+            Discretised binary operator
+
+        """
+        # Post-processing to make sure discretised dimensions match
+        # If neither child has gradients, or both children have gradients
+        # no need to do any averaging
+        if (
+            left.has_gradient_and_not_divergence()
+            == right.has_gradient_and_not_divergence()
+        ):
+            pass
+        # If only left child has gradient, compute diffusivity for right child
+        elif (
+            left.has_gradient_and_not_divergence()
+            and not right.has_gradient_and_not_divergence()
+        ):
+            # Extrapolate at either end depending on the ghost cells (from gradient)
+            extrapolate_left = any(
+                [x.has_left_ghost_cell for x in disc_left.pre_order()]
+            )
+            extrapolate_right = any(
+                [x.has_right_ghost_cell for x in disc_left.pre_order()]
+            )
+            disc_right = self.compute_diffusivity(
+                disc_right, extrapolate_left, extrapolate_right
+            )
+        # If only right child has gradient, compute diffusivity for left child
+        elif (
+            right.has_gradient_and_not_divergence()
+            and not left.has_gradient_and_not_divergence()
+        ):
+            # Extrapolate at either end depending on the ghost cells (from gradient)
+            extrapolate_left = any(
+                [x.has_left_ghost_cell for x in disc_right.pre_order()]
+            )
+            extrapolate_right = any(
+                [x.has_right_ghost_cell for x in disc_right.pre_order()]
+            )
+            disc_left = self.compute_diffusivity(
+                disc_left, extrapolate_left, extrapolate_right
+            )
+        # Return new binary operator with appropriate class
+        return bin_op.__class__(disc_left, disc_right)
+
     def compute_diffusivity(
         self, discretised_symbol, extrapolate_left=False, extrapolate_right=False
     ):
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 575013dcb..18396e975 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -2,7 +2,7 @@
 # A general spatial method class
 #
 import pybamm
-from scipy.sparse import eye, kron
+from scipy.sparse import eye, kron, coo_matrix
 
 
 class SpatialMethod:
@@ -20,6 +20,10 @@ class SpatialMethod:
     """
 
     def __init__(self, mesh):
+        # add npts_for_broadcast to mesh domains for this particular discretisation
+        for dom in mesh.keys():
+            for i in range(len(mesh[dom])):
+                mesh[dom][i].npts_for_broadcast = mesh[dom][i].npts
         self._mesh = mesh
 
     @property
@@ -29,7 +33,7 @@ class SpatialMethod:
     def spatial_variable(self, symbol):
         """
         Convert a :class:`pybamm.SpatialVariable` node to a linear algebra object that
-        can be evaluated (e.g. a :class:`pybamm.Vector`).
+        can be evaluated (here, a :class:`pybamm.Vector` on the nodes).
 
         Parameters
         -----------
@@ -41,7 +45,8 @@ class SpatialMethod:
         :class:`pybamm.Vector`
             Contains the discretised spatial variable
         """
-        raise NotImplementedError
+        symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
+        return pybamm.Vector(symbol_mesh[0].nodes)
 
     def broadcast(self, symbol, domain):
         """
@@ -153,7 +158,9 @@ class SpatialMethod:
     def boundary_value(self, symbol, discretised_symbol, side):
         """
         Returns the boundary value using the approriate expression for the
-        spatial method.
+        spatial method. To do this, we create a sparse vector 'bv_vector' that extracts
+        either the first (for side="left") or last (for side="right") point from
+        'discretised_symbol'.
 
         Parameters
         -----------
@@ -169,7 +176,21 @@ class SpatialMethod:
         :class:`pybamm.Variable`
             The variable representing the surface value.
         """
-        raise NotImplementedError
+        n = sum(self.mesh[dom][0].npts for dom in symbol.domain)
+        if side == "left":
+            # coo_matrix takes inputs (data, (row, col)) and puts data[i] at the point
+            # (row[i], col[i]) for each index of data. Here we just want a single point
+            # with value 1 at (0,0).
+            left_vector = coo_matrix(([1], ([0], [0])), shape=(1, n))
+            bv_vector = pybamm.Matrix(left_vector)
+        elif side == "right":
+            # as above, but now we want a single point with value 1 at (0, n-1)
+            right_vector = coo_matrix(([1], ([0], [n - 1])), shape=(1, n))
+            bv_vector = pybamm.Matrix(right_vector)
+        out = bv_vector @ discretised_symbol
+        # boundary value removes domain
+        out.domain = []
+        return out
 
     def mass_matrix(self, symbol, boundary_conditions):
         """
@@ -208,12 +229,27 @@ class SpatialMethod:
         mass = kron(eye(sec_pts), prim_mass)
         return pybamm.Matrix(mass)
 
-    def compute_diffusivity(
-        self, symbol, extrapolate_left=None, extrapolate_right=None
-    ):
-        """Compute the diffusivity at edges of cells.
-        Could interpret this as: find diffusivity as
-        off grid locations
+    def process_binary_operators(self, bin_op, left, right, disc_left, disc_right):
+        """Discretise binary operators in model equations. Default behaviour is to
+        return a new binary operator with the discretised children.
+
+        Parameters
+        ----------
+        bin_op : :class:`pybamm.BinaryOperator`
+            Binary operator to discretise
+        left : :class:`pybamm.Symbol`
+            The left child of `bin_op`
+        right : :class:`pybamm.Symbol`
+            The right child of `bin_op`
+        disc_left : :class:`pybamm.Symbol`
+            The discretised left child of `bin_op`
+        disc_right : :class:`pybamm.Symbol`
+            The discretised right child of `bin_op`
+
+        Returns
+        -------
+        :class:`pybamm.BinaryOperator`
+            Discretised binary operator
+
         """
-        # Default behaviour (identity operator): return symbol
-        return symbol
+        return bin_op.__class__(disc_left, disc_right)


Test patch:
diff --git a/tests/shared.py b/tests/shared.py
index f991f4c60..a4225955e 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -3,23 +3,15 @@
 #
 import pybamm
 
-from scipy.sparse import eye, coo_matrix
+from scipy.sparse import eye
 
 
 class SpatialMethodForTesting(pybamm.SpatialMethod):
     """Identity operators, no boundary conditions."""
 
     def __init__(self, mesh):
-        for dom in mesh.keys():
-            for i in range(len(mesh[dom])):
-                mesh[dom][i].npts_for_broadcast = mesh[dom][i].npts
         super().__init__(mesh)
 
-    def spatial_variable(self, symbol):
-        # use the cell centres
-        symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
-        return pybamm.Vector(symbol_mesh[0].nodes)
-
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
         n = 0
         for domain in symbol.domain:
@@ -41,21 +33,6 @@ class SpatialMethodForTesting(pybamm.SpatialMethod):
         mass_matrix = pybamm.Matrix(eye(n))
         return mass_matrix
 
-    def boundary_value(self, symbol, discretised_symbol, side):
-        n = 0
-        for domain in symbol.domain:
-            n += self.mesh[domain][0].npts
-        if side == "left":
-            left_vector = coo_matrix(([1], ([0], [0])), shape=(1, n))
-            bv_vector = pybamm.Matrix(left_vector)
-        elif side == "right":
-            right_vector = coo_matrix(([1], ([0], [n - 1])), shape=(1, n))
-            bv_vector = pybamm.Matrix(right_vector)
-        out = bv_vector @ discretised_symbol
-        # boundary value removes domain
-        out.domain = []
-        return out
-
 
 def get_mesh_for_testing(npts=None):
     param = pybamm.ParameterValues(
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 7c96caf64..3708532d3 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -98,12 +98,17 @@ class TestDiscretise(unittest.TestCase):
     def test_process_symbol_base(self):
         # create discretisation
         mesh = get_mesh_for_testing()
-        spatial_methods = {}
+        spatial_methods = {
+            "macroscale": pybamm.SpatialMethod,
+            "negative particle": pybamm.SpatialMethod,
+            "positive particle": pybamm.SpatialMethod,
+        }
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
         # variable
         var = pybamm.Variable("var")
-        disc._y_slices = {var.id: slice(53)}
+        var_vec = pybamm.Variable("var vec", domain=["negative electrode"])
+        disc._y_slices = {var.id: slice(53), var_vec.id: slice(53, 102)}
         var_disc = disc.process_symbol(var)
         self.assertIsInstance(var_disc, pybamm.StateVector)
         self.assertEqual(var_disc._y_slice, disc._y_slices[var.id])
@@ -147,6 +152,18 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(un2_disc, pybamm.AbsoluteValue)
         self.assertIsInstance(un2_disc.children[0], pybamm.Scalar)
 
+        # boundary value
+        bv_left = pybamm.BoundaryValue(var_vec, "left")
+        bv_left_disc = disc.process_symbol(bv_left)
+        self.assertIsInstance(bv_left_disc, pybamm.MatrixMultiplication)
+        self.assertIsInstance(bv_left_disc.left, pybamm.Matrix)
+        self.assertIsInstance(bv_left_disc.right, pybamm.StateVector)
+        bv_right = pybamm.BoundaryValue(var_vec, "left")
+        bv_right_disc = disc.process_symbol(bv_right)
+        self.assertIsInstance(bv_right_disc, pybamm.MatrixMultiplication)
+        self.assertIsInstance(bv_right_disc.left, pybamm.Matrix)
+        self.assertIsInstance(bv_right_disc.right, pybamm.StateVector)
+
         # not implemented
         sym = pybamm.Symbol("sym")
         with self.assertRaises(NotImplementedError):
@@ -240,19 +257,6 @@ class TestDiscretise(unittest.TestCase):
                 eqn_disc.evaluate(None, y), var_disc.evaluate(None, y) ** 2
             )
 
-    def test_core_NotImplementedErrors(self):
-        # create spatial method
-        spatial_method = pybamm.SpatialMethod(None)
-
-        with self.assertRaises(NotImplementedError):
-            spatial_method.gradient(None, None, {})
-        with self.assertRaises(NotImplementedError):
-            spatial_method.divergence(None, None, {})
-        with self.assertRaises(NotImplementedError):
-            spatial_method.integral(None, None, None)
-        with self.assertRaises(NotImplementedError):
-            spatial_method.indefinite_integral(None, None, None)
-
     def test_process_dict(self):
         # one equation
         whole_cell = ["negative electrode", "separator", "positive electrode"]
diff --git a/tests/unit/test_spatial_methods/test_spatial_methods.py b/tests/unit/test_spatial_methods/test_spatial_methods.py
index 7adae1b7a..90ac95c39 100644
--- a/tests/unit/test_spatial_methods/test_spatial_methods.py
+++ b/tests/unit/test_spatial_methods/test_spatial_methods.py
@@ -12,14 +12,14 @@ class TestSpatialMethod(unittest.TestCase):
         mesh = get_mesh_for_testing()
         spatial_method = pybamm.SpatialMethod(mesh)
         self.assertEqual(spatial_method.mesh, mesh)
-        with self.assertRaises(NotImplementedError):
-            spatial_method.spatial_variable(None)
         with self.assertRaises(NotImplementedError):
             spatial_method.gradient(None, None, None)
         with self.assertRaises(NotImplementedError):
             spatial_method.divergence(None, None, None)
         with self.assertRaises(NotImplementedError):
-            spatial_method.boundary_value(None, None, None)
+            spatial_method.integral(None, None, None)
+        with self.assertRaises(NotImplementedError):
+            spatial_method.indefinite_integral(None, None, None)
 
 
 if __name__ == "__main__":