add cylindrical geometry

pybamm/expression_tree/independent_variable.py

@@ -5,7 +5,7 @@
 
 import pybamm
 
-KNOWN_COORD_SYS = ["cartesian", "spherical polar"]
+KNOWN_COORD_SYS = ["cartesian", "cylindrical polar", "spherical polar"]
 
 
 class IndependentVariable(pybamm.Symbol):
@@ -101,9 +101,7 @@ def __init__(self, name, domain=None, auxiliary_domains=None, coord_sys=None):
             raise ValueError("domain must be provided")
 
         # Check symbol name vs domain name
-        if name == "r" and not (len(domain) == 1 and "particle" in domain[0]):
-            raise pybamm.DomainError("domain must be particle if name is 'r'")
-        elif name == "r_n" and domain != ["negative particle"]:
+        if name == "r_n" and domain != ["negative particle"]:
             raise pybamm.DomainError(
                 "domain must be negative particle if name is 'r_n'"
             )

pybamm/geometry/battery_geometry.py

@@ -8,19 +8,23 @@ def battery_geometry(
     include_particles=True,
     options=None,
     current_collector_dimension=0,
+    form_factor="pouch",
 ):
     """
     A convenience function to create battery geometries.
 
     Parameters
     ----------
-    include_particles : bool
-        Whether to include particle domains
-    options : dict
+    include_particles : bool, optional
+        Whether to include particle domains. Can be True (default) or False.
+    options : dict, optional
         Dictionary of model options. Necessary for "particle-size geometry",
         relevant for lithium-ion chemistries.
-    current_collector_dimensions : int, default
-        The dimensions of the current collector. Should be 0 (default), 1 or 2
+    current_collector_dimensions : int, optional
+        The dimensions of the current collector. Can be 0 (default), 1 or 2. For
+        a "cylindrical" form factor the current collector dimension must be 0 or 1.
+    form_factor : str, optional
+        The form factor of the cell. Can be "pouch" (default) or "cylindrical".
 
     Returns
     -------
@@ -33,10 +37,10 @@ def battery_geometry(
     l_n = geo.l_n
     l_s = geo.l_s
     l_n_l_s = l_n + l_s
-
     # Override print_name
     l_n_l_s.print_name = "l_n + l_s"
 
+    # Set up electrode/separator/electrode geometry
     geometry = {
         "negative electrode": {var.x_n: {"min": 0, "max": l_n}},
         "separator": {var.x_s: {"min": l_n, "max": l_n_l_s}},
@@ -51,56 +55,68 @@ def battery_geometry(
             }
         )
     # Add particle size domains
-    if (
-        options is not None and
-        options["particle size"] == "distribution"
-    ):
+    if options is not None and options["particle size"] == "distribution":
         R_min_n = geo.R_min_n
         R_min_p = geo.R_min_p
         R_max_n = geo.R_max_n
         R_max_p = geo.R_max_p
         geometry.update(
             {
-                "negative particle size": {
-                    var.R_n: {"min": R_min_n, "max": R_max_n}
-                },
-                "positive particle size": {
-                    var.R_p: {"min": R_min_p, "max": R_max_p}
-                },
+                "negative particle size": {var.R_n: {"min": R_min_n, "max": R_max_n}},
+                "positive particle size": {var.R_p: {"min": R_min_p, "max": R_max_p}},
             }
         )
-
-    if current_collector_dimension == 0:
-        geometry["current collector"] = {var.z: {"position": 1}}
-    elif current_collector_dimension == 1:
-        geometry["current collector"] = {
-            var.z: {"min": 0, "max": 1},
-            "tabs": {
-                "negative": {"z_centre": geo.centre_z_tab_n},
-                "positive": {"z_centre": geo.centre_z_tab_p},
-            },
-        }
-    elif current_collector_dimension == 2:
-        geometry["current collector"] = {
-            var.y: {"min": 0, "max": geo.l_y},
-            var.z: {"min": 0, "max": geo.l_z},
-            "tabs": {
-                "negative": {
-                    "y_centre": geo.centre_y_tab_n,
-                    "z_centre": geo.centre_z_tab_n,
-                    "width": geo.l_tab_n,
+    # Add current collector domains
+    if form_factor == "pouch":
+        if current_collector_dimension == 0:
+            geometry["current collector"] = {var.z: {"position": 1}}
+        elif current_collector_dimension == 1:
+            geometry["current collector"] = {
+                var.z: {"min": 0, "max": 1},
+                "tabs": {
+                    "negative": {"z_centre": geo.centre_z_tab_n},
+                    "positive": {"z_centre": geo.centre_z_tab_p},
                 },
-                "positive": {
-                    "y_centre": geo.centre_y_tab_p,
-                    "z_centre": geo.centre_z_tab_p,
-                    "width": geo.l_tab_p,
+            }
+        elif current_collector_dimension == 2:
+            geometry["current collector"] = {
+                var.y: {"min": 0, "max": geo.l_y},
+                var.z: {"min": 0, "max": geo.l_z},
+                "tabs": {
+                    "negative": {
+                        "y_centre": geo.centre_y_tab_n,
+                        "z_centre": geo.centre_z_tab_n,
+                        "width": geo.l_tab_n,
+                    },
+                    "positive": {
+                        "y_centre": geo.centre_y_tab_p,
+                        "z_centre": geo.centre_z_tab_p,
+                        "width": geo.l_tab_p,
+                    },
                 },
-            },
-        }
+            }
+        else:
+            raise pybamm.GeometryError(
+                "Invalid current collector dimension '{}' (should be 0, 1 or 2)".format(
+                    current_collector_dimension
+                )
+            )
+    elif form_factor == "cylindrical":
+        if current_collector_dimension == 0:
+            geometry["current collector"] = {var.r: {"position": 1}}
+        elif current_collector_dimension == 1:
+            geometry["current collector"] = {
+                var.r: {"min": geo.r_inner, "max": 1},
+            }
+        else:
+            raise pybamm.GeometryError(
+                "Invalid current collector dimension '{}' (should be 0 or 1 for "
+                "a 'cylindrical' battery geometry)".format(current_collector_dimension)
+            )
     else:
         raise pybamm.GeometryError(
-            "Invalid current collector dimension '{}' (should be 0, 1 or 2)".format(
-                current_collector_dimension
+            "Invalid form factor '{}' (should be 'pouch' or 'cylindrical'".format(
+                form_factor
             )
         )
 

pybamm/parameters/geometric_parameters.py

@@ -37,9 +37,11 @@ def _set_dimensional_parameters(self):
             self.L_n + self.L_s + self.L_p
         )  # Total distance between current collectors
         self.L = self.L_cn + self.L_x + self.L_cp  # Total cell thickness
+        self.L_Li = pybamm.Parameter("Lithium counter electrode thickness [m]")
         self.L_y = pybamm.Parameter("Electrode width [m]")
         self.L_z = pybamm.Parameter("Electrode height [m]")
-        self.L_Li = pybamm.Parameter("Lithium counter electrode thickness [m]")
+        self.r_inner_dimensional = pybamm.Parameter("Inner cell radius [m]")
+        self.r_outer_dimensional = pybamm.Parameter("Outer cell radius [m]")
         self.A_cc = self.L_y * self.L_z  # Area of current collector
         self.A_cooling = pybamm.Parameter("Cell cooling surface area [m2]")
         self.V_cell = pybamm.Parameter("Cell volume [m3]")
@@ -108,7 +110,8 @@ def f_a_dist_n_dimensional(self, R):
             "Negative particle-size variable [m]": R,
         }
         return pybamm.FunctionParameter(
-            "Negative area-weighted particle-size distribution [m-1]", inputs,
+            "Negative area-weighted particle-size distribution [m-1]",
+            inputs,
         )
 
     def f_a_dist_p_dimensional(self, R):
@@ -119,7 +122,8 @@ def f_a_dist_p_dimensional(self, R):
             "Positive particle-size variable [m]": R,
         }
         return pybamm.FunctionParameter(
-            "Positive area-weighted particle-size distribution [m-1]", inputs,
+            "Positive area-weighted particle-size distribution [m-1]",
+            inputs,
         )
 
     def _set_scales(self):
@@ -142,9 +146,11 @@ def _set_dimensionless_parameters(self):
         self.l_p = self.L_p / self.L_x
         self.l_cp = self.L_cp / self.L_x
         self.l_x = self.L_x / self.L_x
+        self.l_Li = self.L_Li / self.L_x
         self.l_y = self.L_y / self.L_z
         self.l_z = self.L_z / self.L_z
-        self.l_Li = self.L_Li / self.L_x
+        self.r_inner = self.r_inner_dimensional / self.r_outer_dimensional
+        self.r_outer = self.r_outer_dimensional / self.r_outer_dimensional
         self.a_cc = self.l_y * self.l_z
         self.a_cooling = self.A_cooling / (self.L_z ** 2)
         self.v_cell = self.V_cell / (self.L_x * self.L_z ** 2)

pybamm/spatial_methods/finite_volume.py

@@ -67,32 +67,6 @@ def spatial_variable(self, symbol):
             entries, domain=symbol.domain, auxiliary_domains=symbol.auxiliary_domains
         )
 
-    def preprocess_external_variables(self, var):
-        """
-        For finite volumes, we need the boundary fluxes for discretising
-        properly. Here, we extrapolate and then add them to the boundary
-        conditions.
-
-        Parameters
-        ----------
-        var : :class:`pybamm.Variable` or :class:`pybamm.Concatenation`
-            The external variable that is to be processed
-
-        Returns
-        -------
-        new_bcs: dict
-            A dictionary containing the new boundary conditions
-        """
-
-        new_bcs = {
-            var: {
-                "left": (pybamm.BoundaryGradient(var, "left"), "Neumann"),
-                "right": (pybamm.BoundaryGradient(var, "right"), "Neumann"),
-            }
-        }
-
-        return new_bcs
-
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the gradient operator.
         See :meth:`pybamm.SpatialMethod.gradient`
@@ -123,6 +97,32 @@ def gradient(self, symbol, discretised_symbol, boundary_conditions):
 
         return out
 
+    def preprocess_external_variables(self, var):
+        """
+        For finite volumes, we need the boundary fluxes for discretising
+        properly. Here, we extrapolate and then add them to the boundary
+        conditions.
+
+        Parameters
+        ----------
+        var : :class:`pybamm.Variable` or :class:`pybamm.Concatenation`
+            The external variable that is to be processed
+
+        Returns
+        -------
+        new_bcs: dict
+            A dictionary containing the new boundary conditions
+        """
+
+        new_bcs = {
+            var: {
+                "left": (pybamm.BoundaryGradient(var, "left"), "Neumann"),
+                "right": (pybamm.BoundaryGradient(var, "right"), "Neumann"),
+            }
+        }
+
+        return new_bcs
+
     def gradient_matrix(self, domain, auxiliary_domains):
         """
         Gradient matrix for finite volumes in the appropriate domain.
