Repository: pybamm-team/PyBaMM
Pull request #4710 (2024-12-23 15:43:40)

Add tertiary integral to FV method
# Description
Adds tertiary integral to the Finite Volume method

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 #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [ ] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)

# Key checklist:

- [ ] No style issues: `$ pre-commit run` (or `$ nox -s pre-commit`) (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) for how to set this up to run automatically when committing locally, in just two lines of code)
- [ ] All tests pass: `$ python -m pytest` (or `$ nox -s tests`)
- [ ] The documentation builds: `$ python -m pytest --doctest-plus src` (or `$ nox -s doctests`)

You can run integration tests, unit tests, and doctests together at once, using `$ nox -s quick`.

## Further checks:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


rtimms: @aabills 
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4710?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
Attention: Patch coverage is `94.11765%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 99.22%. Comparing base [(`a7253b8`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/a7253b8cb6338e46e1468be1769602771c25e0c7?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`bf301b3`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/bf301b32d2d25509ea647416bd8b84b35216a21e?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).
> Report is 6 commits behind head on develop.

| [Files with missing lines](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4710?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) | Patch % | Lines |
|---|---|---|
| [src/pybamm/spatial\_methods/finite\_volume.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4710?src=pr&el=tree&filepath=src%2Fpybamm%2Fspatial_methods%2Ffinite_volume.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-c3JjL3B5YmFtbS9zcGF0aWFsX21ldGhvZHMvZmluaXRlX3ZvbHVtZS5weQ==) | 94.11% | [1 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4710?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |

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


```diff
@@             Coverage Diff             @@
##           develop    #4710      +/-   ##
===========================================
- Coverage    99.22%   99.22%   -0.01%     
===========================================
  Files          303      303              
  Lines        23070    23085      +15     
===========================================
+ Hits         22891    22905      +14     
- Misses         179      180       +1     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4710?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).

aabills: Thanks @rtimms. I can add a quick test to fix coverage if you want
aabills: > Thanks @rtimms. I can add a quick test to fix coverage if you want

Turns out it's not possible to reach the uncovered line, so deleted it.
aabills: @rtimms can we merge 4665 before this, there will be a conflict that will be easier to resolve here.
rtimms: @aabills sounds good

