diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index b8edddcb..d313da6c 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/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 diff --git a/scripts/test.sh b/scripts/test.sh index f58ca36f..6dc40415 100755 --- a/scripts/test.sh +++ b/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 diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 23146df5..9f29b09b 100644 --- a/src/underworld3/systems/ddt.py +++ b/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): diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index 27f18757..1aefe935 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -1,17 +1,23 @@ """Semi-Lagrangian history at the integration points. -Two properties: the value each slot carries is the snapshot evaluated +Three properties. The value each slot carries is the snapshot evaluated exactly at the traced departure point (the floor: for a P2 field and a uniform velocity the sample is exact to round-off, for one and for two -segments), and on a rotating Gaussian the scheme is at least as accurate as -the nodal SLCN it replaces and keeps the peak better. +segments). On a rotating Gaussian the scheme is at least as accurate as the +nodal SLCN it replaces and keeps the peak better. And a history may be a +vector or a tensor, which stores one dof per INDEPENDENT component -- a +symmetric tensor in 2-D is 2x2 symbolically and three columns in storage, and +getting that packing wrong transposes a stress silently. """ +import warnings + import numpy as np import pytest import sympy import underworld3 as uw +from underworld3.systems.ddt import _storage_components pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] @@ -264,3 +270,294 @@ def test_private_trace_when_a_manager_stands_alone(): assert tr.level_valid(1) assert tr.midtime_expr() != tr.V_matrix() # 1.5 v^n - 0.5 v^{n-1} assert tr.n_velocity_evaluations == 2 + + +# --------------------------------------------------------------------------- +# Vector and tensor histories +# --------------------------------------------------------------------------- + +V_UNIFORM = np.array([1.0, 0.5]) +DT = 0.1 + + +def _velocity(): + return sympy.Matrix([[V_UNIFORM[0], V_UNIFORM[1]]]) + + +# Fields chosen to lie in the P2 space, so the sample at the departure point is +# exact to round-off, and with every component DISTINCT so that a packing or +# transposition error cannot hide. +def _scalar_field(X): + return 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + 0.5 * X[:, 0] ** 2 + + +def _vector_field(X): + return np.column_stack([ + 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + 0.5 * X[:, 0] ** 2, + -2.0 + X[:, 0] * X[:, 1] + X[:, 1] ** 2, + ]) + + +def _tensor_entries(X): + """The 2x2 symbolic entries, keyed by (i, j).""" + xx = 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + yy = -2.0 + X[:, 1] + 0.5 * X[:, 1] ** 2 + xy = 0.25 + X[:, 0] * X[:, 1] - 0.5 * X[:, 0] ** 2 + return {(0, 0): xx, (1, 1): yy, (0, 1): xy, (1, 0): xy} + + +def _pack(entries, columns): + return np.column_stack([entries[ij] for ij in columns]) + + +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "vtype,dim,expected", + [ + (uw.VarType.SCALAR, 2, [(0, 0)]), + (uw.VarType.VECTOR, 2, [(0, 0), (0, 1)]), + (uw.VarType.SYM_TENSOR, 2, [(0, 0), (1, 1), (0, 1)]), + (uw.VarType.SYM_TENSOR, 3, [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]), + ], +) +def test_the_storage_order_is_what_the_symbol_reconstructs(vtype, dim, expected): + """Pin the column -> (i, j) convention against the variable itself. + + If the storage order ever changes, this fails here rather than quietly + transposing a stress history somewhere downstream. + """ + if dim == 2: + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, qdegree=3) + else: + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), + cellSize=0.5, qdegree=3, + ) + var = uw.discretisation.MeshVariable(f"sv{vtype.value}{dim}", mesh, + vtype=vtype, degree=1) + columns = _storage_components(vtype, tuple(var.sym.shape)) + assert columns == expected + assert len(columns) == var.num_components + + # a distinct marker per column, read back through the symbol + with uw.synchronised_array_update(): + for c in range(var.num_components): + var.data[:, c] = 10.0 * (c + 1) + point = np.full((1, mesh.dim), 0.5) + got = np.asarray(uw.function.evaluate(var.sym, point)).reshape(var.sym.shape) + for c, (i, j) in enumerate(columns): + assert got[i, j] == pytest.approx(10.0 * (c + 1)), (c, i, j, got) + + +def test_a_vector_history_holds_the_departure_point_values(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=3) + U = uw.discretisation.MeshVariable("Uv", mesh, vtype=uw.VarType.VECTOR, degree=2) + with uw.synchronised_array_update(): + U.data[...] = _vector_field(np.asarray(U.coords)) + + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, U, _velocity(), vtype=uw.VarType.VECTOR, degree=2, order=2) + assert ddt.num_components == 2 + assert all(ps.is_integration_point for ps in ddt.psi_star) + assert ddt.bdf().shape == (1, 2) + + ddt.update_pre_solve(DT) + X = np.asarray(ddt.psi_star[0].coords) + foot = X - V_UNIFORM * DT + inside = (foot > 0.0).all(1) & (foot < 1.0).all(1) + assert inside.sum() > 100 + got = np.asarray(ddt.psi_star[0].data)[inside] + assert np.abs(got - _vector_field(foot[inside])).max() < 1e-12 + + # two segments back, from the older snapshot + ddt.update_post_solve(DT) + ddt.update_pre_solve(DT) + foot2 = X - V_UNIFORM * 2 * DT + inside2 = (foot2 > 0.0).all(1) & (foot2 < 1.0).all(1) + got2 = np.asarray(ddt.psi_star[1].data)[inside2] + assert np.abs(got2 - _vector_field(foot2[inside2])).max() < 1e-12 + + +def test_a_symmetric_tensor_history_transports_every_component(): + """Three independent components, all different, so a packing error or a + transposed off-diagonal cannot pass.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=3) + S = uw.discretisation.MeshVariable("St", mesh, vtype=uw.VarType.SYM_TENSOR, + degree=2) + columns = _storage_components(uw.VarType.SYM_TENSOR, (2, 2)) + with uw.synchronised_array_update(): + S.data[...] = _pack(_tensor_entries(np.asarray(S.coords)), columns) + + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, S, _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) + assert ddt.num_components == 3 + assert ddt._components == columns + assert ddt.bdf().shape == (2, 2) + + ddt.update_pre_solve(DT) + X = np.asarray(ddt.psi_star[0].coords) + foot = X - V_UNIFORM * DT + inside = (foot > 0.0).all(1) & (foot < 1.0).all(1) + assert inside.sum() > 100 + expected = _pack(_tensor_entries(foot[inside]), columns) + got = np.asarray(ddt.psi_star[0].data)[inside] + assert np.abs(got - expected).max() < 1e-12 + + # and the symbol reads back as the right 2x2 matrix, off-diagonal included + point = foot[inside][0].reshape(1, -1) + entries = _tensor_entries(point) + sym = np.asarray( + uw.function.evaluate(ddt.psi_star[0].sym, X[inside][0].reshape(1, -1)) + ).reshape(2, 2) + assert sym[0, 1] == pytest.approx(entries[(0, 1)][0], abs=1e-10) + assert sym[1, 0] == pytest.approx(sym[0, 1]) + assert sym[0, 0] == pytest.approx(entries[(0, 0)][0], abs=1e-10) + assert sym[1, 1] == pytest.approx(entries[(1, 1)][0], abs=1e-10) + assert abs(sym[0, 0] - sym[1, 1]) > 0.1 # the components are distinct + + +def test_a_scalar_history_is_unchanged(): + """The generalisation must not move the scalar answer.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=3) + T = uw.discretisation.MeshVariable("Ts", mesh, 1, degree=2) + with uw.synchronised_array_update(): + T.data[:, 0] = _scalar_field(np.asarray(T.coords)) + + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, T, _velocity(), degree=2, order=1) + assert ddt.num_components == 1 + assert ddt.bdf().shape == (1, 1) + + ddt.update_pre_solve(DT) + X = np.asarray(ddt.psi_star[0].coords) + foot = X - V_UNIFORM * DT + inside = (foot > 0.0).all(1) & (foot < 1.0).all(1) + assert np.abs( + np.asarray(ddt.psi_star[0].data)[inside, 0] - _scalar_field(foot[inside]) + ).max() < 1e-12 + + +@pytest.mark.parametrize("vtype", [uw.VarType.VECTOR, uw.VarType.SYM_TENSOR]) +def test_the_history_symbol_participates_in_expressions(vtype): + """A shaped history has to be usable, not merely storable: its symbol goes + where a mesh variable's symbol goes.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=3) + var = uw.discretisation.MeshVariable(f"Pe{vtype.value}", mesh, vtype=vtype, + degree=2) + columns = _storage_components(vtype, tuple(var.sym.shape)) + with uw.synchronised_array_update(): + var.data[...] = 2.0 + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, var, _velocity(), vtype=vtype, degree=2, order=1) + ddt.update_pre_solve(DT) + + star = ddt.psi_star[0].sym + trace = sum(star[i, i] for i in range(star.shape[0])) + assert float(uw.maths.Integral(mesh, trace).evaluate()) == pytest.approx( + 2.0 * star.shape[0], rel=1e-8) + + # the second invariant of the difference is a legitimate weak-form term + expr = (star - ddt.bdf()).T * (star - ddt.bdf()) + assert expr.shape[0] == star.shape[1] + assert len(columns) == ddt.num_components + + +def test_the_refusal_is_gone_but_the_rule_check_is_not(): + """An undersampled quadrature rule is still refused, whatever the shape: + a delta field cannot carry more values than the rule has points.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=2) + U = uw.discretisation.MeshVariable("Ur", mesh, vtype=uw.VarType.VECTOR, degree=2) + with pytest.raises(RuntimeError, match="qdegree|rule|oversample"): + uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, U, _velocity(), vtype=uw.VarType.VECTOR, degree=2, order=1) + + +def test_the_storage_map_follows_the_shape_not_the_mesh_dimension(): + """On a manifold the topological and embedding dimensions differ — a + spherical surface is dim 2, cdim 3 — and the variable sizes its storage by + the embedding one, which is what ``.sym`` is shaped by. The map therefore + reads the tensor dimension off the shape and never touches the mesh. + """ + assert len(_storage_components(uw.VarType.SYM_TENSOR, (2, 2))) == 3 + assert len(_storage_components(uw.VarType.SYM_TENSOR, (3, 3))) == 6 + assert _storage_components(uw.VarType.VECTOR, (1, 3)) == [(0, 0), (0, 1), (0, 2)] + + +def test_a_vtype_that_does_not_match_psi_fn_is_refused(): + """And refused BEFORE any variable is allocated: a mesh variable created + and then abandoned leaves its field on the DM (#1058).""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25, qdegree=3) + before = len(mesh.vars) + + with pytest.raises(ValueError, match="psi_fn has shape"): + uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, sympy.Matrix([[1.0]]), _velocity(), + vtype=uw.VarType.VECTOR, degree=2, order=1) + + with pytest.raises(ValueError, match="psi_fn has shape"): + uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, sympy.Matrix([[1.0, 2.0]]), _velocity(), + vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) + + assert len(mesh.vars) == before, "a refused history left variables behind" + + +def test_the_shape_guard_is_on_the_setter_not_only_the_constructor(): + """A solver reassigns ``DFDt.psi_fn = flux.T`` on every setup, so a guard + that lives only in ``__init__`` is absent from the one path that is + actually driven. A wrong shape there does not raise: the component writer + reads ``psi_fn[i, j]`` for the slots it already has, so a larger matrix is + silently TRUNCATED.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25, qdegree=3) + S = uw.discretisation.MeshVariable("Sg", mesh, vtype=uw.VarType.SYM_TENSOR, + degree=2) + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, S, _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) + + with pytest.raises(ValueError, match="psi_fn has shape"): + ddt.psi_fn = sympy.eye(3) # would truncate to 3 columns + with pytest.raises(ValueError, match="psi_fn has shape"): + ddt.psi_fn = sympy.Matrix([[1.0]]) # would IndexError later + + ddt.psi_fn = sympy.Matrix([[1.0, 2.0], [2.0, 3.0]]) # the right shape + assert tuple(ddt.psi_fn.shape) == (2, 2) + + +def test_a_full_tensor_is_not_accepted_as_a_symmetric_one(): + """SYM_TENSOR and TENSOR share a symbolic shape and differ in storage + width (3 against 4 in 2-D), so the shape check alone cannot separate them + and the mismatch used to surface as a bare broadcast error naming neither + vtype.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25, qdegree=3) + full = uw.discretisation.MeshVariable("Tg", mesh, vtype=uw.VarType.TENSOR, + degree=2) + with pytest.raises(ValueError, match="stores 4 components"): + uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, full, _velocity(), vtype=uw.VarType.SYM_TENSOR, + degree=2, order=1) + + +def test_an_asymmetric_psi_fn_under_sym_tensor_says_so(): + """A symmetric history stores the upper triangle, so an asymmetric psi_fn + loses its lower entries. It used to do that in silence, and the nodal + SemiLagrangian silently keeps the OTHER triangle (see the TODO(BUG) on + that class), so a user swapping one for the other would get a different + answer with no diagnostic either way.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.3, qdegree=3) + x, y = mesh.X + asymmetric = sympy.Matrix([[1 + x, 2 + y], [100.0, 3 + x * y]]) + + with pytest.warns(UserWarning, match="not symmetric"): + uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, asymmetric, _velocity(), vtype=uw.VarType.SYM_TENSOR, + degree=2, order=1) + + # a symmetric one is silent + symmetric = sympy.Matrix([[1 + x, 2 + y], [2 + y, 3 + x * y]]) + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + uw.systems.ddt.IntegrationPointSemiLagrangian( + mesh, symmetric, _velocity(), vtype=uw.VarType.SYM_TENSOR, + degree=2, order=1)