From e4fff00dbed567efdb041ac918268989ffed048d Mon Sep 17 00:00:00 2001
From: Alexander Bills <alec@ionworks.com>
Date: Wed, 11 Dec 2024 11:19:22 -0800
Subject: [PATCH 01/13] symbolic mesh

---
 src/pybamm/__init__.py                        |  1 +
 src/pybamm/meshes/meshes.py                   | 47 ++++++++---
 .../meshes/one_dimensional_submeshes.py       | 20 +++++
 src/pybamm/parameters/parameter_values.py     |  4 -
 src/pybamm/spatial_methods/finite_volume.py   | 84 ++++++++++++++++---
 5 files changed, 132 insertions(+), 24 deletions(-)

diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py
index b466c3896b..e4ab341c1e 100644
--- a/src/pybamm/__init__.py
+++ b/src/pybamm/__init__.py
@@ -138,6 +138,7 @@
     Chebyshev1DSubMesh,
     UserSupplied1DSubMesh,
     SpectralVolume1DSubMesh,
+    SymbolicUniform1DSubMesh,
 )
 from .meshes.scikit_fem_submeshes import (
     ScikitSubMesh2D,
diff --git a/src/pybamm/meshes/meshes.py b/src/pybamm/meshes/meshes.py
index 3ec291b1b8..bde0430196 100644
--- a/src/pybamm/meshes/meshes.py
+++ b/src/pybamm/meshes/meshes.py
@@ -95,6 +95,8 @@ def __init__(self, geometry, submesh_types, var_pts):
                         if isinstance(sym, pybamm.Symbol):
                             try:
                                 sym_eval = sym.evaluate()
+                            except KeyError:
+                                sym_eval = sym
                             except NotImplementedError as error:
                                 if sym.has_symbol_of_classes(pybamm.Parameter):
                                     raise pybamm.DiscretisationError(
@@ -169,7 +171,12 @@ def combine_submeshes(self, *submeshnames):
         # Check that the final edge of each submesh is the same as the first edge of the
         # next submesh
         for i in range(len(submeshnames) - 1):
-            if self[submeshnames[i]].edges[-1] != self[submeshnames[i + 1]].edges[0]:
+            if self[submeshnames[i]].edges[-1] == self[submeshnames[i + 1]].edges[0]:
+                pass
+            elif hasattr(self[submeshnames[i + 1]], "min"):
+                # Circle back to this, but this should let the case of submesh i+1 be a symbolic submesh
+                pass
+            else:
                 raise pybamm.DomainError("submesh edges are not aligned")
 
             coord_sys = self[submeshnames[i]].coord_sys
@@ -184,10 +191,26 @@ def combine_submeshes(self, *submeshnames):
         )
         coord_sys = self[submeshnames[0]].coord_sys
         submesh = pybamm.SubMesh1D(combined_submesh_edges, coord_sys)
+        if getattr(self[submeshnames[0]], "length", None) is not None:
+            # Assume that the ghost cells have the same length as the first submesh
+            if any("ghost" in submeshname for submeshname in submeshnames):
+                submesh_min = self[submeshnames[0]].min
+                submesh_length = self[submeshnames[0]].length
+            # If not ghost cells, then the length is the sum of the lengths of the submeshes
+            else:
+                submesh_min = self[submeshnames[0]].min
+                submesh_length = sum(
+                    [self[submeshname].length for submeshname in submeshnames]
+                )
+            submesh.length = submesh_length
+            submesh.min = submesh_min
         # add in internal boundaries
-        submesh.internal_boundaries = [
-            self[submeshname].edges[0] for submeshname in submeshnames[1:]
-        ]
+        for submeshname in submeshnames[1:]:
+            if getattr(self[submeshname], "length", None) is not None:
+                min = self[submeshname].min
+            else:
+                min = 0
+            submesh.internal_boundaries.append(self[submeshname].edges[0] + min)
         return submesh
 
     def add_ghost_meshes(self):
@@ -210,16 +233,20 @@ def add_ghost_meshes(self):
 
             # left ghost cell: two edges, one node, to the left of existing submesh
             lgs_edges = np.array([2 * edges[0] - edges[1], edges[0]])
-            self[domain[0] + "_left ghost cell"] = pybamm.SubMesh1D(
-                lgs_edges, submesh.coord_sys
-            )
+            lgs_submesh = pybamm.SubMesh1D(lgs_edges, submesh.coord_sys)
+            if getattr(submesh, "length", None) is not None:
+                lgs_submesh.length = submesh.length
+                lgs_submesh.min = submesh.min
+            self[domain[0] + "_left ghost cell"] = lgs_submesh
 
             # right ghost cell: two edges, one node, to the right of
             # existing submesh
             rgs_edges = np.array([edges[-1], 2 * edges[-1] - edges[-2]])
-            self[domain[0] + "_right ghost cell"] = pybamm.SubMesh1D(
-                rgs_edges, submesh.coord_sys
-            )
+            rgs_submesh = pybamm.SubMesh1D(rgs_edges, submesh.coord_sys)
+            if getattr(submesh, "length", None) is not None:
+                rgs_submesh.length = submesh.length
+                rgs_submesh.min = submesh.min
+            self[domain[0] + "_right ghost cell"] = rgs_submesh
 
     @property
     def geometry(self):
diff --git a/src/pybamm/meshes/one_dimensional_submeshes.py b/src/pybamm/meshes/one_dimensional_submeshes.py
index 027c6d0421..612bc98cea 100644
--- a/src/pybamm/meshes/one_dimensional_submeshes.py
+++ b/src/pybamm/meshes/one_dimensional_submeshes.py
@@ -85,6 +85,26 @@ def to_json(self):
         return json_dict
 
 
+class SymbolicUniform1DSubMesh(SubMesh1D):
+    def __init__(self, lims, npts, tabs=None):
+        spatial_var, spatial_lims, tabs = self.read_lims(lims)
+        if tabs is not None:
+            raise ValueError("Tabs not supported for symbolic uniform submesh")
+        npts = npts[spatial_var.name]
+        length = spatial_lims["max"] - spatial_lims["min"]
+        x0 = spatial_lims["min"]
+        self.edges = np.linspace(0, 1, npts + 1)
+        self.length = length
+        self.min = x0
+        coord_sys = spatial_var.coord_sys
+        self.nodes = (self.edges[1:] + self.edges[:-1]) / 2
+        self.d_edges = np.diff(self.edges)
+        self.d_nodes = np.diff(self.nodes)
+        self.npts = self.nodes.size
+        self.coord_sys = coord_sys
+        self.internal_boundaries = []
+
+
 class Uniform1DSubMesh(SubMesh1D):
     """
     A class to generate a uniform submesh on a 1D domain
diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py
index 9ddbfd2fd3..c23cb0f6f7 100644
--- a/src/pybamm/parameters/parameter_values.py
+++ b/src/pybamm/parameters/parameter_values.py
@@ -604,10 +604,6 @@ def process_geometry(self, geometry):
 
         def process_and_check(sym):
             new_sym = self.process_symbol(sym)
-            if not isinstance(new_sym, pybamm.Scalar):
-                raise ValueError(
-                    "Geometry parameters must be Scalars after parameter processing"
-                )
             return new_sym
 
         for domain in geometry:
diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py
index 5b38b842ba..904b934114 100644
--- a/src/pybamm/spatial_methods/finite_volume.py
+++ b/src/pybamm/spatial_methods/finite_volume.py
@@ -64,7 +64,14 @@ def spatial_variable(self, symbol):
             entries = np.tile(symbol_mesh.edges, repeats)
         else:
             entries = np.tile(symbol_mesh.nodes, repeats)
-        return pybamm.Vector(entries, domains=symbol.domains)
+
+        if hasattr(symbol_mesh, "length"):
+            return (
+                pybamm.Vector(entries, domains=symbol.domains) * symbol_mesh.length
+                + symbol_mesh.min
+            )
+        else:
+            return pybamm.Vector(entries, domains=symbol.domains)
 
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the gradient operator.
@@ -128,8 +135,10 @@ def gradient_matrix(self, domain, domains):
         # Note that this makes column-slicing inefficient, but this should not be an
         # issue
         matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix))
-
-        return pybamm.Matrix(matrix)
+        if getattr(submesh, "length", None) is not None:
+            return pybamm.Matrix(matrix) * 1 / submesh.length
+        else:
+            return pybamm.Matrix(matrix)
 
     def divergence(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the divergence operator.
@@ -145,6 +154,8 @@ def divergence(self, symbol, discretised_symbol, boundary_conditions):
             # create np.array of repeated submesh.edges
             r_edges_numpy = np.kron(np.ones(second_dim_repeats), submesh.edges)
             r_edges = pybamm.Vector(r_edges_numpy)
+            if hasattr(submesh, "length"):
+                r_edges = r_edges * submesh.length
             if submesh.coord_sys == "spherical polar":
                 out = divergence_matrix @ ((r_edges**2) * discretised_symbol)
             elif submesh.coord_sys == "cylindrical polar":
@@ -197,7 +208,25 @@ def divergence_matrix(self, domains):
         # Note that this makes column-slicing inefficient, but this should not be an
         # issue
         matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix))
-        return pybamm.Matrix(matrix)
+        if getattr(submesh, "length", None) is not None:
+            if submesh.coord_sys == "spherical polar":
+                if submesh.min != 0:
+                    raise AssertionError(
+                        "Spherical polar symbolic meshes must have min r = 0"
+                    )
+                else:
+                    return pybamm.Matrix(matrix) * (1 / submesh.length**3)
+            elif submesh.coord_sys == "cylindrical polar":
+                if submesh.min != 0:
+                    raise AssertionError(
+                        "Cylindrical polar symbolic meshes must have min r = 0"
+                    )
+                else:
+                    return pybamm.Matrix(matrix) * (1 / (submesh.length**2))
+            else:
+                return pybamm.Matrix(matrix) * (1 / submesh.length)
+        else:
+            return pybamm.Matrix(matrix)
 
     def laplacian(self, symbol, discretised_symbol, boundary_conditions):
         """
@@ -569,7 +598,15 @@ def internal_neumann_condition(
 
         # Finite volume derivative
         # Remove domains to avoid clash
-        dx = right_mesh.nodes[0] - left_mesh.nodes[-1]
+        if hasattr(right_mesh, "length"):
+            right_mesh_x = right_mesh.min + (right_mesh.nodes[0] * right_mesh.length)
+        else:
+            right_mesh_x = right_mesh.nodes[0]
+        if hasattr(left_mesh, "length"):
+            left_mesh_x = left_mesh.min + (left_mesh.nodes[-1] * left_mesh.length)
+        else:
+            left_mesh_x = left_mesh.nodes[-1]
+        dx = right_mesh_x - left_mesh_x
         dy_r = (right_matrix / dx) @ right_symbol_disc
         dy_r.clear_domains()
         dy_l = (left_matrix / dx) @ left_symbol_disc
@@ -871,6 +908,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             shape=(1, prim_pts),
                         )
                         additive = pybamm.Scalar(0)
+                        multiplicative = pybamm.Scalar(1)
 
                 elif extrap_order_value == "quadratic":
                     if use_bcs and pybamm.has_bc_of_form(
@@ -884,6 +922,7 @@ 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]
+                        multiplicative = pybamm.Scalar(1)
 
                     else:
                         a = (dx0 + dx1) * (dx0 + dx1 + dx2) / (dx1 * (dx1 + dx2))
@@ -895,6 +934,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         )
 
                         additive = pybamm.Scalar(0)
+                        multiplicative = pybamm.Scalar(1)
+
                 else:
                     raise NotImplementedError
 
@@ -909,7 +950,10 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             ([1], ([0], [prim_pts - 1])), shape=(1, prim_pts)
                         )
                         additive = dxN * bcs[child][symbol.side][0]
-
+                        if hasattr(submesh, "length"):
+                            multiplicative = submesh.length
+                        else:
+                            multiplicative = pybamm.Scalar(1)
                     else:
                         # to find value at x* use formula:
                         # f(x*) = f_N - (dxN / dxNm1) (f_N - f_Nm1)
@@ -921,6 +965,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             shape=(1, prim_pts),
                         )
                         additive = pybamm.Scalar(0)
+                        multiplicative = pybamm.Scalar(1)
                 elif extrap_order_value == "quadratic":
                     if use_bcs and pybamm.has_bc_of_form(
                         child, symbol.side, bcs, "Neumann"
@@ -934,7 +979,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         )
 
                         additive = alpha * bcs[child][symbol.side][0]
-
+                        multiplicative = pybamm.Scalar(1)
                     else:
                         a = (
                             (dxN + dxNm1)
@@ -952,6 +997,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             shape=(1, prim_pts),
                         )
                         additive = pybamm.Scalar(0)
+                        multiplicative = pybamm.Scalar(1)
                 else:
                     raise NotImplementedError
 
@@ -960,6 +1006,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]
+                multiplicative = pybamm.Scalar(1)
 
             elif symbol.side == "left":
                 if extrap_order_gradient == "linear":
@@ -968,6 +1015,10 @@ 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)
+                    if hasattr(submesh, "length"):
+                        multiplicative = 1 / submesh.length
+                    else:
+                        multiplicative = pybamm.Scalar(1)
 
                 elif extrap_order_gradient == "quadratic":
                     a = -(2 * dx0 + 2 * dx1 + dx2) / (dx1**2 + dx1 * dx2)
@@ -978,6 +1029,11 @@ 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)
+                    if hasattr(submesh, "length"):
+                        multiplicative = 1 / submesh.length
+                    else:
+                        multiplicative = pybamm.Scalar(1)
+
                 else:
                     raise NotImplementedError
 
@@ -990,6 +1046,10 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         shape=(1, prim_pts),
                     )
                     additive = pybamm.Scalar(0)
+                    if hasattr(submesh, "length"):
+                        multiplicative = 1 / submesh.length
+                    else:
+                        multiplicative = pybamm.Scalar(1)
 
                 elif extrap_order_gradient == "quadratic":
                     a = (2 * dxN + 2 * dxNm1 + dxNm2) / (dxNm1**2 + dxNm1 * dxNm2)
@@ -1004,7 +1064,10 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         shape=(1, prim_pts),
                     )
                     additive = pybamm.Scalar(0)
