Repository: pybamm-team/PyBaMM (python)

Issue:
backward indefinite integral
Add new node `BackwardIndefiniteIntegral` for integrating from x to end.

The relationship
```
IndefiniteIntegral(var, x) + BackwardIndefiniteIntegral(var, x) = Integral(var, x)
```
should hold for all x

Pull request description:
Issue 1005 backward indefinite integral
# Description

Add a symbol for doing the indefinite integral from the end of the domain backwards.

Fixes #1005 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)


# 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


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0218499aa..58ce82fd8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added `BackwardIndefiniteIntegral` symbol ([#1014](https://github.com/pybamm-team/PyBaMM/pull/1014))
 -   Added SEI film resistance as an option ([#994](https://github.com/pybamm-team/PyBaMM/pull/994))
 -   Added tab, edge, and surface cooling ([#965](https://github.com/pybamm-team/PyBaMM/pull/965))
 -   Added functionality to solver to automatically discretise a 0D model ([#947](https://github.com/pybamm-team/PyBaMM/pull/947))
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index e9cae42dd..0d6d2f622 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -836,7 +836,13 @@ class Discretisation(object):
                 return child_spatial_method.boundary_mass_matrix(child, self.bcs)
 
             elif isinstance(symbol, pybamm.IndefiniteIntegral):
-                return child_spatial_method.indefinite_integral(child, disc_child)
+                return child_spatial_method.indefinite_integral(
+                    child, disc_child, "forward"
+                )
+            elif isinstance(symbol, pybamm.BackwardIndefiniteIntegral):
+                return child_spatial_method.indefinite_integral(
+                    child, disc_child, "backward"
+                )
 
             elif isinstance(symbol, pybamm.Integral):
                 out = child_spatial_method.integral(child, disc_child)
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index b1d785c39..78a08181a 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -510,14 +510,8 @@ class Integral(SpatialOperator):
         return False
 
 
-class IndefiniteIntegral(Integral):
-    """A node in the expression tree representing an indefinite integral operator
-
-    .. math::
-        I = \\int_{x_\text{min}}^{x}\\!f(u)\\,du
-
-    where :math:`u\\in\\text{domain}` which can represent either a
-    spatial or temporal variable.
+class BaseIndefiniteIntegral(Integral):
+    """Base class for indefinite integrals (forward or backward).
 
     Parameters
     ----------
@@ -540,15 +534,73 @@ class IndefiniteIntegral(Integral):
         super().__init__(child, integration_variable)
         # overwrite domains with child domains
         self.copy_domains(child)
+
+    def _evaluate_for_shape(self):
+        return self.children[0].evaluate_for_shape()
+
+    def evaluates_on_edges(self):
+        # If child evaluates on edges, indefinite integral doesn't
+        # If child doesn't evaluate on edges, indefinite integral does
+        return not self.child.evaluates_on_edges()
+
+
+class IndefiniteIntegral(BaseIndefiniteIntegral):
+    """A node in the expression tree representing an indefinite integral operator
+
+    .. math::
+        I = \\int_{x_\text{min}}^{x}\\!f(u)\\,du
+
+    where :math:`u\\in\\text{domain}` which can represent either a
+    spatial or temporal variable.
+
+    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:`BaseIndefiniteIntegral`
+    """
+
+    def __init__(self, child, integration_variable):
+        super().__init__(child, integration_variable)
         # Overwrite the name
         self.name = "{} integrated w.r.t {}".format(
-            child.name, integration_variable.name
+            child.name, self.integration_variable[0].name
         )
         if isinstance(integration_variable, pybamm.SpatialVariable):
-            self.name += " on {}".format(integration_variable.domain)
+            self.name += " on {}".format(self.integration_variable[0].domain)
 
-    def _evaluate_for_shape(self):
-        return self.children[0].evaluate_for_shape()
+
+class BackwardIndefiniteIntegral(BaseIndefiniteIntegral):
+    """A node in the expression tree representing a backward indefinite integral
+    operator
+
+    .. math::
+        I = \\int_{x}^{x_\text{max}}\\!f(u)\\,du
+
+    where :math:`u\\in\\text{domain}` which can represent either a
+    spatial or temporal variable.
+
+    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:`BaseIndefiniteIntegral`
+    """
+
+    def __init__(self, child, integration_variable):
+        super().__init__(child, integration_variable)
+        # Overwrite the name
+        self.name = "{} integrated backward w.r.t {}".format(
+            child.name, self.integration_variable[0].name
+        )
+        if isinstance(integration_variable, pybamm.SpatialVariable):
+            self.name += " on {}".format(self.integration_variable[0].domain)
 
 
 class DefiniteIntegralVector(SpatialOperator):
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index d748f5590..671b388c0 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -5,6 +5,7 @@ import pybamm
 
 from scipy.sparse import (
     diags,
+    spdiags,
     eye,
     kron,
     csr_matrix,
@@ -248,6 +249,8 @@ class FiniteVolume(pybamm.SpatialMethod):
         else:
             out = integration_vector @ discretised_child
 
+        out.copy_domains(child)
+
         return out
 
     def definite_integral_matrix(self, domain, vector_type="row"):
@@ -296,15 +299,19 @@ class FiniteVolume(pybamm.SpatialMethod):
         matrix = csr_matrix(kron(eye(second_dim_len), vector))
         return pybamm.Matrix(matrix)
 
-    def indefinite_integral(self, child, discretised_child):
+    def indefinite_integral(self, child, discretised_child, direction):
         """Implementation of the indefinite integral operator. """
 
         # Different integral matrix depending on whether the integrand evaluates on
         # edges or nodes
         if child.evaluates_on_edges():
-            integration_matrix = self.indefinite_integral_matrix_edges(child.domain)
+            integration_matrix = self.indefinite_integral_matrix_edges(
+                child.domain, direction
+            )
         else:
-            integration_matrix = self.indefinite_integral_matrix_nodes(child.domain)
+            integration_matrix = self.indefinite_integral_matrix_nodes(
+                child.domain, direction
+            )
 
         # Don't need to check for spherical domains as spherical polars
         # only change the diveregence (childs here have grad and no div)
@@ -314,10 +321,28 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return out
 
-    def indefinite_integral_matrix_edges(self, domain):
+    def indefinite_integral_matrix_edges(self, domain, direction):
         """
         Matrix for finite-volume implementation of the indefinite integral where the
-        integrand is evaluated on mesh edges
+        integrand is evaluated on mesh edges (shape (n+1, 1)).
+        The integral will then be evaluated on mesh nodes (shape (n, 1)).
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) of integration
+        direction : str
+            The direction of integration (forward or backward). See notes.
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The finite volume integral matrix for the domain
+
+        Notes
+        -----
+
+        **Forward integral**
 
         .. math::
             F(x) = \\int_0^x\\!f(u)\\,du
@@ -335,16 +360,81 @@ class FiniteVolume(pybamm.SpatialMethod):
         Hence we must have
 
         - :math:`F_0 = du_{1/2} * f_{1/2} / 2`
-        - :math:`F_{i+1} = F_i + du * f_{i+1/2}`
+        - :math:`F_{i+1} = F_i + du_{i+1/2} * f_{i+1/2}`
+
+        Note that :math:`f_{-1/2}` and :math:`f_{end+1/2}` are included in the discrete
+        integrand vector `f`, so we add a column of zeros at each end of the
+        indefinite integral matrix to ignore these.
+
+        **Backward integral**
+
+        .. math::
+            F(x) = \\int_x^end\\!f(u)\\,du
 
-        Note that :math:`f_{-1/2}` and :math:`f_{n+1/2}` are included in the discrete
+        The indefinite integral must satisfy the following conditions:
+
+        - :math:`F(end) = 0`
+        - :math:`f(x) = -\\frac{dF}{dx}`
+
+        or, in discrete form,
+
+        - `BoundaryValue(F, "right") = 0`, i.e. :math:`3*F_{end} - F_{end-1} = 0`
+        - :math:`f_{i+1/2} = -(F_{i+1} - F_i) / dx_{i+1/2}`
+
+        Hence we must have
+
+        - :math:`F_{end} = du_{end+1/2} * f_{end-1/2} / 2`
+        - :math:`F_{i-1} = F_i + du_{i-1/2} * f_{i-1/2}`
+
+        Note that :math:`f_{-1/2}` and :math:`f_{end+1/2}` are included in the discrete
         integrand vector `f`, so we add a column of zeros at each end of the
         indefinite integral matrix to ignore these.
+        """
+
+        # Create appropriate submesh by combining submeshes in domain
+        submesh_list = self.mesh.combine_submeshes(*domain)
+        submesh = submesh_list[0]
+        n = submesh.npts
+        sec_pts = len(submesh_list)
+
+        du_n = submesh.d_nodes
+        if direction == "forward":
+            du_entries = [du_n] * (n - 1)
+            offset = -np.arange(1, n, 1)
+            main_integral_matrix = spdiags(du_entries, offset, n, n - 1)
+            bc_offset_matrix = lil_matrix((n, n - 1))
+            bc_offset_matrix[:, 0] = du_n[0] / 2
+        elif direction == "backward":
+            du_entries = [du_n] * (n + 1)
+            offset = np.arange(n, -1, -1)
+            main_integral_matrix = spdiags(du_entries, offset, n, n - 1)
+            bc_offset_matrix = lil_matrix((n, n - 1))
+            bc_offset_matrix[:, -1] = du_n[-1] / 2
+        sub_matrix = main_integral_matrix + bc_offset_matrix
+        # add a column of zeros at each end
+        zero_col = csr_matrix((n, 1))
+        sub_matrix = hstack([zero_col, sub_matrix, zero_col])
+        # 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 = csr_matrix(kron(eye(sec_pts), sub_matrix))
+
+        return pybamm.Matrix(matrix)
+
+    def indefinite_integral_matrix_nodes(self, domain, direction):
+        """
+        Matrix for finite-volume implementation of the (backward) indefinite integral
+        where the integrand is evaluated on mesh nodes (shape (n, 1)).
+        The integral will then be evaluated on mesh edges (shape (n+1, 1)).
+        This is just a straightforward (backward) cumulative sum of the integrand
 
         Parameters
         ----------
         domain : list
             The domain(s) of integration
+        direction : str
+            The direction of integration (forward or backward)
 
         Returns
         -------
@@ -358,16 +448,13 @@ class FiniteVolume(pybamm.SpatialMethod):
         n = submesh.npts
         sec_pts = len(submesh_list)
 
-        du_n = submesh.d_nodes
-        du_entries = [du_n] * (n - 1)
-        offset = -np.arange(1, n, 1)
-        main_integral_matrix = diags(du_entries, offset, shape=(n, n - 1))
-        bc_offset_matrix = lil_matrix((n, n - 1))
-        bc_offset_matrix[:, 0] = du_n[0] / 2
-        sub_matrix = main_integral_matrix + bc_offset_matrix
-        # add a column of zeros at each end
-        zero_col = csr_matrix((n, 1))
-        sub_matrix = hstack([zero_col, sub_matrix, zero_col])
+        du_n = submesh.d_edges
+        du_entries = [du_n] * n
+        if direction == "forward":
+            offset = -np.arange(1, n + 1, 1)  # from -1 down to -n
+        elif direction == "backward":
+            offset = np.arange(n - 1, -1, -1)  # from n-1 down to 0
+        sub_matrix = spdiags(du_entries, offset, n + 1, n)
         # 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
@@ -473,41 +560,6 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return dy / dx
 
-    def indefinite_integral_matrix_nodes(self, domain):
-        """
-        Matrix for finite-volume implementation of the indefinite integral where the
-        integrand is evaluated on mesh nodes.
-        This is just a straightforward cumulative sum of the integrand
-
-        Parameters
-        ----------
-        domain : list
-            The domain(s) of integration
-
-        Returns
-        -------
-        :class:`pybamm.Matrix`
-            The finite volume integral matrix for the domain
-        """
-
-        # Create appropriate submesh by combining submeshes in domain
-        submesh_list = self.mesh.combine_submeshes(*domain)
-        submesh = submesh_list[0]
-        n = submesh.npts
-        sec_pts = len(submesh_list)
-
-        du_n = submesh.d_edges
-        du_entries = [du_n] * (n)
-        offset = -np.arange(1, n + 1, 1)
-        sub_matrix = diags(du_entries, offset, shape=(n + 1, n))
-        # 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 = csr_matrix(kron(eye(sec_pts), sub_matrix))
-
-        return pybamm.Matrix(matrix)
-
     def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
         """
         Add ghost nodes to a symbol.
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 371d50ebc..41d26539b 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -242,7 +242,7 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
-    def indefinite_integral(self, child, discretised_child):
+    def indefinite_integral(self, child, discretised_child, direction):
         """
         Implements the indefinite integral for a spatial method.
 
@@ -252,6 +252,8 @@ class SpatialMethod:
             The symbol to which is being integrated
         discretised_child: :class:`pybamm.Symbol`
             The discretised symbol of the correct size
+        direction : str
+            The direction of integration
 
         Returns
         -------
diff --git a/pybamm/spatial_methods/zero_dimensional_method.py b/pybamm/spatial_methods/zero_dimensional_method.py
index aada27398..845270099 100644
--- a/pybamm/spatial_methods/zero_dimensional_method.py
+++ b/pybamm/spatial_methods/zero_dimensional_method.py
@@ -37,11 +37,16 @@ class ZeroDimensionalSpatialMethod(pybamm.SpatialMethod):
         """
         return pybamm.Matrix(np.ones((1, 1)))
 
-    def indefinite_integral(self, child, discretised_child):
+    def indefinite_integral(self, child, discretised_child, direction):
         """
-        Calculates the zero-dimensional indefinite integral, i.e. the identity operator
+        Calculates the zero-dimensional indefinite integral.
+        If 'direction' is forward, this is the identity operator.
+        If 'direction' is backward, this is the negation operator.
         """
-        return discretised_child
+        if direction == "forward":
+            return discretised_child
+        elif direction == "backward":
+            return -discretised_child
 
     def integral(self, child, discretised_child):
         """


Test patch:
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index a72faf649..e7ddbf305 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -12,6 +12,7 @@ class TestCopy(unittest.TestCase):
         a = pybamm.Scalar(0)
         b = pybamm.Scalar(1)
         v_n = pybamm.Variable("v", "negative electrode")
+        x_n = pybamm.standard_spatial_vars.x_n
         v_s = pybamm.Variable("v", "separator")
         vec = pybamm.Vector(np.array([1, 2, 3, 4, 5]))
         mesh = get_mesh_for_testing()
@@ -29,6 +30,8 @@ class TestCopy(unittest.TestCase):
             pybamm.grad(v_n),
             pybamm.div(pybamm.grad(v_n)),
             pybamm.Integral(a, pybamm.t),
+            pybamm.IndefiniteIntegral(v_n, x_n),
+            pybamm.BackwardIndefiniteIntegral(v_n, x_n),
             pybamm.BoundaryValue(v_n, "right"),
             pybamm.BoundaryGradient(v_n, "right"),
             pybamm.PrimaryBroadcast(a, "domain"),
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 28d206b92..030deaa73 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -173,6 +173,11 @@ class TestUnaryOperators(unittest.TestCase):
         self.assertEqual(
             inta_sec.auxiliary_domains, {"secondary": ["current collector"]}
         )
+        # backward indefinite integral
+        inta = pybamm.BackwardIndefiniteIntegral(a, x)
+        self.assertEqual(
+            inta.name, "a integrated backward w.r.t x on ['negative electrode']"
+        )
 
         # expected errors
         a = pybamm.Symbol("a", domain=["negative electrode"])
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 475c3aa91..211b4586c 100644
--- a/tests/unit/test_spatial_methods/test_base_spatial_method.py
+++ b/tests/unit/test_spatial_methods/test_base_spatial_method.py
@@ -24,7 +24,7 @@ class TestSpatialMethod(unittest.TestCase):
         with self.assertRaises(NotImplementedError):
             spatial_method.integral(None, None)
         with self.assertRaises(NotImplementedError):
-            spatial_method.indefinite_integral(None, None)
+            spatial_method.indefinite_integral(None, None, None)
         with self.assertRaises(NotImplementedError):
             spatial_method.boundary_integral(None, None, None)
         with self.assertRaises(NotImplementedError):
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 8c261c1f9..c023029b7 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -872,6 +872,55 @@ class TestFiniteVolume(unittest.TestCase):
             left_boundary_value_disc.evaluate(y=c_exact), 0
         )
 
+    def test_backward_indefinite_integral(self):
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # --------------------------------------------------------------------
+        # region which doesn't start at zero
+        phi = pybamm.Variable("phi", domain=["negative electrode", "separator"])
+        i = pybamm.grad(phi)  # create test current (variable on edges)
+        x = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
+        int_grad_phi = pybamm.BackwardIndefiniteIntegral(i, x)
+        disc.set_variable_slices([phi])  # i is not a fundamental variable
+        disc._bcs = {
+            phi.id: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        int_grad_phi_disc = disc.process_symbol(int_grad_phi)
+        right_boundary_value = pybamm.BoundaryValue(int_grad_phi, "right")
+        right_boundary_value_disc = disc.process_symbol(right_boundary_value)
+        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+
+        # Test that the backward_integral(grad(phi)) = -phi
+        # constant case
+        phi_exact = np.ones((combined_submesh[0].npts, 1))
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
+        phi_approx += 1  # add constant of integration
+        np.testing.assert_array_equal(phi_exact, phi_approx)
+        self.assertEqual(right_boundary_value_disc.evaluate(y=phi_exact), 0)
+
+        # linear case
+        phi_exact = combined_submesh[0].nodes - combined_submesh[0].edges[-1]
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(phi_exact, -phi_approx)
+        np.testing.assert_array_almost_equal(
+            right_boundary_value_disc.evaluate(y=phi_exact), 0
+        )
+
+        # sine case
+        phi_exact = np.sin(combined_submesh[0].nodes - combined_submesh[0].edges[-1])
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(phi_exact, -phi_approx)
+        np.testing.assert_array_almost_equal(
+            right_boundary_value_disc.evaluate(y=phi_exact), 0
+        )
+
     def test_indefinite_integral_of_broadcasted_to_cell_edges(self):
         # create discretisation
         mesh = get_mesh_for_testing()
@@ -925,19 +974,86 @@ class TestFiniteVolume(unittest.TestCase):
         # constant case
         phi_exact = np.ones((combined_submesh[0].npts, 1))
         int_phi_exact = combined_submesh[0].edges
-        int_phi_approx = int_phi_disc.evaluate(None, phi_exact)[:, 0]
+        int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_equal(int_phi_exact, int_phi_approx)
         # linear case
-        phi_exact = combined_submesh[0].nodes[:, np.newaxis]
-        int_phi_exact = combined_submesh[0].edges[:, np.newaxis] ** 2 / 2
-        int_phi_approx = int_phi_disc.evaluate(None, phi_exact)
+        phi_exact = combined_submesh[0].nodes
+        int_phi_exact = combined_submesh[0].edges ** 2 / 2
+        int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(int_phi_exact, int_phi_approx)
         # cos case
-        phi_exact = np.cos(combined_submesh[0].nodes[:, np.newaxis])
-        int_phi_exact = np.sin(combined_submesh[0].edges[:, np.newaxis])
-        int_phi_approx = int_phi_disc.evaluate(None, phi_exact)
+        phi_exact = np.cos(combined_submesh[0].nodes)
+        int_phi_exact = np.sin(combined_submesh[0].edges)
+        int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(int_phi_exact, int_phi_approx, decimal=5)
 
+    def test_backward_indefinite_integral_on_nodes(self):
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        phi = pybamm.Variable("phi", domain=["negative electrode", "separator"])
+        x = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
+
+        back_int_phi = pybamm.BackwardIndefiniteIntegral(phi, x)
+        disc.set_variable_slices([phi])
+        back_int_phi_disc = disc.process_symbol(back_int_phi)
+
+        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        edges = combined_submesh[0].edges
+
+        # constant case
+        phi_exact = np.ones((combined_submesh[0].npts, 1))
+        back_int_phi_exact = edges[-1] - edges
+        back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(back_int_phi_exact, back_int_phi_approx)
+        # linear case
+        phi_exact = combined_submesh[0].nodes
+        back_int_phi_exact = edges[-1] ** 2 / 2 - edges ** 2 / 2
+        back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(back_int_phi_exact, back_int_phi_approx)
+        # cos case
+        phi_exact = np.cos(combined_submesh[0].nodes)
+        back_int_phi_exact = np.sin(edges[-1]) - np.sin(edges)
+        back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(
+            back_int_phi_exact, back_int_phi_approx, decimal=5
+        )
+
+    def test_forward_plus_backward_integral(self):
+        # Test that forward integral + backward integral = integral
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # On nodes
+        phi = pybamm.Variable("phi", domain=["separator", "positive electrode"])
+        x = pybamm.SpatialVariable("x", ["separator", "positive electrode"])
+
+        disc.set_variable_slices([phi])
+
+        full_int_phi = pybamm.PrimaryBroadcastToEdges(
+            pybamm.Integral(phi, x), ["separator", "positive electrode"]
+        )
+        full_int_phi_disc = disc.process_symbol(full_int_phi)
+        int_plus_back_int_phi = pybamm.IndefiniteIntegral(
+            phi, x
+        ) + pybamm.BackwardIndefiniteIntegral(phi, x)
+        int_plus_back_int_phi_disc = disc.process_symbol(int_plus_back_int_phi)
+
+        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+
+        # test
+        for phi_exact in [
+            np.ones((combined_submesh[0].npts, 1)),
+            combined_submesh[0].nodes,
+            np.cos(combined_submesh[0].nodes),
+        ]:
+            np.testing.assert_array_almost_equal(
+                full_int_phi_disc.evaluate(y=phi_exact).flatten(),
+                int_plus_back_int_phi_disc.evaluate(y=phi_exact).flatten(),
+            )
+
     def test_discretise_spatial_variable(self):
         # create discretisation
         mesh = get_mesh_for_testing()
diff --git a/tests/unit/test_spatial_methods/test_zero_dimensional_method.py b/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
index f3ffdd80f..de893d981 100644
--- a/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
+++ b/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
@@ -16,8 +16,11 @@ class TestZeroDimensionalSpatialMethod(unittest.TestCase):
 
         a = pybamm.Symbol("a")
         self.assertEqual(a, spatial_method.integral(None, a))
-        self.assertEqual(a, spatial_method.indefinite_integral(None, a))
+        self.assertEqual(a, spatial_method.indefinite_integral(None, a, "forward"))
         self.assertEqual(a, spatial_method.boundary_value_or_flux(None, a))
+        self.assertEqual(
+            (-a).id, spatial_method.indefinite_integral(None, a, "backward").id
+        )
 
         mass_matrix = spatial_method.mass_matrix(None, None)
         self.assertIsInstance(mass_matrix, pybamm.Matrix)