diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 58138b8d0..9ecf6cf73 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -105,6 +105,7 @@ from .expression_tree.exceptions import (
     DomainError,
     ModelError,
     SolverError,
+    ShapeError,
     ModelWarning,
 )
 from .expression_tree.simplify import (
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index e4c5246f6..c737a1d0a 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -23,6 +23,14 @@ class SolverError(Exception):
     pass
 
 
+class ShapeError(Exception):
+    """
+    Shape error: cannot evaluate an object to find its shape
+    """
+
+    pass
+
+
 class ModelWarning(UserWarning):
     """
     Model warning: the model does not contain all of the standard output variables
diff --git a/pybamm/models/submodels/electrolyte_current.py b/pybamm/models/submodels/electrolyte_current.py
index 8b2a5795a..a68b71bbc 100644
--- a/pybamm/models/submodels/electrolyte_current.py
+++ b/pybamm/models/submodels/electrolyte_current.py
@@ -484,9 +484,9 @@ class MacInnesCapacitance(ElectrolyteCurrentBaseModel):
         # import parameters and spatial variables
         param = self.set_of_parameters
         i_cell = param.current_with_time
-        x_n = pybamm.standard_spatial_vars.x_n
-        x_s = pybamm.standard_spatial_vars.x_s
-        x_p = pybamm.standard_spatial_vars.x_p
+        # x_n = pybamm.standard_spatial_vars.x_n
+        # x_s = pybamm.standard_spatial_vars.x_s
+        # x_p = pybamm.standard_spatial_vars.x_p
 
         # Combine currents
         i_e_s = pybamm.Broadcast(i_cell, ["separator"])
@@ -498,32 +498,34 @@ class MacInnesCapacitance(ElectrolyteCurrentBaseModel):
         eps_n, eps_s, eps_p = eps.orphans
 
         # Compute potentials
-        phi_e_children = [None] * 3
-        for i, (c_e, eps, i_e, x) in enumerate(
-            [
-                (c_e_n, eps_n, i_e_n, x_n),
-                (c_e_s, eps_s, i_e_s, x_s),
-                (c_e_p, eps_p, i_e_p, x_p),
-            ]
-        ):
-            chi_e = param.chi(c_e)
-            kappa_eff = param.kappa_e(c_e) * (eps ** param.b)
-            d_phi_e__d_x = chi_e / c_e * pybamm.grad(c_e) - param.C_e * i_e / kappa_eff
-            phi_e_children[i] = pybamm.IndefiniteIntegral(d_phi_e__d_x, x)
-
-        # Adjust for boundary conditions and continuity
-        phi_e_n, phi_e_s, phi_e_p = phi_e_children
-        phi_e_n = phi_e_n + pybamm.boundary_value(-delta_phi_n - phi_e_n, "left")
-        phi_e_s = (
-            phi_e_s
-            - pybamm.boundary_value(phi_e_s, "left")
-            + pybamm.boundary_value(phi_e_n, "right")
-        )
-        phi_e_p = (
-            phi_e_p
-            - pybamm.boundary_value(phi_e_p, "left")
-            + pybamm.boundary_value(phi_e_s, "right")
-        )
+        phi_e_n = pybamm.Broadcast(0, "negative electrode")
+        phi_e_s = pybamm.Broadcast(0, "separator")
+        phi_e_p = pybamm.Broadcast(0, "positive electrode")
+        # for i, (c_e, eps, i_e, x) in enumerate(
+        #     [
+        #         (c_e_n, eps_n, i_e_n, x_n),
+        #         (c_e_s, eps_s, i_e_s, x_s),
+        #         (c_e_p, eps_p, i_e_p, x_p),
+        #     ]
+        # ):
+        #     chi_e = param.chi(c_e)
+        #     kappa_eff = param.kappa_e(c_e) * (eps ** param.b)
+        #     d_phi_e__d_x= chi_e / c_e * pybamm.grad(c_e) - param.C_e * i_e / kappa_eff
+        #     phi_e_children[i] = pybamm.IndefiniteIntegral(d_phi_e__d_x, x)
+        #
+        # # Adjust for boundary conditions and continuity
+        # phi_e_n, phi_e_s, phi_e_p = phi_e_children
+        # phi_e_n = phi_e_n + pybamm.boundary_value(-delta_phi_n - phi_e_n, "left")
+        # phi_e_s = (
+        #     phi_e_s
+        #     - pybamm.boundary_value(phi_e_s, "left")
+        #     + pybamm.boundary_value(phi_e_n, "right")
+        # )
+        # phi_e_p = (
+        #     phi_e_p
+        #     - pybamm.boundary_value(phi_e_p, "left")
+        #     + pybamm.boundary_value(phi_e_s, "right")
+        # )
 
         # Concatenate
         phi_e = pybamm.Concatenation(phi_e_n, phi_e_s, phi_e_p)
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 0686b9b5d..2c60ca360 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -74,7 +74,9 @@ class FiniteVolume(pybamm.SpatialMethod):
         # note in 1D spherical grad and normal grad are the same
         gradient_matrix = self.gradient_matrix(domain)
 
-        return gradient_matrix @ discretised_symbol
+        out = gradient_matrix @ discretised_symbol
+        self.test_shape(out)
+        return out
 
     def gradient_matrix(self, domain):
         """
@@ -137,6 +139,7 @@ class FiniteVolume(pybamm.SpatialMethod):
         else:
             out = divergence_matrix @ discretised_symbol
 
+        self.test_shape(out)
         return out
 
     def divergence_matrix(self, domain):
@@ -188,6 +191,8 @@ class FiniteVolume(pybamm.SpatialMethod):
         else:
             out = integration_vector @ discretised_symbol
         out.domain = []
+
+        self.test_shape(out)
         return out
 
     def definite_integral_vector(self, domain):
@@ -240,6 +245,7 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         out.domain = domain
 
+        self.test_shape(out)
         return out
 
     def indefinite_integral_matrix(self, domain):
@@ -461,6 +467,7 @@ class FiniteVolume(pybamm.SpatialMethod):
         else:
             boundary_value.domain = []
 
+        self.test_shape(boundary_value)
         return boundary_value
 
     def process_binary_operators(self, bin_op, left, right, disc_left, disc_right):
@@ -500,7 +507,9 @@ class FiniteVolume(pybamm.SpatialMethod):
         elif right_has_grad_not_div and not left_has_grad_not_div:
             disc_left = self.compute_diffusivity(disc_left)
         # Return new binary operator with appropriate class
-        return bin_op.__class__(disc_left, disc_right)
+        out = bin_op.__class__(disc_left, disc_right)
+        self.test_shape(out)
+        return out
 
     def compute_diffusivity(self, discretised_symbol):
         """
@@ -553,6 +562,9 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         # If discretised_symbol evaluates to number there is no need to average
         if discretised_symbol.evaluates_to_number():
-            return discretised_symbol
+            out = discretised_symbol
         else:
-            return arithmetic_mean(discretised_symbol)
+            out = arithmetic_mean(discretised_symbol)
+
+        self.test_shape(out)
+        return out
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 3f1c7f89c..a6fc38798 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -257,3 +257,29 @@ class SpatialMethod:
 
         """
         return bin_op.__class__(disc_left, disc_right)
+
+    def test_shape(self, symbol):
+        """
+        Check that the discretised symbol has a pybamm `shape`, i.e. can be evaluated
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.Symbol`
+            The symbol to be tested
+
+        Raises
+        ------
+        pybamm.ShapeError
+            If the shape of the object cannot be found
+        """
+        try:
+            symbol.shape
+        except ValueError as e:
+            raise pybamm.ShapeError(
+                """
+                Cannot evaluate the discretised symbol to find its shape
+                (original error: {})
+                """.format(
+                    e
+                )
+            )

Shape test
Add a method `test_shape` to the Finite Volume class that checks, when discretising, that we can evaluate the shape of a discretised object without raising any errors.

Should help to catch bugs earlier.
diff --git a/tests/unit/test_spatial_methods/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume.py
index fcfd6f15f..fa503a1e7 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume.py
@@ -959,7 +959,12 @@ class TestFiniteVolume(unittest.TestCase):
         grad_and_div = pybamm.div(pybamm.grad(phi))
         int_grad_phi = pybamm.IndefiniteIntegral(grad_and_div, x)
         disc.set_variable_slices([phi])
-
+        disc._bcs = {
+            phi.id: {
+                "left": (pybamm.Scalar(0), "Dirichlet"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
         with self.assertRaisesRegex(pybamm.ModelError, "integrated"):
             disc.process_symbol(int_grad_phi)
 
diff --git a/tests/unit/test_spatial_methods/test_spatial_methods.py b/tests/unit/test_spatial_methods/test_spatial_methods.py
index ad61333fd..4a93bcf4f 100644
--- a/tests/unit/test_spatial_methods/test_spatial_methods.py
+++ b/tests/unit/test_spatial_methods/test_spatial_methods.py
@@ -25,6 +25,17 @@ class TestSpatialMethod(unittest.TestCase):
         with self.assertRaisesRegex(TypeError, "Cannot process BoundaryFlux"):
             spatial_method.boundary_value_or_flux(symbol, child)
 
+    def test_test_shape(self):
+        mesh = get_mesh_for_testing()
+        spatial_method = pybamm.SpatialMethod(mesh)
+        # right shape, passes
+        y1 = pybamm.StateVector(slice(0, 10))
+        spatial_method.test_shape(y1)
+        # bad shape, fails
+        y2 = pybamm.StateVector(slice(0, 5))
+        with self.assertRaises(pybamm.ShapeError):
+            spatial_method.test_shape(y1 + y2)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gfortran gcc libopenblas-dev liblapack-dev
3.6
requirements-docs.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs==22.2.0
autograd==1.6.2
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
flake8==5.0.4
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.2.0
iniconfig==1.1.1
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.7.0
mistune==0.8.4
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging==21.3
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy==1.0.0
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py==1.11.0
-e git+https://github.com/pybamm-team/PyBaMM.git@e96eda81d981ad1059fcf58faacdc9f0febc9e02#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing==3.1.4
pyrsistent==0.18.0
pytest==7.0.1
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
tomli==1.2.3
tornado==6.1
traitlets==4.3.3
typing_extensions==4.1.1
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp==3.6.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=21.2.2=py36h06a4308_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - attrs==22.2.0
      - autograd==1.6.2
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - flake8==5.0.4
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.2.0
      - iniconfig==1.1.1
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.7.0
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - packaging==21.3
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - pluggy==1.0.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - py==1.11.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyparsing==3.1.4
      - pyrsistent==0.18.0
      - pytest==7.0.1
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tomli==1.2.3
      - tornado==6.1
      - traitlets==4.3.3
      - typing-extensions==4.1.1
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
      - zipp==3.6.0
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_spatial_methods/test_spatial_methods.py::TestSpatialMethod::test_test_shape
tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_add_ghost_nodes tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_add_ghost_nodes_concatenation tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_compute_diffusivity tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_definite_integral tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_discretise_diffusivity_times_spatial_operator tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_discretise_spatial_variable tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_extrapolate_left_right tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Dirichlet_and_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_mixed_domain tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_indefinite_integral tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_jacobian tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_mass_matrix_shape tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_add_ghost_nodes tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_mass_matrix_shape tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_spherical_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_spherical_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_spatial_methods.py::TestSpatialMethod::test_basics
pybamm/__init__.py pybamm/expression_tree/exceptions.py pybamm/models/submodels/electrolyte_current.py pybamm/spatial_methods/finite_volume.py pybamm/spatial_methods/spatial_method.py
pybamm/__init__.py pybamm/expression_tree/exceptions.py pybamm/models/submodels/electrolyte_current.py pybamm/spatial_methods/finite_volume.py pybamm/spatial_methods/spatial_method.py