-
+                    if hasattr(submesh, "length"):
+                        multiplicative = 1 / submesh.length
+                    else:
+                        multiplicative = pybamm.Scalar(1)
                 else:
                     raise NotImplementedError
 
@@ -1016,7 +1079,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         matrix = csr_matrix(kron(eye(repeats), sub_matrix))
 
         # Return boundary value with domain given by symbol
-        boundary_value = pybamm.Matrix(matrix) @ discretised_child
+        matrix = pybamm.Matrix(matrix) * multiplicative
+        boundary_value = matrix @ discretised_child
         boundary_value.copy_domains(symbol)
 
         additive.copy_domains(symbol)
@@ -1315,6 +1379,7 @@ def harmonic_mean(array):
                 sub_beta = (dx[:-1] / (dx[1:] + dx[:-1]))[:, np.newaxis]
                 beta = pybamm.Array(np.kron(np.ones((second_dim_repeats, 1)), sub_beta))
 
+                # dx_real = dx * length, therefore, beta is unchanged
                 # Compute harmonic mean on internal edges
                 # Note: add small number to denominator to regularise D_eff
                 D_eff = D1 * D2 / (D2 * beta + D1 * (1 - beta) + 1e-16)
@@ -1397,7 +1462,6 @@ def upwind_or_downwind(self, symbol, discretised_symbol, bcs, direction):
         direction : str
             Direction in which to apply the operator (upwind or downwind)
         """
-
         if symbol not in bcs:
             raise pybamm.ModelError(
                 "Boundary conditions must be provided for " f"{direction}ing '{symbol}'"

From 831f97eed5e900942c411196fd5788723403bb2c Mon Sep 17 00:00:00 2001
From: Alexander Bills <alec@ionworks.com>
Date: Wed, 11 Dec 2024 12:33:17 -0800
Subject: [PATCH 02/13] add a few more corrections, be smarter about errors

---
 .../meshes/one_dimensional_submeshes.py       | 12 +++++---
 src/pybamm/spatial_methods/finite_volume.py   | 28 ++++++++++---------
 2 files changed, 23 insertions(+), 17 deletions(-)

diff --git a/src/pybamm/meshes/one_dimensional_submeshes.py b/src/pybamm/meshes/one_dimensional_submeshes.py
index 612bc98cea..706031a5bf 100644
--- a/src/pybamm/meshes/one_dimensional_submeshes.py
+++ b/src/pybamm/meshes/one_dimensional_submeshes.py
@@ -88,15 +88,19 @@ def to_json(self):
 class SymbolicUniform1DSubMesh(SubMesh1D):
     def __init__(self, lims, npts, tabs=None):
         spatial_var, spatial_lims, tabs = self.read_lims(lims)
+        coord_sys = spatial_var.coord_sys
+        x0 = spatial_lims["min"]
         if tabs is not None:
-            raise ValueError("Tabs not supported for symbolic uniform submesh")
+            raise NotImplementedError("Tabs not supported for symbolic uniform submesh")
+        if coord_sys != "cartesian" and spatial_lims["min"] != 0:
+            raise NotImplementedError(
+                "Symbolic uniform submesh with non-cartesian coordinates and non-zero minimum not supported"
+            )
         npts = npts[spatial_var.name]
-        length = spatial_lims["max"] - spatial_lims["min"]
-        x0 = spatial_lims["min"]
+        length = spatial_lims["max"] - x0
         self.edges = np.linspace(0, 1, npts + 1)
         self.length = length
         self.min = x0
-        coord_sys = spatial_var.coord_sys
         self.nodes = (self.edges[1:] + self.edges[:-1]) / 2
         self.d_edges = np.diff(self.edges)
         self.d_nodes = np.diff(self.nodes)
diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py
index 904b934114..9891a191dd 100644
--- a/src/pybamm/spatial_methods/finite_volume.py
+++ b/src/pybamm/spatial_methods/finite_volume.py
@@ -210,19 +210,9 @@ def divergence_matrix(self, domains):
         matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix))
         if getattr(submesh, "length", None) is not None:
             if submesh.coord_sys == "spherical polar":
-                if submesh.min != 0:
-                    raise AssertionError(
-                        "Spherical polar symbolic meshes must have min r = 0"
-                    )
-                else:
-                    return pybamm.Matrix(matrix) * (1 / submesh.length**3)
+                return pybamm.Matrix(matrix) * (1 / submesh.length**3)
             elif submesh.coord_sys == "cylindrical polar":
-                if submesh.min != 0:
-                    raise AssertionError(
-                        "Cylindrical polar symbolic meshes must have min r = 0"
-                    )
-                else:
-                    return pybamm.Matrix(matrix) * (1 / (submesh.length**2))
+                return pybamm.Matrix(matrix) * (1 / (submesh.length**2))
             else:
                 return pybamm.Matrix(matrix) * (1 / submesh.length)
         else:
@@ -336,7 +326,15 @@ def definite_integral_matrix(
         # not supported by the default kron format
         # Note that this makes column-slicing inefficient, but this should not be an
         # issue
-        return pybamm.Matrix(csr_matrix(matrix))
+        matrix = pybamm.Matrix(csr_matrix(matrix))
+        if hasattr(submesh, "length"):
+            if submesh.coord_sys == "spherical polar":
+                matrix = matrix * submesh.length**3
+            elif submesh.coord_sys == "cylindrical polar":
+                matrix = matrix * submesh.length**2
+            else:
+                matrix = matrix * submesh.length
+        return matrix
 
     def indefinite_integral(self, child, discretised_child, direction):
         """Implementation of the indefinite integral operator."""
@@ -463,6 +461,8 @@ def indefinite_integral_matrix_edges(self, domains, direction):
         # add a column of zeros at each end
         zero_col = csr_matrix((n, 1))
         sub_matrix = hstack([zero_col, sub_matrix, zero_col])
+        if hasattr(submesh, "length"):
+            sub_matrix = sub_matrix * submesh.length
         # Convert to csr_matrix so that we can take the index (row-slicing), which is
         # not supported by the default kron format
         # Note that this makes column-slicing inefficient, but this should not be an
@@ -508,6 +508,8 @@ def indefinite_integral_matrix_nodes(self, domains, direction):
         # Note that this makes column-slicing inefficient, but this should not be an
         # issue
         matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix))
+        if hasattr(submesh, "length"):
+            matrix = matrix * submesh.length
 
         return pybamm.Matrix(matrix)
 

From d468be4f90022c652b00fad3a0d4dc2af26082e7 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alec@ionworks.com>
Date: Wed, 11 Dec 2024 13:14:00 -0800
Subject: [PATCH 03/13] fix bug, reimplement value error, make sure
 multiplicative is always defined

---
 src/pybamm/parameters/parameter_values.py     |  6 ++++
 src/pybamm/spatial_methods/finite_volume.py   | 28 ++++++++++++++++++-
 .../test_one_dimensional_submesh.py           | 25 +++++++++++++++++
 3 files changed, 58 insertions(+), 1 deletion(-)

diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py
index c23cb0f6f7..ce26335347 100644
--- a/src/pybamm/parameters/parameter_values.py
+++ b/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:
diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py
index 9891a191dd..bd22409fc3 100644
--- a/src/pybamm/spatial_methods/finite_volume.py
+++ b/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,6 +905,11 @@ 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(
@@ -910,6 +917,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":
@@ -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
 
diff --git a/tests/unit/test_meshes/test_one_dimensional_submesh.py b/tests/unit/test_meshes/test_one_dimensional_submesh.py
index 94895b6c4b..65462fa085 100644
--- a/tests/unit/test_meshes/test_one_dimensional_submesh.py
+++ b/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"}

From 1c82b6c2ada1e6f1c305314e227e6b1e19e12c48 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alec@ionworks.com>
Date: Thu, 12 Dec 2024 08:42:05 -0800
Subject: [PATCH 04/13] begin working on tests

---
 src/pybamm/parameters/parameter_values.py            | 12 ++++++------
 .../unit/test_meshes/test_one_dimensional_submesh.py |  4 ++++
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py
index ce26335347..8a0c1d62a1 100644
--- a/src/pybamm/parameters/parameter_values.py
+++ b/src/pybamm/parameters/parameter_values.py
@@ -604,12 +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"
-                )
+            #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:
diff --git a/tests/unit/test_meshes/test_one_dimensional_submesh.py b/tests/unit/test_meshes/test_one_dimensional_submesh.py
index 65462fa085..7349327c85 100644
--- a/tests/unit/test_meshes/test_one_dimensional_submesh.py
+++ b/tests/unit/test_meshes/test_one_dimensional_submesh.py
@@ -94,6 +94,10 @@ def test_exceptions(self):
         lims = {"a": 1, "b": 2}
         with pytest.raises(pybamm.GeometryError):
             pybamm.SymbolicUniform1DSubMesh(lims, None)
+        tabs = {"negative": {"z_centre": 0}, "positive": {"z_centre": 1}}
+
+        with pytest.raises(NotImplementedError):
+            pybamm.SymbolicUniform1DSubMesh(lims, 20, tabs=tabs)
 
     def test_symmetric_mesh_creation_no_parameters(self, r, geometry):
         submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}

From 76a77f0c6defccc4a1b7d7bcc53d279122558c2b Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 12 Dec 2024 16:42:35 +0000
Subject: [PATCH 05/13] style: pre-commit fixes

---
 src/pybamm/parameters/parameter_values.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py
index 8a0c1d62a1..5ae13cc32f 100644
--- a/src/pybamm/parameters/parameter_values.py
+++ b/src/pybamm/parameters/parameter_values.py
@@ -604,9 +604,9 @@ 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(
+            # 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"
             #    )

From 037d4b4699950c240d75819d8337261b60b4bfbb Mon Sep 17 00:00:00 2001
From: Alexander Bills <alec@ionworks.com>
Date: Fri, 20 Dec 2024 13:57:51 -0800
Subject: [PATCH 06/13] test mesh

---
 .../meshes/one_dimensional_submeshes.py       |  2 +-
 .../test_one_dimensional_submesh.py           | 49 +++++++++++++++++--
 .../test_parameters/test_parameter_values.py  |  6 ---
 3 files changed, 46 insertions(+), 11 deletions(-)

diff --git a/src/pybamm/meshes/one_dimensional_submeshes.py b/src/pybamm/meshes/one_dimensional_submeshes.py
index 706031a5bf..0e95b2ac59 100644
--- a/src/pybamm/meshes/one_dimensional_submeshes.py
+++ b/src/pybamm/meshes/one_dimensional_submeshes.py
@@ -93,7 +93,7 @@ def __init__(self, lims, npts, tabs=None):
         if tabs is not None:
             raise NotImplementedError("Tabs not supported for symbolic uniform submesh")
         if coord_sys != "cartesian" and spatial_lims["min"] != 0:
-            raise NotImplementedError(
+            raise pybamm.GeometryError(
                 "Symbolic uniform submesh with non-cartesian coordinates and non-zero minimum not supported"
             )
         npts = npts[spatial_var.name]
diff --git a/tests/unit/test_meshes/test_one_dimensional_submesh.py b/tests/unit/test_meshes/test_one_dimensional_submesh.py
index 7349327c85..e50c911150 100644
--- a/tests/unit/test_meshes/test_one_dimensional_submesh.py
+++ b/tests/unit/test_meshes/test_one_dimensional_submesh.py
@@ -11,6 +11,13 @@ def r():
     return r
 
 
+@pytest.fixture()
+def x():
+    return pybamm.SpatialVariable(
+        "x", domain=["negative electrode"], coord_sys="cartesian"
+    )
+
+
 @pytest.fixture()
 def geometry(r):
     geometry = {
@@ -90,18 +97,34 @@ def test_symmetric_mesh_creation_no_parameters(self, r, geometry):
 
 
 class TestSymbolicUniform1DSubMesh:
-    def test_exceptions(self):
+    def test_exceptions(self, r):
         lims = {"a": 1, "b": 2}
         with pytest.raises(pybamm.GeometryError):
             pybamm.SymbolicUniform1DSubMesh(lims, None)
+        lims = {"x_n": {"min": 0, "max": 1}}
+        npts = {"x_n": 10}
         tabs = {"negative": {"z_centre": 0}, "positive": {"z_centre": 1}}
+        lims["tabs"] = tabs
 
         with pytest.raises(NotImplementedError):
-            pybamm.SymbolicUniform1DSubMesh(lims, 20, tabs=tabs)
+            pybamm.SymbolicUniform1DSubMesh(lims, npts, tabs=tabs)
 
-    def test_symmetric_mesh_creation_no_parameters(self, r, geometry):
-        submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
+        submesh_types = {"negative particle": pybamm.SymbolicUniform1DSubMesh}
         var_pts = {r: 20}
+        geometry = {
+            "negative particle": {
+                r: {"min": pybamm.InputParameter("min"), "max": pybamm.Scalar(2)}
+            }
+        }
+        with pytest.raises(pybamm.GeometryError):
+            pybamm.Mesh(geometry, submesh_types, var_pts)
+
+    def test_mesh_creation(self, r, x):
+        submesh_types = {"negative particle": pybamm.SymbolicUniform1DSubMesh}
+        var_pts = {r: 20}
+        geometry = {
+            "negative particle": {r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(2)}}
+        }
 
         # create mesh
         mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
@@ -110,6 +133,10 @@ def test_symmetric_mesh_creation_no_parameters(self, r, geometry):
         assert mesh["negative particle"].edges[0] == 0
         assert mesh["negative particle"].edges[-1] == 1
 
+        # check scaling and min/max
+        assert mesh["negative particle"].length == 2
+        assert mesh["negative particle"].min == 0
+
         # check number of edges and nodes
         assert len(mesh["negative particle"].nodes) == var_pts[r]
         assert (
@@ -117,6 +144,20 @@ def test_symmetric_mesh_creation_no_parameters(self, r, geometry):
             == len(mesh["negative particle"].nodes) + 1
         )
 
+        # Check that length and min are scaled correctly
+        submesh_types = {"negative electrode": pybamm.SymbolicUniform1DSubMesh}
+        var_pts = {x: 20}
+        geometry = {
+            "negative electrode": {
+                x: {"min": pybamm.InputParameter("min"), "max": pybamm.Scalar(2)}
+            }
+        }
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+        assert mesh["negative electrode"].length == pybamm.Scalar(
+            2
+        ) - pybamm.InputParameter("min")
+        assert mesh["negative electrode"].min == pybamm.InputParameter("min")
+
 
 class TestExponential1DSubMesh:
     def test_symmetric_mesh_creation_no_parameters_even(self, r, geometry):
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 0e32f60476..52cc00628c 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -914,12 +914,6 @@ def test_process_complex_expression(self):
         exp_param = param.process_symbol(expression)
         assert exp_param == 3.0 * (2.0**var2) / ((-4.0 + var1) + var2)
 
-    def test_process_geometry(self):
-        var = pybamm.Variable("var")
-        geometry = {"negative electrode": {"x": {"min": 0, "max": var}}}
-        with pytest.raises(ValueError, match="Geometry parameters must be Scalars"):
-            pybamm.ParameterValues({}).process_geometry(geometry)
-
     def test_process_model(self):
         model = pybamm.BaseModel()
         a = pybamm.Parameter("a")

From 1a606b86e18bf46c3d1e4a897a94d010c2a99b86 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alec@ionworks.com>
Date: Mon, 23 Dec 2024 15:23:41 -0600
Subject: [PATCH 07/13] test meshes

---
 tests/unit/test_meshes/test_meshes.py | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py
index eecc1a8911..50822fa91e 100644
--- a/tests/unit/test_meshes/test_meshes.py
+++ b/tests/unit/test_meshes/test_meshes.py
@@ -233,6 +233,32 @@ def test_ghost_cells(self, submesh_types):
             mesh["positive electrode"].edges[-1],
         )
 
+    def test_symbolic_mesh_ghost_cells(self, submesh_types):
+        param = get_param()
+
+        submesh_types.update({"negative electrode": pybamm.SymbolicUniform1DSubMesh})
+
+        geometry = pybamm.battery_geometry()
+        param.process_geometry(geometry)
+
+        var_pts = {"x_n": 10, "x_s": 10, "x_p": 12, "r_n": 5, "r_p": 6}
+
+        # create mesh
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        lgs_submesh = mesh["negative electrode_left ghost cell"]
+        rgs_submesh = mesh["negative electrode_right ghost cell"]
+        submesh = mesh["negative electrode"]
+
+        np.testing.assert_array_equal(
+            lgs_submesh.edges[1] * lgs_submesh.length + lgs_submesh.min,
+            submesh.edges[0],
+        )
+        np.testing.assert_array_equal(
+            rgs_submesh.edges[0] * rgs_submesh.length + rgs_submesh.min,
+            submesh.edges[-1] * submesh.length + submesh.min,
+        )
+
     def test_mesh_coord_sys(self, submesh_types):
         param = get_param()
 

From e8e628f1fcfa6a42a4c10133c1bfa19153501069 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alec@ionworks.com>
Date: Mon, 30 Dec 2024 11:37:53 -0800
Subject: [PATCH 08/13] add test graddiv

---
 .../test_grad_div_shapes.py                   | 44 +++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py b/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
index 9e7f993e2a..78a8772707 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
@@ -12,6 +12,16 @@
 import numpy as np
 
 
+def get_mesh_for_testing_symbolic():
+    submesh_types = {"domain": pybamm.SymbolicUniform1DSubMesh}
+    geometry = {
+        "domain": {"x": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(2)}},
+    }
+    var_pts = {"x": 15}
+    mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+    return mesh
+
+
 class TestFiniteVolumeGradDiv:
     def test_grad_div_shapes_Dirichlet_bcs(self):
         """
