Skip to content

Fix FiniteVolume2D boundary extrapolation on meshes with fewer than 3 nodes #5839

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
aabills wants to merge 3 commits into
base: main
Choose a base branch
from
+139 −39
Open

Fix FiniteVolume2D boundary extrapolation on meshes with fewer than 3 nodes #5839

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
833dc72
Compute FiniteVolume2D extrapolation spacings only when needed
aabills Sep 28, 2026
965bea9
Add PR link to CHANGELOG entry
aabills Oct 2, 2026
99fec58
Merge remote-tracking branch 'origin/main' into fix/fv2d-coarse-mesh-…
aabills Oct 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,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`. ([#5839](https://github.com/pybamm-team/PyBaMM/pull/5839))
- A per-state `atol` given to `IDAKLUSolver` (or as `model.atol`) can now be a list, tuple or array, so it survives a `to_config`/`from_config` round trip, which returns a JSON list. One without exactly one value per state is rejected with a `SolverError` instead of reaching the integrator, as are booleans, even one among numbers such as `[True, 1e-6]`, and non-real values. ([#5782](https://github.com/pybamm-team/PyBaMM/pull/5782))
- `IDAKLUSolver`'s `print_stats` with an iterative linear solver and `preconditioner="none"` no longer reports a garbage "Number of calls to residual function in preconditioner": without a BBD preconditioner, SUNDIALS read that count out of the model's own data. ([#5782](https://github.com/pybamm-team/PyBaMM/pull/5782))
- Discretising a model whose differential variables are all scalars or on zero-dimensional domains (such as a lumped current collector) no longer raises SciPy 1.18's `block_diag` `DeprecationWarning`, and its mass matrix stays a sparse matrix rather than becoming a sparse array once SciPy changes `block_diag`'s return type. ([#5782](https://github.com/pybamm-team/PyBaMM/pull/5782))
Expand Down
88 changes: 49 additions & 39 deletions packages/pybamm/src/pybamm/spatial_methods/finite_volume_2d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1103,20 +1103,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

Expand Down Expand Up @@ -1168,7 +1178,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
Expand Down Expand Up @@ -1211,8 +1221,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))
Expand Down Expand Up @@ -1265,7 +1275,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
Expand Down Expand Up @@ -1307,8 +1317,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)
Expand Down Expand Up @@ -1373,7 +1383,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)
Expand All @@ -1396,8 +1406,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))
Expand Down Expand Up @@ -1438,7 +1448,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)
Expand All @@ -1450,7 +1460,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)
Expand All @@ -1473,8 +1483,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)
Expand Down Expand Up @@ -1549,7 +1559,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)]
Expand Down Expand Up @@ -1602,7 +1612,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)]
Expand Down Expand Up @@ -1632,7 +1642,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
Expand All @@ -1652,8 +1662,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)
Expand Down Expand Up @@ -1685,8 +1695,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
Expand All @@ -1706,8 +1716,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)
Expand Down Expand Up @@ -1737,7 +1747,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)
Expand All @@ -1756,8 +1766,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)
Expand All @@ -1783,7 +1793,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)
Expand All @@ -1802,8 +1812,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)
Expand Down
89 changes: 89 additions & 0 deletions packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_extrapolation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Loading