@@ -291,7 +291,7 @@ def definite_integral_matrix(
             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
-        if submesh.coord_sys == "cylindrical polar":
+        elif submesh.coord_sys == "cylindrical polar":
             r_edges_left = submesh.edges[:-1]
             r_edges_right = submesh.edges[1:]
             d_edges = 2 * np.pi * (r_edges_right ** 2 - r_edges_left ** 2) / 3

pybamm/spatial_methods/spectral_volume.py

@@ -282,7 +282,7 @@ def gradient(self, symbol, discretised_symbol, boundary_conditions):
                     symbol, reconstructed_symbol, bcs
                 )
 
-        # note in 1D spherical grad and normal grad are the same
+        # note in 1D cartesian, cylindrical and spherical grad are the same
         gradient_matrix = self.gradient_matrix(domain, symbol.auxiliary_domains)
         penalty_matrix = self.penalty_matrix(domain, symbol.auxiliary_domains)
 

tests/shared.py

@@ -42,7 +42,14 @@ def mass_matrix(self, symbol, boundary_conditions):
 
 
 def get_mesh_for_testing(
-    xpts=None, rpts=10, Rpts=10, ypts=15, zpts=15, geometry=None, cc_submesh=None
+    xpts=None,
+    rpts=10,
+    Rpts=10,
+    ypts=15,
+    zpts=15,
+    rcellpts=15,
+    geometry=None,
+    cc_submesh=None,
 ):
     param = pybamm.ParameterValues(
         values={
@@ -59,6 +66,8 @@ def get_mesh_for_testing(
             "Positive electrode thickness [m]": 0.3,
             "Negative particle radius [m]": 0.5,
             "Positive particle radius [m]": 0.5,
+            "Inner cell radius [m]": 0.2,
+            "Outer cell radius [m]": 1.0,
             "Negative minimum particle radius [m]": 0.0,
             "Negative maximum particle radius [m]": 1.0,
             "Positive minimum particle radius [m]": 0.0,
@@ -96,10 +105,10 @@ def get_mesh_for_testing(
         var.r_p: rpts,
         var.y: ypts,
         var.z: zpts,
+        var.r: rcellpts,
         var.R_n: Rpts,
         var.R_p: Rpts,
     }
-
     return pybamm.Mesh(geometry, submesh_types, var_pts)
 
 
@@ -123,7 +132,7 @@ def get_size_distribution_mesh_for_testing(
         Rpts=Rpts,
         zpts=zpts,
         geometry=geometry,
-        cc_submesh=cc_submesh
+        cc_submesh=cc_submesh,
     )
 
 
@@ -195,6 +204,23 @@ def get_unit_2p1D_mesh_for_testing(ypts=15, zpts=15, include_particles=True):
     return pybamm.Mesh(geometry, submesh_types, var_pts)
 
 
+def get_cylindrical_mesh_for_testing(
+    xpts=10, rpts=10, rcellpts=15, include_particles=False
+):
+    geometry = pybamm.battery_geometry(
+        include_particles=include_particles,
+        current_collector_dimension=1,
+        form_factor="cylindrical",
+    )
+    return get_mesh_for_testing(
+        xpts=xpts,
+        rpts=rpts,
+        rcellpts=rcellpts,
+        geometry=geometry,
+        cc_submesh=pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
+    )
+
+
 def get_discretisation_for_testing(
     xpts=None, rpts=10, mesh=None, cc_method=SpatialMethodForTesting
 ):
@@ -223,7 +249,8 @@ def get_size_distribution_disc_for_testing(xpts=None, rpts=10, Rpts=10, zpts=15)
 
 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)
+        mesh=get_1p1d_mesh_for_testing(xpts, rpts, zpts),
+        cc_method=pybamm.FiniteVolume,
     )
 
 
@@ -234,3 +261,12 @@ def get_2p1d_discretisation_for_testing(
         mesh=get_2p1d_mesh_for_testing(xpts, rpts, ypts, zpts, include_particles),
         cc_method=pybamm.ScikitFiniteElement,
     )
+
+
+def get_cylindrical_discretisation_for_testing(
+    xpts=10, rpts=10, rcellpts=15, include_particles=False
+):
+    return get_discretisation_for_testing(
+        mesh=get_cylindrical_mesh_for_testing(xpts, rpts, rcellpts, include_particles),
+        cc_method=pybamm.FiniteVolume,
+    )

tests/unit/test_expression_tree/test_independent_variable.py

@@ -50,8 +50,6 @@ def test_spatial_variable(self):
 
         with self.assertRaisesRegex(ValueError, "domain must be"):
             pybamm.SpatialVariable("x", [])
-        with self.assertRaises(pybamm.DomainError):
-            pybamm.SpatialVariable("r", ["negative electrode"])
         with self.assertRaises(pybamm.DomainError):
             pybamm.SpatialVariable("r_n", ["positive particle"])
         with self.assertRaises(pybamm.DomainError):