@@ -636,3 +646,37 @@ def test_grad_1plus1d(self):
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y), expected
         )
+
+    def test_grad_div_shapes_symbolic_mesh(self):
+        mesh = get_mesh_for_testing_symbolic()
+        spatial_methods = {"domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain="domain")
+        grad_eqn = pybamm.grad(var)
+        div_eqn = pybamm.div(grad_eqn)
+        boundary_conditions = {
+            var: {
+                "left": (pybamm.Scalar(1), "Neumann"),
+                "right": (pybamm.Scalar(1), "Neumann"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        disc.set_variable_slices([var])
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+        div_eqn_disc = disc.process_symbol(div_eqn)
+
+        # Evaluate grad
+        dom = ("domain_left ghost cell", "domain", "domain_right ghost cell")
+        linear_y = mesh[dom].nodes * mesh[dom].length + mesh[dom].min
+        expected = np.ones((16, 1))
+        np.testing.assert_array_almost_equal(
+            grad_eqn_disc.evaluate(None, linear_y), expected
+        )
+
+        # Evaluate div
+        div_eqn = pybamm.div(pybamm.grad(var))
+        div_eqn_disc = disc.process_symbol(div_eqn)
+        div_eval = div_eqn_disc.evaluate(None, linear_y)
+        div_eval = np.reshape(div_eval, [15, 1])
+        np.testing.assert_array_almost_equal(div_eval, np.zeros([15, 1]))

From 64bb7b595af4a4e77c6d422444f8f8d57ee52483 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Mon, 30 Dec 2024 12:44:56 -0800
Subject: [PATCH 09/13] finish testing meshes

---
 tests/unit/test_meshes/test_meshes.py | 60 +++++++++++++++++++++++++++
 1 file changed, 60 insertions(+)

diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py
index 50822fa91e..b51f902fd0 100644
--- a/tests/unit/test_meshes/test_meshes.py
+++ b/tests/unit/test_meshes/test_meshes.py
@@ -208,6 +208,28 @@ def test_combine_submeshes(self, submesh_types):
         ):
             mesh.combine_submeshes()
 
+        # test symbolic submesh
+        new_submesh_types = submesh_types.copy()
+        geometry = {
+            "negative electrode": {
+                "x_n": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
+            },
+            "separator": {
+                "x_s": {"min": pybamm.Scalar(1), "max": pybamm.InputParameter("L")}
+            },
+        }
+        param.process_geometry(geometry)
+        for k in new_submesh_types:
+            new_submesh_types[k] = pybamm.SymbolicUniform1DSubMesh
+        mesh = pybamm.Mesh(geometry, new_submesh_types, var_pts)
+        submesh = mesh[("negative electrode", "separator")]
+        mesh.combine_submeshes("negative electrode", "separator")
+        assert (
+            submesh.length
+            == mesh["separator"].length + mesh["negative electrode"].length
+        )
+        assert submesh.min == mesh["negative electrode"].min
+
     def test_ghost_cells(self, submesh_types):
         param = get_param()
 
@@ -233,6 +255,44 @@ def test_ghost_cells(self, submesh_types):
             mesh["positive electrode"].edges[-1],
         )
 
+        # test symbolic mesh
+        geometry = {
+            "negative electrode": {
+                "x_n": {"min": pybamm.Scalar(0), "max": pybamm.InputParameter("L_n")}
+            },
+            "separator": {
+                "x_s": {
+                    "min": pybamm.InputParameter("L_n"),
+                    "max": pybamm.InputParameter("L"),
+                }
+            },
+        }
+        submesh_types = {
+            "negative electrode": pybamm.SymbolicUniform1DSubMesh,
+            "separator": pybamm.SymbolicUniform1DSubMesh,
+        }
+        var_pts = {
+            "x_n": 10,
+            "x_s": 10,
+        }
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges[1],
+            mesh["negative electrode"].edges[0],
+        )
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges[0],
+            -mesh["negative electrode"].edges[1],
+        )
+        np.testing.assert_array_equal(
+            mesh["separator_left ghost cell"].edges[1],
+            mesh["separator"].edges[0],
+        )
+        np.testing.assert_array_equal(
+            mesh["separator_left ghost cell"].edges[0],
+            -mesh["separator"].edges[1],
+        )
+
     def test_symbolic_mesh_ghost_cells(self, submesh_types):
         param = get_param()
 

