Repository: tinosulzer/PyBaMM

Issue:
integration expression tree node
need to be able to define integrations over space of the form:
$$
f(x) = \int_0^x  f(x')dx'
$$

These will be replaced by a matrix-vector multiply by the discretisation

Discussion:

Patch:
diff --git a/docs/source/expression_tree/unary_operator.rst b/docs/source/expression_tree/unary_operator.rst
index c547d1af..45b50e6b 100644
--- a/docs/source/expression_tree/unary_operator.rst
+++ b/docs/source/expression_tree/unary_operator.rst
@@ -22,6 +22,9 @@ Unary Operators
 .. autoclass:: pybamm.SurfaceValue
   :members:
 
+.. autoclass:: pybamm.Integral
+  :members:
+
 .. autofunction:: pybamm.grad
 
 .. autofunction:: pybamm.div
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 3a879f1b..06979c74 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -77,6 +77,7 @@ from .expression_tree.unary_operators import (
     Gradient,
     Divergence,
     SurfaceValue,
+    Integral,
     grad,
     div,
     surf,
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index dd090985..3c38e6f3 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -218,6 +218,13 @@ class Discretisation(object):
                 child, discretised_child, self._bcs
             )
 
+        elif isinstance(symbol, pybamm.Integral):
+            child = symbol.children[0]
+            discretised_child = self.process_symbol(child)
+            return self._spatial_methods[symbol.domain[0]].integral(
+                child, discretised_child
+            )
+
         elif isinstance(symbol, pybamm.Broadcast):
             # Process child first
             new_child = self.process_symbol(symbol.children[0])
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 88a939ce..55ffc092 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -17,7 +17,6 @@ class UnaryOperator(pybamm.Symbol):
 
     Parameters
     ----------
-
     name : str
         name of the node
     child : :class:`Symbol`
@@ -128,6 +127,54 @@ class Divergence(SpatialOperator):
         super().__init__("div", child)
 
 
+class Integral(SpatialOperator):
+    """A node in the expression tree representing an integral operator (definite or
+    indefinite)
+
+    .. math::
+        \\text{definite}: \\quad I = \\int_{a}^{b}\\!f(u)\\,du,
+
+        \\text{indefinite}: \\quad I(s) = \\int_{a}^{s}\\!f(u)\\,du,
+
+    where :math:`a` and :math:`b` are the left-hand and right-hand boundaries of
+    the domain respectively, and :math:`s\\in\\text{domain}`.
+    Can be integration with respect to time or space.
+
+    Parameters
+    ----------
+    function : :class:`pybamm.Symbol`
+        The function to be integrated (will become self.children[0])
+    integration_variable : :class:`pybamm.IndependentVariable`
+        The variable over which to integrate
+
+    **Extends:** :class:`SpatialOperator`
+    """
+
+    def __init__(self, child, integration_variable):
+        if isinstance(integration_variable, pybamm.SpatialVariable):
+            # Check that child and integration_variable domains agree
+            if child.domain != integration_variable.domain:
+                raise pybamm.DomainError(
+                    """child and integration_variable must have the same domain"""
+                )
+        elif not isinstance(integration_variable, pybamm.IndependentVariable):
+            raise ValueError(
+                """integration_variable must be of type pybamm.IndependentVariable,
+                   not {}""".format(
+                    type(integration_variable)
+                )
+            )
+        name = "integral d{}".format(integration_variable.name)
+        if isinstance(integration_variable, pybamm.SpatialVariable):
+            name += " {}".format(integration_variable.domain)
+        super().__init__(name, child)
+        self._integration_variable = integration_variable
+
+    @property
+    def integration_variable(self):
+        return self._integration_variable
+
+
 class SurfaceValue(SpatialOperator):
     """A node in the expression tree which gets the surface value of a variable.
 
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index c0553017..298daa92 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -126,7 +126,7 @@ class FiniteVolume(pybamm.SpatialMethod):
 
     def divergence(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the divergence operator.
-        See :meth:`pybamm.SpatialMethod.gradient`
+        See :meth:`pybamm.SpatialMethod.divergence`
         """
         # Check that boundary condition keys are hashes (ids)
         for key in boundary_conditions.keys():
@@ -142,20 +142,15 @@ class FiniteVolume(pybamm.SpatialMethod):
             discretised_symbol = pybamm.NumpyConcatenation(lbc, discretised_symbol, rbc)
 
         domain = symbol.domain
-        # check for spherical domains
+        # check for particle domain
+        divergence_matrix = self.divergence_matrix(domain)
         if ("negative particle" or "positive particle") in domain:
-
-            # implement spherical operator
-            divergence_matrix = self.divergence_matrix(domain)
-
-            submesh = self.mesh[domain[0]]
+            submesh = self.mesh.combine_submeshes(*domain)
             r = pybamm.Vector(submesh.nodes)
             r_edges = pybamm.Vector(submesh.edges)
-
             out = (1 / (r ** 2)) * (
                 divergence_matrix @ ((r_edges ** 2) * discretised_symbol)
             )
-
         else:
             divergence_matrix = self.divergence_matrix(domain)
             out = divergence_matrix @ discretised_symbol
@@ -189,6 +184,49 @@ class FiniteVolume(pybamm.SpatialMethod):
         matrix = spdiags(data, diags, n - 1, n)
         return pybamm.Matrix(matrix)
 
+    def integral(self, symbol, discretised_symbol):
+        """Vector-vector dot product to implement the integral operator.
+        See :meth:`pybamm.BaseDiscretisation.integral`
+        """
+        # Calculate integration vector
+        integration_vector = self.definite_integral_vector(symbol.domain)
+        # Check for particle domain
+        if ("negative particle" or "positive particle") in symbol.domain:
+            submesh = self.mesh.combine_submeshes(*symbol.domain)
+            r = pybamm.Vector(submesh.nodes)
+            out = 2 * np.pi * integration_vector @ (discretised_symbol * r)
+        else:
+            out = integration_vector @ discretised_symbol
+        return out
+
+    def definite_integral_vector(self, domain):
+        """
+        Vector for finite-volume implementation of the definite integral
+
+        .. math::
+            I = \\int_{a}^{b}\\!f(s)\\,ds
+
+        for where :math:`a` and :math:`b` are the left-hand and right-hand boundaries of
+        the domain respectively
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) of integration
+
+        Returns
+        -------
+        :class:`pybamm.Vector`
+            The finite volume integral vector for the domain
+        """
+        # Create appropriate submesh by combining submeshes in domain
+        submesh = self.mesh.combine_submeshes(*domain)
+
+        # Create vector of ones using submesh
+        vector = submesh.d_edges * np.ones_like(submesh.nodes)
+
+        return pybamm.Vector(vector)
+
     def add_ghost_nodes(self, discretised_symbol, lbc, rbc):
         """
         Add Dirichlet boundary conditions via ghost nodes.
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 18927420..79a5e9bb 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -67,7 +67,7 @@ class SpatialMethod:
         ----------
         symbol: :class:`pybamm.Symbol`
             The symbol that we will take the gradient of.
-        broadcasted_symbol: class: pybamm.Array
+        discretised_symbol: class: pybamm.Array
             The discretised symbol of the correct size
 
         boundary_conditions : dict
@@ -82,7 +82,7 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
-    def divergence(self, symbol, broadcasted_symbol, boundary_conditions):
+    def divergence(self, symbol, discretised_symbol, boundary_conditions):
         """
         Implements the divergence for a spatial method.
 
@@ -90,12 +90,12 @@ class SpatialMethod:
         ----------
         symbol: :class:`pybamm.Symbol`
             The symbol that we will take the gradient of.
-        broadcasted_symbol: class: pybamm.Array
+        discretised_symbol: class: pybamm.Array
             The discretised symbol of the correct size
-
         boundary_conditions : dict
             The boundary conditions of the model
             ({symbol.id: {"left": left bc, "right": right bc}})
+
         Returns
         -------
         :class: `pybamm.Array`
@@ -104,6 +104,25 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
+    def integral(self, symbol, discretised_symbol):
+        """
+        Implements the integral for a spatial method.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol that we will take the gradient of.
+        discretised_symbol: class: pybamm.Array
+            The discretised symbol of the correct size
+
+        Returns
+        -------
+        :class: `pybamm.Array`
+            Contains the result of acting the discretised integral on
+            the child discretised_symbol
+        """
+        raise NotImplementedError
+
     def surface_value(self, discretised_symbol):
         """
         Returns the surface value using the approriate expression for the


Test patch:
diff --git a/tests/test_discretisations/test_discretisation.py b/tests/test_discretisations/test_discretisation.py
index a08c90d2..dfd82b0d 100644
--- a/tests/test_discretisations/test_discretisation.py
+++ b/tests/test_discretisations/test_discretisation.py
@@ -235,14 +235,15 @@ class TestDiscretise(unittest.TestCase):
             )
 
     def test_core_NotImplementedErrors(self):
-        # create discretisation
-        pybamm.Discretisation(None, {})
-
-        # TODO: implement and equivelent test in spatial_methods
-        # with self.assertRaises(NotImplementedError):
-        #     disc.gradient(None, None, {})
-        # with self.assertRaises(NotImplementedError):
-        #     disc.divergence(None, None, {})
+        # 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)
 
     def test_process_dict(self):
         # one equation
diff --git a/tests/test_expression_tree/test_unary_operators.py b/tests/test_expression_tree/test_unary_operators.py
index afe3bf9e..d685b2fc 100644
--- a/tests/test_expression_tree/test_unary_operators.py
+++ b/tests/test_expression_tree/test_unary_operators.py
@@ -65,6 +65,48 @@ class TestUnaryOperators(unittest.TestCase):
         grad = pybamm.Gradient(a)
         self.assertEqual(grad.children[0].name, a.name)
 
+    def test_integral(self):
+        # time integral
+        a = pybamm.Symbol("a")
+        t = pybamm.t
+        inta = pybamm.Integral(a, t)
+        self.assertEqual(inta.name, "integral dtime")
+        # self.assertTrue(inta.definite)
+        self.assertEqual(inta.children[0].name, a.name)
+        self.assertEqual(inta.integration_variable, t)
+        self.assertEqual(inta.domain, [])
+
+        # space integral
+        a = pybamm.Symbol("a", domain=["negative electrode"])
+        x = pybamm.SpatialVariable("x", ["negative electrode"])
+        inta = pybamm.Integral(a, x)
+        self.assertEqual(inta.name, "integral dx ['negative electrode']")
+        self.assertEqual(inta.children[0].name, a.name)
+        self.assertEqual(inta.integration_variable, x)
+        self.assertEqual(inta.domain, ["negative electrode"])
+
+        # # Indefinite
+        # for inta in [
+        #     pybamm.Integral(a, x, definite=False),
+        #     pybamm.IndefiniteIntegral(a, x),
+        # ]:
+        #     self.assertEqual(
+        #         inta.name, "indefinite integral dspace (['negative electrode'])"
+        #     )
+        #     self.assertFalse(inta.definite)
+        #     self.assertEqual(inta.children[0].name, a.name)
+        #     self.assertEqual(inta.integration_variable, x)
+        #     self.assertEqual(inta.domain, ["negative electrode"])
+
+        # expected errors
+        a = pybamm.Symbol("a", domain=["negative electrode"])
+        x = pybamm.SpatialVariable("x", ["separator"])
+        y = pybamm.Variable("y")
+        with self.assertRaises(pybamm.DomainError):
+            pybamm.Integral(a, x)
+        with self.assertRaises(ValueError):
+            pybamm.Integral(a, y)
+
     def test_printing(self):
         a = pybamm.Symbol("a")
         self.assertEqual(str(-a), "-a")
diff --git a/tests/test_spatial_methods/test_finite_volume.py b/tests/test_spatial_methods/test_finite_volume.py
index 637c1ff7..c8c24e09 100644
--- a/tests/test_spatial_methods/test_finite_volume.py
+++ b/tests/test_spatial_methods/test_finite_volume.py
@@ -402,6 +402,152 @@ class TestFiniteVolume(unittest.TestCase):
             div_eqn_disc.evaluate(None, linear_y), np.zeros_like(combined_submesh.nodes)
         )
 
+    def test_definite_integral(self):
+        # create discretisation
+        mesh = get_mesh_for_testing(200)
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+            "positive particle": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        # lengths
+        ln = mesh["negative electrode"].edges[-1]
+        ls = mesh["separator"].edges[-1] - ln
+        lp = 1 - (ln + ls)
+
+        # macroscale variable
+        var = pybamm.Variable("var", domain=["negative electrode", "separator"])
+        x = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
+        integral_eqn = pybamm.Integral(var, x)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+
+        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        constant_y = np.ones_like(combined_submesh.nodes)
+        self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), ln + ls)
+        linear_y = combined_submesh.nodes
+        self.assertAlmostEqual(
+            integral_eqn_disc.evaluate(None, linear_y), (ln + ls) ** 2 / 2
+        )
+        cos_y = np.cos(combined_submesh.nodes)
+        self.assertAlmostEqual(
+            integral_eqn_disc.evaluate(None, cos_y), np.sin(ln + ls), places=4
+        )
+
+        # domain not starting at zero
+        var = pybamm.Variable("var", domain=["separator", "positive electrode"])
+        x = pybamm.SpatialVariable("x", ["separator", "positive electrode"])
+        integral_eqn = pybamm.Integral(var, x)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+
+        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+        constant_y = np.ones_like(combined_submesh.nodes)
+        self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), ls + lp)
+        linear_y = combined_submesh.nodes
+        self.assertAlmostEqual(
+            integral_eqn_disc.evaluate(None, linear_y), (1 - (ln) ** 2) / 2
+        )
+        cos_y = np.cos(combined_submesh.nodes)
+        self.assertAlmostEqual(
+            integral_eqn_disc.evaluate(None, cos_y), np.sin(1) - np.sin(ln), places=4
+        )
+
+        # microscale variable
+        var = pybamm.Variable("var", domain=["negative particle"])
+        r = pybamm.SpatialVariable("r", ["negative particle"])
+        integral_eqn = pybamm.Integral(var, r)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+
+        constant_y = np.ones_like(mesh["negative particle"].nodes)
+        self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), np.pi)
+        linear_y = mesh["negative particle"].nodes
+        self.assertAlmostEqual(
+            integral_eqn_disc.evaluate(None, linear_y), 2 * np.pi / 3, places=4
+        )
+        one_over_y = 1 / mesh["negative particle"].nodes
+        self.assertEqual(integral_eqn_disc.evaluate(None, one_over_y), 2 * np.pi)
+
+    @unittest.skip("indefinite integral not yet implemented")
+    def test_indefinite_integral(self):
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+            "positive particle": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        # lengths
+        ln = mesh["negative electrode"].edges[-1]
+
+        # macroscale variable
+        var = pybamm.Variable("var", domain=["negative electrode", "separator"])
+        x = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
+        integral_eqn = pybamm.IndefiniteIntegral(var, x)
+        disc.set_variable_slices([var])
+        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
+        np.testing.assert_array_equal(
+            integral_eqn_disc.evaluate(None, constant_y_edges), linear_y
+        )
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc.evaluate(None, linear_y_edges), linear_y ** 2 / 2
+        )
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc.evaluate(None, np.cos(linear_y_edges)),
+            np.sin(linear_y),
+            places=4,
+        )
+
+        # domain not starting at zero
+        var = pybamm.Variable("var", domain=["separator", "positive electrode"])
+        x = pybamm.SpatialVariable("x", ["separator", "positive electrode"])
+        integral_eqn = pybamm.IndefiniteIntegral(var, x)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+
+        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+        np.testing.assert_array_equal(
+            integral_eqn_disc.evaluate(None, constant_y), linear_y - ln
+        )
+        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)
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc.evaluate(None, cos_y),
+            np.sin(linear_y) - np.sin(ln),
+            places=4,
+        )
+
+        # microscale variable
+        var = pybamm.Variable("var", domain=["negative particle"])
+        r = pybamm.SpatialVariable("r", ["negative particle"])
+        integral_eqn = pybamm.IndefiniteIntegral(var, r)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+
+        constant_y = np.ones_like(mesh["negative particle"].nodes)
+        np.testing.assert_array_equal(
+            integral_eqn_disc.evaluate(None, constant_y), np.pi
+        )
+        linear_y = mesh["negative particle"].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
+        np.testing.assert_array_equal(
+            integral_eqn_disc.evaluate(None, one_over_y), 2 * np.pi
+        )
+
     def test_grad_convergence_without_bcs(self):
         # Convergence
         whole_cell = ["negative electrode", "separator", "positive electrode"]