From 6e9f6f7399ad0737b15827ccd72d160c97cca89c Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Tue, 2 Mar 2021 16:34:04 +0000 Subject: [PATCH 1/2] #1399 fix polynomial and average bugs --- pybamm/expression_tree/unary_operators.py | 27 +++++++++---- .../full_battery_models/base_battery_model.py | 18 ++++----- .../particle/polynomial_single_particle.py | 40 +++++++++++++------ .../test_unary_operators.py | 14 +++++-- .../test_zero_dimensional_method.py | 19 ++++++++- 5 files changed, 83 insertions(+), 35 deletions(-) diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index 35efac3b9f..cebd2ae2f1 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/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 are 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): diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py index 04e3149312..28162e5e9e 100644 --- a/pybamm/models/full_battery_models/base_battery_model.py +++ b/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` """ diff --git a/pybamm/models/submodels/particle/polynomial_single_particle.py b/pybamm/models/submodels/particle/polynomial_single_particle.py index 83d2ce7cc4..75b7c43273 100644 --- a/pybamm/models/submodels/particle/polynomial_single_particle.py +++ b/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, + ) } ) diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index 474df9ff7a..07661c825f 100644 --- a/tests/unit/test_expression_tree/test_unary_operators.py +++ b/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, []) 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 bf827a6204..2ba6f4d22f 100644 --- a/tests/unit/test_spatial_methods/test_zero_dimensional_method.py +++ b/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") From 43e4789a213a7517b37c407f9d455e1afacbb9c3 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 3 Mar 2021 13:18:21 +0000 Subject: [PATCH 2/2] #1399 changelog --- CHANGELOG.md | 2 ++ pybamm/expression_tree/unary_operators.py | 2 +- pybamm/models/submodels/particle/polynomial_single_particle.py | 2 +- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 996e4b4232..d168ee73f0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -31,6 +31,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)) diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index cebd2ae2f1..3fe4db0b8b 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -1265,7 +1265,7 @@ def yz_average(symbol): return symbol.orphans[0] # Otherwise, use Integral to calculate average value else: - # We compute the are as Integral(1, [y,z]) as this will be easier to identify + # 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 diff --git a/pybamm/models/submodels/particle/polynomial_single_particle.py b/pybamm/models/submodels/particle/polynomial_single_particle.py index 75b7c43273..1dc4bfad19 100644 --- a/pybamm/models/submodels/particle/polynomial_single_particle.py +++ b/pybamm/models/submodels/particle/polynomial_single_particle.py @@ -278,7 +278,7 @@ 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 + # 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)