#632 working on finite volume discretisation

pybamm/discretisations/discretisation.py

@@ -845,7 +845,12 @@ def _process_symbol(self, symbol):
                 )
 
             elif isinstance(symbol, pybamm.Integral):
-                out = child_spatial_method.integral(child, disc_child)
+                integral_spatial_method = self.spatial_methods[
+                    symbol.integration_variable[0].domain[0]
+                ]
+                out = integral_spatial_method.integral(
+                    child, disc_child, symbol._integration_dimension
+                )
                 out.copy_domains(symbol)
                 return out
 

pybamm/expression_tree/unary_operators.py

@@ -441,17 +441,17 @@ def __init__(self, child, integration_variable):
             if isinstance(var, pybamm.SpatialVariable):
                 # Check that child and integration_variable domains agree
                 if var.domain == child.domain:
-                    self._integration_domain = "primary"
+                    self._integration_dimension = "primary"
                 elif (
                     "secondary" in child.auxiliary_domains
                     and var.domain == child.auxiliary_domains["secondary"]
                 ):
-                    self._integration_domain = "secondary"
+                    self._integration_dimension = "secondary"
                 elif (
                     "tertiary" in child.auxiliary_domains
                     and var.domain == child.auxiliary_domains["tertiary"]
                 ):
-                    self._integration_domain = "tertiary"
+                    self._integration_dimension = "tertiary"
                 else:
                     raise pybamm.DomainError(
                         "integration_variable must be the same as child domain or "
@@ -464,7 +464,7 @@ def __init__(self, child, integration_variable):
                 )
             name += " d{}".format(var.name)
 
-        if self._integration_domain == "primary":
+        if self._integration_dimension == "primary":
             # integral of a child takes the domain from auxiliary domain of the child
             if child.auxiliary_domains != {}:
                 domain = child.auxiliary_domains["secondary"]
@@ -478,15 +478,15 @@ def __init__(self, child, integration_variable):
             else:
                 domain = []
                 auxiliary_domains = {}
-        elif self._integration_domain == "secondary":
+        elif self._integration_dimension == "secondary":
             # integral in the secondary dimension keeps the same domain, moves tertiary
             # domain to secondary domain
             domain = child.domain
             if "tertiary" in child.auxiliary_domains:
                 auxiliary_domains = {"secondary": child.auxiliary_domains["tertiary"]}
             else:
                 auxiliary_domains = {}
-        elif self._integration_domain == "tertiary":
+        elif self._integration_dimension == "tertiary":
             # integral in the tertiary dimension keeps the domain and secondary domain
             domain = child.domain
             auxiliary_domains = {"secondary": child.auxiliary_domains["secondary"]}

pybamm/spatial_methods/finite_volume.py

@@ -234,13 +234,19 @@ def laplacian(self, symbol, discretised_symbol, boundary_conditions):
         grad = self.gradient(symbol, discretised_symbol, boundary_conditions)
         return self.divergence(grad, grad, boundary_conditions)
 
-    def integral(self, child, discretised_child):
+    def integral(self, child, discretised_child, integration_dimension):
         """Vector-vector dot product to implement the integral operator. """
+        if integration_dimension == "primary":
+            domain = child.domain
+        else:
+            domain = child.auxiliary_domains[integration_dimension]
         # Calculate integration vector
-        integration_vector = self.definite_integral_matrix(child.domain)
+        integration_vector = self.definite_integral_matrix(
+            domain, integration_dimension=integration_dimension
+        )
 
         # Check for spherical domains
-        submesh_list = self.mesh.combine_submeshes(*child.domain)
+        submesh_list = self.mesh.combine_submeshes(*domain)
         if submesh_list[0].coord_sys == "spherical polar":
             second_dim = len(submesh_list)
             r_numpy = np.kron(np.ones(second_dim), submesh_list[0].nodes)
