Swarm proxy as a polynomial per cell, PIC transport, and one characteristic trace per solver

docs/developer/design/eulerian-supg-transport.md

@@ -742,6 +742,49 @@ semi-Lagrangian Stokes stress history (`DFDt` on a viscoelastic Stokes solve) is
 untouched: it advects a stress that is not an unknown of the solve, which is a different
 job from the one the contract describes.
 
+### The histories share one characteristic trace
+
+The composable terms are the right design, and they exposed a cost: two histories on
+the same nodes each traced their own characteristics. A profile of the nodal
+`AdvDiffusionSLCN` step (2026-09-08) counted, per step, two managers (the value and the
+flux history) each recording its field, tracing the same nodes with the same velocity
+(two evaluations), sampling its history and caching the same velocity level: ten
+evaluate-class calls and a projection where the work is one trace, two samples and a
+copy. A hidden `simplify` in the parallel evaluator, run on every evaluation of an
+expression holding a mesh variable, had been adding 14 of 25 seconds on top.
+
+The mathematics stays with the term. What the terms now share is the evaluation behind
+`update_pre_solve`: a `CharacteristicTrace` (`systems/ddt.py`) owned by the solver
+holds the departure points per node set and segment structure for the current step, and
+the velocity levels $v^{n-1}, v^{n-2}, \dots$ cached by evaluation at the true nodes. Each
+history asks it for "the feet of my nodes through these segments" and samples its own
+field; the second history on the same nodes, and the older slots of a one-segment
+history, are served from the cache. The solver delimits the step (`begin_step`,
+`finish_step`, which records the velocity used this step); a manager used on its own
+owns a private trace and delimits its own steps, so nothing changes for standalone use.
+`share_characteristics(DuDt, DFDt)` attaches one trace to every manager that follows
+the same velocity on the same mesh; the Navier-Stokes SLCN solver shares the velocity
+levels this way even though its stress history lives on different nodes.
+
+A history is skipped only when it is a derived quantity the weak form does not read:
+the diffusive flux history under BDF or theta = 1, whose old-level weight is zero. It is
+never a state history. The flux history of a viscoelastic solve is state (the stress at
+the old time cannot be rebuilt from the present velocity gradient and rheology), and
+Crank-Nicolson is what keeps the elastic response undamped, so that history is always
+carried.
+
+Measured on the rotating Gaussian (h = 0.1, C = 0.25, one revolution, answer identical
+to every printed digit, L2 5.142e-2, peak 0.6859):
+
+| | nodal SLCN | integration-point SLCN |
+|---|---|---|
+| before | 1350 ms per step | 2250 |
+| `simplify` forwarded | 510 | 1250 |
+| shared trace, unread flux skipped | 250 | 660 |
+
+The cost ratio against SUPG in the convection comparison above was measured before
+both fixes and should be re-read with that in mind.
+
 ## What the timestep estimate means
 
 The cell-crossing time is not a stability limit for either scheme and says

docs/developer/subsystems/integration-point-variables.md

@@ -260,3 +260,142 @@ integration-point scheme loses 45 times less of it.
 The trace-back samples twelve points per cell rather than the P2 nodes, and
 the snapshot evaluation at the moving feet misses the locator cache every
 step, so the update costs about three times the nodal one.