From c63a44570ff84f92e834fb245210dd322dc680a2 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Mon, 30 Dec 2024 15:48:29 -0800
Subject: [PATCH 10/13] add grad/div tests for spherical and cylindrical

---
 tests/__init__.py                             |   3 +
 tests/shared.py                               |  35 ++++++
 .../test_finite_volume/test_finite_volume.py  |  13 ++
 .../test_grad_div_shapes.py                   | 114 ++++++++++++++++--
 4 files changed, 155 insertions(+), 10 deletions(-)

diff --git a/tests/__init__.py b/tests/__init__.py
index 9ca19981bb..8227cd0650 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -24,11 +24,14 @@
 
 from .shared import (
     get_mesh_for_testing,
+    get_mesh_for_testing_symbolic,
     get_p2d_mesh_for_testing,
     get_size_distribution_mesh_for_testing,
     get_1p1d_mesh_for_testing,
     get_2p1d_mesh_for_testing,
     get_cylindrical_mesh_for_testing,
+    get_spherical_mesh_for_testing_symbolic,
+    get_cylindrical_mesh_for_testing_symbolic,
     get_discretisation_for_testing,
     get_p2d_discretisation_for_testing,
     get_size_distribution_disc_for_testing,
diff --git a/tests/shared.py b/tests/shared.py
index 15c8d428a6..1af9152cd3 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -345,3 +345,38 @@ def assert_domain_equal(a, b):
     a_dict = {k: v for k, v in a.items() if v != []}
     b_dict = {k: v for k, v in b.items() if v != []}
     assert a_dict == b_dict
+
+
+def get_mesh_for_testing_symbolic():
+    submesh_types = {"domain": pybamm.SymbolicUniform1DSubMesh}
+    geometry = {
+        "domain": {"x": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(2)}},
+    }
+    var_pts = {"x": 15}
+    mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+    return mesh
+
+
+def get_spherical_mesh_for_testing_symbolic():
+    submesh_types = {"spherical domain": pybamm.SymbolicUniform1DSubMesh}
+    geometry = {
+        "spherical domain": {"r_n": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(2)}},
+    }
+    var_pts = {"r_n": 15}
+    mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+    return mesh
+
+
+def get_cylindrical_mesh_for_testing_symbolic():
+    submesh_types = {"cylindrical domain": pybamm.SymbolicUniform1DSubMesh}
+    cylindrical_r = pybamm.SpatialVariable(
+        "r", ["cylindrical domain"], coord_sys="cylindrical polar"
+    )
+    geometry = {
+        "cylindrical domain": {
+            cylindrical_r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(2)}
+        },
+    }
+    var_pts = {cylindrical_r: 15}
+    mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+    return mesh
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 204e831855..7266843e39 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -7,6 +7,7 @@
     get_mesh_for_testing,
     get_p2d_mesh_for_testing,
     get_1p1d_mesh_for_testing,
