Add unstructured finite volume spatial method

[aabills]

Summary

• Add FiniteVolumeUnstructured spatial method, discretisation hooks (div_D_grad, internal BC dispatch), and unstructured ProcessedVariable* classes.
• Depends on Generalise VectorField to N components #5686 (VectorField) and Add unstructured mesh infrastructure #5687 (Meshing); both are merged into this branch.

Stack

1. Generalise VectorField to N components #5686 — VectorField N-comp
2. Add unstructured mesh infrastructure #5687 — Meshing
3. This PR (Add unstructured finite volume spatial method #5688) — Spatial method + ProcessedVariable
4. Add VTK plotting for unstructured meshes #5689 — VTK plotting
5. Add unstructured 2D/3D DFN battery models #5690 — Unstructured DFN models

Review tip: the spatial-only commit is feat: add unstructured finite volume spatial method; earlier commits are the two dependencies.

Test plan

• test_finite_volume_unstructured.py + mesh + tensor-field tests (158 passed locally)
• CI unit suite

Also in this stack: #5691 deprecates pybamm.Magnitude.

[codecov]

Codecov Report

❌ Patch coverage is 98.67752% with 13 lines in your changes missing coverage. Please review.
✅ Project coverage is 98.15%. Comparing base (78fe10a) to head (f80e2ec).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
...es/pybamm/src/pybamm/solvers/processed_variable.py	95.27%	6 Missing ⚠️
...ybamm/src/pybamm/discretisations/discretisation.py	90.69%	4 Missing ⚠️
...bamm/spatial_methods/finite_volume_unstructured.py	99.62%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #5688      +/-   ##
==========================================
- Coverage   98.21%   98.15%   -0.06%     
==========================================
  Files         340      341       +1     
  Lines       32743    33712     +969     
==========================================
+ Hits        32158    33090     +932     
- Misses        585      622      +37     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[MarcBerliner]

I think these should be defined in the visualization, not in the processed variables

[aabills]

Agreed, and done. The visualisation grid (N_VIS/N_VIS_3D and the *_dim_pts it populated), the 3D mid-plane slices, the quiver sampling and the slice-plane positions are gone from this PR: the unstructured processed variables are now pure point interpolators (pv(t, x=, y=, z=)), and QuickPlot here raises a clear NotImplementedError for unstructured variables. The sampling lives on the plotting side in #5689 as pybamm.plotting.unstructured_plot_grid (plot_grid, default_slice_positions, midplane_slices, quiver_data), and QuickPlot there owns the grid and slice positions per plotted variable instead of mutating the variable.

[MarcBerliner]

I think this functionality is fine but should live in a mesh processing file, not the processed variables

[aabills]

Moved: point-in-domain now lives on the mesh as UnstructuredSubMesh.contains_points (2D even-odd rule over the cached boundary loops, 3D delegating to contains_points_3d), with a direct test. The processed variable just asks the mesh.

[MarcBerliner]

same

[aabills]

Same treatment for the vector-field processed variable: N_QUIVER/get_quiver_data and the copied grid attributes are removed here; it only wraps the per-component interpolators. Quiver sampling is quiver_data in #5689's plotting module.

[MarcBerliner]

The docstring calls this "the non-orthogonality correction," but multiplying by |n·e| isn't one — it drops the tangential (cross-diffusion) part of ∇u·n entirely, and that term is O(1) on a fixed-angle non-orthogonal grid, not O(h). Plain TPFA on cell centroids is only consistent when the mesh is orthogonal (quad/hex here); on triangles/tets it doesn't converge at all — see the numbers in the review summary.

A real correction needs an over-relaxed decomposition with an explicit gradient-reconstructed cross term, or MPFA. Alternatively, restrict the method to orthogonal quad/hex meshes — but then the module docstring ("unstructured simplex meshes (2D triangles / 3D tets)"), the class docstring ("Supports triangles (2D) and tetrahedra (3D)"), and the generator's 2D default all need to change to match.

Same expression at line 542 in _div_D_grad_matrices.

[aabills]

You were right, and it was worse than the review summary suggested: on the structured tri/tet meshes there is no error cancellation at all (Poisson L2 rate ≈ 0 under refinement, Laplacian of a linear field = 2.9 on the default 2D mesh).

Implemented the full decomposition n = α ê + k, fully implicit since everything is linear: the compact part carries α (u_j − u_i)/d and the cross term k·∇u_f is a matrix built from a cell-gradient reconstruction. α is an option ("non-orthogonal correction": "over-relaxed" = 1/cosθ default, "minimum" = cosθ); both are exact on linear fields, the choice only sets the implicit/explicit split. The gradient had to become a weighted least-squares fit rather than Green-Gauss, because Green-Gauss with distance weights is itself O(1) wrong on tets (0.74 error on a linear field on the Kuhn mesh — face centroids are not on the centroid lines). The same correction is applied across concatenated-domain interfaces (cross-mesh LS rows), and faces orthogonal to within 1e-8 are skipped so hex/quad meshes are unchanged.

Result: linear patch test 1e-13 on tri, perturbed tri, Kuhn tet and perturbed tet; Poisson rates 1.95–2.04 (tri) and 1.74–1.88 (tet) for both α; quad/hex still 2.0. Docstrings now describe what is actually computed.

[MarcBerliner]

d_perp is the Euclidean distance |face_centroid − cell_centroid|, not the perpendicular one — the name asserts the opposite of what the code does. The correct boundary transmissibility is A / (Δ·n), i.e. this divided by cos θ. Same at line 651 in div_D_grad. Related to the interior issue on _tpfa_matrix, but independently fixable.

[aabills]

Fixed. The boundary term now uses the perpendicular distance Δ·n̂ (i.e. A/(Δ·n̂) for the over-relaxed α, consistent with the interior decomposition) plus its own k·∇u cross term from the owner cell's least-squares gradient. The Neumann-side gradient rows use the normal direction directly.

[MarcBerliner]

This is identical to _green_gauss_matrices (line 798) apart from local variable names — I confirmed the assembled matrices match exactly. ~55 duplicated lines, and they're cached under two separate keys so both get built. Could be _divergence_matrices = _green_gauss_matrices (or one delegating to the other) with a comment on why divergence and Green-Gauss gradient share an assembly.

[aabills]

Aliased: _divergence_matrices now returns _green_gauss_matrices(submesh) with a note on why the two assemblies coincide. (Green-Gauss stays as the divergence assembly; the gradient operator itself moved to least squares, see the thread on the correction.)

[MarcBerliner]

This returns a raw scipy.sparse.csr_matrix rather than a pybamm.Matrix. Discretisation._process_symbol hands it straight back for DefiniteIntegralVector, and process_symbol then calls .test_shape() on it:

AttributeError: 'csr_matrix' object has no attribute 'test_shape'

Reproduces for both vector_type="row" and "column" — so DefiniteIntegralVector is unusable on any unstructured mesh. Untested.

Separately, vector_type is accepted and silently ignored: Discretisation passes vector_type=symbol.vector_type, so a "column" request quietly gets a (1, n) row. FiniteVolume raises NotImplementedError for the unsupported combination.

[aabills]

Fixed: definite_integral_matrix returns a pybamm.Matrix (lifted by the auxiliary-domain repeats), vector_type is honoured for "row"/"column" and rejected otherwise, and there is a Discretisation-level test that DefiniteIntegralVector processes in both orientations.

[MarcBerliner]

integration_dimension is accepted and dropped — Discretisation passes symbol._integration_dimension here, so an Integral over a secondary dimension silently integrates the primary volume instead. Worth raising rather than being silently wrong, as the rest of this file does.

[aabills]

Fixed: integral now passes integration_dimension through and definite_integral_matrix raises NotImplementedError for anything but "primary" instead of silently integrating the cell volume; tested.

[MarcBerliner]

The nearest fill is mostly wasted work: out-of-hull points are filled here and then overwritten with fill_value two lines down. On the default 200×200 grid over a convex domain I measured 25,220 of 40,000 query points nearest-filled, 24,380 of them (97%) immediately NaN'd. Computing outside first and filling only mask & ~outside avoids building the NearestNDInterpolator in the common case.

[aabills]

Fixed: the outside mask is computed first and excluded from the nearest-neighbour fill, so NearestNDInterpolator is only built when some in-domain point actually fell outside the hull.

[MarcBerliner]

The 3D branch actually returns 10 values — (X1, Z1, u_xz, w_xz, y_mid, X2, Y2, u_xy, v_xy, z_mid) — not (X, Y, Z, U, V, W).

[aabills]

Docstring corrected — and the function has since moved out of the processed variable entirely (see the thread at the top of the class); the 3D return is documented as the two mid-plane slices there.

[MarcBerliner]

Thanks!