Repository: tinosulzer/PyBaMM

Issue:
Evaluate "left" and "right"
Adapt `SurfaceValue` (from #159) to work with any object, not just `StateVector`, and to extrapolate at either end of the vector

**Motivation**
We need this to get the voltage from DFN 

Discussion:

Patch:
diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst
index 45b50e6b..6884637f 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.SurfaceValue
+.. autoclass:: pybamm.BoundaryValue
   :members:
 
 .. autoclass:: pybamm.Integral
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 464d9cb5..e294ff64 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -75,7 +75,7 @@ from .expression_tree.unary_operators import (
     SpatialOperator,
     Gradient,
     Divergence,
-    SurfaceValue,
+    BoundaryValue,
     Integral,
     grad,
     div,
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index b20b3694..bf7991a1 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -235,11 +235,11 @@ class Discretisation(object):
                 )
             return symbol
 
-        elif isinstance(symbol, pybamm.SurfaceValue):
+        elif isinstance(symbol, pybamm.BoundaryValue):
             child = symbol.children[0]
             discretised_child = self.process_symbol(child)
-            return self._spatial_methods[child.domain[0]].surface_value(
-                discretised_child
+            return self._spatial_methods[child.domain[0]].boundary_value(
+                discretised_child, symbol.side
             )
 
         elif isinstance(symbol, pybamm.BinaryOperator):
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 4d12e5b2..a6ccb55e 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -207,17 +207,24 @@ class Integral(SpatialOperator):
         return self._integration_variable
 
 
-class SurfaceValue(SpatialOperator):
-    """A node in the expression tree which gets the surface value of a 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")
 
     **Extends:** :class:`SpatialOperator`
     """
 
-    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
+    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
         self.domain = []
 
 
@@ -285,4 +292,4 @@ def surf(variable):
         the surface value of ``variable``
     """
 
-    return SurfaceValue(variable)
+    return BoundaryValue(variable, "right")
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index e00c226e..b3d66c27 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -179,6 +179,10 @@ 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/ode_solver.py b/pybamm/solvers/ode_solver.py
index be96eaef..46dc8609 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -45,6 +45,7 @@ class OdeSolver(pybamm.BaseSolver):
         events = [event_fun(event) 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):
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 6714f70f..2352b69e 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -387,38 +387,38 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return new_discretised_symbol
 
-    def surface_value(self, discretised_symbol):
+    def boundary_value(self, discretised_symbol, side):
         """
-        Uses linear extrapolation to get the surface value of a variable in the
+        Uses linear extrapolation to get the boundary value of a variable in the
         Finite Volume Method.
 
         Parameters
         -----------
         discretised_symbol : :class:`pybamm.StateVector`
-            The discretised variable (a state vector) from which to calculate
-            the surface value.
+            The discretised variable from which to calculate the boundary value
+        side : string
+            Which side to take the boundary value on ("left" or "right")
 
         Returns
         -------
-        :class:`pybamm.Variable`
-            The variable representing the surface value.
+        :class:`pybamm.Symbol`
+            The variable representing the boundary 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 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
+
+        return BoundaryValueEvaluated(discretised_symbol, linear_extrapolation)
 
     #######################################################
     # Can probably be moved outside of the spatial method
     ######################################################
 
-    def compute_diffusivity(self, symbol):
+    def compute_diffusivity(self, discretised_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 +429,7 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         Parameters
         ----------
-        symbol : :class:`pybamm.Symbol`
+        discretised_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,7 +445,37 @@ class FiniteVolume(pybamm.SpatialMethod):
             """Calculate the arithemetic mean of an array"""
             return (array[1:] + array[:-1]) / 2
 
-        return pybamm.NodeToEdge(symbol, arithmetic_mean)
+        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)
 
 
 class NodeToEdge(pybamm.SpatialOperator):
@@ -455,8 +485,6 @@ class NodeToEdge(pybamm.SpatialOperator):
     Parameters
     ----------
 
-    name : str
-        name of the node
     child : :class:`Symbol`
         child node
     node_to_edge_function : method
@@ -478,7 +506,6 @@ 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


Test patch:
diff --git a/tests/test_models/test_li_ion.py b/tests/test_models/test_li_ion.py
index dfe48abd..904ab6ae 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_concentrartion(self):
+    def test_surface_concentration(self):
         model = pybamm.li_ion.SPM()
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
@@ -26,12 +26,12 @@ 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))
diff --git a/tests/test_spatial_methods/test_finite_volume.py b/tests/test_spatial_methods/test_finite_volume.py
index 42e92b14..17300738 100644
--- a/tests/test_spatial_methods/test_finite_volume.py
+++ b/tests/test_spatial_methods/test_finite_volume.py
@@ -32,14 +32,40 @@ 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_surface_value(self):
+    def test_extrapolate_left_right(self):
         # create discretisation
         mesh = get_mesh_for_testing()
-        spatial_methods = {"negative particle": pybamm.FiniteVolume}
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+        }
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
-        combined_submesh = mesh.combine_submeshes("negative particle")
+        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)
+
+        # Microscale
         # create variable
         var = pybamm.Variable("var", domain="negative particle")
         surf_eqn = pybamm.surf(var)
@@ -47,12 +73,12 @@ class TestFiniteVolume(unittest.TestCase):
         surf_eqn_disc = disc.process_symbol(surf_eqn)
 
         # check constant extrapolates to constant
-        constant_y = np.ones_like(combined_submesh[0].nodes)
+        constant_y = np.ones_like(micro_submesh[0].nodes)
         self.assertEqual(surf_eqn_disc.evaluate(None, constant_y), 1)
 
         # check linear variable extrapolates correctly
-        linear_y = combined_submesh[0].nodes
-        y_surf = combined_submesh[0].nodes[-1] + combined_submesh[0].d_nodes[-1] / 2
+        linear_y = micro_submesh[0].nodes
+        y_surf = micro_submesh[0].nodes[-1] + micro_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):