+    get_mesh_for_testing_symbolic,
 )
 import numpy as np
 from scipy.sparse import kron, eye
@@ -174,6 +175,18 @@ def test_discretise_spatial_variable(self):
             r_disc.evaluate(), 3 * disc.mesh["negative particle"].nodes[:, np.newaxis]
         )
 
+    def test_spatial_variable_symbolic(self):
+        mesh = get_mesh_for_testing_symbolic()
+        spatial_methods = {"domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        x = pybamm.SpatialVariable("x", ["domain"])
+        x_disc = disc.process_symbol(x)
+        assert isinstance(x_disc, pybamm.Vector)
+        np.testing.assert_array_equal(
+            x_disc.evaluate(),
+            mesh["domain"].nodes[:, np.newaxis] * mesh["domain"].length,
+        )
+
     def test_mass_matrix_shape(self):
         # Create model
         whole_cell = ["negative electrode", "separator", "positive electrode"]
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py b/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
index 78a8772707..41a4ddee44 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
@@ -8,20 +8,13 @@
     get_p2d_mesh_for_testing,
     get_1p1d_mesh_for_testing,
     get_cylindrical_mesh_for_testing,
+    get_mesh_for_testing_symbolic,
+    get_spherical_mesh_for_testing_symbolic,
+    get_cylindrical_mesh_for_testing_symbolic,
 )
 import numpy as np
 
 
-def get_mesh_for_testing_symbolic():
-    submesh_types = {"domain": pybamm.SymbolicUniform1DSubMesh}
-    geometry = {
-        "domain": {"x": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(2)}},
-    }
-    var_pts = {"x": 15}
-    mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-    return mesh
-
-
 class TestFiniteVolumeGradDiv:
     def test_grad_div_shapes_Dirichlet_bcs(self):
         """
@@ -440,6 +433,55 @@ def test_cylindrical_grad_div_shapes_Neumann_bcs(self):
             4 * np.ones((npts, 1)),
         )
 
+    def test_cylindrical_grad_div_shapes_Neumann_bcs_symbolic(self):
+        mesh = get_cylindrical_mesh_for_testing_symbolic()
+        spatial_methods = {"cylindrical domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # Test gradient
+        var = pybamm.Variable("var", domain="cylindrical domain")
+        disc.set_variable_slices([var])
+        grad_eqn = pybamm.grad(var)
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+        constant_y = np.ones_like(mesh["cylindrical domain"].nodes[:, np.newaxis])
+        np.testing.assert_array_equal(
+            grad_eqn_disc.evaluate(None, constant_y),
+            np.zeros((14, 1)),
+        )
+
+        # Test divergence
+        boundary_conditions = {
+            var: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        div_eqn_disc = disc.process_symbol(pybamm.div(grad_eqn))
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(None, constant_y),
+            np.zeros((15, 1)),
+        )
+
+        # Test divergence of gradient
+        # div(grad(r^2)) = 4, N_left = 2*r_inner, N_right = 2
+        submesh = mesh["cylindrical domain"]
+        y_squared = (submesh.nodes * submesh.length) ** 2
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            var: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(4), "Neumann"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        div_eqn_disc = disc.process_symbol(div_eqn)
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(None, y_squared),
+            4 * np.ones((15, 1)),
+        )
+
     def test_spherical_grad_div_shapes_Neumann_bcs(self):
         """
         Test grad and div with Neumann boundary conditions spherical polar
@@ -502,6 +544,58 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
             6 * np.ones((submesh.npts, 1)),
         )
 
