fix bug, reimplement value error, make sure multiplicative is always …

…defined

src/pybamm/parameters/parameter_values.py

@@ -604,6 +604,12 @@ def process_geometry(self, geometry):
 
         def process_and_check(sym):
             new_sym = self.process_symbol(sym)
+            if not isinstance(new_sym, pybamm.Scalar) and not isinstance(
+                new_sym, pybamm.InputParameter
+            ):
+                raise ValueError(
+                    "Geometry parameters must be Scalars or InputParameters after parameter processing"
+                )
             return new_sym
 
         for domain in geometry:

src/pybamm/spatial_methods/finite_volume.py

@@ -891,6 +891,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                 # just use the value from the bc: f(x*)
                 sub_matrix = csr_matrix((1, prim_pts))
                 additive = bcs[child][symbol.side][0]
+                additive_multiplicative = pybamm.Scalar(1)
+                multiplicative = pybamm.Scalar(1)
 
             elif symbol.side == "left":
                 if extrap_order_value == "linear":
@@ -903,13 +905,19 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         sub_matrix = csr_matrix(([1], ([0], [0])), shape=(1, prim_pts))
 
                         additive = -dx0 * bcs[child][symbol.side][0]
+                        if hasattr(submesh, "length"):
+                            additive_multiplicative = submesh.length
+                        else:
+                            additive_multiplicative = pybamm.Scalar(1)
+                        multiplicative = pybamm.Scalar(1)
 
                     else:
                         sub_matrix = csr_matrix(
                             ([1 + (dx0 / dx1), -(dx0 / dx1)], ([0, 0], [0, 1])),
                             shape=(1, prim_pts),
                         )
                         additive = pybamm.Scalar(0)
+                        additive_multiplicative = pybamm.Scalar(1)
                         multiplicative = pybamm.Scalar(1)
 
                 elif extrap_order_value == "quadratic":
@@ -924,6 +932,10 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             ([a, b], ([0, 0], [0, 1])), shape=(1, prim_pts)
                         )
                         additive = alpha * bcs[child][symbol.side][0]
+                        if hasattr(submesh, "length"):
+                            additive_multiplicative = submesh.length
+                        else:
+                            additive_multiplicative = pybamm.Scalar(1)
                         multiplicative = pybamm.Scalar(1)
 
                     else:
@@ -936,6 +948,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         )
 
                         additive = pybamm.Scalar(0)
+                        additive_multiplicative = pybamm.Scalar(1)
                         multiplicative = pybamm.Scalar(1)
 
                 else:
@@ -954,8 +967,10 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         additive = dxN * bcs[child][symbol.side][0]
                         if hasattr(submesh, "length"):
                             multiplicative = submesh.length
+                            additive_multiplicative = submesh.length
                         else:
                             multiplicative = pybamm.Scalar(1)
+                            additive_multiplicative = pybamm.Scalar(1)
                     else:
                         # to find value at x* use formula:
                         # f(x*) = f_N - (dxN / dxNm1) (f_N - f_Nm1)
@@ -967,6 +982,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             shape=(1, prim_pts),
                         )
                         additive = pybamm.Scalar(0)
+                        additive_multiplicative = pybamm.Scalar(1)
                         multiplicative = pybamm.Scalar(1)
                 elif extrap_order_value == "quadratic":
                     if use_bcs and pybamm.has_bc_of_form(
@@ -981,6 +997,10 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         )
 
                         additive = alpha * bcs[child][symbol.side][0]
+                        if hasattr(submesh, "length"):
+                            additive_multiplicative = submesh.length
+                        else:
+                            additive_multiplicative = pybamm.Scalar(1)
                         multiplicative = pybamm.Scalar(1)
                     else:
                         a = (
@@ -999,6 +1019,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             shape=(1, prim_pts),
                         )
                         additive = pybamm.Scalar(0)
+                        additive_multiplicative = pybamm.Scalar(1)
                         multiplicative = pybamm.Scalar(1)
                 else:
                     raise NotImplementedError
@@ -1008,6 +1029,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                 # just use the value from the bc: f'(x*)
                 sub_matrix = csr_matrix((1, prim_pts))
                 additive = bcs[child][symbol.side][0]
+                additive_multiplicative = pybamm.Scalar(1)
                 multiplicative = pybamm.Scalar(1)
 
             elif symbol.side == "left":
@@ -1017,6 +1039,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         ([-1, 1], ([0, 0], [0, 1])), shape=(1, prim_pts)
                     )
                     additive = pybamm.Scalar(0)
+                    additive_multiplicative = pybamm.Scalar(1)
                     if hasattr(submesh, "length"):
                         multiplicative = 1 / submesh.length
                     else:
@@ -1031,6 +1054,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         ([a, b, c], ([0, 0, 0], [0, 1, 2])), shape=(1, prim_pts)
                     )
                     additive = pybamm.Scalar(0)
+                    additive_multiplicative = pybamm.Scalar(1)
                     if hasattr(submesh, "length"):
                         multiplicative = 1 / submesh.length
                     else:
@@ -1048,6 +1072,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         shape=(1, prim_pts),
                     )
                     additive = pybamm.Scalar(0)
+                    additive_multiplicative = pybamm.Scalar(1)
                     if hasattr(submesh, "length"):
                         multiplicative = 1 / submesh.length
                     else:
@@ -1066,6 +1091,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         shape=(1, prim_pts),
                     )
                     additive = pybamm.Scalar(0)
+                    additive_multiplicative = pybamm.Scalar(1)
                     if hasattr(submesh, "length"):
                         multiplicative = 1 / submesh.length
                     else:
@@ -1086,7 +1112,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         boundary_value.copy_domains(symbol)
 
         additive.copy_domains(symbol)
-        boundary_value += additive
+        boundary_value += additive * additive_multiplicative
 
         return boundary_value
 

tests/unit/test_meshes/test_one_dimensional_submesh.py

@@ -89,6 +89,31 @@ def test_symmetric_mesh_creation_no_parameters(self, r, geometry):
         )
 
 
+class TestSymbolicUniform1DSubMesh:
+    def test_exceptions(self):
+        lims = {"a": 1, "b": 2}
+        with pytest.raises(pybamm.GeometryError):
+            pybamm.SymbolicUniform1DSubMesh(lims, None)
+
+    def test_symmetric_mesh_creation_no_parameters(self, r, geometry):
+        submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
+        var_pts = {r: 20}
+
+        # create mesh
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        # check boundary locations
+        assert mesh["negative particle"].edges[0] == 0
+        assert mesh["negative particle"].edges[-1] == 1
+
+        # check number of edges and nodes
+        assert len(mesh["negative particle"].nodes) == var_pts[r]
+        assert (
+            len(mesh["negative particle"].edges)
+            == len(mesh["negative particle"].nodes) + 1
+        )
+
+
 class TestExponential1DSubMesh:
     def test_symmetric_mesh_creation_no_parameters_even(self, r, geometry):
         submesh_params = {"side": "symmetric"}