Merge pull request #1782 from pybamm-team/spherical-finite-volume

Spherical finite volume

CHANGELOG.md

@@ -1,7 +1,11 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Features
+
+-   Half-cell SPM and SPMe have been implemented (  [#1731](https://github.com/pybamm-team/PyBaMM/pull/1731))
 ## Bug fixes
 
+-   Fixed finite volume discretization in spherical polar coordinates ([#1782](https://github.com/pybamm-team/PyBaMM/pull/1782))
 -   Fixed `sympy` operators for `Arctan` and `Exponential` ([#1786](https://github.com/pybamm-team/PyBaMM/pull/1786))
 
 # [v21.10](https://github.com/pybamm-team/PyBaMM/tree/v21.9) - 2021-10-31

pybamm/models/full_battery_models/base_battery_model.py

@@ -414,8 +414,8 @@ def default_var_pts(self):
             var.x_n: 20,
             var.x_s: 20,
             var.x_p: 20,
-            var.r_n: 30,
-            var.r_p: 30,
+            var.r_n: 20,
+            var.r_p: 20,
             var.y: 10,
             var.z: 10,
             var.R_n: 30,

pybamm/spatial_methods/finite_volume.py

@@ -171,18 +171,12 @@ def divergence(self, symbol, discretised_symbol, boundary_conditions):
         # check for particle domain
         if submesh.coord_sys == "spherical polar":
             second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains)
-            edges = submesh.edges
 
-            # create np.array of repeated submesh.nodes
-            r_numpy = np.kron(np.ones(second_dim_repeats), submesh.nodes)
-            r_edges_numpy = np.kron(np.ones(second_dim_repeats), edges)
-
-            r = pybamm.Vector(r_numpy)
+            # create np.array of repeated submesh.edges
+            r_edges_numpy = np.kron(np.ones(second_dim_repeats), submesh.edges)
             r_edges = pybamm.Vector(r_edges_numpy)
 
-            out = (1 / (r ** 2)) * (
-                divergence_matrix @ ((r_edges ** 2) * discretised_symbol)
-            )
+            out = divergence_matrix @ ((r_edges ** 2) * discretised_symbol)
         else:
             out = divergence_matrix @ discretised_symbol
 
@@ -205,7 +199,14 @@ def divergence_matrix(self, domains):
         """
         # Create appropriate submesh by combining submeshes in domain
         submesh = self.mesh.combine_submeshes(*domains["primary"])
-        e = 1 / submesh.d_edges
+        if submesh.coord_sys == "spherical polar":
+            r_edges_left = submesh.edges[:-1]
+            r_edges_right = submesh.edges[1:]
+            d_edges = (r_edges_right ** 3 - r_edges_left ** 3) / 3
+        else:
+            d_edges = submesh.d_edges
+
+        e = 1 / d_edges
 
         # Create matrix using submesh
         n = submesh.npts + 1
@@ -234,17 +235,7 @@ def integral(self, child, discretised_child, integration_dimension):
         integration_vector = self.definite_integral_matrix(
             child, integration_dimension=integration_dimension
         )
-
-        # Check for spherical domains
-        domain = child.domains[integration_dimension]
-        submesh = self.mesh.combine_submeshes(*domain)
-        if submesh.coord_sys == "spherical polar":
-            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 * integration_vector @ (discretised_child * r ** 2)
-        else:
-            out = integration_vector @ discretised_child
+        out = integration_vector @ discretised_child
 
         return out
 
@@ -277,33 +268,40 @@ def definite_integral_matrix(
             The finite volume integral matrix for the domain
         """
         domains = child.domains
+        if vector_type != "row" and integration_dimension == "secondary":
+            raise NotImplementedError(
+                "Integral in secondary vector only implemented in 'row' form"
+            )
+
+        domain = child.domains[integration_dimension]
+        submesh = self.mesh.combine_submeshes(*domain)
+        if submesh.coord_sys == "spherical polar":
+            r_edges_left = submesh.edges[:-1]
+            r_edges_right = submesh.edges[1:]
+            d_edges = 4 * np.pi * (r_edges_right ** 3 - r_edges_left ** 3) / 3
+        else:
+            d_edges = submesh.d_edges
+
         if integration_dimension == "primary":
             # Create appropriate submesh by combining submeshes in domain
             submesh = self.mesh.combine_submeshes(*domains["primary"])
 
             # Create vector of ones for primary domain submesh
-            vector = submesh.d_edges
 
             if vector_type == "row":
-                vector = vector[np.newaxis, :]
+                d_edges = d_edges[np.newaxis, :]
             elif vector_type == "column":
-                vector = vector[:, np.newaxis]
+                d_edges = d_edges[:, np.newaxis]
 
             # repeat matrix for each node in secondary dimensions
             second_dim_repeats = self._get_auxiliary_domain_repeats(domains)
             # generate full matrix from the submatrix
-            matrix = kron(eye(second_dim_repeats), vector)
+            matrix = kron(eye(second_dim_repeats), d_edges)
         elif integration_dimension == "secondary":
-            if vector_type != "row":
-                raise NotImplementedError(
-                    "Integral in secondary vector only implemented in 'row' form"
-                )
             # Create appropriate submesh by combining submeshes in domain
             primary_submesh = self.mesh.combine_submeshes(*domains["primary"])
-            secondary_submesh = self.mesh.combine_submeshes(*domains["secondary"])
 
             # Create matrix which integrates in the secondary dimension
-            d_edges = secondary_submesh.d_edges
             # Different number of edges depending on whether child evaluates on edges
             # in the primary dimensions
             if child.evaluates_on_edges("primary"):

pybamm/util.py

@@ -276,6 +276,7 @@ def load_function(filename):
     # Strip absolute path to pybamm/input/example.py
     if "pybamm" in filename:
         root_path = filename[filename.rfind("pybamm") :]
+    # Commenting not removing these lines in case we get problems later
     elif os.getcwd() in filename:
         root_path = filename.replace(os.getcwd(), "")
         root_path = root_path[1:]

tests/integration/test_models/standard_output_tests.py

@@ -377,23 +377,19 @@ def test_concentration_limits(self):
     def test_conservation(self):
         """Test amount of lithium stored across all particles and in SEI layers is
         constant."""
-        self.c_s_tot = (
+        c_s_tot = (
             self.c_s_n_tot(self.solution.t)
             + self.c_s_p_tot(self.solution.t)
             + self.c_SEI_tot(self.solution.t)
             + self.c_Li_tot(self.solution.t)
         )
-        diff = (self.c_s_tot[1:] - self.c_s_tot[:-1]) / self.c_s_tot[:-1]
-        if "profile" in self.model.options["particle"]:
-            np.testing.assert_array_almost_equal(diff, 0, decimal=10)
-        elif self.model.options["particle size"] == "distribution":
+        diff = (c_s_tot[1:] - c_s_tot[:-1]) / c_s_tot[:-1]
+        if self.model.options["particle"] == "quartic profile":
             np.testing.assert_array_almost_equal(diff, 0, decimal=10)
+        # elif self.model.options["particle size"] == "distribution":
+        #     np.testing.assert_array_almost_equal(diff, 0, decimal=10)
         elif self.model.options["surface form"] == "differential":
             np.testing.assert_array_almost_equal(diff, 0, decimal=10)
-        elif self.model.options["SEI"] == "ec reaction limited":
-            np.testing.assert_array_almost_equal(diff, 0, decimal=11)
-        elif self.model.options["lithium plating"] == "irreversible":
-            np.testing.assert_array_almost_equal(diff, 0, decimal=13)
         else:
             np.testing.assert_array_almost_equal(diff, 0, decimal=15)
 
@@ -779,9 +775,9 @@ def __init__(self, model, param, disc, solution, operating_condition):
 
     def test_degradation_modes(self):
         """Test degradation modes are between 0 and 100%"""
-        np.testing.assert_array_less(-1e-3, self.LLI(self.t))
-        np.testing.assert_array_less(-1e-3, self.LAM_ne(self.t))
-        np.testing.assert_array_less(-1e-3, self.LAM_pe(self.t))
+        np.testing.assert_array_less(-3e-3, self.LLI(self.t))
+        np.testing.assert_array_less(-1e-13, self.LAM_ne(self.t))
+        np.testing.assert_array_less(-1e-13, self.LAM_pe(self.t))
         np.testing.assert_array_less(self.LLI(self.t), 100)
         np.testing.assert_array_less(self.LAM_ne(self.t), 100)
         np.testing.assert_array_less(self.LAM_pe(self.t), 100)

tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py

@@ -81,8 +81,8 @@ def test_conservation_each_electrode(self):
             pos_Li.append(pos)
 
         # compare
-        np.testing.assert_array_almost_equal(neg_Li[0], neg_Li[1], decimal=14)
-        np.testing.assert_array_almost_equal(pos_Li[0], pos_Li[1], decimal=14)
+        np.testing.assert_array_almost_equal(neg_Li[0], neg_Li[1], decimal=13)
+        np.testing.assert_array_almost_equal(pos_Li[0], pos_Li[1], decimal=13)
 
 
 if __name__ == "__main__":

tests/integration/test_spatial_methods/test_finite_volume.py

@@ -126,11 +126,10 @@ def get_error(n):
             r_edge = pybamm.SpatialVariableEdge("r_n", domain=["negative particle"])
 
             # Define flux and bcs
-            N = r_edge ** 2 * pybamm.sin(r_edge)
+            N = pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
-            # N = r**3 --> div(N) = 5 * r**2
-            div_exact = 4 * r * np.sin(r) + r ** 2 * np.cos(r)
+            div_exact = 2 / r * np.sin(r) + np.cos(r)
 
             # Discretise and evaluate
             div_eqn_disc = disc.process_symbol(div_eqn)
@@ -140,7 +139,7 @@ def get_error(n):
             return div_approx[:, 0] - div_exact
 
         # Get errors
-        ns = 10 * 2 ** np.arange(6)
+        ns = 10 * 2 ** np.arange(1, 7)
         errs = {n: get_error(int(n)) for n in ns}
         # expect quadratic convergence everywhere
         err_norm = np.array([np.linalg.norm(errs[n], np.inf) for n in ns])
@@ -195,11 +194,12 @@ def get_error(m):
             r = submesh.nodes
             r_edge = pybamm.standard_spatial_vars.r_n_edge
 
-            N = r_edge ** 2 * pybamm.sin(r_edge)
+            # N = r_edge ** 2 * pybamm.sin(r_edge)
+            N = pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
-            # N = r**2*sin(r) --> div(N) = 4*r*sin(r) - r**2*cos(r)
-            div_exact = 4 * r * np.sin(r) + r ** 2 * np.cos(r)
+            # N = sin(r) --> div(N) = 1/r2 * d/dr(r2*N) = 2/r*sin(r) + cos(r)
+            div_exact = 2 / r * np.sin(r) + np.cos(r)
             div_exact = np.kron(np.ones(mesh["negative electrode"].npts), div_exact)
 
             # Discretise and evaluate
@@ -209,7 +209,7 @@ def get_error(m):
             return div_approx[:, 0] - div_exact
 
         # Get errors
-        ns = 10 * 2 ** np.arange(6)
+        ns = 10 * 2 ** np.arange(1, 7)
         errs = {n: get_error(int(n)) for n in ns}
         # expect quadratic convergence everywhere
         err_norm = np.array([np.linalg.norm(errs[n], np.inf) for n in ns])
@@ -233,13 +233,11 @@ def get_error(m):
             r_edge = pybamm.standard_spatial_vars.r_n_edge
             x = pybamm.standard_spatial_vars.x_n
 
-            N = pybamm.PrimaryBroadcast(x, "negative particle") * (
-                r_edge ** 2 * pybamm.sin(r_edge)
-            )
+            N = pybamm.PrimaryBroadcast(x, "negative particle") * pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
             # N = r**2*sin(r) --> div(N) = 4*r*sin(r) - r**2*cos(r)
-            div_exact = 4 * r * np.sin(r) + r ** 2 * np.cos(r)
+            div_exact = 2 / r * np.sin(r) + np.cos(r)
             div_exact = np.kron(mesh["negative electrode"].nodes, div_exact)
 
             # Discretise and evaluate
@@ -249,7 +247,7 @@ def get_error(m):
             return div_approx[:, 0] - div_exact
 
         # Get errors
-        ns = 10 * 2 ** np.arange(6)
+        ns = 10 * 2 ** np.arange(1, 7)
         errs = {n: get_error(int(n)) for n in ns}
         # expect quadratic convergence everywhere
         err_norm = np.array([np.linalg.norm(errs[n], np.inf) for n in ns])

tests/integration/test_spatial_methods/test_spectral_volume.py

@@ -182,7 +182,7 @@ def get_error(n):
         np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
 
     def test_spherical_div_convergence_quadratic(self):
-        # test div( r**2 * sin(r) ) == 4*r*sin(r) - r**2*cos(r)
+        # test div( r**2 * sin(r) ) == 2/r*sin(r) + cos(r)
         spatial_methods = {"negative particle": pybamm.SpectralVolume()}
 
         # Function for convergence testing
@@ -195,11 +195,11 @@ def get_error(n):
             r_edge = pybamm.SpatialVariableEdge("r_n", domain=["negative particle"])
 
             # Define flux and bcs
-            N = r_edge ** 2 * pybamm.sin(r_edge)
+            N = pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
             # N = r**3 --> div(N) = 5 * r**2
-            div_exact = 4 * r * np.sin(r) + r ** 2 * np.cos(r)
+            div_exact = 2 / r * np.sin(r) + np.cos(r)
 
             # Discretise and evaluate
             div_eqn_disc = disc.process_symbol(div_eqn)
@@ -252,7 +252,7 @@ def get_error(n):
         np.testing.assert_array_less(0.99 * np.ones_like(rates), rates)
 
     def test_p2d_spherical_convergence_quadratic(self):
-        # test div( r**2 * sin(r) ) == 4*r*sin(r) - r**2*cos(r)
+        # test div( r**2 * sin(r) ) == 2/r*sin(r) + cos(r)
         spatial_methods = {"negative particle": pybamm.SpectralVolume()}
 
         # Function for convergence testing
@@ -264,11 +264,11 @@ def get_error(m):
             r = submesh.nodes
             r_edge = pybamm.standard_spatial_vars.r_n_edge
 
-            N = r_edge ** 2 * pybamm.sin(r_edge)
+            N = pybamm.sin(r_edge)
             div_eqn = pybamm.div(N)
             # Define exact solutions
-            # N = r**2*sin(r) --> div(N) = 4*r*sin(r) - r**2*cos(r)
-            div_exact = 4 * r * np.sin(r) + r ** 2 * np.cos(r)
+            # N = r**2*sin(r) --> div(N) = 2/r*sin(r) + cos(r)
+            div_exact = 2 / r * np.sin(r) + np.cos(r)
             div_exact = np.kron(np.ones(mesh["negative electrode"].npts), div_exact)
 
             # Discretise and evaluate
@@ -286,7 +286,7 @@ def get_error(m):
         np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
 
     def test_p2d_with_x_dep_bcs_spherical_convergence(self):
-        # test div_r( (r**2 * sin(r)) * x ) == (4*r*sin(r) - r**2*cos(r)) * x
+        # test div_r( sin(r) * x ) == (2/r*sin(r) + cos(r)) * x
         spatial_methods = {
             "negative particle": pybamm.SpectralVolume(),
             "negative electrode": pybamm.SpectralVolume(),
@@ -302,13 +302,11 @@ def get_error(m):
             r_edge = pybamm.standard_spatial_vars.r_n_edge
             x = pybamm.standard_spatial_vars.x_n
 
-            N = pybamm.PrimaryBroadcast(x, "negative particle") * (
-                r_edge ** 2 * pybamm.sin(r_edge)
-            )
+            N = pybamm.PrimaryBroadcast(x, "negative particle") * (pybamm.sin(r_edge))
             div_eqn = pybamm.div(N)
             # Define exact solutions
-            # N = r**2*sin(r) --> div(N) = 4*r*sin(r) - r**2*cos(r)
-            div_exact = 4 * r * np.sin(r) + r ** 2 * np.cos(r)
+            # N = sin(r) --> div(N) = 2/r*sin(r) + cos(r)
+            div_exact = 2 / r * np.sin(r) + np.cos(r)
             div_exact = np.kron(mesh["negative electrode"].nodes, div_exact)
 
             # Discretise and evaluate

tests/unit/test_models/test_full_battery_models/test_base_battery_model.py

@@ -124,8 +124,8 @@ def test_default_var_pts(self):
             var.x_n: 20,
             var.x_s: 20,
             var.x_p: 20,
-            var.r_n: 30,
-            var.r_p: 30,
+            var.r_n: 20,
+            var.r_p: 20,
             var.y: 10,
             var.z: 10,
             var.R_n: 30,

tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.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), 4 * np.pi / 3, decimal=4
+            integral_eqn_disc.evaluate(None, constant_y), 4 * np.pi / 3
         )
         linear_y = mesh["negative particle"].nodes
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, linear_y), np.pi, decimal=3
+            integral_eqn_disc.evaluate(None, linear_y), np.pi, decimal=4
         )
-        one_over_y_squared = 1 / mesh["negative particle"].nodes ** 2
+        one_over_y = 1 / mesh["negative particle"].nodes
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, one_over_y_squared), 4 * np.pi
+            integral_eqn_disc.evaluate(None, one_over_y), 2 * np.pi, decimal=3
         )
 
         # test failure for secondary dimension column form