Merge pull request #1118 from pybamm-team/issue-1114-r-average

Issue 1114 r average

Author: valentinsulzer

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added "R-averaged particle concentration" variables` ([#1118](https://github.com/pybamm-team/PyBaMM/pull/1118))
 -   Added support for sensitivity calculations to the casadi solver ([#1109](https://github.com/pybamm-team/PyBaMM/pull/1109))
 -   Added support for index 1 semi-explicit dae equations and sensitivity calculations to JAX BDF solver ([#1107](https://github.com/pybamm-team/PyBaMM/pull/1107))
 -   Allowed keyword arguments to be passed to `Simulation.plot()` ([#1099](https://github.com/pybamm-team/PyBaMM/pull/1099))
@@ -10,6 +11,8 @@
 
 ## Bug fixes
 
+-   Fixed `r_average` to work with `SecondaryBroadcast` ([#1118](https://github.com/pybamm-team/PyBaMM/pull/1118))
+-   Fixed finite volume discretisation of spherical integrals ([#1118](https://github.com/pybamm-team/PyBaMM/pull/1118))
 -   `t_eval` now gets changed to a `linspace` if a list of length 2 is passed ([#1113](https://github.com/pybamm-team/PyBaMM/pull/1113))
 -   Fixed bug when setting a function with an `InputParameter` ([#1111](https://github.com/pybamm-team/PyBaMM/pull/1111))
 

examples/notebooks/spatial_methods/finite-volumes.ipynb

@@ -184,7 +184,11 @@ def check_and_set_domains(
         self, child, broadcast_type, broadcast_domain, broadcast_auxiliary_domains
     ):
         "See :meth:`Broadcast.check_and_set_domains`"
-
+        if child.domain == []:
+            raise TypeError(
+                "Cannot take SecondaryBroadcast of an object with empty domain. "
+                "Use PrimaryBroadcast instead."
+            )
         # Can only do secondary broadcast from particle to electrode or from
         # electrode to current collector
         if child.domain[0] in [

pybamm/expression_tree/broadcasts.py

@@ -1167,13 +1167,25 @@ def r_average(symbol):
     # Can't take average if the symbol evaluates on edges
     if symbol.evaluates_on_edges("primary"):
         raise ValueError("Can't take the r-average of a symbol that evaluates on edges")
-    # If symbol doesn't have a particle domain, its r-averaged value is itself
-    if symbol.domain not in [["positive particle"], ["negative particle"]]:
+    # Otherwise, if symbol doesn't have a particle domain,
+    # its r-averaged value is itself
+    elif symbol.domain not in [["positive particle"], ["negative particle"]]:
         new_symbol = symbol.new_copy()
         new_symbol.parent = None
         return new_symbol
-    # If symbol is a Broadcast, its average value is its child
-    elif isinstance(symbol, pybamm.Broadcast):
+    # If symbol is a secondary broadcast onto "negative electrode" or
+    # "positive electrode", take the r-average of the child then broadcast back
+    elif isinstance(symbol, pybamm.SecondaryBroadcast) and symbol.domains[
+        "secondary"
+    ] in [["positive electrode"], ["negative electrode"]]:
+        child = symbol.orphans[0]
+        child_av = pybamm.r_average(child)
+        return pybamm.PrimaryBroadcast(child_av, symbol.domains["secondary"])
+    # If symbol is a Broadcast onto a particle domain, its average value is its child
+    elif isinstance(symbol, pybamm.PrimaryBroadcast) and symbol.domain in [
+        ["positive particle"],
+        ["negative particle"],
+    ]:
         return symbol.orphans[0]
     else:
         r = pybamm.SpatialVariable("r", symbol.domain)

pybamm/expression_tree/unary_operators.py

@@ -34,15 +34,22 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
         elif self.domain == "Positive":
             c_scale = self.param.c_p_max
             active_volume = geo_param.a_p_dim * geo_param.R_p / 3
+
+        c_s_rav = pybamm.r_average(c_s)
         c_s_av = pybamm.r_average(c_s_xav)
         c_s_av_vol = active_volume * c_s_av
+
         variables = {
             self.domain + " particle concentration": c_s,
             self.domain + " particle concentration [mol.m-3]": c_s * c_scale,
             "X-averaged " + self.domain.lower() + " particle concentration": c_s_xav,
             "X-averaged "
             + self.domain.lower()
             + " particle concentration [mol.m-3]": c_s_xav * c_scale,
+            "R-averaged " + self.domain.lower() + " particle concentration": c_s_rav,
+            "R-averaged "
+            + self.domain.lower()
+            + " particle concentration [mol.m-3]": c_s_rav * c_scale,
             self.domain + " particle surface concentration": c_s_surf,
             self.domain
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf,

pybamm/models/submodels/particle/base_particle.py

@@ -242,7 +242,7 @@ def integral(self, child, discretised_child, integration_dimension):
             second_dim_repeats = self._get_auxiliary_domain_repeats(child.domains)
             r_numpy = np.kron(np.ones(second_dim_repeats), submesh.nodes)
             r = pybamm.Vector(r_numpy)
-            out = 4 * np.pi ** 2 * integration_vector @ (discretised_child * r)
+            out = 4 * np.pi * integration_vector @ (discretised_child * r ** 2)
         else:
             out = integration_vector @ discretised_child
 

pybamm/spatial_methods/finite_volume.py

@@ -55,6 +55,9 @@ def test_secondary_broadcast(self):
             {"secondary": ["negative electrode"], "tertiary": ["current collector"]},
         )
 
+        a = pybamm.Symbol("a")
+        with self.assertRaisesRegex(TypeError, "empty domain"):
+            pybamm.SecondaryBroadcast(a, "current collector")
         a = pybamm.Symbol("a", domain="negative particle")
         with self.assertRaisesRegex(
             pybamm.DomainError, "Secondary broadcast from particle"

tests/unit/test_expression_tree/test_broadcasts.py

@@ -439,6 +439,15 @@ def test_r_average(self):
             # electrode domains go to current collector when averaged
             self.assertEqual(av_a.domain, [])
 
+        # r-average of a symbol that is broadcast to x
+        # takes the average of the child then broadcasts it
+        a = pybamm.Scalar(1, domain="positive particle")
+        broad_a = pybamm.SecondaryBroadcast(a, "positive electrode")
+        average_broad_a = pybamm.r_average(broad_a)
+        self.assertIsInstance(average_broad_a, pybamm.PrimaryBroadcast)
+        self.assertEqual(average_broad_a.domain, ["positive electrode"])
+        self.assertEqual(average_broad_a.children[0].id, pybamm.r_average(a).id)
+
         # r-average of symbol that evaluates on edges raises error
         symbol_on_edges = pybamm.PrimaryBroadcastToEdges(1, "domain")
         with self.assertRaisesRegex(

tests/unit/test_expression_tree/test_unary_operators.py

@@ -725,15 +725,15 @@ def test_definite_integral(self):
 
         constant_y = np.ones_like(mesh["negative particle"].nodes[:, np.newaxis])
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, constant_y), 2 * np.pi ** 2
+            integral_eqn_disc.evaluate(None, constant_y), 4 * np.pi / 3, decimal=4
         )
         linear_y = mesh["negative particle"].nodes
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, linear_y), 4 * np.pi ** 2 / 3, decimal=3
+            integral_eqn_disc.evaluate(None, linear_y), np.pi, decimal=3
         )
-        one_over_y = 1 / mesh["negative particle"].nodes
+        one_over_y_squared = 1 / mesh["negative particle"].nodes ** 2
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, one_over_y), 4 * np.pi ** 2
+            integral_eqn_disc.evaluate(None, one_over_y_squared), 4 * np.pi
         )
 
         # test failure for secondary dimension column form