+
+## Swarm proxy at the integration points
+
+A swarm variable normally reaches the weak form through a nodal proxy:
+the particle field is reconstructed at the proxy's nodes from the nearest
+particles and the assembler interpolates it to the integration points with
+the basis. With `proxy_location="integration_points"` the proxy is an
+integration-point variable, reconstructed from the nearest particles at
+every integration point and read there directly. That is the
+Ellipsis / Underworld PIC-LIP mapping: material properties are sampled
+from the particles around each integration point, not smoothed to the
+nodes and back.
+
+```python
+swarm = uw.swarm.Swarm(mesh)
+M = uw.swarm.SwarmVariable("M", swarm, 1, proxy_location="integration_points")
+swarm.populate(fill_param=3)
+M.data[:, 0] = ...                                   # per particle
+stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_0 * M.sym[0] + eta_1 * (1 - M.sym[0])
+```
+
+The reconstruction itself is unchanged (a linear-exact RBF over the nearest
+particles, `rbf_interpolate`); only its target moved. 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 or 1 one cell away from
+the interface (`tests/test_0067_integration_point_proxy.py`).
+`proxy_degree` and `proxy_continuous` are ignored for this proxy; the proxy
+has no gradient, so a derivative of the swarm variable's symbol is refused.
+Vector and tensor swarm variables get a multi-component proxy.
+
+`Lagrangian_Swarm(..., proxy_location="integration_points")` applies the
+same to the fully Lagrangian history: the slots carried on the particles
+are reconstructed at the integration points and the weak form reads them
+there, with no nodal history field. This is the Lagrangian option for large
+particle swarms, where the particles carry the state and the mesh only
+integrates it.
+
+## Swarm proxy as a polynomial per cell
+
+`proxy_location="cells"` is the third target. The proxy is a discontinuous
+mesh variable of `proxy_degree`, and every cell holds the least-squares
+polynomial through the particles that cell holds
+(`utilities/cell_polynomial_projection.py`). The assembler reads it at the
+integration points through the ordinary basis, so:
+
+- a polynomial particle field up to `proxy_degree` is reproduced exactly;
+- the value at the rule is a polynomial on the mesh cell, so the default
+  rule integrates it exactly and the oversampling guard of the
+  integration-point history does not apply;
+- a material step on a cell edge is exactly 0 or 1 on either side, with no
+  overshoot (the RBF reconstruction overshoots a step by up to 14%);
+- the proxy has a gradient, so Crank-Nicolson and the Adams-Moulton flux
+  of the history work;
+- each rank fits its own cells from its own particles: no neighbour search
+  across ranks, no halo particles.
+
+A cell with fewer particles than the basis size plus two takes a linear
+fit to the particles nearest its centroid, which is the RBF's
+neighbourhood; linear, because a higher-degree polynomial extrapolated
+from a distant neighbourhood is unbounded (a P2 extrapolation into the
+emptied corner cells of a rotating box reached twice the field maximum).
+That cell is consistent to first order but no longer a cell-local fit, and
+a light swarm degrades the same way the RBF proxy does. A cell with no
+particles keeps its previous proxy value: no particles is no information,
+and that holds until the swarm is repopulated. The threshold and patch
+size are `nmin` and `patch_nnn` on `CellPolynomialProjector.fit`.
+
+```python
+M = uw.swarm.SwarmVariable("M", swarm, 1, proxy_location="cells", proxy_degree=2)
+```
+
+### Why a least-squares fit and not a conservative transfer
+
+The conservative particle-to-mesh transfer solves the rule mass matrix
+against the particle moments, `M u = sum_p V_p phi(x_p) psi_p` (PETSc's
+`DMSwarmProjectFields` on a Plex does exactly this, with unit weights). It
+hands the mesh the particle sums exactly, but its nodal values carry the
+error of the particle "quadrature", which scales with the field value over
+the square root of the particle count, not with the field's variation over
+the cell. Measured on `UnstructuredSimplexBox(cellSize=0.05)` with the
+`populate()` lattice jittered by 30% of the particle spacing, L2 error at
+the integration points (`~/+Simulations/integration_point_proxy`,
+2026-09-08):
+
+| field, fill (particles per cell) | RBF at the points | conservative P1 | conservative P2 | least squares P1 | least squares P2 |
+|---|---|---|---|---|---|
+| linear, 3 | 7e-16 | 3.5 | 5.8 | 2e-10 | 6e-10 |
+| linear, 21 | 6e-16 | 0.89 | 2.6 | 9e-16 | 1e-15 |
+| Gaussian (width 0.1), 3 | 1.1e-3 | 1.6e-1 | 2.7e-1 | 2.8e-3 | 4.4e-4 |
+| Gaussian, 10 | 3.7e-4 | 6.8e-2 | 1.5e-1 | 1.7e-3 | 1.2e-4 |
+| Gaussian, 21 | 1.2e-4 | 3.8e-2 | 1.3e-1 | 1.6e-3 | 5.5e-5 |
+| Gaussian, 1 per cell on average (34% of cells empty) | 3.9e-3 | 2.5e-1 | 3.8e-1 | 8.9e-3 | 3.1e-3 |
+
+Moment-matched particle weights (chosen so constants transfer exactly)
+repair the conservative transfer's constant mode and cut the linear error
+fifty-fold, but go negative on thin cells and still trail the fit by two
+orders of magnitude on the Gaussian. Conservation and light sampling are in
+tension: what makes a light swarm usable is polynomial reproduction with a
+support that widens when the cell is thin, and the fit degree has to reach
+the mesh degree to profit from particle density (the P1 fit is limited by
+cell size and is worse than the RBF; the P2 fit beats the RBF three times
+over at ten particles per cell and at every density tested). The cell-mean
+of the fit matches the particle mean of the cell exactly; the integral of
+the fitted field differs from the particle sum by the particle-quadrature
+error, 1e-4 relative at ten particles per cell.
+
+The refresh costs about the same as the RBF path: at ten particles per
+cell on 944 cells, 8 ms (locate 6 ms, fit 2 ms) against 22 ms for the RBF
+proxy with its kd-tree rebuilt.
+
+### As the transport term of an advection-diffusion solve
+
+`Lagrangian_Swarm(..., proxy_location="cells")` composed into
+`uw.systems.AdvDiffusion` gives a particle-in-cell transport scheme: the
+particles are advected (`swarm.advection`), each history slot is fitted
+to the cells, the mesh solves the diffusion against that history, and the
+particles re-read the solution (`particle_update="pic"`, the default;
+`"flip"` adds the mesh increment instead and is kept for the MPM line of
+work, it accumulates the projection increments). The history is sampled
+at the particles the first time the swarm moves, through the swarm's
+pre-advection hook; sampled at the first solve instead, it would see the
+landed positions and lose a step.
+
+Rotating Gaussian (sigma 0.1 at radius 0.4, one revolution, h = 0.1, P2,
+C = 0.25, `~/+Simulations/integration_point_proxy/scripts/transport_gaussian.py`),
+L2 error of the mesh field against the exact solution:
+
+| Pe_h | nodal SLCN | integration-point SLCN | PIC, cells P2, 10 particles per cell | PIC, 21 per cell |
+|---|---|---|---|---|
+| infinite | 5.1e-2 (peak 0.69) | 9.1e-3 (peak 0.99) | 1.6e-2 (peak 0.86) | 5.2e-3 (peak 0.995) |
+| 400 | 4.3e-2 | 6.8e-3 | 1.3e-2 | 4.1e-3 |
+| 100 | 2.7e-2 | 3.4e-3 | 8.2e-3 | 2.4e-3 |
+| ms per step | 510 | 1250 | 140 to 170 | 180 to 250 |
+
+The particle scheme's error is the per-step re-projection (fit, then
+Galerkin projection, then read-back) and falls with particle density; at
+21 particles per cell it is below the integration-point history at a
+fifth of the cost. At C = 2 all three schemes are limited by the midpoint
+RK2 trajectory (half a radian per step, 4% phase lead), not by transport.

