From d226492b5878789fc7fd6bc569ea7a33d7cda61e Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 20 May 2020 12:01:15 -0400
Subject: [PATCH 1/8] #1005 add backward indefinite integral operator

---
 pybamm/discretisations/discretisation.py      |   8 +-
 pybamm/expression_tree/unary_operators.py     |  72 +++++++--
 pybamm/spatial_methods/finite_volume.py       | 140 ++++++++++++------
 .../test_unary_operators.py                   |   5 +
 .../test_finite_volume/test_finite_volume.py  |  52 +++++++
 5 files changed, 218 insertions(+), 59 deletions(-)

diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index e9cae42dd8..0d6d2f6229 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -836,7 +836,13 @@ def _process_symbol(self, symbol):
                 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 b1d785c39a..3450d5cc3f 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -510,14 +510,8 @@ def evaluates_on_edges(self):
         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,6 +534,37 @@ def __init__(self, child, integration_variable):
         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
@@ -547,8 +572,35 @@ def __init__(self, child, integration_variable):
         if isinstance(integration_variable, pybamm.SpatialVariable):
             self.name += " on {}".format(integration_variable.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, integration_variable.name
+        )
+        if isinstance(integration_variable, pybamm.SpatialVariable):
+            self.name += " on {}".format(integration_variable.domain)
 
 
 class DefiniteIntegralVector(SpatialOperator):
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index d748f55909..41bba884a8 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -248,6 +248,8 @@ def integral(self, child, discretised_child):
         else:
             out = integration_vector @ discretised_child
 
+        out.copy_domains(child)
+
         return out
 
     def definite_integral_matrix(self, domain, vector_type="row"):
@@ -296,15 +298,19 @@ def definite_integral_matrix(self, domain, vector_type="row"):
         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 +320,28 @@ def indefinite_integral(self, child, discretised_child):
 
         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.
+
+        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
@@ -341,15 +365,30 @@ def indefinite_integral_matrix_edges(self, domain):
         integrand vector `f`, so we add a column of zeros at each end of the
         indefinite integral matrix to ignore these.
 
-        Parameters
-        ----------
-        domain : list
-            The domain(s) of integration
+        Backward integral
+        ^^^^^^^^^^^^^^^^^
 
