Repository: pybamm-team/PyBaMM
Pull request #5108 (2025-07-09 20:16:49)

Boundary mesh size
# Description

follow up to #5107 providing the distance from the center of the control volume to the edge at the boundary

Fixes # (issue)

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #)

# Important checks:

Please confirm the following before marking the PR as ready for review:
- No style issues: `nox -s pre-commit`
- All tests pass: `nox -s tests`
- The documentation builds: `nox -s doctests`
- Code is commented for hard-to-understand areas
- Tests added that prove fix is effective or that feature works


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5108?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.10%. Comparing base [(`f04a1f8`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/f04a1f8ed421fc4653071eaa089be9b85b0855de?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`95f5f14`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/95f5f1452360ba1f9c909f788000687dba84b78c?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).
> Report is 1 commits behind head on develop.

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff            @@
##           develop    #5108   +/-   ##
========================================
  Coverage    99.09%   99.10%           
========================================
  Files          306      306           
  Lines        23659    23688   +29     
========================================
+ Hits         23446    23475   +29     
  Misses         213      213           
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/5108?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

<details><summary> :rocket: New features to boost your workflow: </summary>

- :snowflake: [Test Analytics](https://docs.codecov.com/docs/test-analytics): Detect flaky tests, report on failures, and find test suite problems.
</details>

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b1179b6c1..98e65c7f6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,7 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
 ## Features
-
+- Creates a boundary mesh size object that returns the distance from the center of the leftmost/rightmost control volume to the boundary of the domain ([#5108](https://github.com/pybamm-team/PyBaMM/pull/5108))
 - Introduced entry points for models, similar to parameter sets, and moved entry point handling to `pybamm.dispatch.entry_points`. There is now experimental support for loading third-party models outside of the PyBaMM framework via `pybamm.Model("model_name")`. This API is currently unstable until further notice and may be subject to change without warning. ([#4490](https://github.com/pybamm-team/PyBaMM/pull/4490))
 - Allow for overriding the spatial method's extrapolation and for using constant extrapolation of boundary values ([#5107](https://github.com/pybamm-team/PyBaMM/pull/5107))
 - Creates `BaseProcessedVariable` to enable object combination when adding solutions together ([#5076](https://github.com/pybamm-team/PyBaMM/pull/5076))
@@ -356,7 +356,7 @@ package to install PyBaMM with only the required dependencies. ([conda-forge/pyb
 - Fixed a bug where supplying an initial soc did not work with half cell models ([#3456](https://github.com/pybamm-team/PyBaMM/pull/3456))
 - Fixed a bug where empty lists passed to QuickPlot resulted in an IndexError and did not return a meaningful error message ([#3359](https://github.com/pybamm-team/PyBaMM/pull/3359))
 - Fixed a bug where there was a missing thermal conductivity in the thermal pouch cell models ([#3330](https://github.com/pybamm-team/PyBaMM/pull/3330))
-- Fixed a bug that caused incorrect results of “{Domain} electrode thickness change [m]” due to the absence of dimension for the variable `electrode_thickness_change`([#3329](https://github.com/pybamm-team/PyBaMM/pull/3329)).
+- Fixed a bug that caused incorrect results of "{Domain} electrode thickness change [m]" due to the absence of dimension for the variable `electrode_thickness_change`([#3329](https://github.com/pybamm-team/PyBaMM/pull/3329)).
 - Fixed a bug that occured in `check_ys_are_not_too_large` when trying to reference `y-slice` where the referenced variable was not a `pybamm.StateVector` ([#3313](https://github.com/pybamm-team/PyBaMM/pull/3313)
 - Fixed a bug with `_Heaviside._evaluate_for_shape` which meant some expressions involving heaviside function and subtractions did not work ([#3306](https://github.com/pybamm-team/PyBaMM/pull/3306))
 - Attributes of `pybamm.Simulation` objects (models, parameter values, geometries, choice of solver, and output variables) are now private and as such cannot be edited in-place after the simulation has been created ([#3267](https://github.com/pybamm-team/PyBaMM/pull/3267)
diff --git a/src/pybamm/expression_tree/unary_operators.py b/src/pybamm/expression_tree/unary_operators.py
index 1471e943f..8a9f21359 100644
--- a/src/pybamm/expression_tree/unary_operators.py
+++ b/src/pybamm/expression_tree/unary_operators.py
@@ -1085,6 +1085,22 @@ class BoundaryValue(BoundaryOperator):
             return sympy.Symbol(latex_child)
 
 
+class BoundaryMeshSize(BoundaryOperator):
+    """
+    A node in the expression tree which gets the width of the control volume at the boundary of a variable's domain.
+
+    Parameters
+    ----------
+    child : :class:`pybamm.Symbol`
+        The variable whose boundary value to take
+    side : str
+        Which side to take the boundary value on ("left" or "right")
+    """
+
+    def __init__(self, child, side):
+        super().__init__("boundary mesh size", child, side)
+
+
 class ExplicitTimeIntegral(UnaryOperator):
     def __init__(self, children, initial_condition):
         super().__init__("explicit time integral", children)
diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py
index 717039223..04407593e 100644
--- a/src/pybamm/spatial_methods/finite_volume.py
+++ b/src/pybamm/spatial_methods/finite_volume.py
@@ -851,6 +851,25 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return new_gradient
 
+    def _boundary_mesh_size(self, child, side):
+        """
+        Get the mesh size at the boundary of a variable's domain.
+        """
+        submesh = self.mesh[child.domain]
+        if side == "left":
+            if hasattr(submesh, "length"):
+                val = submesh.length * pybamm.Scalar(submesh.d_nodes[0])
+            else:
+                val = pybamm.Scalar(submesh.d_nodes[0])
+        elif side == "right":
+            if hasattr(submesh, "length"):
+                val = submesh.length * pybamm.Scalar(submesh.d_nodes[-1])
+            else:
+                val = pybamm.Scalar(submesh.d_nodes[-1])
+        else:
+            raise ValueError(f"Invalid side: {side}")
+        return val
+
     def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         """
         Uses extrapolation to get the boundary value or flux of a variable in the
@@ -869,10 +888,12 @@ class FiniteVolume(pybamm.SpatialMethod):
             bcs = {}
 
         extrap_order_gradient = (
-            symbol.order or self.options["extrapolation"]["order"]["gradient"]
+            getattr(symbol, "order", None)
+            or self.options["extrapolation"]["order"]["gradient"]
         )
         extrap_order_value = (
-            symbol.order or self.options["extrapolation"]["order"]["value"]
+            getattr(symbol, "order", None)
+            or self.options["extrapolation"]["order"]["value"]
         )
         use_bcs = self.options["extrapolation"]["use bcs"]
 
@@ -892,7 +913,9 @@ class FiniteVolume(pybamm.SpatialMethod):
         # Create submatrix to compute boundary values or fluxes
         # Derivation of extrapolation formula can be found at:
         # https://github.com/Scottmar93/extrapolation-coefficents/tree/master
-        if isinstance(symbol, pybamm.BoundaryValue):
+        if isinstance(symbol, pybamm.BoundaryMeshSize):
+            return self._boundary_mesh_size(child, symbol.side)
+        elif isinstance(symbol, pybamm.BoundaryValue):
             if use_bcs and pybamm.has_bc_of_form(child, symbol.side, bcs, "Dirichlet"):
                 # just use the value from the bc: f(x*)
                 sub_matrix = csr_matrix((1, prim_pts))


Test patch:
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 2dc17bb0c..d0a512c52 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
@@ -808,3 +808,185 @@ class TestExtrapolation:
             extrap_pos_disc.evaluate(None, constant_y),
             extrap_right_disc.evaluate(None, constant_y),
         )
+
+    def test_boundary_mesh_size(self):
+        """Test BoundaryMeshSize functionality."""
+        # Create a simple 1D geometry
+        geometry = {"domain": {"x": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
+        submesh_types = {"domain": pybamm.Uniform1DSubMesh}
+        var_pts = {"x": 5}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": False,
+            }
+        }
+        spatial_methods = {"domain": pybamm.FiniteVolume(method_options)}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # Create a variable
+        var = pybamm.Variable("var", domain="domain")
+
+        # Test left boundary mesh size
+        left_mesh_size = pybamm.BoundaryMeshSize(var, "left")
+        right_mesh_size = pybamm.BoundaryMeshSize(var, "right")
+
+        disc.set_variable_slices([var])
+        left_mesh_size_disc = disc.process_symbol(left_mesh_size)
+        right_mesh_size_disc = disc.process_symbol(right_mesh_size)
+
+        # Get the submesh for direct comparison
+        submesh = mesh["domain"]
+
+        # Calculate expected values
+        expected_left = submesh.d_nodes[0]
+        expected_right = submesh.d_nodes[-1]
+
+        # Test that the discretized symbols return the correct values
+        np.testing.assert_allclose(
+            left_mesh_size_disc.evaluate(), expected_left, rtol=1e-15, atol=1e-15
+        )
+        np.testing.assert_allclose(
+            right_mesh_size_disc.evaluate(), expected_right, rtol=1e-15, atol=1e-15
+        )
+
+        # Test with different mesh sizes
+        for npts in [10, 20, 50]:
+            var_pts = {"x": npts}
+            mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+            spatial_methods = {"domain": pybamm.FiniteVolume(method_options)}
+            disc = pybamm.Discretisation(mesh, spatial_methods)
+
+            var = pybamm.Variable("var", domain="domain")
+            left_mesh_size = pybamm.BoundaryMeshSize(var, "left")
+            right_mesh_size = pybamm.BoundaryMeshSize(var, "right")
+
+            disc.set_variable_slices([var])
+            left_mesh_size_disc = disc.process_symbol(left_mesh_size)
+            right_mesh_size_disc = disc.process_symbol(right_mesh_size)
+
+            submesh = mesh["domain"]
+            expected_left = submesh.d_nodes[0]
+            expected_right = submesh.d_nodes[-1]
+
+            np.testing.assert_allclose(
+                left_mesh_size_disc.evaluate(), expected_left, rtol=1e-15, atol=1e-15
+            )
+            np.testing.assert_allclose(
+                right_mesh_size_disc.evaluate(), expected_right, rtol=1e-15, atol=1e-15
+            )
+
+        # Test with non-uniform mesh
+        geometry_nonuniform = {
+            "negative particle": {"r_n": {"min": 0, "max": 1}},
+        }
+        submesh_types_nonuniform = {
+            "negative particle": pybamm.MeshGenerator(pybamm.Exponential1DSubMesh),
+        }
+        var_pts_nonuniform = {"r_n": 10}
+        mesh_nonuniform = pybamm.Mesh(
+            geometry_nonuniform, submesh_types_nonuniform, var_pts_nonuniform
+        )
+
+        spatial_methods_nonuniform = {
+            "negative particle": pybamm.FiniteVolume(method_options)
+        }
+        disc_nonuniform = pybamm.Discretisation(
+            mesh_nonuniform, spatial_methods_nonuniform
+        )
+
+        var_nonuniform = pybamm.Variable("var", domain="negative particle")
+        left_mesh_size_nonuniform = pybamm.BoundaryMeshSize(var_nonuniform, "left")
+        right_mesh_size_nonuniform = pybamm.BoundaryMeshSize(var_nonuniform, "right")
+
+        disc_nonuniform.set_variable_slices([var_nonuniform])
+        left_mesh_size_disc_nonuniform = disc_nonuniform.process_symbol(
+            left_mesh_size_nonuniform
+        )
+        right_mesh_size_disc_nonuniform = disc_nonuniform.process_symbol(
+            right_mesh_size_nonuniform
+        )
+
+        submesh_nonuniform = mesh_nonuniform["negative particle"]
+        expected_left_nonuniform = submesh_nonuniform.d_nodes[0]
+        expected_right_nonuniform = submesh_nonuniform.d_nodes[-1]
+
+        np.testing.assert_allclose(
+            left_mesh_size_disc_nonuniform.evaluate(),
+            expected_left_nonuniform,
+            rtol=1e-15,
+            atol=1e-15,
+        )
+        np.testing.assert_allclose(
+            right_mesh_size_disc_nonuniform.evaluate(),
+            expected_right_nonuniform,
+            rtol=1e-15,
+            atol=1e-15,
+        )
+
+        # Test with symbolic submesh
+        geometry_symbolic = {
+            "domain": {"x": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}
+        }
+        submesh_types_symbolic = {"domain": pybamm.SymbolicUniform1DSubMesh}
+        var_pts_symbolic = {"x": 8}
+        mesh_symbolic = pybamm.Mesh(
+            geometry_symbolic, submesh_types_symbolic, var_pts_symbolic
+        )
+
+        spatial_methods_symbolic = {"domain": pybamm.FiniteVolume(method_options)}
+        disc_symbolic = pybamm.Discretisation(mesh_symbolic, spatial_methods_symbolic)
+
+        var_symbolic = pybamm.Variable("var", domain="domain")
+        left_mesh_size_symbolic = pybamm.BoundaryMeshSize(var_symbolic, "left")
+        right_mesh_size_symbolic = pybamm.BoundaryMeshSize(var_symbolic, "right")
+
+        disc_symbolic.set_variable_slices([var_symbolic])
+        left_mesh_size_disc_symbolic = disc_symbolic.process_symbol(
+            left_mesh_size_symbolic
+        )
+        right_mesh_size_disc_symbolic = disc_symbolic.process_symbol(
+            right_mesh_size_symbolic
+        )
+
+        submesh_symbolic = mesh_symbolic["domain"]
+        expected_left_symbolic = submesh_symbolic.d_nodes[0]
+        expected_right_symbolic = submesh_symbolic.d_nodes[-1]
+
+        np.testing.assert_allclose(
+            left_mesh_size_disc_symbolic.evaluate(),
+            expected_left_symbolic,
+            rtol=1e-15,
+            atol=1e-15,
+        )
+        np.testing.assert_allclose(
+            right_mesh_size_disc_symbolic.evaluate(),
+            expected_right_symbolic,
+            rtol=1e-15,
+            atol=1e-15,
+        )
+
+        # Test error handling for invalid side
+        var_error = pybamm.Variable("var", domain="domain")
+        with pytest.raises(ValueError, match="Invalid side"):
+            invalid_mesh_size = pybamm.BoundaryMeshSize(var_error, "invalid")
+            disc.set_variable_slices([var_error])
+            disc.process_symbol(invalid_mesh_size)
+
+        # Test that the BoundaryMeshSize symbol has the correct properties
+        left_mesh_size_symbol = pybamm.BoundaryMeshSize(var, "left")
+        right_mesh_size_symbol = pybamm.BoundaryMeshSize(var, "right")
+
+        # Check that it's a BoundaryMeshSize instance
+        assert isinstance(left_mesh_size_symbol, pybamm.BoundaryMeshSize)
+        assert isinstance(right_mesh_size_symbol, pybamm.BoundaryMeshSize)
+
+        # Check side property
+        assert left_mesh_size_symbol.side == "left"
+        assert right_mesh_size_symbol.side == "right"
+
+        # Check child property
+        assert left_mesh_size_symbol.child == var
+        assert right_mesh_size_symbol.child == var