src/underworld3/cython/petsc_quadrature_fe.pyx

@@ -30,7 +30,7 @@ from petsc4py import PETSc
 from petsc4py.PETSc cimport FE, PetscFE, Quad, PetscQuadrature, DM, PetscDM
 from petsc4py.PETSc cimport PetscSpace, PetscDualSpace, PetscObject, MPI_Comm
 from petsc4py.PETSc cimport CHKERR as CHKERRQ
-from underworld3.cython.petsc_types cimport PetscInt, PetscReal, PetscErrorCode
+from underworld3.cython.petsc_types cimport PetscInt, PetscReal, PetscErrorCode, PetscBool
 
 import numpy as np
 
@@ -62,6 +62,8 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode PetscQuadratureDestroy(PetscQuadrature*)
 
     PetscErrorCode PetscFECreateFromSpaces(PetscSpace, PetscDualSpace, PetscQuadrature, PetscQuadrature, PetscFE*)
+    PetscErrorCode PetscFECreateVector(PetscFE, PetscInt, PetscBool, PetscBool, PetscFE*)
+    PetscErrorCode PetscFEDestroy(PetscFE*)
     PetscErrorCode PetscObjectReference(PetscObject)
     PetscErrorCode PetscObjectSetName(PetscObject, const char*)
     PetscErrorCode PetscMalloc(size_t, void**)
@@ -90,8 +92,12 @@ cdef extern from "uw_delta_space.h" nogil:
 CHKERRQ(UWDeltaSpaceRegister())
 
 