@@ -249,11 +255,11 @@ 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"):
+    def definite_integral_matrix(
+        self, domain, vector_type="row", integration_dimension="primary"
+    ):
         """
         Matrix for finite-volume implementation of the definite integral in the
         primary dimension
@@ -268,21 +274,23 @@ def definite_integral_matrix(self, domain, vector_type="row"):
         ----------
         domain : list
             The domain(s) of integration
+        vector_type : str, optional
+            Whether to return a row or column vector in the primary dimension
+            (default is row)
+        integration_dimension : str, optional
+            The dimension in which to integrate (default is "primary")
 
         Returns
         -------
         :class:`pybamm.Matrix`
             The finite volume integral matrix for the domain
-        vector_type : str, optional
-            Whether to return a row or column vector in the primary dimension
-            (default is row)
         """
         # Create appropriate submesh by combining submeshes in domain
         submesh_list = self.mesh.combine_submeshes(*domain)
 
         # Create vector of ones for primary domain submesh
         submesh = submesh_list[0]
-        vector = submesh.d_edges * np.ones_like(submesh.nodes)
+        vector = submesh.d_edges
 
         if vector_type == "row":
             vector = vector[np.newaxis, :]
@@ -309,12 +317,22 @@ def indefinite_integral(self, child, discretised_child, direction):
                 child.domain, direction
             )
         else:
+            # Check coordinate system is not spherical polar for the case where child
+            # evaluates on edges
+            # If it becomes necessary to implement this, will need to think about what
+            # the spherical polar indefinite integral should be
+            submesh_list = self.mesh.combine_submeshes(*child.domain)
+            if submesh_list[0].coord_sys == "spherical polar":
+                raise NotImplementedError(
+                    "Indefinite integral on a spherical polar domain is not implemented"
+                )
             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)
+        # Don't need to check for spherical domains as we have ruled out spherical
+        # polars in the case that involves integrating a divergence
+        # (child evaluates on nodes)
         out = integration_matrix @ discretised_child
 
         out.copy_domains(child)

pybamm/spatial_methods/scikit_finite_element.py

@@ -293,7 +293,7 @@ def stiffness_form(u, du, v, dv, w):
 
         return pybamm.Matrix(stiffness)
 
-    def integral(self, child, discretised_child):
+    def integral(self, child, discretised_child, integration_dimension):
         """Vector-vector dot product to implement the integral operator.
         See :meth:`pybamm.SpatialMethod.integral`
         """
@@ -344,7 +344,7 @@ def integral_form(v, dv, w):
         elif vector_type == "column":
             return pybamm.Matrix(vector[:, np.newaxis])
 
-    def indefinite_integral(self, child, discretised_child):
+    def indefinite_integral(self, child, discretised_child, direction):
         """Implementation of the indefinite integral operator. The
         input discretised child must be defined on the internal mesh edges.
         See :meth:`pybamm.SpatialMethod.indefinite_integral`

pybamm/spatial_methods/spatial_method.py

@@ -223,7 +223,7 @@ def gradient_squared(self, symbol, discretised_symbol, boundary_conditions):
         """
         raise NotImplementedError
 
-    def integral(self, child, discretised_child):
+    def integral(self, child, discretised_child, integration_dimension):
         """
         Implements the integral for a spatial method.
 
@@ -233,6 +233,8 @@ def integral(self, child, discretised_child):
             The symbol to which is being integrated
         discretised_child: :class:`pybamm.Symbol`
             The discretised symbol of the correct size
+        integration_dimension : str, optional
+            The dimension in which to integrate (default is "primary")
 
         Returns
         -------

pybamm/spatial_methods/zero_dimensional_method.py

@@ -48,7 +48,7 @@ def indefinite_integral(self, child, discretised_child, direction):
         elif direction == "backward":
             return -discretised_child
 
-    def integral(self, child, discretised_child):
+    def integral(self, child, discretised_child, integration_dimension):
         """
         Calculates the zero-dimensional integral, i.e. the identity operator
         """

tests/shared.py

@@ -99,23 +99,32 @@ def get_p2d_mesh_for_testing(xpts=None, rpts=10):
 
 
 def get_1p1d_mesh_for_testing(
-    xpts=None, zpts=15, cc_submesh=pybamm.MeshGenerator(pybamm.Uniform1DSubMesh)
+    xpts=None,
+    rpts=10,
+    zpts=15,
+    cc_submesh=pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
 ):
-    geometry = pybamm.Geometry("1+1D macro")
+    geometry = pybamm.Geometry("1+1D macro", "1+1D micro")
     return get_mesh_for_testing(
-        xpts=xpts, zpts=zpts, geometry=geometry, cc_submesh=cc_submesh
+        xpts=xpts, rpts=rpts, zpts=zpts, geometry=geometry, cc_submesh=cc_submesh
     )
 
 
 def get_2p1d_mesh_for_testing(
     xpts=None,
+    rpts=10,
     ypts=15,
     zpts=15,
     cc_submesh=pybamm.MeshGenerator(pybamm.ScikitUniform2DSubMesh),
 ):
-    geometry = pybamm.Geometry("2+1D macro")
+    geometry = pybamm.Geometry("2+1D macro", "1+1D micro")
     return get_mesh_for_testing(
-        xpts=xpts, zpts=zpts, geometry=geometry, cc_submesh=cc_submesh
+        xpts=xpts,
+        rpts=rpts,
+        ypts=ypts,
+        zpts=zpts,
+        geometry=geometry,
+        cc_submesh=cc_submesh,
     )
 
 
@@ -170,12 +179,14 @@ def get_p2d_discretisation_for_testing(xpts=None, rpts=10):
     return get_discretisation_for_testing(mesh=get_p2d_mesh_for_testing(xpts, rpts))
 
 
-def get_1p1d_discretisation_for_testing(xpts=None, zpts=15):
-    return get_discretisation_for_testing(mesh=get_1p1d_mesh_for_testing(xpts, zpts))
+def get_1p1d_discretisation_for_testing(xpts=None, rpts=10, zpts=15):
+    return get_discretisation_for_testing(
+        mesh=get_1p1d_mesh_for_testing(xpts, rpts, zpts)
+    )
 
 
-def get_2p1d_discretisation_for_testing(xpts=None, ypts=15, zpts=15):
+def get_2p1d_discretisation_for_testing(xpts=None, rpts=10, ypts=15, zpts=15):
     return get_discretisation_for_testing(
-        mesh=get_2p1d_mesh_for_testing(xpts, ypts, zpts),
+        mesh=get_2p1d_mesh_for_testing(xpts, rpts, ypts, zpts),
         cc_method=pybamm.ScikitFiniteElement,
     )

tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py

@@ -691,6 +691,86 @@ def test_definite_integral(self):
             integral_eqn_disc.evaluate(None, one_over_y), 4 * np.pi ** 2
         )
 
+    def test_integral_secondary_domain(self):
+        # create discretisation
+        mesh = get_1p1d_mesh_for_testing(xpts=200, rpts=200, zpts=100)
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume(),
+            "negative particle": pybamm.FiniteVolume(),
+            "positive particle": pybamm.FiniteVolume(),
+            "current collector": pybamm.FiniteVolume(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        # lengths
+        ln = mesh["negative electrode"][0].edges[-1]
+        ls = mesh["separator"][0].edges[-1] - ln
+        lp = 1 - (ln + ls)
+
+        var = pybamm.Variable(
+            "var",
+            domain="positive particle",
+            auxiliary_domains={
+                "secondary": "positive electrode",
+                "tertiary": "current collector",
+            },
+        )
+        x = pybamm.SpatialVariable("x", "positive electrode")
+        integral_eqn = pybamm.Integral(var, x)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+
+        submesh = mesh["positive particle"]
+        constant_y = np.ones_like(submesh[0].nodes[:, np.newaxis])
+        self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), ls + lp)
+        linear_y = submesh[0].nodes
+        self.assertAlmostEqual(
+            integral_eqn_disc.evaluate(None, linear_y)[0][0], (1 - (ln) ** 2) / 2
+        )
+        cos_y = np.cos(submesh[0].nodes[:, np.newaxis])
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc.evaluate(None, cos_y), np.sin(1) - np.sin(ln), decimal=4
+        )
+
+    def test_integral_primary_then_secondary_same_result(self):
+        # Test that integrating in r then in x gives the same result as integrating in
+        # x then in r
+        # create discretisation
+        mesh = get_1p1d_mesh_for_testing(xpts=200, rpts=200, zpts=100)
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume(),
+            "negative particle": pybamm.FiniteVolume(),
+            "positive particle": pybamm.FiniteVolume(),
+            "current collector": pybamm.FiniteVolume(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable(
+            "var",
+            domain="positive particle",
+            auxiliary_domains={
+                "secondary": "positive electrode",
+                "tertiary": "current collector",
+            },
+        )
+        x = pybamm.SpatialVariable("x", "positive electrode")
+        r = pybamm.SpatialVariable("r", "positive particle")
+        integral_eqn_x_then_r = pybamm.Integral(pybamm.Integral(var, x), r)
+        integral_eqn_r_then_x = pybamm.Integral(pybamm.Integral(var, r), x)
+
+        # discretise
+        disc.set_variable_slices([var])
+        integral_eqn_x_then_r_disc = disc.process_symbol(integral_eqn_x_then_r)
+        integral_eqn_r_then_x_disc = disc.process_symbol(integral_eqn_r_then_x)
+
+        # test
+        submesh = mesh["positive particle"]
+        cos_y = np.cos(submesh[0].nodes[:, np.newaxis])
+        np.testing.assert_array_almost_equal(
+            integral_eqn_x_then_r_disc.evaluate(None, cos_y),
+            integral_eqn_r_then_x_disc.evaluate(None, cos_y),
+            decimal=4,
+        )
+
     def test_definite_integral_vector(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {
@@ -833,7 +913,7 @@ def test_indefinite_integral(self):
         # --------------------------------------------------------------------
         # micrsoscale case
         c = pybamm.Variable("c", domain=["negative particle"])
-        N = pybamm.grad(c)  # create test current (variable on edges)
+        N = pybamm.grad(c)  # create test flux (variable on edges)
         r_n = pybamm.SpatialVariable("r_n", ["negative particle"])
         c_integral = pybamm.IndefiniteIntegral(N, r_n)
         disc.set_variable_slices([c])  # N is not a fundamental variable
@@ -987,6 +1067,22 @@ def test_indefinite_integral_on_nodes(self):
         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)
 
+        # microscale case should fail
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"negative particle": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        c = pybamm.Variable("c", domain=["negative particle"])
+        r = pybamm.SpatialVariable("r", ["negative particle"])
+
+        int_c = pybamm.IndefiniteIntegral(c, r)
+        disc.set_variable_slices([c])
+        with self.assertRaisesRegex(
+            NotImplementedError,
+            "Indefinite integral on a spherical polar domain is not implemented",
+        ):
+            disc.process_symbol(int_c)
+
     def test_backward_indefinite_integral_on_nodes(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}

tests/unit/test_spatial_methods/test_scikit_finite_element.py

@@ -16,7 +16,7 @@ def test_not_implemented(self):
         with self.assertRaises(NotImplementedError):
             spatial_method.divergence(None, None, None)
         with self.assertRaises(NotImplementedError):
-            spatial_method.indefinite_integral(None, None)
+            spatial_method.indefinite_integral(None, None, None)
 
     def test_discretise_equations(self):
         # get mesh