Merge pull request #1014 from pybamm-team/issue-1005-backward-indefinite-integral

Issue 1005 backward indefinite integral

Author: valentinsulzer

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))

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)

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,15 +534,73 @@ 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
+            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):

pybamm/spatial_methods/finite_volume.py

@@ -5,6 +5,7 @@
 
 from scipy.sparse import (
     diags,
+    spdiags,
     eye,
     kron,
     csr_matrix,
@@ -248,6 +249,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 +299,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 +321,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 (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 @@ def indefinite_integral_matrix_edges(self, domain):
         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 @@ def indefinite_integral_matrix_edges(self, domain):
         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 @@ 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.

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
         -------

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):
         """

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"),