Skip to content

SLCN at the integration points: vector and tensor histories - #720

Merged
lmoresi merged 4 commits into
developmentfrom
feature/slcn-ip-tensor
Sep 11, 2026
Merged

SLCN at the integration points: vector and tensor histories#720
lmoresi merged 4 commits into
developmentfrom
feature/slcn-ip-tensor

Conversation

@lmoresi

@lmoresi lmoresi commented Sep 10, 2026

Copy link
Copy Markdown
Member

IntegrationPointSemiLagrangian refused anything but a scalar:

if vtype != VarType.SCALAR:
    raise NotImplementedError(
        "IntegrationPointSemiLagrangian: scalar histories only for now")

so a Navier-Stokes momentum history or a viscoelastic stress history could not use the integration-point route at all. Lagrangian_Swarm has carried both since the VE stress history (test_0070); this brings the two level. The choice between them is where the state lives, not what shape it can take.

What it took

The trace-back, the characteristic cache and the weighted sums were already shape-agnostic — only the fills were scalar-bound, at four sites. The one real subtlety is that a shaped field stores one dof per independent component rather than one per matrix entry:

vtype (2-D) symbolic shape stored columns
SCALAR 1×1 1
VECTOR 1×2 2
SYM_TENSOR 2×2 3
DFDt = uw.systems.ddt.IntegrationPointSemiLagrangian(
    mesh, stress, v.sym, vtype=uw.VarType.SYM_TENSOR, degree=2, order=1)

DFDt.psi_star[0].sym          # 2x2 symbolic matrix
DFDt.psi_star[0].data.shape   # (npoints, 3)
DFDt.bdf()                    # 2x2, as psi_fn is

The column order was measured, not assumed: diagonal first, then the off-diagonals in row-major upper-triangular order — (0,0), (1,1), (0,1) in 2-D and (0,0), (1,1), (2,2), (0,1), (0,2), (1,2) in 3-D. Getting it wrong transposes a stress in silence, so _storage_components is pinned by a test against what the variable's own .sym reconstructs, in both dimensions.

Fills are now component-wise (_write_components), because a symmetric tensor's symbolic form repeats its off-diagonals and only the independent columns exist in storage.

Accuracy

The scalar property, per component: with a uniform velocity and a field in the P2 space, every slot holds the snapshot evaluated at the exact departure point to < 1e-12, for one segment and for two. The scalar path is unchanged, and a test asserts that.

The tensor test uses three distinct independent components so a packing error or a transposed off-diagonal cannot pass.

Docs

The documentation contradicted itself, which is likely why this never reached anyone: docs/developer/subsystems/integration-point-variables.md said "Scalar components only for now; use one variable per component" while enhanced_variables.py said vector and tensor were supported. Both now agree, and there is a new Vector and tensor histories section with the table, the storage convention, and a pointer to Lagrangian_Swarm for the same history on particles.

Tests

tests/test_0068_integration_point_slcn_tensor.py (10). Full level_1 and tier_a: 1197 passed, 3 skipped, 1 xfailed (development is 1187 + these 10).

Underworld development team with AI support from Claude Code

🤖 Generated with Claude Code

https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G

IntegrationPointSemiLagrangian refused anything but a scalar:

    if vtype != VarType.SCALAR:
        raise NotImplementedError(
            "IntegrationPointSemiLagrangian: scalar histories only for now")

so a Navier-Stokes momentum history or a viscoelastic stress history could
not use it. Lagrangian_Swarm has carried both since the VE stress history
(test_0070); this brings the integration-point route level with it. The
choice between them is where the state lives, not what shape it can take.

The trace-back, the characteristic cache and the weighted sums were already
shape-agnostic. Only the fills were scalar-bound, at four sites, and the one
real subtlety is that a shaped field stores one dof per INDEPENDENT
component rather than one per matrix entry: a symmetric tensor in 2-D is 2x2
symbolically and THREE columns in storage.

  vtype (2-D)   symbolic   stored columns
  SCALAR        1x1        1
  VECTOR        1x2        2
  SYM_TENSOR    2x2        3

The column order was measured, not assumed: diagonal first, then the
off-diagonals in row-major upper-triangular order, in 2-D ((0,0), (1,1),
(0,1)) and in 3-D ((0,0), (1,1), (2,2), (0,1), (0,2), (1,2)). Getting it
wrong transposes a stress in silence, so _storage_components is pinned by a
test against what the variable's own .sym reconstructs, in both dimensions.

Fills are now component-wise (_write_components) because a symmetric
tensor's symbolic form repeats its off-diagonals and only the independent
columns exist in storage.

Accuracy is the scalar property, per component: with a uniform velocity and
a field in the P2 space every slot holds the snapshot evaluated at the exact
departure point to < 1e-12, for one segment and for two. The scalar path is
unchanged and asserted so.

The docs contradicted themselves, which is likely why this never reached
anyone: the subsystem page said "Scalar components only for now; use one
variable per component" while enhanced_variables.py said vector and tensor
were supported. Both now agree, and there is a "Vector and tensor histories"
section with the table, the storage convention and the Lagrangian_Swarm
pointer.

Tests: test_0068_integration_point_slcn_tensor.py (10).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
Copilot AI lite review requested due to automatic review settings September 10, 2026 18:23

Copilot AI 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.

🟡 Changes recommended

IntegrationPointSemiLagrangian currently derives SYM_TENSOR storage components using mesh.dim instead of mesh.cdim (breaking manifold meshes) and should add an early psi_fn shape vs vtype validation to prevent incorrect/misleading behaviour for expression inputs.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

This PR extends uw.systems.ddt.IntegrationPointSemiLagrangian so it can store and transport vector and tensor histories at integration points (previously scalar-only), aligning its capabilities with Lagrangian_Swarm and enabling integration-point histories for Navier–Stokes momentum and viscoelastic stress.

Changes:

  • Add a component-mapping helper (_storage_components) and a component-wise writer (_write_components) to support vector/tensor (including symmetric tensor) storage semantics.
  • Update integration-point history slots (psi_star) and nodal snapshots (psi_snap) to be created with the requested vtype, and fill them component-by-component.
  • Add focused tests for storage order, vector/tensor transport accuracy, scalar regression, and expression usability; update developer docs to describe the storage convention.
File summaries
File Description
src/underworld3/systems/ddt.py Generalises IntegrationPointSemiLagrangian to allocate/fill vector/tensor histories and introduces component mapping + per-component fill logic.
tests/test_0068_integration_point_slcn_tensor.py Adds tests that pin storage packing order and verify transport accuracy for vector and symmetric-tensor histories (plus scalar regression).
docs/developer/subsystems/integration-point-variables.md Updates subsystem documentation to describe vector/tensor support and the independent-component storage convention.
Review details
  • Files reviewed: 3/3 changed files
  • Comments generated: 1
  • Review effort level: Lite

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment on lines +4369 to +4378
self.num_components = int(self.psi_star[0].num_components)
self._components = _storage_components(
vtype, tuple(self.psi_star[0].sym.shape), mesh.dim
)
if len(self._components) != self.num_components:
raise RuntimeError(
f"IntegrationPointSemiLagrangian: {vtype} maps "
f"{len(self._components)} components onto "
f"{self.num_components} stored columns"
)
Two findings, both real.

`_storage_components` took `mesh.dim`, which is wrong on a manifold: a
spherical surface is dim 2, cdim 3, and the variable sizes its vector/tensor
storage by the EMBEDDING dimension. A symmetric tensor there is 3x3
symbolically with six columns, and the dim-2 map would have built three and
tripped the length check. The tensor dimension is now read off the symbolic
shape, which is what the variable is shaped by, so the map never touches the
mesh at all.

There was also no check that psi_fn's shape matches the caller's vtype: a
scalar expression with vtype=VECTOR would have failed later with an
IndexError from the component writer, or stored the wrong thing. It now
raises with the shape it got and the shape it needs — and it raises BEFORE
any variable is allocated, since a mesh variable created and then abandoned
leaves its field on the DM (#1058). The test asserts the variable count is
unchanged after a refusal.

Tests: 12 (was 10).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
@lmoresi

lmoresi commented Sep 10, 2026

Copy link
Copy Markdown
Member Author

Review from the 2026-09 style/clutter audit. One finding, scoped to lines this PR adds.

ddt.py:4526 — flat .data where §7 mandates .array.

var.data[:, column] = np.asarray(...)

The value written has been through np.asarray conversion, so this is not the
one sanctioned flat-.data use (a raw variable-to-variable copy inside the
non-dimensionalisation boundary) and the units evaluation it bypasses matters.

For contrast, ddt.py:4567 in this same PR —

self.psi_star[k].data[...] = self.psi_star[0].data[...]

is the sanctioned raw copy and needs no change. Listing both so the line is
visibly drawn rather than guessed at.

Filed as #722 with the instances in the merged PRs, since the pattern is wider
than this branch.

Also note #721: the test_006* band this work is tested in matches no glob in
scripts/test.sh, so the integration-point suite has not been running in CI.

Underworld development team with AI support from Claude Code

lmoresi and others added 2 commits September 10, 2026 14:38
…ht up the tests

Three findings from adversarial review, all real.

**The shape guard was in the wrong place.** It lived only in __init__, but a
solver reassigns DFDt.psi_fn = flux.T on every setup (six call sites in
solvers.py), so the guard was absent from the one path that is actually
driven. A larger matrix silently TRUNCATED -- the component writer reads
psi_fn[i, j] for the slots it already has -- and a smaller one died later with
"IndexError: Index out of range: a[1]". The check now lives in the setter.

**SYM_TENSOR and TENSOR share a symbolic shape** and differ only in storage
width (3 against 4 in 2-D), so the shape check could not separate them: a full
tensor handed to a symmetric history passed and failed later with a bare
broadcast error naming neither vtype. When psi_fn is a variable it knows its
own width, so compare that.

**A non-symmetric psi_fn under SYM_TENSOR loses its lower entries**, and the
two implementations disagree about which triangle survives. On
[[1+x, 2+y], [100.0, 3+x*y]]:

    nodal SemiLagrangian        -> [[1.4499, 100.0 ], [100.0 , 3.2137]]
    IntegrationPointSemiLagrangian -> [[1.4603, 2.4667], [2.4667, 3.2148]]

Neither averages and neither warned. This class now warns and names the
entries it drops; the nodal divergence is pre-existing, so it carries a
TODO(BUG) with the measurement rather than a silent change (charter section 9).

**The tests never ran.** test_006[2-9] and test_0070 matched no batch glob in
scripts/test.sh -- ten files, including the whole integration-point suite
(0064-0067), swarm repopulation, mid-time velocity and the VE stress history.
They are not covered by the disabled test_06*py line either; that one is
0600-0699. Verified all 83 passing, then wired the range in.

The tensor tests move into tests/test_0066_integration_point_slcn.py: they are
SLCN-at-the-integration-points tests, that file is their family, and the file
they were in duplicated the number of the existing
test_0068_swarm_repopulation.py.

Full level_1 and tier_a: 1202 passed, 3 skipped, 1 xfailed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
test_007x was still unglobbed after the previous commit: tests/test_0071,
0072, 0073 (the material-index and materials suites on feature/particle-demos)
would have landed dark exactly as test_0068 did. Enumerating the gap invites
the next one; the band is now taken whole. 103 tests, all passing.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
@lmoresi
lmoresi merged commit fae3fbd into development Sep 11, 2026
2 checks passed
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.

2 participants