Skip to content

Generalise VectorField to N components - #5686

Merged
aabills merged 6 commits into
mainfrom
ufv-0-vector-field
Aug 5, 2026
Merged

Generalise VectorField to N components#5686
aabills merged 6 commits into
mainfrom
ufv-0-vector-field

Conversation

@aabills

@aabills aabills commented Jul 31, 2026

Copy link
Copy Markdown
Contributor

Summary

  • Generalise VectorField from fixed 2D (lr, tb) to N components (needed for 3D unstructured FV).
  • Add Component / Norm operators and matching discretisation + solution handling.
  • First of a stacked split of Unstructured finite volume #5397 (plan B).

Stack

  1. This PR (Generalise VectorField to N components #5686) — VectorField N-comp
  2. Add unstructured mesh infrastructure #5687 — Meshing
  3. 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

Full pre-split branch preserved as backup/unstructured-finite-volume-full and original #5397 (unstructured-finite-volume).

Test plan

  • test_finite_volume_2d/test_tensor_field.py (includes new VectorField N-comp cases)
  • CI unit suite

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

Support 3D vector fields via N-component VectorField, Component/Norm
operators, and matching discretisation/solution handling. Extracted from
the unstructured finite-volume work for review in isolation.

Co-authored-by: Cursor <cursoragent@cursor.com>
@aabills
aabills requested a review from a team as a code owner July 31, 2026 21:39
@codecov

codecov Bot commented Jul 31, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.11%. Comparing base (1827bda) to head (bef253d).

Additional details and impacted files
@@            Coverage Diff             @@
##             main    #5686      +/-   ##
==========================================
+ Coverage   98.09%   98.11%   +0.02%     
==========================================
  Files         340      340              
  Lines       32671    32732      +61     
==========================================
+ Hits        32049    32116      +67     
+ Misses        622      616       -6     

☔ 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.

Add unit tests for the uncovered patch lines: Component/Norm
discretisation (success and error paths), _unary_new_copy and the
component/norm convenience functions, _disc_state_vector propagation
through binary/unary operators, and the per-component casadi handling
for VectorField variables in Solution.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

@rtimms rtimms left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reviewed this PR, then checked the findings against the rest of the stack (ufv-1ufv-5) to avoid raising things that are already handled downstream. Verified behaviour by running the code at both this branch and the stack tip; targeted unit suites are green at both (1733 serialisation + expression-tree tests, 755 spatial-method/solver/discretisation tests).

Three things I think need fixing before this lands, then some dead surface and housekeeping.

Already covered downstream — ignore these

For the record, so they don't get re-raised:

  • Magnitude vs Component duplication. Magnitude(vf, "lr"|"tb") is exactly Component(vf, 0|1), and Magnitude is misnamed (it extracts a component; the new Norm is the actual magnitude). #5691 handles this properly — deprecation warning, both BasicDFN2D call sites migrated, FiniteVolume2D._edge_direction so edge-averaging still resolves a direction. Deprecating is the right call for a public symbol.
  • _disc_state_vector has a real writer. It's set at finite_volume_unstructured.py:617 in #5688, with genuine coverage at test_finite_volume_unstructured.py:874. See the note below about where it belongs, but the concept isn't dead.

Blockers

1. Mismatched component counts silently drop data

n is taken from whichever side happens to be a VectorField, with no check that both agree (discretisation.py:973-975). Once N can vary, that's a correctness hole:

2-comp + 3-comp -> returns 2 components; the third is silently dropped, no error
3-comp + 2-comp -> IndexError: list index out of range

Both reproduce on this branch. Suggest an explicit check before the loop:

if disc_left.n_components != disc_right.n_components:
    raise pybamm.DiscretisationError(
        f"Cannot combine VectorFields with {disc_left.n_components} and "
        f"{disc_right.n_components} components"
    )

Worth fixing here rather than downstream, because finite_volume_unstructured.py:948-977 in #5688 is a near-verbatim second copy of this block — same n logic, same missing guard, same hasattr/break. Either fix both, or factor the "broadcast the scalar side, then zip components" step into one shared helper that both call.

2. solution["<vector field>"] is unreadable on structured 2D meshes

The solution.py change stores a list of casadi functions, but process_variable routes structured 2D FV meshes to ProcessedVariable2DFVM at processed_variable.py:1750, before the VectorField dispatch — and that dispatch is gated on UnstructuredSubMesh, so it never catches this case even at the tip of the stack. Result, verified on both this branch and ufv-5:

processed class: ProcessedVariable2DFVM
call(0.5)  -> TypeError: unhashable type: 'list'
entries    -> TypeError: unhashable type: 'list'
data       -> TypeError: unhashable type: 'list'

Not a functional regression — on main this path raises TypeError: Cannot convert symbol of type VectorField to CasADi — but it swaps a clear message for an opaque one, on a path that stays broken through all six PRs. Either hoist an isinstance(base_variables[0], pybamm.VectorField) check above line 1750, or raise NotImplementedError there so the failure names the actual limitation.

Also worth noting: the added test only asserts isinstance(casadi_components, list), which is why it passes despite every read path raising. If the plumbing stays in this PR, the test should read a value.

Separately, if you'd rather drop the solution.py change from this PR and land it with its consumer in #5688, drop VectorField._to_casadi with it — with _to_casadi present and the special case gone, solution["flux"] silently returns wrong shapes instead of erroring.

3. hasattr sniffing of a private attribute across objects

result = pybamm.VectorField(*new_comps)
for src in (disc_left, disc_right):
    if hasattr(src, "_disc_state_vector"):
        result._disc_state_vector = src._disc_state_vector
        break
return result

By the tip of the stack this pattern exists in three places (discretisation.py:993, discretisation.py:1194, finite_volume_unstructured.py:974). Declaring the attribute removes the duck-typing from all three:

class VectorField(TensorField):
    _disc_state_vector = None
result._disc_state_vector = disc_left._disc_state_vector or disc_right._disc_state_vector
# and in the unary branch
result._disc_state_vector = disc_child._disc_state_vector

Given the only writer is in #5688, the propagation plus test_disc_state_vector_propagation would sit more naturally there — the test here has to fabricate the attribute by hand (disc_vf._disc_state_vector = marker), so it doesn't exercise anything this PR ships.

Dead surface

Grepped all six branches; none of the following has a single use in src/:

  • fb_field — unused even by the unstructured 3D models in #5690, and vf[2] already works via TensorField.__getitem__. It's also documented as a "backward-compatible alias" when it's brand new. Suggest dropping it and its two tests until something needs a third component.
  • pybamm.component() / pybamm.norm() — zero call sites anywhere in the stack; every real caller writes pybamm.Component(N_e, 0). They're also bare constructor aliases with no simplify_if_constant, unlike neighbours such as sign().
  • Norm — zero uses in src/ across the stack, and finite_volume_unstructured.gradient_squared (#5688, lines 936-942) hand-rolls exactly the sum-of-squares that the Norm discretisation branch builds. Either give Norm that caller or defer it.

Reuse and style

  • TensorField already exposes components and __getitem__, but private _components access grows from 7 sites in this PR to 13 in src/ by the tip (processed_variable.py:1232,1238, finite_volume_unstructured.py:770,940,967). Worth switching to the public accessors while the count is still small.
  • VectorField._to_casadi reimplements the inherited helper; casadi.vertcat(*self._children_to_casadi(t, y, y_dot, inputs, casadi_symbols)) is equivalent.
  • The Norm discretisation branch reads more directly as return sum(c**2 for c in disc_child.components) ** 0.5.
  • n_components is a third spelling of len(components) / shape[0]. Fine to keep, but the codebase should settle on one.
  • New disc branches raise bare ValueError; AGENTS.md asks for DiscretisationError. (The adjacent Magnitude branch sets the precedent, but it'd be good not to extend it.)

Housekeeping

  • The CHANGELOG bullet has no PR link — AGENTS.md requires one, and the surrounding bug-fix bullets have them. Same applies to the other four feature bullets added across the stack.
  • New public Component / Norm have no docs/source/api/expression_tree/unary_operator.rst entry. (VectorField / TensorField were already undocumented, so this is arguably pre-existing.)

@aabills

aabills commented Aug 3, 2026

Copy link
Copy Markdown
Contributor Author

Thanks for the thorough review — and for checking against the rest of the stack so we don't churn on things already handled downstream.

Blockers

  1. Mismatched component counts — Agreed. Added a DiscretisationError guard and factored the "broadcast scalar side, then zip components" step into _process_vector_field_binary so the same logic can be reused downstream.

  2. solution["<vector field>"] on structured 2D — Raising a clear NotImplementedError naming the limitation (instead of the opaque TypeError: unhashable type: 'list'). Updated the test to assert that failure mode.

  3. hasattr / _disc_state_vector — Declared _disc_state_vector = None on VectorField and assign directly. Propagation coverage for the real writer can land with Add unstructured finite volume spatial method #5688.

Dead surface

  • Dropped fb_field (and its tests); vf[2] is enough for now.
  • Dropped pybamm.component() / pybamm.norm() — callers use the constructors.
  • Keeping Norm.

Reuse / style / housekeeping

  • Switched to public components / __getitem__ instead of _components.
  • VectorField._to_casadi now uses _children_to_casadi.
  • Simplified the Norm disc branch to sum(c**2 for c in disc_child.components) ** 0.5.
  • New disc branches raise DiscretisationError instead of bare ValueError.
  • CHANGELOG bullet now links to this PR; Component / Norm added to the unary-operator API docs.

Will push the follow-up shortly.

aabills and others added 2 commits August 3, 2026 11:54
Guard mismatched component counts, raise a clear error for structured-2D
VectorField solution reads, and drop unused fb_field / convenience helpers.

Co-authored-by: Cursor <cursoragent@cursor.com>
`or` bool-evaluates StateVector and raises; use an explicit None check.

Co-authored-by: Cursor <cursoragent@cursor.com>
rtimms
rtimms previously approved these changes Aug 4, 2026

@rtimms rtimms left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks!

@aabills
aabills enabled auto-merge (squash) August 5, 2026 18:16
@aabills
aabills merged commit 863922c into main Aug 5, 2026
32 checks passed
@aabills
aabills deleted the ufv-0-vector-field branch August 5, 2026 18:17
aabills added a commit that referenced this pull request Aug 6, 2026
Adds FiniteVolumeUnstructured (TPFA Laplacian, fused div_D_grad,
Green-Gauss gradient), unstructured processed variables, and the
discretisation dispatch for div(D*grad(u)) and graph-topology internal
boundary conditions, on top of the unstructured meshing (#5687) and
N-component VectorField (#5686) already on main.

Includes the fixes from the review of the previous revision:
- Neumann values on named axis sides are coordinate-direction
  derivatives (matching FiniteVolume); custom tags stay outward-normal
- unknown BC sides and bc_types raise DiscretisationError instead of
  being silently dropped; boundary_integral supports "entire"
- auxiliary domains are handled in laplacian/gradient BC assembly and
  in secondary/tertiary broadcasts
- divergence of a BC-bearing flux raises rather than silently dropping
  the boundary flux; the Green-Gauss gradient warns for BC-less buckets
- scalar reductions (Max/Min) on unstructured domains post-process as 0D
- operator matrices are cached per submesh (invalidated when the cell
  ordering changes) and BC assembly is vectorised
- ParameterSubstitutor.process_boundary_conditions processes every
  boundary side present, not a fixed whitelist, so tab and named-region
  Dirichlet conditions are no longer silently dropped

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
aabills added a commit that referenced this pull request Sep 5, 2026
* feat: add unstructured finite volume spatial method

Adds FiniteVolumeUnstructured (TPFA Laplacian, fused div_D_grad,
Green-Gauss gradient), unstructured processed variables, and the
discretisation dispatch for div(D*grad(u)) and graph-topology internal
boundary conditions, on top of the unstructured meshing (#5687) and
N-component VectorField (#5686) already on main.

Includes the fixes from the review of the previous revision:
- Neumann values on named axis sides are coordinate-direction
  derivatives (matching FiniteVolume); custom tags stay outward-normal
- unknown BC sides and bc_types raise DiscretisationError instead of
  being silently dropped; boundary_integral supports "entire"
- auxiliary domains are handled in laplacian/gradient BC assembly and
  in secondary/tertiary broadcasts
- divergence of a BC-bearing flux raises rather than silently dropping
  the boundary flux; the Green-Gauss gradient warns for BC-less buckets
- scalar reductions (Max/Min) on unstructured domains post-process as 0D
- operator matrices are cached per submesh (invalidated when the cell
  ordering changes) and BC assembly is vectorised
- ParameterSubstitutor.process_boundary_conditions processes every
  boundary side present, not a fixed whitelist, so tab and named-region
  Dirichlet conditions are no longer silently dropped

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* fix: tile per-face BC vectors across auxiliary-domain repeats in div_D_grad

A BC value with one entry per boundary face broke with a ShapeError when
the variable had auxiliary domains; lift it to n_bnd * repeats entries,
matching _bc_contribution's convention.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* fix: raise on y= queries of 2D unstructured processed variables

y was silently dropped when building the query grid, returning
z-midplane values for a query the user thinks is at y; raise like the
existing r/R check. 3D unchanged.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* perf: interpolate all time steps of unstructured variables in one pass

scipy's interpolators accept (n_points, n_t) value arrays, so the
nearest-neighbour KDTree and the domain boundary mask are now built once
per query instead of once per time step (~40x on a 100-step animation).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* fix: raise DomainError for unresolvable unstructured spatial variables

spatial_variable silently fell back to the x column for unknown names
and directions, matched loose prefixes (zeta -> z), and accepted y/fb on
2D x-z meshes. Match the leading name token exactly and raise otherwise.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* refactor: define set_internal_bcs_for_concat on the base SpatialMethod

Replace the hasattr duck-typing dispatch with a documented base-class
hook returning None (use the default 1D-stack routine); the unstructured
method's implementation becomes a plain override.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* fix: match legacy tab BC sides exactly instead of by substring

With arbitrary boundary tags now allowed, a region name containing
'tab' (e.g. 'tab_weld') was misrouted into check_tab_conditions and
raised ModelError; share one exact-name set across both check sites.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* style: pre-commit fixes

* fix: fail loudly when unstructured meshes lack tags or interface data

build() now warns for a submesh with exterior faces but no boundary
tags (BCs cannot apply, interface discovery cannot pair it), and
_internal_neumann_unstructured raises instead of returning zeros when
no interface pairs two meshes, which silently decoupled the domains.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Add implicit non-orthogonal correction to the unstructured TPFA operators

The two-point flux only carried alpha (u_j - u_i)/d with alpha = cos(theta)
and dropped the k . grad(u) remainder of n = alpha e + k, so the Laplacian
of a linear field was nonzero on any non-orthogonal mesh and the scheme did
not converge on triangles or tetrahedra (L2 rate ~0 under refinement).

- Split the face normal as n = alpha e_ij + k and add the k . grad(u)_f
  cross flux implicitly, for internal faces (laplacian, div_D_grad) and
  Dirichlet boundary faces (perpendicular distance delta . n).
- "non-orthogonal correction" option: "over-relaxed" (default,
  alpha = 1/cos theta, floored) or "minimum" (alpha = cos theta).
- Replace the Green-Gauss gradient with a batched weighted least-squares
  reconstruction that is exact on linear fields on skewed meshes; internal,
  Dirichlet and Neumann faces each contribute one directional-derivative
  row. Green-Gauss remains the divergence assembly.
- Raise GeometryError for inverted cells, warn at build above 70 degrees
  of non-orthogonality, skip the correction entirely on orthogonal meshes.
- _divergence_matrices now aliases _green_gauss_matrices.

Linear patch test is at 1e-13 on tri, perturbed tri and Kuhn tet meshes;
Poisson convergence rates are 1.95-2.04 (tri) and 1.74-1.88 (tet) for both
options, unchanged (2.0) on quad/hex.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Test the build-time non-orthogonality warning; make zip strict

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Return a pybamm.Matrix from the unstructured definite integral

definite_integral_matrix returned a raw csr_matrix, so DefiniteIntegralVector
failed with AttributeError in process_symbol's shape check; it also ignored
vector_type (a "column" request came back as a row) and integral() dropped
integration_dimension, silently integrating the primary volume for a
secondary-dimension Integral. Support row/column, raise NotImplementedError
for non-primary dimensions, and lift by the auxiliary-domain repeats.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Move point-in-domain onto UnstructuredSubMesh; skip nearest fill outside it

UnstructuredSubMesh.contains_points now owns the 2D even-odd loop test
(loops cached on the mesh) and delegates to contains_points_3d, so the
processed variable only asks the mesh. In _interpolate_spatial the outside
mask is computed first and excluded from the nearest-neighbour fill, which
previously filled ~97% of out-of-hull points only to overwrite them with
fill_value. Also correct the get_quiver_data docstring for the 3D return.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Apply the non-orthogonal correction across unstructured domain interfaces

The interface Neumann value between concatenated domains was the plain
two-point difference (u_r - u_l)/d, the same alpha = 1 truncation that made
the Laplacian inconsistent inside a domain, so a linear field was not a
steady state of a two-domain tet model (residual 4.6 on the Kuhn mesh).

- internal_neumann_condition now returns alpha (u_r - u_l)/d + k . grad(u)_f
  with the same n = alpha e + k split as the interior; the face gradient is
  the distance-weighted mean of both sides' least-squares gradients, which
  take the interface faces as cross-mesh rows towards the paired cell and
  each side's external boundary conditions (passed by
  set_internal_bcs_for_concat). Orthogonal interfaces are unchanged.
- interface_data records left_faces/right_faces on both pairing paths.
- Least-squares rows for boundary faces without a condition now fit a zero
  normal derivative, matching the operators' zero-flux treatment, so cells
  on untagged boundaries stay fully determined.
- UnstructuredSubMesh.contains_points gets a direct test.

Two-domain linear field on tets: residual 4.6 -> 5e-14; interface flux is
conservative to 1e-10 for a random field.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Test the corrected interface gradient on a linear field

The old expectation hard-coded the two-point (u_r - u_l)/d formula, which
is only the full interface gradient on orthogonal pairs; check u = x gives
exactly 1 on the skewed triangle interface and keep the two-point check
for a quad interface where it still holds.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Rename the axis token variable that Bandit mistakes for a password

Codacy's Bandit B105 flags `token == "y"` as a hardcoded password purely
because of the variable name.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* docs: drop duplicate Component/Norm entries after merging main

main already documents both in unary_operator.rst; the second copy made
Sphinx fail with "duplicate object description" under -W.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Cover the inconsistent-face-count error and the no-loops containment path

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Treat every scalar-shaped boundary value as a broadcast scalar

Time- or input-dependent scalars (e.g. a current-density Neumann value)
evaluate for shape to () rather than (1, 1), so _bc_contribution built
Matrix @ scalar and failed once div_D_grad started assembling the
least-squares gradient with boundary rows. Centralise the check.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Use the harmonic mean of D at unstructured faces in div(D * grad(u))

The face coefficient was a linear interpolation of the two cell values.
For a flux the physically consistent face value is the resistances-in-
series (distance-weighted harmonic) mean, which FiniteVolume already uses
for coefficients of a gradient: it reproduces the exact two-cell flux
across a material interface and gives zero flux into an impermeable cell,
whereas the arithmetic mean lets half the conductive side leak through.
On combined concatenation meshes this is the current-collector/electrode
face, where sigma jumps by ~1e5. Linear weights are kept for the face
gradient of the cross term; D must be strictly positive.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Treat faces orthogonal to within 1e-8 as orthogonal, not 1e-12

Centroid rounding on high-aspect-ratio cells (10 um thick, cm wide, as in
a pouch cell) puts ~1e-11 into the face direction, so with a 1e-12 cut-off
the cross term and boundary cross terms were assembled on every hex mesh
with coefficients of 1e-11. That doubled the Jacobian stencil (57k -> 101k
nonzeros on the pouch demo) and the IDAKLU time for no change in the
answer. k is now zeroed below a dimensionless 1e-8 (an angle far below any
real skew) and explicit zeros are dropped, so orthogonal faces never widen
the stencil, also on mixed meshes.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Keep plotting out of the unstructured processed variables

The visualisation grid (N_VIS/N_VIS_3D and the first/second/third_dim_pts
it populated), mid-plane 3D slices, quiver sampling and slice positions
only exist to feed QuickPlot, whose unstructured support lives in the
plotting PR. Processed variables here now only interpolate at requested
points; QuickPlot refuses unstructured variables with a clear message until
that support lands, instead of failing deep inside its 2D branch.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants