Issue 1174 raise error bcs sphere other than no-flux

CHANGELOG.md

@@ -1,3 +1,19 @@
+# [Unreleased](https://github.com/pybamm-team/PyBaMM)
+
+## Features
+
+
+## Optimizations
+
+
+## Bug fixes
+
+-   Raise error if the boundary condition at the origin in a spherical domain is other than no-flux ([#1175](https://github.com/pybamm-team/PyBaMM/pull/1175))
+-   Fix boundary conditions at r = 0 for Creating Models notebooks ([#1173](https://github.com/pybamm-team/PyBaMM/pull/1173))
+
+## Breaking changes
+
+
 # [v0.2.4](https://github.com/pybamm-team/PyBaMM/tree/v0.2.4) - 2020-09-07
 
 This release adds new operators for more complex models, some basic sensitivity analysis, and a spectral volumes spatial method, as well as some small bug fixes.

pybamm/discretisations/discretisation.py

@@ -492,6 +492,19 @@ def process_boundary_conditions(self, model):
         for key, bcs in model.boundary_conditions.items():
             processed_bcs[key.id] = {}
 
+            # check if the boundary condition at the origin for sphere domains is other
+            # than no flux
+            if key not in model.external_variables:
+                for subdomain in key.domain:
+                    if self.mesh[subdomain].coord_sys == "spherical polar":
+                        if bcs["left"][0].value != 0 or bcs["left"][1] != "Neumann":
+                            raise pybamm.ModelError(
+                                "Boundary condition at r = 0 must be a homogeneous "
+                                "Neumann condition for {} coordinates".format(
+                                    self.mesh[subdomain].coord_sys
+                                )
+                            )
+
             # Handle any boundary conditions applied on the tabs
             if any("tab" in side for side in list(bcs.keys())):
                 bcs = self.check_tab_conditions(key, bcs)

tests/unit/test_discretisations/test_discretisation.py

@@ -1237,6 +1237,58 @@ def test_process_input_variable(self):
         n = disc.mesh.combine_submeshes(*a.domain).npts
         self.assertEqual(a_disc._expected_size, n)
 
+    def test_bc_symmetry(self):
+        # define model
+        model = pybamm.BaseModel()
+        c = pybamm.Variable("Concentration", domain="negative particle")
+        N = -pybamm.grad(c)
+        dcdt = -pybamm.div(N)
+        model.rhs = {c: dcdt}
+
+        # initial conditions
+        model.initial_conditions = {c: pybamm.Scalar(1)}
+
+        # define geometry
+        r = pybamm.SpatialVariable(
+            "r", domain=["negative particle"], coord_sys="spherical polar"
+        )
+        geometry = {
+            "negative particle": {r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}
+        }
+
+        # mesh
+        submesh_types = {
+            "negative particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh)
+        }
+        var_pts = {r: 20}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"negative particle": pybamm.FiniteVolume()}
+
+        # boundary conditions (Dirichlet)
+        lbc = pybamm.Scalar(0)
+        rbc = pybamm.Scalar(2)
+        model.boundary_conditions = {
+            c: {"left": (lbc, "Dirichlet"), "right": (rbc, "Neumann")}
+        }
+
+        # discretise
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        with self.assertRaisesRegex(pybamm.ModelError, "Boundary condition at r = 0"):
+            disc.process_model(model)
+
+        # boundary conditions (non-homog Neumann)
+        lbc = pybamm.Scalar(0)
+        rbc = pybamm.Scalar(2)
+        model.boundary_conditions = {
+            c: {"left": (rbc, "Neumann"), "right": (rbc, "Neumann")}
+        }
+
+        # discretise
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        with self.assertRaisesRegex(pybamm.ModelError, "Boundary condition at r = 0"):
+            disc.process_model(model)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py

@@ -1348,7 +1348,7 @@ def test_p2d_mass_matrix_shape(self):
         model.rhs = {c: pybamm.div(N)}
         model.initial_conditions = {c: pybamm.Scalar(0)}
         model.boundary_conditions = {
-            c: {"left": (0, "Dirichlet"), "right": (0, "Dirichlet")}
+            c: {"left": (0, "Neumann"), "right": (0, "Dirichlet")}
         }
         model.variables = {"c": c, "N": N}
         mesh = get_p2d_mesh_for_testing()