-def create_delta_fe(Quad quad, int polytope, name="quadrature_point_fe"):
-    r"""Build the scalar quadrature-point element on ``quad``.
+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,
+    component-minor, so a local vector reshapes to ``(ncells * Nq, Nc)``.
 
     Parameters
     ----------
@@ -107,7 +113,8 @@ def create_delta_fe(Quad quad, int polytope, name="quadrature_point_fe"):
     Returns
     -------
     petsc4py.PETSc.FE
-        Element of dimension ``Nq`` (points in the rule), one component,
+        Element of ``Nq * num_components`` basis functions (``Nq`` points in
+        the rule, ``num_components`` interleaved components, one by default),
         with ``quad`` as its cell quadrature and no face quadrature.
     """
     cdef PetscInt qdim = 0, qNc = 0, Nq = 0, i, d
@@ -120,7 +127,10 @@ def create_delta_fe(Quad quad, int polytope, name="quadrature_point_fe"):
     cdef PetscDualSpace Q = NULL
     cdef PetscDM refcell = NULL
     cdef PetscFE cfe = NULL
+    cdef PetscFE vfe = NULL
     cdef FE pyfe
+    if num_components < 1:
+        raise ValueError("num_components must be >= 1")
 
     CHKERRQ(PetscQuadratureGetData(quad.quad, &qdim, &qNc, &Nq, &points, &weights))
     if qNc != 1:
@@ -161,6 +171,10 @@ def create_delta_fe(Quad quad, int polytope, name="quadrature_point_fe"):
     # caller's Quad alive by taking a reference first. No face quadrature.
     CHKERRQ(PetscObjectReference(<PetscObject>quad.quad))
     CHKERRQ(PetscFECreateFromSpaces(P, Q, quad.quad, NULL, &cfe))
+    if num_components > 1:
+        CHKERRQ(PetscFECreateVector(cfe, num_components, <PetscBool>1, <PetscBool>1, &vfe))
+        CHKERRQ(PetscFEDestroy(&cfe))      # the vector element holds its own reference
+        cfe = vfe
     CHKERRQ(PetscObjectSetName(<PetscObject>cfe, name.encode()))
 
     pyfe = FE()
@@ -231,3 +245,79 @@ def cell_quadrature_points(DM dm, Quad quad):
             for d in range(cdim):
                 ov[c - cStart, q, d] = v[q * cdim + d]
     return out
+
+
+def tabulate_with_derivatives(FE fe, points):
+    r"""Tabulate ``fe``'s basis and its reference gradient at ``points``.
+
+    Returns ``(B, D)`` shaped ``(Np, Nb, Nc)`` and ``(Np, Nb, Nc, dim)``;
+    ``D`` is the gradient with respect to the reference coordinates, so a
+    physical gradient is ``invJ^T D``.
+    """
+    cdef PetscTabulation T = NULL
+    cdef PetscInt Np, Nb, Nc, cdim, p, b, c, d
+    pts = np.ascontiguousarray(points, dtype=np.float64)
+    if pts.ndim != 2:
+        raise ValueError("points must be (Np, dim)")
+    cdef double[:, ::1] pv = pts
+    Np = pts.shape[0]
+    CHKERRQ(PetscFECreateTabulation(fe.fe, 1, Np, &pv[0, 0], 1, &T))
+    Nb = T.Nb
+    Nc = T.Nc
+    cdim = T.cdim
+    B = np.empty((Np, Nb, Nc), dtype=np.float64)
+    D = np.empty((Np, Nb, Nc, cdim), dtype=np.float64)
+    cdef double[:, :, ::1] bv = B
+    cdef double[:, :, :, ::1] dv = D
+    for p in range(Np):
+        for b in range(Nb):
+            for c in range(Nc):
+                bv[p, b, c] = T.T[0][(p * Nb + b) * Nc + c]
+                for d in range(cdim):
+                    dv[p, b, c, d] = T.T[1][((p * Nb + b) * Nc + c) * cdim + d]
+    CHKERRQ(PetscTabulationDestroy(&T))
+    return B, D
+
+
+def cell_affine_maps(DM dm):
+    r"""Affine reference map of every local cell.
+
+    Returns ``(v0, invJ, detJ)`` shaped ``(ncells, cdim)``, ``(ncells, cdim,
+    cdim)`` and ``(ncells,)`` in local cell order, from
+    ``DMPlexComputeCellGeometryFEM`` with no rule (the affine map). The
+    reference coordinate of a physical point ``x`` in cell ``c`` is
+    ``invJ[c] @ (x - v0[c]) - 1`` in PETSc's ``[-1, 1]`` reference frame
+    (``v0`` is the image of the reference corner ``(-1, ..., -1)``), which is
+    the frame :func:`tabulate` expects.
+    """
+    cdef PetscInt cStart = 0, cEnd = 0, cdim = 0, c, d, e
+    cdef PetscReal *v = NULL
+    cdef PetscReal *J = NULL
+    cdef PetscReal *invJ = NULL
+    cdef PetscReal detJ = 0.0
+    CHKERRQ(DMPlexGetHeightStratum(dm.dm, 0, &cStart, &cEnd))
+    CHKERRQ(DMGetCoordinateDim(dm.dm, &cdim))
+    ncells = cEnd - cStart
+    v0 = np.empty((ncells, cdim), dtype=np.float64)
+    iJ = np.empty((ncells, cdim, cdim), dtype=np.float64)
+    dJ = np.empty((ncells,), dtype=np.float64)
+    cdef double[:, ::1] v0v = v0
+    cdef double[:, :, ::1] iJv = iJ
+    cdef double[::1] dJv = dJ
+    vbuf = np.empty(cdim, dtype=np.float64)
+    Jbuf = np.empty(cdim * cdim, dtype=np.float64)
+    iJbuf = np.empty(cdim * cdim, dtype=np.float64)
+    cdef double[::1] vv = vbuf
+    cdef double[::1] Jv = Jbuf
+    cdef double[::1] iJvb = iJbuf
+    if ncells == 0:
+        return v0, iJ, dJ
+    v = &vv[0]; J = &Jv[0]; invJ = &iJvb[0]
+    for c in range(cStart, cEnd):
+        CHKERRQ(DMPlexComputeCellGeometryFEM(dm.dm, c, NULL, v, J, invJ, &detJ))
+        dJv[c - cStart] = detJ
+        for d in range(cdim):
+            v0v[c - cStart, d] = v[d]
+            for e in range(cdim):
+                iJv[c - cStart, d, e] = invJ[d * cdim + e]
+    return v0, iJ, dJ

src/underworld3/discretisation/discretisation_mesh_variables.py

@@ -3525,7 +3525,9 @@ class _BaseIntegrationPointVariable(_BaseMeshVariable):
     ``coords`` are the physical integration points in the same order.
 
     Derivatives of the symbol are meaningless (the tabulated gradient is zero)
-    and the JIT refuses them. Scalar components only for now.
+    and the JIT refuses them. Any number of components: the element is the
+    scalar delta element wrapped as a vector element, dofs point-major and
+    component-minor within a cell.
     """
 
     is_integration_point = True
@@ -3551,15 +3553,11 @@ def _basis_key(self):
 
     def _create_petsc_fe(self, dim, prefix):
         from underworld3.cython.petsc_quadrature_fe import create_delta_fe
-        if self.num_components != 1:
-            raise NotImplementedError(
-                "IntegrationPointVariable: scalar components only for now - "
-                "use one variable per component"
-            )
         cStart, _ = self.mesh.dm.getHeightStratum(0)
-        fe = create_delta_fe(self.mesh.integration_rule, self.mesh.dm.getCellType(cStart),
-                             name=f"{prefix}integration_point_fe")
-        return fe
+        return create_delta_fe(
+            self.mesh.integration_rule, self.mesh.dm.getCellType(cStart),
+            name=f"{prefix}integration_point_fe", num_components=self.num_components,
+        )
 
     # -- geometry ---------------------------------------------------------------
 

src/underworld3/discretisation/enhanced_variables.py

@@ -66,6 +66,7 @@ class EnhancedMeshVariable(DimensionalityMixin, MathematicalMixin):
     # The storage class this wrapper delegates to; IntegrationPointVariable
     # swaps in the quadrature-point element.
     _base_variable_class = _BaseMeshVariable
+    is_integration_point = False
 
     def __new__(cls, varname, mesh, *args, **kwargs):
         """Custom __new__ to ensure proper initialization and registration."""
@@ -947,7 +948,8 @@ class IntegrationPointVariable(EnhancedMeshVariable):
     nearest-integration-point partition of each cell; ``uw.function.evaluate``
     returns that, so a query agrees with what the assembler used at the same
     point. Derivatives of the symbol are refused by the JIT (the gradient is
-    identically zero). Scalar components only for now.
+    identically zero). Vector and tensor variables are supported (one dof
+    per component per point).
 
     Examples
     --------