Patch:
diff --git a/src/pybamm/spatial_methods/finite_volume.py b/src/pybamm/spatial_methods/finite_volume.py
index 9ea609192..47ca8a1ac 100644
--- a/src/pybamm/spatial_methods/finite_volume.py
+++ b/src/pybamm/spatial_methods/finite_volume.py
@@ -264,9 +264,9 @@ class FiniteVolume(pybamm.SpatialMethod):
             The finite volume integral matrix for the domain
         """
         domains = child.domains
-        if vector_type != "row" and integration_dimension == "secondary":
+        if vector_type != "row" and integration_dimension != "primary":
             raise NotImplementedError(
-                "Integral in secondary vector only implemented in 'row' form"
+                f"Integral in {integration_dimension} vector only implemented in 'row' form"
             )
 
         domain = child.domains[integration_dimension]
@@ -282,7 +282,7 @@ class FiniteVolume(pybamm.SpatialMethod):
                 d_edges = 2 * np.pi * (r_edges_right**2 - r_edges_left**2) / 2
         else:
             d_edges = submesh.d_edges
-
+        possible_dimensions = ["primary", "secondary", "tertiary", "quaternary"]
         if integration_dimension == "primary":
             # Create appropriate submesh by combining submeshes in domain
             submesh = self.mesh[domains["primary"]]
@@ -298,29 +298,29 @@ class FiniteVolume(pybamm.SpatialMethod):
             second_dim_repeats = self._get_auxiliary_domain_repeats(domains)
             # generate full matrix from the submatrix
             matrix = kron(eye(second_dim_repeats), d_edges)
-        elif integration_dimension == "secondary":
-            # Create appropriate submesh by combining submeshes in domain
-            primary_submesh = self.mesh[domains["primary"]]
-
-            # Create matrix which integrates in the secondary dimension
-            # Different number of edges depending on whether child evaluates on edges
-            # in the primary dimensions
-            if child.evaluates_on_edges("primary"):
-                n_primary_pts = primary_submesh.npts + 1
-            else:
-                n_primary_pts = primary_submesh.npts
-            int_matrix = hstack([d_edge * eye(n_primary_pts) for d_edge in d_edges])
-
-            # repeat matrix for each node in higher dimensions
-            third_dim_repeats = self._get_auxiliary_domain_repeats(
-                {
-                    k: v
-                    for k, v in domains.items()
-                    if (k == "tertiary" or k == "quaternary")
-                }
+        elif integration_dimension in possible_dimensions[1:]:
+            this_dimension_index = possible_dimensions.index(integration_dimension)
+            # get lower dimensions and the corresponding domains, i.e. if integration_dimension is "secondary",
+            # lower_dimensions is ["primary"] and lower_domains is [child.domains["primary"]]
+            lower_dimensions = possible_dimensions[:this_dimension_index]
+            lower_domains = [child.domains[dimension] for dimension in lower_dimensions]
+            # get higher dimensions, i.e. if integration_dimension is "secondary",
+            # higher_dimensions is ["tertiary", "quaternary"]
+            higher_dimensions = possible_dimensions[this_dimension_index + 1 :]
+            n_lower_pts = 1
+            #  Lower dimensions should be repeated, so add them to the eye matrix
+            for lower_domain, lower_dimension in zip(lower_domains, lower_dimensions):
+                lower_submesh = self.mesh[lower_domain]
+                if child.evaluates_on_edges(lower_dimension):
+                    n_lower_pts *= lower_submesh.npts + 1
+                else:
+                    n_lower_pts *= lower_submesh.npts
+            int_matrix = hstack([d_edge * eye(n_lower_pts) for d_edge in d_edges])
+            # Higher dimensions should be tiled, so repeat the matrix for each higher dimension.
+            higher_repeats = self._get_auxiliary_domain_repeats(
+                {k: v for k, v in domains.items() if (k in higher_dimensions)}
             )
-            # generate full matrix from the submatrix
-            matrix = kron(eye(third_dim_repeats), int_matrix)
+            matrix = kron(eye(higher_repeats), int_matrix)
         # generate full matrix from the submatrix
         # Convert to csr_matrix so that we can take the index (row-slicing), which is
         # not supported by the default kron format


Test patch:
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 3abfa6466..c40f7fdf9 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
@@ -173,7 +173,7 @@ class TestFiniteVolumeIntegration:
             integral_eqn_disc.evaluate(None, const_y), 4 * np.pi * 2.0**3 / 3, decimal=4
         )
 
-    def test_integral_secondary_domain(self):
+    def test_integral_secondary_tertiary_domain(self):
         # create discretisation
         mesh = get_1p1d_mesh_for_testing()
         spatial_methods = {
@@ -186,6 +186,7 @@ class TestFiniteVolumeIntegration:
         # lengths
         ln = mesh["negative electrode"].edges[-1]
         ls = mesh["separator"].edges[-1] - ln
+        lc = mesh["current collector"].edges[-1]
         lp = 1 - (ln + ls)
 
         var = pybamm.Variable(
@@ -197,9 +198,12 @@ class TestFiniteVolumeIntegration:
             },
         )
         x = pybamm.SpatialVariable("x", "positive electrode")
+        z = pybamm.SpatialVariable("z", "current collector")
         integral_eqn = pybamm.Integral(var, x)
+        integral_eqn_tertiary = pybamm.Integral(var, z)
         disc.set_variable_slices([var])
-        integral_eqn_disc = disc.process_symbol(integral_eqn)
+        integral_eqn_disc_secondary = disc.process_symbol(integral_eqn)
+        integral_eqn_disc_tertiary = disc.process_symbol(integral_eqn_tertiary)
 
         submesh = mesh["positive particle"]
         constant_y = np.ones(
@@ -210,16 +214,22 @@ class TestFiniteVolumeIntegration:
                 1,
             )
         )
+        # Test with constant y
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, constant_y),
+            integral_eqn_disc_secondary.evaluate(None, constant_y),
             lp * np.ones((submesh.npts * mesh["current collector"].npts, 1)),
         )
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc_tertiary.evaluate(None, constant_y),
+            lc * np.ones((submesh.npts * mesh["positive electrode"].npts, 1)),
+        )
+        # Test with linear
         linear_in_x = np.tile(
             np.repeat(mesh["positive electrode"].nodes, submesh.npts),
             mesh["current collector"].npts,
         )
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, linear_in_x),
+            integral_eqn_disc_secondary.evaluate(None, linear_in_x),
             (1 - (ln + ls) ** 2)
             / 2
             * np.ones((submesh.npts * mesh["current collector"].npts, 1)),
@@ -229,16 +239,25 @@ class TestFiniteVolumeIntegration:
             mesh["positive electrode"].npts * mesh["current collector"].npts,
         )
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, linear_in_r).flatten(),
+            integral_eqn_disc_secondary.evaluate(None, linear_in_r).flatten(),
             lp * np.tile(submesh.nodes, mesh["current collector"].npts),
         )
         cos_y = np.cos(linear_in_x)
         np.testing.assert_array_almost_equal(
-            integral_eqn_disc.evaluate(None, cos_y),
+            integral_eqn_disc_secondary.evaluate(None, cos_y),
             (np.sin(1) - np.sin(ln + ls))
             * np.ones((submesh.npts * mesh["current collector"].npts, 1)),
             decimal=4,
         )
+        # test tertiary integration with linear function in y
+        linear_in_z = np.repeat(
+            mesh["current collector"].nodes,
+            submesh.npts * mesh["positive electrode"].npts,
+        )
+        np.testing.assert_array_almost_equal(
+            integral_eqn_disc_tertiary.evaluate(None, linear_in_z),
+            (lc**2 / 2) * np.ones((submesh.npts * mesh["positive electrode"].npts, 1)),
+        )
 
     def test_integral_primary_then_secondary_same_result(self):
         # Test that integrating in r then in x gives the same result as integrating in