+    def test_spherical_grad_div_shapes_Neumann_bcs_symbolic(self):
+        mesh = get_spherical_mesh_for_testing_symbolic()
+        spatial_methods = {"spherical domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # Test gradient
+        var = pybamm.Variable("var", domain="spherical domain")
+        disc.set_variable_slices([var])
+        grad_eqn = pybamm.grad(var)
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+        constant_y = np.ones_like(mesh["spherical domain"].nodes[:, np.newaxis])
+        np.testing.assert_array_equal(
+            grad_eqn_disc.evaluate(None, constant_y),
+            np.zeros((14, 1)),
+        )
+
+        # Test divergence
+        boundary_conditions = {
+            var: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        disc.set_variable_slices([var])
+        div_eqn_disc = disc.process_symbol(pybamm.div(grad_eqn))
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(
+                None, np.ones_like(mesh["spherical domain"].nodes[:, np.newaxis])
+            ),
+            np.zeros_like(mesh["spherical domain"].nodes[:, np.newaxis]),
+        )
+
+        # Test divergence of gradient
+        # div(grad(r^2)) = 6, N_left = 0, N_right = 2
+        submesh = mesh["spherical domain"]
+        quadratic_y = (submesh.nodes * submesh.length) ** 2
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            var: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(4), "Neumann"),
+            }
+        }
+        disc.bcs = boundary_conditions
+        div_eqn_disc = disc.process_symbol(div_eqn)
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(None, quadratic_y),
+            6 * np.ones_like(mesh["spherical domain"].nodes[:, np.newaxis]),
+        )
+
     def test_p2d_spherical_grad_div_shapes_Neumann_bcs(self):
         """
         Test grad and div with Neumann boundary conditions in the pseudo

From f1ebbf5700b7afdfe42573bdd8111f1adcfadd44 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Tue, 31 Dec 2024 12:07:37 -0800
Subject: [PATCH 11/13] add tests for integrals

---
 src/pybamm/spatial_methods/finite_volume.py   |  4 +-
 .../test_discretisation.py                    | 38 ++++++++
 .../test_finite_volume/test_integration.py    | 87 +++++++++++++++++++
 3 files changed, 127 insertions(+), 2 deletions(-)

diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py
index bd22409fc3..9ea6091920 100644
--- a/src/pybamm/spatial_methods/finite_volume.py
+++ b/src/pybamm/spatial_methods/finite_volume.py
@@ -461,13 +461,13 @@ def indefinite_integral_matrix_edges(self, domains, direction):
         # add a column of zeros at each end
         zero_col = csr_matrix((n, 1))
         sub_matrix = hstack([zero_col, sub_matrix, zero_col])
-        if hasattr(submesh, "length"):
-            sub_matrix = sub_matrix * submesh.length
         # Convert to csr_matrix so that we can take the index (row-slicing), which is
         # not supported by the default kron format
         # Note that this makes column-slicing inefficient, but this should not be an
         # issue
         matrix = csr_matrix(kron(eye(second_dim_repeats), sub_matrix))
+        if hasattr(submesh, "length"):
+            matrix = matrix * submesh.length
 
         return pybamm.Matrix(matrix)
 
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 78b6cc0d59..55d6546bd8 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -68,6 +68,44 @@ def test_add_internal_boundary_conditions(self):
         for child in c_e.children:
             assert child in disc.bcs.keys()
 
+    def test_add_internal_boundary_conditions_symbolic(self):
+        submesh_types = {
+            "left domain": pybamm.SymbolicUniform1DSubMesh,
+            "middle domain": pybamm.SymbolicUniform1DSubMesh,
+            "right domain": pybamm.SymbolicUniform1DSubMesh,
+        }
+        spatial_methods = {
+            "left domain": pybamm.FiniteVolume(),
+            "middle domain": pybamm.FiniteVolume(),
+            "right domain": pybamm.FiniteVolume(),
+        }
+        x_l = pybamm.SpatialVariable("x_l", ["left domain"], coord_sys="cartesian")
+        x_m = pybamm.SpatialVariable("x_m", ["middle domain"], coord_sys="cartesian")
+        x_r = pybamm.SpatialVariable("x_r", ["right domain"], coord_sys="cartesian")
+        geometry = {
+            "left domain": {x_l: {"min": 0, "max": 1}},
+            "middle domain": {x_m: {"min": 1, "max": 2}},
+            "right domain": {x_r: {"min": 2, "max": 3}},
+        }
+        var_pts = {x_l: 10, x_m: 10, x_r: 10}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        model = pybamm.BaseModel()
+        c_e_n = pybamm.Variable("c_e_n", domain=["left domain"])
+        c_e_s = pybamm.Variable("c_e_s", domain=["middle domain"])
+        c_e_p = pybamm.Variable("c_e_p", domain=["right domain"])
+        c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p)
+        lbc = (pybamm.Scalar(0), "Neumann")
+        rbc = (pybamm.Scalar(0), "Neumann")
+        model.boundary_conditions = {c_e: {"left": lbc, "right": rbc}}
+        disc.bcs = model.boundary_conditions
+        disc.set_variable_slices([c_e_n, c_e_s, c_e_p])
+        disc.set_internal_boundary_conditions(model)
+
+        for child in c_e.children:
+            assert child in disc.bcs.keys()
+
     def test_discretise_slicing(self):
         # create discretisation
         mesh = get_mesh_for_testing()
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py b/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py
index bf6f44059a..3abfa6466c 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py
@@ -8,6 +8,9 @@
     get_mesh_for_testing,
     get_1p1d_mesh_for_testing,
     get_cylindrical_mesh_for_testing,
+    get_mesh_for_testing_symbolic,
+    get_cylindrical_mesh_for_testing_symbolic,
+    get_spherical_mesh_for_testing_symbolic,
 )
 import numpy as np
 
@@ -128,6 +131,48 @@ def test_definite_integral(self):
         ):
             finite_volume.definite_integral_matrix(var, "column", "secondary")
 
+    def test_definite_integral_symbolic_mesh(self):
+        mesh = get_mesh_for_testing_symbolic()
+        spatial_methods = {"domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="domain")
+        x = pybamm.SpatialVariable("x", "domain")
+        integral_eqn = pybamm.Integral(var, x)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+        const_y = np.ones_like(mesh["domain"].nodes[:, np.newaxis])
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc.evaluate(None, const_y), 2.0, decimal=4
+        )
+
+        mesh = get_cylindrical_mesh_for_testing_symbolic()
+        spatial_methods = {"cylindrical domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="cylindrical domain")
+        r = pybamm.SpatialVariable(
+            "r", "cylindrical domain", coord_sys="cylindrical polar"
+        )
+        integral_eqn = pybamm.Integral(var, r)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+        const_y = np.ones_like(mesh["cylindrical domain"].nodes[:, np.newaxis])
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc.evaluate(None, const_y), 2 * np.pi * 2.0, decimal=4
+        )
+
+        mesh = get_spherical_mesh_for_testing_symbolic()
+        spatial_methods = {"spherical domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="spherical domain")
+        r = pybamm.SpatialVariable("r", "spherical domain", coord_sys="spherical polar")
+        integral_eqn = pybamm.Integral(var, r)
+        disc.set_variable_slices([var])
+        integral_eqn_disc = disc.process_symbol(integral_eqn)
+        const_y = np.ones_like(mesh["spherical domain"].nodes[:, np.newaxis])
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc.evaluate(None, const_y), 4 * np.pi * 2.0**3 / 3, decimal=4
+        )
+
     def test_integral_secondary_domain(self):
         # create discretisation
         mesh = get_1p1d_mesh_for_testing()
@@ -543,6 +588,31 @@ def test_indefinite_integral_of_broadcasted_to_cell_edges(self):
         phi_approx = int_int_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(x_end**3 / 6, phi_approx, decimal=4)
 
+    def test_indefinite_integral_on_edges_symbolic(self):
+        mesh = get_mesh_for_testing_symbolic()
+        spatial_methods = {"domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # make a variable 'phi' and a vector 'i' which is broadcast onto edges
+        # the integral of this should then be put onto the nodes
+        phi = pybamm.Variable("phi", domain=["domain"])
+        i = pybamm.PrimaryBroadcastToEdges(1, phi.domain)
+        x = pybamm.SpatialVariable("x", phi.domain)
+        disc.set_variable_slices([phi])
+        submesh = mesh["domain"]
+        x_end = submesh.edges[-1] * submesh.length + submesh.min
+
+        # take indefinite integral
+        int_phi = pybamm.IndefiniteIntegral(i * phi, x)
+        # take integral again
+        int_int_phi = pybamm.Integral(int_phi, x)
+        int_int_phi_disc = disc.process_symbol(int_int_phi)
+
+        # constant case
+        phi_exact = np.ones_like(submesh.nodes)
+        phi_approx = int_int_phi_disc.evaluate(None, phi_exact)
+        np.testing.assert_array_almost_equal(x_end**2 / 2, phi_approx)
+
     def test_indefinite_integral_on_nodes(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
@@ -589,6 +659,23 @@ def test_indefinite_integral_on_nodes(self):
         ):
             disc.process_symbol(int_c)
 
+    def test_indefinite_integral_on_nodes_symbolic(self):
+        mesh = get_mesh_for_testing_symbolic()
+        spatial_methods = {"domain": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        phi = pybamm.Variable("phi", domain=["domain"])
+        x = pybamm.SpatialVariable("x", ["domain"])
+        int_phi = pybamm.IndefiniteIntegral(phi, x)
+        disc.set_variable_slices([phi])
+        int_phi_disc = disc.process_symbol(int_phi)
+
+        submesh = mesh["domain"]
+        constant_y = np.ones((submesh.npts, 1))
+        phi_exact = submesh.edges * submesh.length + submesh.min
+        phi_approx = int_phi_disc.evaluate(None, constant_y).flatten()
+        np.testing.assert_array_almost_equal(phi_exact, phi_approx)
+
     def test_backward_indefinite_integral_on_nodes(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}

From 97b9b1c97a7d9a36a5e96d6209e01361c57a499c Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Tue, 31 Dec 2024 12:54:07 -0800
Subject: [PATCH 12/13] add extrapolation test

---
 .../test_finite_volume/test_extrapolation.py  | 90 +++++++++++++++++++
 1 file changed, 90 insertions(+)

diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py b/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py
index 5a8253b09d..6461dd4502 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py
@@ -7,6 +7,7 @@
     get_mesh_for_testing,
     get_p2d_mesh_for_testing,
     get_1p1d_mesh_for_testing,
+    get_mesh_for_testing_symbolic,
 )
 import numpy as np
 
@@ -332,6 +333,95 @@ def test_linear_extrapolate_left_right(self):
             surf_eqn_disc.evaluate(None, linear_y), y_surf
         )
 
+    def test_extrapolate_symbolic(self):
+        mesh = get_mesh_for_testing_symbolic()
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": True,
+            }
+        }
+        spatial_methods = {
+            "domain": pybamm.FiniteVolume(method_options),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain="domain")
+        extrap_left = pybamm.BoundaryValue(var, "left")
+        extrap_right = pybamm.BoundaryValue(var, "right")
+        extrap_grad_left = pybamm.BoundaryGradient(var, "left")
+        extrap_grad_right = pybamm.BoundaryGradient(var, "right")
+        disc.set_variable_slices([var])
+        extrap_left_disc = disc.process_symbol(extrap_left)
+        extrap_right_disc = disc.process_symbol(extrap_right)
+        extrap_grad_left_disc = disc.process_symbol(extrap_grad_left)
+        extrap_grad_right_disc = disc.process_symbol(extrap_grad_right)
+
+        # check constant extrapolates to constant
+        constant_y = np.ones_like(mesh["domain"].nodes[:, np.newaxis])
+        assert extrap_left_disc.evaluate(None, constant_y) == 1
+        assert extrap_right_disc.evaluate(None, constant_y) == 1
+
+        # check linear variable extrapolates correctly
+        linear_y = mesh["domain"].nodes
+        np.testing.assert_array_almost_equal(
+            extrap_left_disc.evaluate(None, linear_y), 0
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_right_disc.evaluate(None, linear_y), 1
+        )
+
+        # check gradient extrapolates correctly
+        np.testing.assert_array_almost_equal(
+            extrap_grad_left_disc.evaluate(None, linear_y), 1 / 2
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_grad_right_disc.evaluate(None, linear_y), 1 / 2
+        )
+
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "quadratic", "value": "quadratic"},
+                "use bcs": True,
+            }
+        }
+        spatial_methods = {"domain": pybamm.FiniteVolume(method_options)}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        var = pybamm.Variable("var", domain="domain")
+        extrap_left = pybamm.BoundaryValue(var, "left")
+        extrap_right = pybamm.BoundaryValue(var, "right")
+        extrap_grad_left = pybamm.BoundaryGradient(var, "left")
+        extrap_grad_right = pybamm.BoundaryGradient(var, "right")
+        disc.set_variable_slices([var])
+        extrap_left_disc = disc.process_symbol(extrap_left)
+        extrap_right_disc = disc.process_symbol(extrap_right)
+        extrap_grad_left_disc = disc.process_symbol(extrap_grad_left)
+        extrap_grad_right_disc = disc.process_symbol(extrap_grad_right)
+
+        # check constant extrapolates to constant
+        np.testing.assert_array_almost_equal(
+            extrap_left_disc.evaluate(None, constant_y), 1
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_right_disc.evaluate(None, constant_y), 1
+        )
+
+        # check linear variable extrapolates correctly
+        np.testing.assert_array_almost_equal(
+            extrap_left_disc.evaluate(None, linear_y), 0
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_right_disc.evaluate(None, linear_y), 1
+        )
+
+        # check gradient extrapolates correctly
+        np.testing.assert_array_almost_equal(
+            extrap_grad_left_disc.evaluate(None, linear_y), 1 / 2
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_grad_right_disc.evaluate(None, linear_y), 1 / 2
+        )
+
     def test_quadratic_extrapolate_left_right(self):
         # create discretisation
         mesh = get_mesh_for_testing()

From edb9402321f40c7e1c7f7f9f20c631e3588a63ac Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Tue, 31 Dec 2024 13:07:15 -0800
Subject: [PATCH 13/13] changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c15aba8bbe..014e09be49 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+- Added symbolic mesh which allows for using InputParameters for geometric parameters ([#4665](https://github.com/pybamm-team/PyBaMM/pull/4665))
 - Enhanced the `search` method to accept multiple search terms in the form of a string or a list. ([#4650](https://github.com/pybamm-team/PyBaMM/pull/4650))
 - Made composite electrode model compatible with particle size distribution ([#4687](https://github.com/pybamm-team/PyBaMM/pull/4687))
 - Added `Symbol.post_order()` method to return an iterable that steps through the tree in post-order fashion. ([#4684](https://github.com/pybamm-team/PyBaMM/pull/4684))