-        Returns
-        -------
-        :class:`pybamm.Matrix`
-            The finite volume integral matrix for the domain
+        .. math::
+            F(x) = \\int_x^end\\!f(u)\\,du
+
+        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_{1/2} * f_{end-1/2} / 2`
+        - :math:`F_{i+1} = F_i + du * f_{i+1/2}`
+
+        Note that :math:`f_{-1/2}` and :math:`f_{n+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
@@ -376,6 +415,46 @@ def indefinite_integral_matrix_edges(self, domain):
 
         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.
+        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
+        -------
+        :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)
+        if direction == "forward":
+            offset = -np.arange(1, n + 1, 1)
+        elif direction == "backward":
+            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 delta_function(self, symbol, discretised_symbol):
         """
         Delta function. Implemented as a vector whose only non-zero element is the
@@ -473,41 +552,6 @@ def internal_neumann_condition(
 
         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/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 28d206b92e..030deaa73d 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 @@ def test_integral(self):
         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_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 8c261c1f96..ec5a6eeb00 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,58 @@ def test_indefinite_integral(self):
             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=["separator", "positive electrode"])
+        i = pybamm.grad(phi)  # create test current (variable on edges)
+        x = pybamm.SpatialVariable("x", ["separator", "positive electrode"])
+        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)
+        left_boundary_value = pybamm.BoundaryValue(int_grad_phi, "left")
+        left_boundary_value_disc = disc.process_symbol(left_boundary_value)
+        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+
+        # 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(left_boundary_value_disc.evaluate(y=phi_exact), 0)
+
+        # linear case
+        phi_exact = (
+            combined_submesh[0].nodes[:, np.newaxis] - combined_submesh[0].edges[0]
+        )
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
+        np.testing.assert_array_almost_equal(phi_exact, phi_approx)
+        np.testing.assert_array_almost_equal(
+            left_boundary_value_disc.evaluate(y=phi_exact), 0
+        )
+
+        # sine case
+        phi_exact = np.sin(
+            combined_submesh[0].nodes[:, np.newaxis] - combined_submesh[0].edges[0]
+        )
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
+        np.testing.assert_array_almost_equal(phi_exact, phi_approx)
+        np.testing.assert_array_almost_equal(
+            left_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()

From bd360061f3a8f2b518b261be488798c372c14c9d Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 20 May 2020 16:29:27 -0400
Subject: [PATCH 2/8] #1005 get backward integral on nodes working

---
 pybamm/spatial_methods/finite_volume.py       | 13 +--
 .../test_finite_volume/test_finite_volume.py  | 80 +++++++++++++++++--
 2 files changed, 81 insertions(+), 12 deletions(-)

diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 41bba884a8..aecd32e080 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -5,6 +5,7 @@
 
 from scipy.sparse import (
     diags,
+    spdiags,
     eye,
     kron,
     csr_matrix,
@@ -323,7 +324,8 @@ def indefinite_integral(self, child, discretised_child, direction):
     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. The integral will then be evaluated on
+        mesh nodes.
 
         Parameters
         ----------
@@ -418,7 +420,8 @@ def indefinite_integral_matrix_edges(self, domain, direction):
     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.
+        where the integrand is evaluated on mesh nodes. The integral will then be
+        evaluated on mesh edges.
         This is just a straightforward (backward) cumulative sum of the integrand
 
         Parameters
@@ -441,12 +444,12 @@ def indefinite_integral_matrix_nodes(self, domain, direction):
         sec_pts = len(submesh_list)
 
         du_n = submesh.d_edges
-        du_entries = [du_n] * (n)
+        du_entries = [du_n] * n
         if direction == "forward":
             offset = -np.arange(1, n + 1, 1)
         elif direction == "backward":
-            offset = -np.arange(1, n + 1, 1)
-        sub_matrix = diags(du_entries, offset, shape=(n + 1, n))
+            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
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 ec5a6eeb00..3540c08343 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
@@ -977,19 +977,85 @@ def test_indefinite_integral_on_nodes(self):
         # 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)
+
+        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()

From efafa95c5e9d24748d97da3f88a3c8f6fd1d0da0 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 20 May 2020 17:34:49 -0400
Subject: [PATCH 3/8] #1005 get backward integral on edges working

---
 pybamm/spatial_methods/finite_volume.py       | 41 +++++++++++--------
 .../test_finite_volume/test_finite_volume.py  | 34 ++++++++-------
 2 files changed, 40 insertions(+), 35 deletions(-)

diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index aecd32e080..bba525d420 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -324,8 +324,8 @@ def indefinite_integral(self, child, discretised_child, direction):
     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. The integral will then be evaluated on
-        mesh nodes.
+        integrand is evaluated on mesh edges (shape (n+1, 1)).
+        The integral will then be evaluated on mesh nodes (shape (n, 1)).
 
         Parameters
         ----------
@@ -361,9 +361,9 @@ def indefinite_integral_matrix_edges(self, domain, direction):
         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_{n+1/2}` are included in the discrete
+        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.
 
@@ -376,19 +376,19 @@ def indefinite_integral_matrix_edges(self, domain, direction):
         The indefinite integral must satisfy the following conditions:
 
         - :math:`F(end) = 0`
-        - :math:`f(x) = \\frac{dF}{dx}`
+        - :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}`
+        - :math:`f_{i+1/2} = -(F_{i+1} - F_i) / dx_{i+1/2}`
 
         Hence we must have
 
-        - :math:`F_end = du_{1/2} * f_{end-1/2} / 2`
-        - :math:`F_{i+1} = F_i + du * f_{i+1/2}`
+        - :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_{n+1/2}` are included in the discrete
+        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.
         """
@@ -400,11 +400,18 @@ def indefinite_integral_matrix_edges(self, domain, direction):
         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
+        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))
@@ -420,8 +427,8 @@ def indefinite_integral_matrix_edges(self, domain, direction):
     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. The integral will then be
-        evaluated on mesh edges.
+        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
@@ -446,7 +453,7 @@ def indefinite_integral_matrix_nodes(self, domain, direction):
         du_n = submesh.d_edges
         du_entries = [du_n] * n
         if direction == "forward":
-            offset = -np.arange(1, n + 1, 1)
+            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)
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 3540c08343..c023029b71 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
@@ -881,9 +881,9 @@ def test_backward_indefinite_integral(self):
 
         # --------------------------------------------------------------------
         # region which doesn't start at zero
-        phi = pybamm.Variable("phi", domain=["separator", "positive electrode"])
+        phi = pybamm.Variable("phi", domain=["negative electrode", "separator"])
         i = pybamm.grad(phi)  # create test current (variable on edges)
-        x = pybamm.SpatialVariable("x", ["separator", "positive electrode"])
+        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 = {
@@ -893,35 +893,32 @@ def test_backward_indefinite_integral(self):
             }
         }
         int_grad_phi_disc = disc.process_symbol(int_grad_phi)
-        left_boundary_value = pybamm.BoundaryValue(int_grad_phi, "left")
-        left_boundary_value_disc = disc.process_symbol(left_boundary_value)
-        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+        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(left_boundary_value_disc.evaluate(y=phi_exact), 0)
+        self.assertEqual(right_boundary_value_disc.evaluate(y=phi_exact), 0)
 
         # linear case
