Allow for cylindrical and spherical shell domains

CHANGELOG.md

@@ -1,7 +1,7 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
 ## Features
-
+- Spherical and cylindrical shell domains can now be solved with any boundary conditions ([#3237](https://github.com/pybamm-team/PyBaMM/pull/3237))
 - Processed variables now get the spatial variables automatically, allowing plotting of more generic models ([#3234](https://github.com/pybamm-team/PyBaMM/pull/3234))
 
 ## Breaking changes

pybamm/discretisations/discretisation.py

@@ -458,7 +458,11 @@ def process_boundary_conditions(self, model):
             # check if the boundary condition at the origin for sphere domains is other
             # than no flux
             for subdomain in key.domain:
-                if self.mesh[subdomain].coord_sys == "spherical polar":
+                if (
+                    self.mesh[subdomain].coord_sys
+                    in ["spherical polar", "cylindrical polar"]
+                    and list(self.mesh.geometry[subdomain].values())[0]["min"] == 0
+                ):
                     if bcs["left"][0].value != 0 or bcs["left"][1] != "Neumann":
                         raise pybamm.ModelError(
                             "Boundary condition at r = 0 must be a homogeneous "

tests/integration/test_spatial_methods/test_finite_volume.py

@@ -294,6 +294,51 @@ def get_error(m):
         np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
 
 
+def solve_laplace_equation(coord_sys="cartesian"):
+    model = pybamm.BaseModel()
+    r = pybamm.SpatialVariable("r", domain="domain", coord_sys=coord_sys)
+    u = pybamm.Variable("u", domain="domain")
+    del_u = pybamm.div(pybamm.grad(u))
+    model.boundary_conditions = {
+        u: {
+            "left": (pybamm.Scalar(0), "Dirichlet"),
+            "right": (pybamm.Scalar(1), "Dirichlet"),
+        }
+    }
+    model.algebraic = {u: del_u}
+    model.initial_conditions = {u: pybamm.Scalar(0)}
+    model.variables = {"u": u, "r": r}
+    geometry = {"domain": {r: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}}}
+    submesh_types = {"domain": pybamm.Uniform1DSubMesh}
+    var_pts = {r: 500}
+    mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+    spatial_methods = {"domain": pybamm.FiniteVolume()}
+    disc = pybamm.Discretisation(mesh, spatial_methods)
+    disc.process_model(model)
+    solver = pybamm.CasadiAlgebraicSolver()
+    return solver.solve(model)
+
+
+class TestFiniteVolumeLaplacian(TestCase):
+    def test_laplacian_cartesian(self):
+        solution = solve_laplace_equation(coord_sys="cartesian")
+        np.testing.assert_array_almost_equal(
+            solution["u"].entries, solution["r"].entries - 1, decimal=10
+        )
+
+    def test_laplacian_cylindrical(self):
+        solution = solve_laplace_equation(coord_sys="cylindrical polar")
+        np.testing.assert_array_almost_equal(
+            solution["u"].entries, np.log(solution["r"].entries) / np.log(2), decimal=5
+        )
+
+    def test_laplacian_spherical(self):
+        solution = solve_laplace_equation(coord_sys="spherical polar")
+        np.testing.assert_array_almost_equal(
+            solution["u"].entries, 2 - 2 / solution["r"].entries, decimal=5
+        )
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys