Swarm proxy as a polynomial per cell, PIC transport, and one characteristic trace per solver - #707
Conversation
The scalar delta element is wrapped with PetscFECreateVector (interleaved basis and components), so a cell's dofs are point-major, component-minor and the local vector reshapes to (ncells * Nq, Nc). Tests: the two- component element tabulates to the point-major identity on triangle, tet, quad and hex and to zero off-rule; a two-component variable has the (ncells*Nq, 2) layout, evaluates exactly at its own points, and each component is reproduced by a P2 vector projection to 1e-9. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
…ntegration points directly The proxy is an IntegrationPointVariable reconstructed from the nearest particles at every integration point of the mesh rule and read there by the assembler with no basis interpolation: the Ellipsis / Underworld PIC-LIP mapping. The reconstruction is unchanged (linear-exact RBF at the target's own coordinates); only the target moved. Lagrangian_Swarm takes the same option for its history slots, so the fully Lagrangian history skips the nodal proxy. EnhancedMeshVariable exposes is_integration_point. tests/test_0067: a particle-carried material step is reproduced at the integration points with less than half the L2 error of the nodal proxy and is exactly 0/1 away from the interface; a derivative of the proxy is refused; the Lagrangian history slot's proxy reproduces a linear field at the integration points to 1e-8; a symmetric-tensor swarm variable gets a three-component proxy. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
…ntegration-point proxy Viscosity step 1e3 carried by particles with the interface on mesh edges: the exact piecewise-linear velocity lies in the P2 space, so the only error is the proxy's representation of the step. Integration-point proxy: 4e-7 (solver tolerance); nodal P1 proxy: 0.2, because a node on the interface averages both materials. Study and the cell-cutting sweep in ~/+Simulations/integration_point_proxy. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
A third proxy target: a discontinuous mesh variable of proxy_degree whose every cell holds the least-squares polynomial through the particles it holds (utilities/cell_polynomial_projection.py). Exact for polynomial particle fields up to the degree, integrated exactly by the default rule (no oversampling guard), sharp at cell edges, with a gradient, and rank-local. A thin cell (fewer than basis + 2 particles) is fitted to the particles nearest its centroid, so light swarms degrade the way the RBF proxy does. Lagrangian_Swarm accepts the option for its history slots. The conservative particle-to-mesh transfer (rule mass against particle moments, PETSc's DMSwarmProjectFields) was measured and rejected: with irregular particles its nodal values carry the particle-quadrature error, O(|psi| / sqrt N), 25% of a linear field at 21 particles per cell, while the least-squares P2 fit is exact and beats the RBF three-fold at ten particles per cell. UW3 swarms are DMSWARM_BASIC (no cell DM), so the fit is built from the owning-cell locator and PetscFE tabulation instead; cell_affine_maps and tabulate_with_derivatives are added to cython/petsc_quadrature_fe.pyx for it. Tests: polynomial reproduction through the weak form, gradient, light swarm on the patch fit, edge-aligned step exact, vector variable and Lagrangian_Swarm slots. Docs section with the measurements. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
… thin-cell fit; Lagrangian_Swarm FLIP option global_evaluate_nd called evaluate_nd without its `simplify` argument, so the local evaluator's default (True) ran sympy.simplify on every call for any expression holding a mesh variable. A semi-Lagrangian step whose mid-time velocity carries the cached velocity level paid 14 of 25 s in simplify with a tanh velocity; forwarding the flag (the public default is False) cuts the nodal SLCN step from 1350 to 510 ms and the integration-point SLCN from 2250 to 1250 ms on the rotating Gaussian, with identical answers. Cells proxy: a thin cell now takes a LINEAR fit to the particles nearest its centroid (a P2 extrapolation into the emptied corner cells of a rotating box reached twice the field maximum), and after the first fit a cell with no particles keeps its previous value. Lagrangian_Swarm gains particle_update="pic"|"flip" with residual_retention: FLIP adds the mesh increment (solution minus the proxy the mesh saw) to the particle value. Measured on the rotating Gaussian it accumulates the projection increments (9% overshoot on pure advection) and diverges next to held cells; kept as an option for the MPM line of work, PIC remains the default. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
The history initialised itself on the first update_pre_solve, which comes after the user's swarm.advection(): the first sample saw the landed positions and the first step transported nothing, a one-step lag of the whole field (0.05 of displacement on the rotating Gaussian, L2 3e-2 against 5e-3 once fixed). Swarm.advection() now runs pre-advection hooks before any particle moves and Lagrangian_Swarm registers its first sampling there (weak reference; the first-solve fallback remains). Docs: the composed PIC scheme's results on the rotating Gaussian against nodal and integration-point SLCN. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
CharacteristicTrace (systems/ddt.py) holds, for one advecting velocity on
one mesh, the departure points of the current step per node set and
segment chain, and the velocity levels v^{n-1}, v^{n-2}, ... cached by
evaluation at the true nodes. The nodal SemiLagrangian and the
IntegrationPointSemiLagrangian managers trace through it; a solver
attaches one trace to every manager that follows the same velocity
(share_characteristics) and delimits the steps; a manager used on its
own keeps a private trace. The mathematics of each history term is
unchanged: its symbols, stencils and view are its own; only the
evaluation behind update_pre_solve is shared.
AdvDiffusionSLCN shares the trace between its value and flux histories
and skips the flux history when the weak form does not read it (BDF
orders, theta = 1: old-level weight zero). A state history (the
viscoelastic stress) is never skipped. The Navier-Stokes SLCN solver
shares the velocity levels the same way.
Rotating Gaussian, h = 0.1, C = 0.25: answer identical to every printed
digit (L2 5.142e-2, peak 0.6859); nodal SLCN 510 -> 250 ms per step,
integration-point SLCN 1250 -> 660. Tests count the evaluations: two per
step for two histories, one cache hit, none when the flux is unread.
Underworld development team with AI support from Claude Code
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
The Cython evaluators (global_evaluate_nd, evaluate_nd, petsc_interpolate, rbf_evaluate) defaulted simplify=True while the public evaluate and global_evaluate default to False; any internal call that dropped the flag ran sympy.simplify on every evaluation. Off by default everywhere now; simplify=True remains available on request. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
…ation-point-proxy
There was a problem hiding this comment.
🟡 Changes recommended
The new FLIP update path can compute residuals against a stale proxy (missing proxy refresh in _proxy_values_at_particles), and there is at least one mutable-default argument that should be corrected before approval.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
Pull request overview
This PR extends Underworld3’s transport infrastructure to support higher-fidelity and more efficient semi-Lagrangian / PIC-style advection workflows, without changing default solver mathematics. It introduces (1) shared characteristic tracing per solver/velocity, (2) new swarm proxy targets at integration points and as per-cell polynomials, and (3) particle transport updates (PIC by default, optional FLIP-style update) alongside evaluator performance fixes.
Changes:
- Add
CharacteristicTraceand sharing (share_characteristics) so multiple histories following the same velocity reuse departure points and cached velocity levels, and allow skipping unread flux histories. - Add new swarm proxy locations:
proxy_location="integration_points"(direct proxy at quadrature points) andproxy_location="cells"(rank-local per-cell least-squares polynomial projection viaCellPolynomialProjector). - Adjust evaluator defaults/forwarding to avoid unintended
sympy.simplifycosts; expand integration-point variable support to multi-component layouts; add/extend tests and docs.
File summaries
| File | Description |
|---|---|
| tests/test_0067_integration_point_proxy.py | New tests covering integration-point proxies, cells polynomial proxies, and PIC/FLIP transport behaviors. |
| tests/test_0066_integration_point_slcn.py | Tests for shared characteristic tracing and standalone/private trace semantics. |
| tests/test_0065_integration_point_variable.py | Updates integration-point variable tests for multi-component support and projection/evaluate behavior. |
| tests/test_0064_quadrature_point_fe.py | Extends delta FE tests to vector layout and derivative tabulation helpers. |
| src/underworld3/utilities/cell_polynomial_projection.py | New per-cell least-squares polynomial projection used by proxy_location="cells". |
| src/underworld3/systems/solvers.py | Shares a characteristic trace across value/flux histories; skips flux-history update when weak form won’t read it. |
| src/underworld3/systems/ddt.py | Adds CharacteristicTrace, shared-trace wiring, and Lagrangian swarm transport extensions (PIC/FLIP, proxy_location support). |
| src/underworld3/swarm.py | Adds proxy_location to SwarmVariable, implements integration-point and cell-polynomial proxy updates, and adds pre-advection hooks. |
| src/underworld3/function/_function.pyx | Changes defaults to simplify=False and forwards simplify through global→local evaluator paths. |
| src/underworld3/discretisation/enhanced_variables.py | Declares/clarifies integration-point variable support for vector/tensor components. |
| src/underworld3/discretisation/discretisation_mesh_variables.py | Enables multi-component integration-point variables by creating vector-wrapped delta elements. |
| src/underworld3/cython/petsc_quadrature_fe.pyx | Adds num_components support to create_delta_fe() and derivative tabulation / affine-map helpers. |
| docs/developer/subsystems/integration-point-variables.md | Documents new swarm proxy locations and PIC/FLIP transport composition. |
| docs/developer/design/eulerian-supg-transport.md | Design note describing shared characteristic trace motivation and measured speedups. |
Review details
Suppressed comments (1)
src/underworld3/systems/ddt.py:4161
- Same issue as above: this _proxy_values_at_particles() path should refresh the swarm proxy before evaluating the proxy mesh variable; otherwise FLIP can use stale proxy values if slot.sym wasn't accessed earlier in the step.
def _proxy_values_at_particles(self, slot, coords, evalf):
"""The slot's proxy evaluated at the particles, shaped like ``slot.data``.
A ``"cells"`` proxy is read through its own fitted polynomials (exact,
no locator round trip); any other proxy through ``evaluate`` of the
proxy mesh variable's symbol.
"""
- Files reviewed: 14/14 changed files
- Comments generated: 3
- Review effort level: Lite
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| def _proxy_values_at_particles(self, slot, coords, evalf): | ||
| """The slot's proxy evaluated at the particles, shaped like ``slot.data``. | ||
|
|
||
| A ``"cells"`` proxy is read through its own fitted polynomials (exact, | ||
| no locator round trip); any other proxy through ``evaluate`` of the | ||
| proxy mesh variable's symbol. | ||
| """ | ||
| projector = getattr(slot, "_cell_projector", None) | ||
| if projector is not None and getattr(slot, "_proxy_location", None) == "cells": | ||
| slot._update_proxy_if_stale() | ||
| vals = projector.interpolate(np.asarray(slot._meshVar.data), coords) | ||
| return np.nan_to_num(vals) | ||
| mv = slot._meshVar | ||
| out = np.empty((coords.shape[0], slot.data.shape[1])) | ||
| for i in range(slot.shape[0]): | ||
| for j in range(slot.shape[1]): | ||
| ij = slot._data_layout(i, j) | ||
| out[:, ij] = np.asarray( | ||
| uw.function.evaluate(mv.sym[i, j], coords, evalf=evalf) | ||
| ).reshape(-1) | ||
| return out |
There was a problem hiding this comment.
Fixed in the follow-up commit: the proxy is refreshed first on every path of _proxy_values_at_particles, the misplaced copy of the helper in the nodal-swarm Lagrangian class is removed, bcs defaults to None, and the create_delta_fe docstring now describes the multi-component element.
| varsymbol: Optional[str] = r"u", | ||
| verbose: Optional[bool] = False, | ||
| bcs=[], | ||
| order=1, | ||
| smoothing=0.0, |
There was a problem hiding this comment.
Fixed in the follow-up commit: the proxy is refreshed first on every path of _proxy_values_at_particles, the misplaced copy of the helper in the nodal-swarm Lagrangian class is removed, bcs defaults to None, and the create_delta_fe docstring now describes the multi-component element.
| def create_delta_fe(Quad quad, int polytope, name="quadrature_point_fe", int num_components=1): | ||
| r"""Build the quadrature-point element on ``quad``. | ||
|
|
||
| ``num_components > 1`` wraps the scalar element with ``PetscFECreateVector`` | ||
| (interleaved basis and components): the dofs of a cell are point-major, |
There was a problem hiding this comment.
Fixed in the follow-up commit: the proxy is refreshed first on every path of _proxy_values_at_particles, the misplaced copy of the helper in the nodal-swarm Lagrangian class is removed, bcs defaults to None, and the create_delta_fe docstring now describes the multi-component element.
…splaced helper copy, bcs=None, delta-element docstring Copilot's three findings: _proxy_values_at_particles read the proxy mesh variable without refreshing a stale proxy on the non-cells path (now refreshed first on every path); the same helper had been inserted into the nodal-swarm Lagrangian class as well, where nothing uses it (removed); Lagrangian_Swarm's bcs default was a mutable list; create_delta_fe's docstring still described a one-component element. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G
Stacked on #703 (merged). Three things, none of which changes a solver's default mathematics.
Efficiency, for every semi-Lagrangian solver
CharacteristicTrace(systems/ddt.py): one trace per solver per advecting velocity. It holds the departure points of the step per node set and segment chain, and the velocity levels cached by evaluation at the true nodes. The nodal and the integration-point managers trace through it; a solver attaches one trace to every history that follows the same velocity (share_characteristics) and delimits the step; a standalone manager keeps a private trace. Each history term keeps its own symbols, stencils and view.AdvDiffusionSLCNskips its flux history when the weak form does not read it (BDF orders, theta = 1). A state history (the viscoelastic stress) is never skipped.simplifyflag on the way to the local evaluator, so any expression holding a mesh variable was sympy-simplified at every evaluation. Forwarded, and the internal defaults now match the public API (off).Rotating Gaussian, h = 0.1, C = 0.25, answer identical to every printed digit: nodal SLCN 1350 to 250 ms per step, integration-point SLCN 2250 to 660. The SUPG-versus-SLCN cost ratio measured earlier carried both traps.
Swarm proxy as a polynomial per cell
SwarmVariable(proxy_location="cells", proxy_degree=2): a discontinuous proxy whose every cell holds the least-squares polynomial through the particles it holds (utilities/cell_polynomial_projection.py). Exact for polynomial particle fields, integrated exactly by the default rule, sharp at cell edges, with a gradient, rank-local. Thin cells take a linear patch fit, empty cells hold their previous value. The conservative particle-to-mesh transfer was measured and rejected: with irregular particles its error scales with the field value over the square root of the particle count (25% of a linear field at 21 particles per cell).Particle transport
Lagrangian_Swarmwith the cells proxy composed intoAdvDiffusionis a particle-in-cell transport scheme. Two traps fixed on the way: the swarm history sampled itself at the first solve, after the first advection, losing a step (now sampled through a swarm pre-advection hook); and a P2 extrapolation into emptied cells was unbounded.particle_update="flip"is added as an option for the MPM line of work; PIC stays the default. On the rotating Gaussian PIC at 21 particles per cell reaches L2 5.2e-3 against the integration-point history's 9.1e-3 at a fifth of the cost.Docs:
docs/developer/subsystems/integration-point-variables.md,docs/developer/design/eulerian-supg-transport.md. Tests: 0064, 0066, 0067. Study:~/+Simulations/integration_point_proxy.Underworld development team with AI support from Claude Code
🤖 Generated with Claude Code
https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G