-        phi_exact = (
-            combined_submesh[0].nodes[:, np.newaxis] - combined_submesh[0].edges[0]
-        )
-        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
-        np.testing.assert_array_almost_equal(phi_exact, phi_approx)
+        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(
-            left_boundary_value_disc.evaluate(y=phi_exact), 0
+            right_boundary_value_disc.evaluate(y=phi_exact), 0
         )
 
         # sine case
-        phi_exact = np.sin(
-            combined_submesh[0].nodes[:, np.newaxis] - combined_submesh[0].edges[0]
-        )
-        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
-        np.testing.assert_array_almost_equal(phi_exact, phi_approx)
+        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(
-            left_boundary_value_disc.evaluate(y=phi_exact), 0
+            right_boundary_value_disc.evaluate(y=phi_exact), 0
         )
 
     def test_indefinite_integral_of_broadcasted_to_cell_edges(self):
@@ -1029,6 +1026,7 @@ def test_forward_plus_backward_integral(self):
         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"])
 

From 05b3f6ce91b68ad3b4cba4e367616ba25bbc1bcd Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 20 May 2020 18:15:22 -0400
Subject: [PATCH 4/8] #1005 fix integration variable

---
 pybamm/expression_tree/unary_operators.py | 8 ++++----
 pybamm/spatial_methods/finite_volume.py   | 6 ++----
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 3450d5cc3f..78a08181a0 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -567,10 +567,10 @@ 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)
 
 
 class BackwardIndefiniteIntegral(BaseIndefiniteIntegral):
@@ -597,10 +597,10 @@ def __init__(self, child, integration_variable):
         super().__init__(child, integration_variable)
         # Overwrite the name
         self.name = "{} integrated backward 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)
 
 
 class DefiniteIntegralVector(SpatialOperator):
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index bba525d420..671b388c0b 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -342,8 +342,7 @@ def indefinite_integral_matrix_edges(self, domain, direction):
         Notes
         -----
 
-        Forward integral
-        ^^^^^^^^^^^^^^^^
+        **Forward integral**
 
         .. math::
             F(x) = \\int_0^x\\!f(u)\\,du
@@ -367,8 +366,7 @@ def indefinite_integral_matrix_edges(self, domain, direction):
         integrand vector `f`, so we add a column of zeros at each end of the
         indefinite integral matrix to ignore these.
 
-        Backward integral
-        ^^^^^^^^^^^^^^^^^
+        **Backward integral**
 
         .. math::
             F(x) = \\int_x^end\\!f(u)\\,du

From b47932220a059d35b4a593aeded6170b4e42be1f Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 21 May 2020 12:33:42 -0400
Subject: [PATCH 5/8] #1005 add integrals to copy test

---
 tests/unit/test_expression_tree/test_operations/test_copy.py | 3 +++
 1 file changed, 3 insertions(+)

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 a72faf6495..e7ddbf3057 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 @@ def test_symbol_new_copy(self):
         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 @@ def test_symbol_new_copy(self):
             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"),

From c952bfdc5bdcea253a1afd513bfa8bc4f8a7028a Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 21 May 2020 12:35:21 -0400
Subject: [PATCH 6/8] #1005 changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0218499aaa..58ce82fd8e 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))

From 863da6336484bf830a1b40f8fe743e346ce41508 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 21 May 2020 13:58:38 -0400
Subject: [PATCH 7/8] #1005 fix 0D indefinite integral

---
 pybamm/spatial_methods/spatial_method.py              |  4 +++-
 pybamm/spatial_methods/zero_dimensional_method.py     | 11 ++++++++---
 .../test_zero_dimensional_method.py                   |  5 ++++-
 3 files changed, 15 insertions(+), 5 deletions(-)

diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 371d50ebc4..41d26539bd 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -242,7 +242,7 @@ def integral(self, child, discretised_child):
         """
         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 @@ def indefinite_integral(self, child, discretised_child):
             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 aada273987..8452700993 100644
--- a/pybamm/spatial_methods/zero_dimensional_method.py
+++ b/pybamm/spatial_methods/zero_dimensional_method.py
@@ -37,11 +37,16 @@ def mass_matrix(self, symbol, boundary_conditions):
         """
         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):
         """
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 f3ffdd80f3..de893d9812 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 @@ def test_identity_ops(self):
 
         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)

From 14bf307d0d0643bf2c75ca096ee3fc62a082c1ac Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 21 May 2020 14:20:06 -0400
Subject: [PATCH 8/8] #1005 fix another test

---
 tests/unit/test_spatial_methods/test_base_spatial_method.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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 475c3aa91f..211b4586cd 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 @@ def test_basics(self):
         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):