Repository: pybamm-team/PyBaMM

Issue:
Allow for cylindrical & spherical shell domains
### Description

At the moment we cannot simulate spherical shells with a boundary condition at the inner boundary other than no flux, as PyBaMM assumes the domain is always `[0, R]` and thus throws an error, see
https://github.com/pybamm-team/PyBaMM/blob/fffe650253596f14b6918d9e57fdd912bc3cfbdc/pybamm/discretisations/discretisation.py#L457-L464

Apparently it works for cylindrical polar (the error does not catch that case) but it makes me suspect that this feature is not well tested, so it would be good to also add some easy tests against analytical solutions to ensure this functionality works correctly.

### Motivation

Current models do not make use of this feature, but one can imagine how it could be relevant (e.g. core shell models).

In any case, for completeness of the functionality it would be good to enable this anyway.

### Possible Implementation

_No response_

### Additional context

_No response_

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ebc259555..6488dbef4 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index ca2d5b400..c04b41c88 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -458,7 +458,11 @@ class Discretisation(object):
             # 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 "


Test patch:
diff --git a/tests/integration/test_spatial_methods/test_finite_volume.py b/tests/integration/test_spatial_methods/test_finite_volume.py
index 0060e4909..2f9f6ef1b 100644
--- a/tests/integration/test_spatial_methods/test_finite_volume.py
+++ b/tests/integration/test_spatial_methods/test_finite_volume.py
@@ -294,6 +294,51 @@ class TestFiniteVolumeConvergence(TestCase):
         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