SLCN at the integration points: vector and tensor histories

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

@@ -86,7 +86,12 @@ raises if they differ. Boundary integrals evaluate the field on the face
 rule and see zeros; that is correct for a history term and worth knowing for
 anything else.
 
-Scalar components only for now; use one variable per component.
+Vector and tensor variables are supported: one dof per **independent**
+component per point, so a symmetric tensor in 2-D is `2x2` symbolically and
+three columns in storage. The column order is the diagonal first, then the
+off-diagonals in row-major upper-triangular order — `(0,0), (1,1), (0,1)` in
+2-D — and `tests/test_0066_integration_point_slcn.py` pins it against
+what the variable's own `.sym` reconstructs.
 
 ## Implementation
 
@@ -132,8 +137,8 @@ adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V_fn, DuDt=DuDt, order=1
 ```
 
 The diffusive flux history (`DFDt`) keeps its nodal projection, since it
-carries derivatives. Scalar histories only; no ALE or old-frame trace-back,
-no checkpoint state yet.
+carries derivatives. Scalar, vector and tensor histories are all carried (see
+below); no ALE or old-frame trace-back, no checkpoint state yet.
 
 It is the transport manager of either advection-diffusion solver. In the
 composed `uw.systems.AdvDiffusion` (#688) it runs at `order=2` (BDF2) and
@@ -144,6 +149,61 @@ old level, which a delta field cannot supply, and the JIT guard refuses
 with a clear message; for that scheme use `AdvDiffusionSLCN`, whose
 diffusive history is a separate nodal `DFDt`.
 
+### Vector and tensor histories
+
+`vtype` selects the shape, and the slots hold one value per **independent**
+component per integration point:
+
+```python
+# a momentum history for Navier-Stokes
+DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian(
+    mesh, v, v.sym, vtype=uw.VarType.VECTOR, degree=2, order=2)
+
+# a viscoelastic stress history
+DFDt = uw.systems.ddt.IntegrationPointSemiLagrangian(
+    mesh, stress, v.sym, vtype=uw.VarType.SYM_TENSOR, degree=2, order=1)
+
+DFDt.psi_star[0].sym          # a 2x2 symbolic matrix
+DFDt.psi_star[0].data.shape   # (npoints, 3) -- three stored columns
+DFDt.bdf()                    # 2x2, as psi_fn is
+```
+
+| `vtype` (2-D) | symbolic shape | stored columns |
+|---|---|---|
+| `SCALAR` | 1×1 | 1 |
+| `VECTOR` | 1×2 | 2 |
+| `SYM_TENSOR` | 2×2 | **3** |
+
+The trace-back, the characteristic cache and the weighted sums are all
+shape-agnostic — only the fills know the shape, and they write component by
+component (`_write_components`) because a symmetric tensor's symbolic form
+repeats its off-diagonals and only the independent columns exist in storage.
+The order is 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. Get that wrong and a
+stress transposes silently, so `_storage_components` is pinned by test against
+what the variable's own `.sym` reconstructs.
+
+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 round-off (`< 1e-12`), for one segment and for two.
+
+For the same history carried on **particles** rather than at the rule, use
+`Lagrangian_Swarm`, which has been vector- and tensor-capable since the
+viscoelastic stress history (see below). The choice between them is where the
+state lives, not what shape it can take.
+
+A symmetric history stores the **upper** triangle, so an asymmetric `psi_fn`
+loses its lower entries — the manager warns rather than transporting half the
+field silently. Note the nodal `SemiLagrangian` keeps the *other* triangle in
+the same situation and does not warn; that divergence is marked with a
+`TODO(BUG)` on that class.
+
+`psi_fn` is re-checked on every assignment, not only at construction, because
+a solver reassigns it (`DFDt.psi_fn = flux.T`) on each setup.
+
+Tests: `tests/test_0066_integration_point_slcn.py`.
+
 ### The mid-point velocity is taken at the mid time
 
 The RK2 trace, `x_mid = x - dt/2 v(x)`, `x_dep = x - dt v(x_mid)`, is second

scripts/test.sh

@@ -94,7 +94,13 @@ if [ $PARALLEL_ONLY -eq 0 ]; then
   # Run simple tests (0000-0299: basic functionality, imports, simple operations)
   $PYTEST tests/test_00[0-4]*py || status=1
   $PYTEST tests/test_0050*py || status=1
-  $PYTEST tests/test_005[1-9]*py tests/test_006[0-1]*py || status=1
+  # test_006[2-9] and test_007x matched NO batch glob and so never ran in CI:
+  # 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 below either - that one is 0600-0699. The band is
+  # taken whole (00[6-7]) rather than enumerated, so a test added next to its
+  # siblings is not dark again. Verified passing (103 tests) before wiring in.
+  $PYTEST tests/test_005[1-9]*py tests/test_00[6-7]*py || status=1
   $PYTEST tests/test_01*py || status=1
   $PYTEST tests/test_02*py || status=1
 

src/underworld3/systems/ddt.py

@@ -2079,6 +2079,13 @@ def _matrix_of(V):
     return V
 
 
+# TODO(BUG): a non-symmetric psi_fn under vtype=SYM_TENSOR is silently
+# reduced here, and this class keeps the LOWER entry where
+# IntegrationPointSemiLagrangian keeps the UPPER one. Measured 2026-09-10 on
+# [[1+x, 2+y], [100.0, 3+x*y]]: nodal psi_star -> [[1.45, 100.0], [100.0, 3.21]],
+# integration-point -> [[1.46, 2.47], [2.47, 3.21]]. Neither averages and
+# neither warns. The integration-point path now warns; this one should too,
+# and the two should agree on which triangle wins.
 class SemiLagrangian(_DDtBase):
     r"""
     Semi-Lagrangian history manager.
@@ -4225,6 +4232,45 @@ def update_post_solve(
 
 
 
+def _psi_shape_for(vtype, cdim):
+    """The symbolic shape a history of ``vtype`` must have."""
+    if vtype == uw.VarType.SCALAR:
+        return (1, 1)
+    if vtype == uw.VarType.VECTOR:
+        return (1, cdim)
+    if vtype in (uw.VarType.SYM_TENSOR, uw.VarType.TENSOR):
+        return (cdim, cdim)
+    return None                     # MATRIX and friends: shape is the caller's
+
+
+def _storage_components(vtype, shape):
+    """(i, j) of the symbolic matrix that each stored column holds.
+
+    A vector or tensor field stores one dof per INDEPENDENT component, which
+    is not the same as one per matrix entry: a symmetric tensor in 2-D has a
+    2x2 symbolic form and three stored columns. The order is the one the
+    variable's own ``.sym`` reconstructs from — diagonal first, then the
+    off-diagonals in row-major upper-triangular order — and
+    ``test_0066_integration_point_slcn.py`` asserts that round trip, so
+    a change of convention fails there rather than silently transposing a
+    stress.
+
+    The tensor dimension is read off ``shape``, not off ``mesh.dim``: on a
+    manifold the two differ (a spherical surface is dim 2, cdim 3) and the
+    variable sizes its storage by the embedding dimension, which is what
+    ``.sym`` is shaped by.
+    """
+    if vtype == uw.VarType.SCALAR or shape == (1, 1):
+        return [(0, 0)]
+    if shape[0] == 1:                                    # VECTOR
+        return [(0, j) for j in range(shape[1])]
+    if vtype == uw.VarType.SYM_TENSOR:
+        dim = shape[0]
+        return ([(i, i) for i in range(dim)]
+                + [(i, j) for i in range(dim) for j in range(i + 1, dim)])
+    return [(i, j) for i in range(shape[0]) for j in range(shape[1])]
+
+
 class IntegrationPointSemiLagrangian(_DDtBase):
     r"""Semi-Lagrangian history stored at the mesh integration points.
 
@@ -4247,16 +4293,25 @@ class IntegrationPointSemiLagrangian(_DDtBase):
     ``n-k`` at the foot. Every slot carries one evaluation error rather than
     one per generation.
 
-    What is not here (yet): vector/tensor histories, units-aware velocity
-    reduction, ALE / old-frame trace-back, forcing history, checkpoint state.
-    Use :class:`SemiLagrangian` for those.
+    Scalar, vector and tensor histories are all carried: pass ``vtype``, and
+    the slots hold one value per INDEPENDENT component per integration point
+    (a symmetric tensor in 2-D is 2x2 symbolically and three columns in
+    storage). The trace-back and the weighted sums are shape-agnostic; only
+    the fills know the shape. A vector history is what a Navier-Stokes
+    momentum term needs, a symmetric tensor what a viscoelastic stress
+    history needs.
+
+    What is not here (yet): units-aware velocity reduction, ALE / old-frame
+    trace-back, forcing history, checkpoint state. Use
+    :class:`SemiLagrangian` for those, or :class:`Lagrangian_Swarm` when the
+    history should ride on particles rather than on the rule.
 
     Parameters
     ----------
     mesh, psi_fn, V_fn, degree, continuous, varsymbol, verbose, bcs, order, theta
-        As for :class:`SemiLagrangian`. ``psi_fn`` may be a scalar
-        ``MeshVariable`` (its nodal data is then copied into the snapshot
-        rather than re-evaluated) or a scalar expression.
+        As for :class:`SemiLagrangian`. ``psi_fn`` may be a ``MeshVariable``
+        (its nodal data is then copied into the snapshot rather than
+        re-evaluated) or an expression, of any ``vtype``.
     ``V_fn`` may be any expression (``-v``, ``v/2``, ``c(t) v``); the
     velocity history caches it by evaluation at each time level.
     """
@@ -4278,10 +4333,7 @@ def __init__(
         **_unsupported,
     ):
         super().__init__()
-        if vtype != VarType.SCALAR:
-            raise NotImplementedError(
-                "IntegrationPointSemiLagrangian: scalar histories only for now"
-            )
+        self.vtype = vtype
         self.monotone_mode = monotone_mode
         self.mesh = mesh
         self.bcs = list(bcs) if bcs is not None else []   # per instance, never a shared default
@@ -4299,6 +4351,8 @@ def __init__(
             self._psi_meshVar = None
             self._psi_fn = psi_fn if isinstance(psi_fn, sympy.Matrix) else sympy.Matrix([[psi_fn]])
 
+        self._check_psi_shape(self._psi_fn)
+
         self._init_history_tracking(order)
         self._check_rule_oversampling(degree)
 
@@ -4311,10 +4365,17 @@ def __init__(
             psi_units = None
         self._psi_units = psi_units
 
+        # A vector or tensor history is one dof per INDEPENDENT component per
+        # point. The trace-back and the weighted sums are shape-agnostic, so
+        # the shape changes only how much each slot stores and how many
+        # columns the fills write. Let the variable derive its own count from
+        # the vtype (a symmetric tensor in 2-D is 2x2 symbolically and three
+        # columns in storage) and read it back.
+
         # History slots at the integration points (injected, never sampled).
         self.psi_star = [
             uw.discretisation.IntegrationPointVariable(
-                f"psi_star_ip_{inst}_{k}", mesh,
+                f"psi_star_ip_{inst}_{k}", mesh, vtype=vtype,
                 varsymbol=rf"{{ {varsymbol}^{{ {'*' * (k + 1)} }} }}",
                 units=psi_units,
             )
@@ -4324,12 +4385,23 @@ def __init__(
         # (sampled at the departure points).
         self.psi_snap = [
             uw.discretisation.MeshVariable(
-                f"psi_snap_ip_{inst}_{k}", mesh, 1, degree=degree, continuous=continuous,
+                f"psi_snap_ip_{inst}_{k}", mesh, vtype=vtype,
+                degree=degree, continuous=continuous,
                 varsymbol=rf"{{ {varsymbol}^{{ (n-{k}) }} }}",
                 units=psi_units,
             )
             for k in range(order)
         ]
+        self.num_components = int(self.psi_star[0].num_components)
+        self._components = _storage_components(
+            vtype, tuple(self.psi_star[0].sym.shape)
+        )
+        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"
+            )
         # At least two velocity levels: the current interval's mid-time
         # velocity is extrapolated from v^n and v^{n-1}. Each level caches
         # V_fn evaluated at the nodes at that time, so V_fn may be any
@@ -4412,8 +4484,70 @@ def psi_fn(self):
 
     @psi_fn.setter
     def psi_fn(self, new_fn):
+        # Re-checked on every assignment, not only at construction: a solver
+        # reassigns this on each setup (``DFDt.psi_fn = flux.T``), so the
+        # constructor's guard would be bypassed on the one path that is
+        # actually driven. A wrong shape here silently truncates -- the
+        # component writer reads psi_fn[i, j] for the slots it already has.
+        new_fn = new_fn if isinstance(new_fn, sympy.Matrix) else sympy.Matrix([[new_fn]])
+        self._check_psi_shape(new_fn)
         self._psi_meshVar = None
-        self._psi_fn = new_fn if isinstance(new_fn, sympy.Matrix) else sympy.Matrix([[new_fn]])
+        self._psi_fn = new_fn
+
+    def _check_psi_shape(self, psi_fn):
+        """Refuse a psi_fn whose shape does not match this history's vtype."""
+        expected = _psi_shape_for(self.vtype, self.mesh.cdim)
+        if expected is not None and tuple(psi_fn.shape) != expected:
+            raise ValueError(
+                f"IntegrationPointSemiLagrangian: psi_fn has shape "
+                f"{tuple(psi_fn.shape)} but vtype={self.vtype} on a cdim="
+                f"{self.mesh.cdim} mesh needs {expected}. Pass the vtype that "
+                "matches the field, or reshape psi_fn."
+            )
+        # SYM_TENSOR and TENSOR have the SAME symbolic shape and different
+        # storage widths (3 and 4 in 2-D), so shape alone cannot tell them
+        # apart. When psi_fn is a variable, it knows its own width; without
+        # this a full tensor handed to a symmetric history passes the check
+        # above and dies later in the component writer with a bare broadcast
+        # error that names neither vtype.
+        supplied_var = getattr(self, "_psi_meshVar", None)
+        wanted = len(_storage_components(self.vtype, tuple(psi_fn.shape)))
+        if supplied_var is not None and expected is not None:
+            if int(supplied_var.num_components) != wanted:
+                raise ValueError(
+                    f"IntegrationPointSemiLagrangian: psi_fn stores "
+                    f"{supplied_var.num_components} components but vtype="
+                    f"{self.vtype} stores {wanted}. A full tensor and a "
+                    "symmetric tensor share a shape; pass the vtype the field "
+                    "was built with."
+                )
+        components = getattr(self, "num_components", None)
+        if components is not None and expected is not None and wanted != components:
+            raise ValueError(
+                f"IntegrationPointSemiLagrangian: psi_fn needs {wanted} stored "
+                f"components but this history has {components}; vtype="
+                f"{self.vtype} is probably not the vtype of the field."
+            )
+
+        # A symmetric history stores the upper triangle, so an asymmetric
+        # psi_fn loses its lower entries without trace. Say so rather than
+        # quietly transporting half the field the user wrote.
+        if self.vtype == uw.VarType.SYM_TENSOR:
+            dropped = [
+                (i, j) for i in range(psi_fn.shape[0])
+                for j in range(i + 1, psi_fn.shape[1])
+                if psi_fn[i, j] != psi_fn[j, i]
+            ]
+            if dropped:
+                import warnings
+
+                warnings.warn(
+                    f"IntegrationPointSemiLagrangian: psi_fn is not symmetric at "
+                    f"{dropped} but vtype=SYM_TENSOR stores only the upper "
+                    "triangle, so the lower entries are discarded (not averaged). "
+                    "Symmetrise psi_fn explicitly, or use VarType.TENSOR.",
+                    stacklevel=3,
+                )
 
     def _object_viewer(self):
         from IPython.display import Latex, Markdown, display
@@ -4441,8 +4575,22 @@ def _record_current(self):
         ):
             ps.data[...] = self._psi_meshVar.data[...]
         else:
-            vals = uw.function.evaluate(self.psi_fn[0], self._nudged_node_coords(ps))
-            ps.data[:, 0] = np.asarray(_to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1)
+            self._write_components(ps, self.psi_fn, self._nudged_node_coords(ps))
+
+    def _write_components(self, var, expr, coords, evaluate=None, **kwargs):
+        """Evaluate ``expr`` at ``coords`` and store it component by component.
+
+        One evaluation per independent component rather than one of the whole
+        matrix: a symmetric tensor's symbolic form repeats its off-diagonals,
+        and only the independent columns exist in storage.
+        """
+        if evaluate is None:
+            evaluate = uw.function.evaluate
+        for column, (i, j) in enumerate(self._components):
+            vals = evaluate(expr[i, j], coords, **kwargs)
+            var.data[:, column] = np.asarray(
+                _to_nondim_ndarray(vals, units=self._psi_units)
+            ).reshape(-1)
 
     def _segment_dt(self, j, dt):
         """Length of segment ``j`` (0 = the current step)."""
@@ -4465,11 +4613,11 @@ def _fill_slots(self, dt, evalf):
             segments.append(("first", 0, self._segment_dt(0, dt)) if k == 0
                             else ("older", k, self._segment_dt(k, dt)))
             X = trace.departure_points(key, X0, tuple(segments), evalf=evalf)
-            vals = uw.function.global_evaluate(
-                self.psi_snap[k].sym[0], X, evalf=evalf, monotone=self.monotone_mode
+            self._write_components(
+                self.psi_star[k], self.psi_snap[k].sym, X,
+                evaluate=uw.function.global_evaluate,
+                evalf=evalf, monotone=self.monotone_mode,
             )
-            self.psi_star[k].data[:, 0] = np.asarray(
-                _to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1)
 
     def initialise_history(self):
         """Start every snapshot and slot from the current field, so
@@ -4479,10 +4627,9 @@ def initialise_history(self):
             self.psi_snap[k].data[...] = self.psi_snap[0].data[...]
         self.characteristics.initialise_levels(self._n_v)
         X = np.asarray(self.psi_star[0].coords_nd)
-        vals = uw.function.evaluate(self.psi_snap[0].sym[0], X)
-        vals = np.asarray(_to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1)
-        for k in range(self.order):
-            self.psi_star[k].data[:, 0] = vals
+        self._write_components(self.psi_star[0], self.psi_snap[0].sym, X)
+        for k in range(1, self.order):
+            self.psi_star[k].data[...] = self.psi_star[0].data[...]
         self._history_initialised = True
 
     def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored):