Merge pull request #1411 from pybamm-team/issue-1399-poly-bug

#1399 fix polynomial and average bugs

CHANGELOG.md

@@ -33,6 +33,8 @@
 
 ## Bug fixes
 
+-   Fixed a bug where the `PolynomialSingleParticle` submodel gave incorrect results with "dimensionality" equal to 2 ([#1411](https://github.com/pybamm-team/PyBaMM/pull/1411))
+-   Fixed a bug where volume averaging in 0D gave the wrong result ([#1411](https://github.com/pybamm-team/PyBaMM/pull/1411))
 -   Fixed a sign error in the positive electrode ohmic losses ([#1407](https://github.com/pybamm-team/PyBaMM/pull/1407))
 -   Fixed the formulation of the EC reaction SEI model ([#1397](https://github.com/pybamm-team/PyBaMM/pull/1397))
 -   Simulations now stop when an experiment becomes infeasible ([#1395](https://github.com/pybamm-team/PyBaMM/pull/1395))

pybamm/expression_tree/unary_operators.py

@@ -1112,8 +1112,14 @@ def surf(symbol):
     return boundary_value(symbol, "right")
 
 
+#
+# Methods for averaging
+#
+
+
 def x_average(symbol):
-    """convenience function for creating an average in the x-direction
+    """
+    convenience function for creating an average in the x-direction
 
     Parameters
     ----------
@@ -1219,10 +1225,13 @@ def z_average(symbol):
         return symbol.orphans[0]
     # Otherwise, use Integral to calculate average value
     else:
-        geo = pybamm.geometric_parameters
+        # We compute the length as Integral(1, z) as this will be easier to identify
+        # for simplifications later on and it gives the correct behaviour when using
+        # ZeroDimensionalSpatialMethod
         z = pybamm.standard_spatial_vars.z
-        l_z = geo.l_z
-        return Integral(symbol, z) / l_z
+        v = pybamm.ones_like(symbol)
+        l = pybamm.Integral(v, z)
+        return Integral(symbol, z) / l
 
 
 def yz_average(symbol):
@@ -1256,12 +1265,14 @@ def yz_average(symbol):
         return symbol.orphans[0]
     # Otherwise, use Integral to calculate average value
     else:
-        geo = pybamm.geometric_parameters
+        # We compute the area as Integral(1, [y,z]) as this will be easier to identify
+        # for simplifications later on and it gives the correct behaviour when using
+        # ZeroDimensionalSpatialMethod
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
-        l_y = geo.l_y
-        l_z = geo.l_z
-        return Integral(symbol, [y, z]) / (l_y * l_z)
+        v = pybamm.ones_like(symbol)
+        A = pybamm.Integral(v, [y, z])
+        return Integral(symbol, [y, z]) / A
 
 
 def r_average(symbol):

pybamm/models/full_battery_models/base_battery_model.py

@@ -46,14 +46,17 @@ class BaseBatteryModel(pybamm.BaseModel):
             * "lithium plating" : str, optional
                 Sets the model for lithium plating. Can be "none" (default),
                 "reversible" or "irreversible".
-            * "loss of active material" : str, optional
-                Sets the model for loss of active material. Can be "none" (default) or
-                "example", which is a placeholder for LAM models.
-            * "particle" : str, optional
             * "loss of active material" : str
                 Sets the model for loss of active material. Can be "none" (default),
                 "positive", "negative" or "both" to enable it for the specific
                 electrode.
+            * "operating mode" : str
+                Sets the operating mode for the model. Can be "current" (default),
+                "voltage" or "power". Alternatively, the operating mode can be
+                controlled with an arbitrary function by passing the function directly
+                as the option. In this case the function must define the residual of
+                an algebraic equation. The applied current will be solved for such
+                that the algebraic constraint is satisfied.
             * "particle" : str
                 Sets the submodel to use to describe behaviour within the particle.
                 Can be "Fickian diffusion" (default), "uniform profile",
@@ -130,13 +133,6 @@ class BaseBatteryModel(pybamm.BaseModel):
                 solve an algebraic equation for it. Default is "false", unless "SEI film
                 resistance" is distributed in which case it is automatically set to
                 "true".
-            * "operating mode" : str
-                Sets the operating mode for the model. Can be "current" (default),
-                "voltage" or "power". Alternatively, the operating mode can be
-                controlled with an arbitrary function by passing the function directly
-                as the option. In this case the function must define the residual of
-                an algebraic equation. The applied current will be solved for such
-                that the algebraic constraint is satisfied.
 
     **Extends:** :class:`pybamm.BaseModel`
     """

pybamm/models/submodels/particle/polynomial_single_particle.py

@@ -278,6 +278,10 @@ def get_coupled_variables(self, variables):
         return variables
 
     def set_rhs(self, variables):
+        # Note: we have to use `pybamm.source(rhs, var)` in the rhs dict so that
+        # the scalar source term gets multplied by the correct mass matrix when
+        # using this model with 2D current collectors with the finite element
+        # method (see #1399)
 
         c_s_rxav = variables[
             "Average " + self.domain.lower() + " particle concentration"
@@ -289,10 +293,16 @@ def set_rhs(self, variables):
         ]
 
         if self.domain == "Negative":
-            self.rhs = {c_s_rxav: -3 * j_xav / self.param.a_R_n}
+            self.rhs = {
+                c_s_rxav: pybamm.source(-3 * j_xav / self.param.a_R_n, c_s_rxav)
+            }
 
         elif self.domain == "Positive":
-            self.rhs = {c_s_rxav: -3 * j_xav / self.param.a_R_p / self.param.gamma_p}
+            self.rhs = {
+                c_s_rxav: pybamm.source(
+                    -3 * j_xav / self.param.a_R_p / self.param.gamma_p, c_s_rxav
+                )
+            }
 
         if self.name == "quartic profile":
             # We solve an extra ODE for the average particle concentration gradient
@@ -308,21 +318,27 @@ def set_rhs(self, variables):
             if self.domain == "Negative":
                 self.rhs.update(
                     {
-                        q_s_rxav: -30
-                        * self.param.D_n(c_s_surf_xav, T_xav)
-                        * q_s_rxav
-                        / self.param.C_n
-                        - 45 * j_xav / self.param.a_R_n / 2
+                        q_s_rxav: pybamm.source(
+                            -30
+                            * self.param.D_n(c_s_surf_xav, T_xav)
+                            * q_s_rxav
+                            / self.param.C_n
+                            - 45 * j_xav / self.param.a_R_n / 2,
+                            q_s_rxav,
+                        )
                     }
                 )
             elif self.domain == "Positive":
                 self.rhs.update(
                     {
-                        q_s_rxav: -30
-                        * self.param.D_p(c_s_surf_xav, T_xav)
-                        * q_s_rxav
-                        / self.param.C_p
-                        - 45 * j_xav / self.param.a_R_p / self.param.gamma_p / 2
+                        q_s_rxav: pybamm.source(
+                            -30
+                            * self.param.D_p(c_s_surf_xav, T_xav)
+                            * q_s_rxav
+                            / self.param.C_p
+                            - 45 * j_xav / self.param.a_R_p / self.param.gamma_p / 2,
+                            q_s_rxav,
+                        )
                     }
                 )
 

tests/unit/test_expression_tree/test_unary_operators.py

@@ -632,18 +632,26 @@ def test_yz_average(self):
         self.assertEqual(z_average_broad_a.evaluate(), np.array([1]))
         self.assertEqual(yz_average_broad_a.evaluate(), np.array([1]))
 
-        a = pybamm.Symbol("a", domain=["current collector"])
+        a = pybamm.Variable("a", domain=["current collector"])
         y = pybamm.SpatialVariable("y", ["current collector"])
         z = pybamm.SpatialVariable("z", ["current collector"])
         z_av_a = pybamm.z_average(a)
         yz_av_a = pybamm.yz_average(a)
 
         self.assertIsInstance(yz_av_a, pybamm.Division)
-        self.assertIsInstance(z_av_a, pybamm.Integral)
+        self.assertIsInstance(z_av_a, pybamm.Division)
+        self.assertIsInstance(z_av_a.children[0], pybamm.Integral)
         self.assertIsInstance(yz_av_a.children[0], pybamm.Integral)
-        self.assertEqual(z_av_a.integration_variable[0].domain, z.domain)
+        self.assertEqual(z_av_a.children[0].integration_variable[0].domain, z.domain)
         self.assertEqual(yz_av_a.children[0].integration_variable[0].domain, y.domain)
         self.assertEqual(yz_av_a.children[0].integration_variable[1].domain, z.domain)
+        self.assertIsInstance(z_av_a.children[1], pybamm.Integral)
+        self.assertIsInstance(yz_av_a.children[1], pybamm.Integral)
+        self.assertEqual(z_av_a.children[1].integration_variable[0].domain, a.domain)
+        self.assertEqual(z_av_a.children[1].children[0].id, pybamm.ones_like(a).id)
+        self.assertEqual(yz_av_a.children[1].integration_variable[0].domain, y.domain)
+        self.assertEqual(yz_av_a.children[1].integration_variable[0].domain, z.domain)
+        self.assertEqual(yz_av_a.children[1].children[0].id, pybamm.ones_like(a).id)
         self.assertEqual(z_av_a.domain, [])
         self.assertEqual(yz_av_a.domain, [])
 

tests/unit/test_spatial_methods/test_zero_dimensional_method.py

@@ -4,7 +4,7 @@
 import numpy as np
 import pybamm
 import unittest
-from tests import get_mesh_for_testing
+from tests import get_mesh_for_testing, get_discretisation_for_testing
 
 
 class TestZeroDimensionalSpatialMethod(unittest.TestCase):
@@ -55,6 +55,23 @@ def test_discretise_spatial_variable(self):
                 var_disc.evaluate()[:, 0], mesh.combine_submeshes(*var.domain).edges
             )
 
+    def test_averages(self):
+        # create discretisation
+        disc = get_discretisation_for_testing(
+            cc_method=pybamm.ZeroDimensionalSpatialMethod
+        )
+        # create and discretise variable
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+        var_disc = disc.process_symbol(var)
+        # check average returns the same value
+        y = np.array([1])
+        for expression in [pybamm.z_average(var), pybamm.yz_average(var)]:
+            expr_disc = disc.process_symbol(expression)
+            np.testing.assert_array_equal(
+                var_disc.evaluate(y=y), expr_disc.evaluate(y=y)
+            )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")