From 833dc72ff69ae7b4c191d538de16b0f39cc52b14 Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Mon, 28 Sep 2026 16:20:02 -0700 Subject: [PATCH 1/2] Compute FiniteVolume2D extrapolation spacings only when needed Boundary extrapolation read the second node spacing in both directions unconditionally, so any 2D mesh with fewer than 3 nodes in either direction raised IndexError. Spacings are now looked up lazily and a DiscretisationError is raised when the mesh is too coarse for the requested extrapolation. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 4 + .../spatial_methods/finite_volume_2d.py | 88 ++++++++++-------- .../test_extrapolation.py | 89 +++++++++++++++++++ 3 files changed, 142 insertions(+), 39 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 02e63cacc2..5af7e2a172 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) +## Bug fixes + +- Fixed `FiniteVolume2D` boundary values, gradients and integrals raising an `IndexError` on 2D meshes with fewer than 3 nodes in either direction, even when the extrapolation only needed the other direction. Meshes too coarse for the requested extrapolation order now raise a `DiscretisationError`. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX)) + # [v26.9.0.0](https://github.com/pybamm-team/PyBaMM/tree/pybamm-v26.9.0.0) - 2026-09-28 ## Breaking changes diff --git a/packages/pybamm/src/pybamm/spatial_methods/finite_volume_2d.py b/packages/pybamm/src/pybamm/spatial_methods/finite_volume_2d.py index 9d47087b21..b3c299ecab 100644 --- a/packages/pybamm/src/pybamm/spatial_methods/finite_volume_2d.py +++ b/packages/pybamm/src/pybamm/spatial_methods/finite_volume_2d.py @@ -1138,20 +1138,30 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): edges_tb = submesh.edges_tb dx0_lr = nodes_lr[0] - edges_lr[0] - dx1_lr = submesh.d_nodes_lr[0] - dx2_lr = submesh.d_nodes_lr[1] - dxN_lr = edges_lr[-1] - nodes_lr[-1] - dxNm1_lr = submesh.d_nodes_lr[-1] - dxNm2_lr = submesh.d_nodes_lr[-2] - dx0_tb = nodes_tb[0] - edges_tb[0] - dx1_tb = submesh.d_nodes_tb[0] - dx2_tb = submesh.d_nodes_tb[1] - dxN_tb = edges_tb[-1] - nodes_tb[-1] - dxNm1_tb = submesh.d_nodes_tb[-1] - dxNm2_tb = submesh.d_nodes_tb[-2] + + if isinstance(symbol, pybamm.BoundaryGradient): + extrap_order = extrap_order_gradient + else: + extrap_order = extrap_order_value + + # Node spacings are looked up lazily so that coarse meshes only fail when + # the requested extrapolation actually needs the missing nodes + def node_spacing(direction, index): + d_nodes = getattr(submesh, f"d_nodes_{direction}") + required_npts = index + 2 if index >= 0 else 1 - index + if len(d_nodes) + 1 < required_npts: + raise pybamm.DiscretisationError( + f"{extrap_order.capitalize()} extrapolation for " + f"'{symbol.name}' requires at least {required_npts} nodes in " + f"the '{direction}' direction of domain " + f"{discretised_child.domain}, but the mesh has " + f"{len(d_nodes) + 1}. Refine the mesh or use a lower " + "extrapolation order." + ) + return d_nodes[index] child = symbol.child @@ -1203,7 +1213,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): else: dx0 = dx0_lr - dx1 = dx1_lr + dx1 = node_spacing("lr", 0) row_indices = np.arange(0, n_tb) col_indices_0 = np.arange(0, n_tb * n_lr, n_lr) col_indices_1 = col_indices_0 + 1 @@ -1246,8 +1256,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): else: dx0 = dx0_lr - dx1 = dx1_lr - dx2 = dx2_lr + dx1 = node_spacing("lr", 0) + dx2 = node_spacing("lr", 1) a = (dx0 + dx1) * (dx0 + dx1 + dx2) / (dx1 * (dx1 + dx2)) b = -dx0 * (dx0 + dx1 + dx2) / (dx1 * dx2) c = dx0 * (dx0 + dx1) / (dx2 * (dx1 + dx2)) @@ -1300,7 +1310,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): # to find value at x* use formula: # f(x*) = f_N - (dxN / dxNm1) (f_N - f_Nm1) dxN = dxN_lr - dxNm1 = dxNm1_lr + dxNm1 = node_spacing("lr", -1) row_indices = np.arange(0, n_tb) col_indices_Nm1 = np.arange(n_lr - 2, n_lr * n_tb, n_lr) col_indices_N = col_indices_Nm1 + 1 @@ -1342,8 +1352,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): raise NotImplementedError else: dxN = dxN_lr - dxNm1 = dxNm1_lr - dxNm2 = dxNm2_lr + dxNm1 = node_spacing("lr", -1) + dxNm2 = node_spacing("lr", -2) a = ( (dxN + dxNm1) * (dxN + dxNm1 + dxNm2) @@ -1408,7 +1418,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): additive = -dx0 * bcs[child][side_first][0] else: dx0 = dx0_tb - dx1 = dx1_tb + dx1 = node_spacing("tb", 0) first_val = (1 + (dx0 / dx1)) * np.ones(n_lr) second_val = -(dx0 / dx1) * np.ones(n_lr) rows_first = np.arange(0, n_lr) @@ -1431,8 +1441,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): raise NotImplementedError else: dx0 = dx0_tb - dx1 = dx1_tb - dx2 = dx2_tb + dx1 = node_spacing("tb", 0) + dx2 = node_spacing("tb", 1) a = (dx0 + dx1) * (dx0 + dx1 + dx2) / (dx1 * (dx1 + dx2)) b = -dx0 * (dx0 + dx1 + dx2) / (dx1 * dx2) c = dx0 * (dx0 + dx1) / (dx2 * (dx1 + dx2)) @@ -1473,7 +1483,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): if use_bcs and pybamm.has_bc_of_form( child, side_first, bcs, "Neumann" ): - dxNm1 = dxNm1_tb + dxNm1 = node_spacing("tb", -1) dxN = dxN_tb val_N = np.ones(n_lr) rows = np.arange(0, n_lr) @@ -1485,7 +1495,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): additive = dxN * bcs[child][side_first][0] else: dx0 = dxN_tb - dx1 = dxNm1_tb + dx1 = node_spacing("tb", -1) first_val = -(dx0 / dx1) * np.ones(n_lr) second_val = (1 + (dx0 / dx1)) * np.ones(n_lr) rows_first = np.arange(0, n_lr) @@ -1508,8 +1518,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): raise NotImplementedError else: dxN = dxN_tb - dxNm1 = dxNm1_tb - dxNm2 = dxNm2_tb + dxNm1 = node_spacing("tb", -1) + dxNm2 = node_spacing("tb", -2) a = ( (dxN + dxNm1) * (dxN + dxNm1 + dxNm2) @@ -1584,7 +1594,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): else: dx0 = dx0_tb - dx1 = dx1_tb + dx1 = node_spacing("tb", 0) row_indices = [0, 0] col_indices = [0, 1] vals = [1.0 + (dx0 / dx1), -(dx0 / dx1)] @@ -1637,7 +1647,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): else: dxN = dxN_tb - dxNm1 = dxNm1_tb + dxNm1 = node_spacing("tb", -1) row_indices = [0, 0] col_indices = [n_tb - 2, n_tb - 1] vals = [-(dxN / dxNm1), 1.0 + (dxN / dxNm1)] @@ -1667,7 +1677,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): elif side_first == "left": if extrap_order_gradient == "linear": # f'(x*) = (f_2 - f_1) / dx1 - dx1 = dx1_lr + dx1 = node_spacing("lr", 0) row_indices = np.arange(0, n_tb) col_indices_0 = np.arange(0, n_tb * n_lr, n_lr) col_indices_1 = col_indices_0 + 1 @@ -1687,8 +1697,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): elif extrap_order_gradient == "quadratic": dx0 = dx0_lr - dx1 = dx1_lr - dx2 = dx2_lr + dx1 = node_spacing("lr", 0) + dx2 = node_spacing("lr", 1) a = -(2 * dx0 + 2 * dx1 + dx2) / (dx1**2 + dx1 * dx2) b = (2 * dx0 + dx1 + dx2) / (dx1 * dx2) c = -(2 * dx0 + dx1) / (dx1 * dx2 + dx2**2) @@ -1720,8 +1730,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): if extrap_order_gradient == "linear": # use formula: # f'(x*) = (f_N - f_Nm1) / dxNm1 - dxN = dxNm1_lr - dxNm1 = dxNm1_lr + dxN = node_spacing("lr", -1) + dxNm1 = node_spacing("lr", -1) row_indices = np.arange(0, n_tb) col_indices_Nm1 = np.arange(n_lr - 2, n_lr * n_tb, n_lr) col_indices_N = col_indices_Nm1 + 1 @@ -1741,8 +1751,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): elif extrap_order_gradient == "quadratic": dxN = dxN_lr - dxNm1 = dxNm1_lr - dxNm2 = dxNm2_lr + dxNm1 = node_spacing("lr", -1) + dxNm2 = node_spacing("lr", -2) a = (2 * dxN + 2 * dxNm1 + dxNm2) / (dxNm1**2 + dxNm1 * dxNm2) b = -(2 * dxN + dxNm1 + dxNm2) / (dxNm1 * dxNm2) c = (2 * dxN + dxNm1) / (dxNm1 * dxNm2 + dxNm2**2) @@ -1772,7 +1782,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): elif side_first == "bottom": if extrap_order_gradient == "linear": - dx1 = dx1_tb + dx1 = node_spacing("tb", 0) row_indices = np.arange(0, n_lr) col_indices_0 = np.arange(0, n_lr) col_indices_1 = np.arange(n_lr, 2 * n_lr) @@ -1791,8 +1801,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): additive = pybamm.Scalar(0) elif extrap_order_gradient == "quadratic": dx0 = dx0_tb - dx1 = dx1_tb - dx2 = dx2_tb + dx1 = node_spacing("tb", 0) + dx2 = node_spacing("tb", 1) a = -(2 * dx0 + 2 * dx1 + dx2) / (dx1**2 + dx1 * dx2) b = (2 * dx0 + dx1 + dx2) / (dx1 * dx2) c = -(2 * dx0 + dx1) / (dx1 * dx2 + dx2**2) @@ -1818,7 +1828,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): additive = pybamm.Scalar(0) elif side_first == "top": if extrap_order_gradient == "linear": - dxNm1 = dxNm1_tb + dxNm1 = node_spacing("tb", -1) row_indices = np.arange(0, n_lr) col_indices_0 = np.arange((n_tb - 2) * n_lr, (n_tb - 1) * n_lr) col_indices_1 = np.arange((n_tb - 1) * n_lr, n_tb * n_lr) @@ -1837,8 +1847,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None): additive = pybamm.Scalar(0) elif extrap_order_gradient == "quadratic": dxN = dxN_tb - dxNm1 = dxNm1_tb - dxNm2 = dxNm2_tb + dxNm1 = node_spacing("tb", -1) + dxNm2 = node_spacing("tb", -2) a = (2 * dxN + 2 * dxNm1 + dxNm2) / (dxNm1**2 + dxNm1 * dxNm2) b = -(2 * dxN + dxNm1 + dxNm2) / (dxNm1 * dxNm2) c = (2 * dxN + dxNm1) / (dxNm1 * dxNm2 + dxNm2**2) diff --git a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_extrapolation.py b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_extrapolation.py index ddff166e24..84ff6390b3 100644 --- a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_extrapolation.py +++ b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_extrapolation.py @@ -528,3 +528,92 @@ def test_boundary_value_constant_extrapolation(self, mesh_2d): discretised_bottom = disc.process_symbol(boundary_value_bottom) result_bottom = discretised_bottom.evaluate(y=tb).flatten() np.testing.assert_allclose(result_bottom, expected_bottom) + + +def get_coarse_mesh_2d(npts_lr, npts_tb): + x = pybamm.SpatialVariable("x", ["negative electrode"], direction="lr") + z = pybamm.SpatialVariable("z", ["negative electrode"], direction="tb") + geometry = { + "negative electrode": { + x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}, + z: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(2)}, + } + } + return pybamm.Mesh( + geometry, + {"negative electrode": pybamm.Uniform2DSubMesh}, + {x: npts_lr, z: npts_tb}, + ) + + +class TestExtrapolationCoarseMesh2D: + @pytest.mark.parametrize("order", ["linear", "quadratic"]) + def test_lr_extrapolation_with_two_tb_nodes(self, order): + mesh = get_coarse_mesh_2d(npts_lr=10, npts_tb=2) + submesh = mesh["negative electrode"] + disc = pybamm.Discretisation( + mesh, + { + "negative electrode": pybamm.FiniteVolume2D( + { + "extrapolation": { + "order": {"gradient": "linear", "value": order}, + "use bcs": False, + } + } + ) + }, + ) + var = pybamm.Variable("var", ["negative electrode"]) + disc.set_variable_slices([var]) + + # f(x, z) = x is reproduced exactly by linear and quadratic extrapolation + LR, _ = np.meshgrid(submesh.nodes_lr, submesh.nodes_tb) + y = LR.flatten() + for side, expected in [("left", 0), ("right", 1)]: + boundary_value = disc.process_symbol(pybamm.BoundaryValue(var, side)) + np.testing.assert_allclose( + boundary_value.evaluate(y=y).flatten(), expected, atol=1e-12 + ) + + # the right boundary has length 2, so integrating f = x over it gives 2 + boundary_integral = disc.process_symbol(pybamm.BoundaryIntegral(var, "right")) + np.testing.assert_allclose(boundary_integral.evaluate(y=y), 2, rtol=1e-12) + + @pytest.mark.parametrize( + "symbol_class,order,side,npts_lr,npts_tb,required_npts", + [ + (pybamm.BoundaryValue, "quadratic", "top", 10, 2, 3), + (pybamm.BoundaryValue, "quadratic", "bottom", 10, 2, 3), + (pybamm.BoundaryValue, "linear", "left", 1, 10, 2), + (pybamm.BoundaryValue, "quadratic", "right", 2, 10, 3), + (pybamm.BoundaryGradient, "quadratic", "left", 2, 10, 3), + (pybamm.BoundaryGradient, "linear", "top", 10, 1, 2), + ], + ) + def test_mesh_too_coarse_for_extrapolation( + self, symbol_class, order, side, npts_lr, npts_tb, required_npts + ): + mesh = get_coarse_mesh_2d(npts_lr=npts_lr, npts_tb=npts_tb) + disc = pybamm.Discretisation( + mesh, + { + "negative electrode": pybamm.FiniteVolume2D( + { + "extrapolation": { + "order": {"gradient": order, "value": order}, + "use bcs": False, + } + } + ) + }, + ) + var = pybamm.Variable("var", ["negative electrode"]) + disc.set_variable_slices([var]) + direction = "lr" if side in ["left", "right"] else "tb" + with pytest.raises( + pybamm.DiscretisationError, + match=rf"{order.capitalize()} extrapolation .* at least {required_npts} " + rf"nodes in the '{direction}' direction", + ): + disc.process_symbol(symbol_class(var, side)) From 965bea9661adb6336d749e3b67a1400219c1bc8d Mon Sep 17 00:00:00 2001 From: Alexander Bills Date: Fri, 2 Oct 2026 15:41:30 -0700 Subject: [PATCH 2/2] Add PR link to CHANGELOG entry Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5af7e2a172..cce69dbd2f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,7 @@ ## Bug fixes -- Fixed `FiniteVolume2D` boundary values, gradients and integrals raising an `IndexError` on 2D meshes with fewer than 3 nodes in either direction, even when the extrapolation only needed the other direction. Meshes too coarse for the requested extrapolation order now raise a `DiscretisationError`. ([#XXXX](https://github.com/pybamm-team/PyBaMM/pull/XXXX)) +- Fixed `FiniteVolume2D` boundary values, gradients and integrals raising an `IndexError` on 2D meshes with fewer than 3 nodes in either direction, even when the extrapolation only needed the other direction. Meshes too coarse for the requested extrapolation order now raise a `DiscretisationError`. ([#5839](https://github.com/pybamm-team/PyBaMM/pull/5839)) # [v26.9.0.0](https://github.com/pybamm-team/PyBaMM/tree/pybamm-v26.9.0.0) - 2026-09-28