From f95687de0cd726f1fe58a9dd015a0fe355d35094 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 10 Sep 2026 11:22:44 -0700 Subject: [PATCH 1/4] SLCN at the integration points: vector and tensor histories 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 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 57 ++++- src/underworld3/systems/ddt.py | 103 ++++++-- ...test_0068_integration_point_slcn_tensor.py | 219 ++++++++++++++++++ 3 files changed, 354 insertions(+), 25 deletions(-) create mode 100644 tests/test_0068_integration_point_slcn_tensor.py diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index b8edddcb..bb52368e 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_0068_integration_point_slcn_tensor.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,52 @@ 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. + +Tests: `tests/test_0068_integration_point_slcn_tensor.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/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 23146df5..c22894e2 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -4225,6 +4225,28 @@ def update_post_solve( +def _storage_components(vtype, shape, dim): + """(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_0068_integration_point_slcn_tensor.py`` asserts that round trip, so + a change of convention fails there rather than silently transposing a + stress. + """ + 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: + 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 +4269,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 +4309,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 @@ -4311,10 +4339,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 +4359,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), 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" + ) # 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 @@ -4441,8 +4487,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 +4525,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 +4539,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_0068_integration_point_slcn_tensor.py b/tests/test_0068_integration_point_slcn_tensor.py new file mode 100644 index 00000000..a51f61b1 --- /dev/null +++ b/tests/test_0068_integration_point_slcn_tensor.py @@ -0,0 +1,219 @@ +"""Vector and tensor histories at the integration points. + +``IntegrationPointSemiLagrangian`` used to refuse anything but a scalar. It +now carries a vector or a tensor, which is what a Navier-Stokes momentum +history or a viscoelastic stress history needs. + +The properties checked here are the same ones the scalar case rests on — each +slot holds the snapshot evaluated exactly at the traced departure point — plus +the one the shape introduces: a field with N independent components is stored +in N columns, not in one per matrix entry, and the packing has to survive the +round trip. A symmetric tensor in 2-D is 2x2 symbolically and 3 columns in +storage; getting that wrong transposes a stress silently. +""" + +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] + +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), mesh.dim) + 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), 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), mesh.dim) + 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) From 7b25a038750c038dc5501916fde3b95d22e244e1 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 10 Sep 2026 11:33:44 -0700 Subject: [PATCH 2/4] Review (#720): the storage map follows the shape, not the mesh dimension MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- src/underworld3/systems/ddt.py | 30 ++++++++++++++-- ...test_0068_integration_point_slcn_tensor.py | 36 +++++++++++++++++-- 2 files changed, 61 insertions(+), 5 deletions(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index c22894e2..4f370fe1 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -4225,7 +4225,18 @@ def update_post_solve( -def _storage_components(vtype, shape, dim): +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 @@ -4236,12 +4247,18 @@ def _storage_components(vtype, shape, dim): ``test_0068_integration_point_slcn_tensor.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])] @@ -4327,6 +4344,15 @@ def __init__( self._psi_meshVar = None self._psi_fn = psi_fn if isinstance(psi_fn, sympy.Matrix) else sympy.Matrix([[psi_fn]]) + expected = _psi_shape_for(vtype, mesh.cdim) + if expected is not None and tuple(self._psi_fn.shape) != expected: + raise ValueError( + f"IntegrationPointSemiLagrangian: psi_fn has shape " + f"{tuple(self._psi_fn.shape)} but vtype={vtype} on a cdim=" + f"{mesh.cdim} mesh needs {expected}. Pass the vtype that " + "matches the field, or reshape psi_fn." + ) + self._init_history_tracking(order) self._check_rule_oversampling(degree) @@ -4368,7 +4394,7 @@ def __init__( ] self.num_components = int(self.psi_star[0].num_components) self._components = _storage_components( - vtype, tuple(self.psi_star[0].sym.shape), mesh.dim + vtype, tuple(self.psi_star[0].sym.shape) ) if len(self._components) != self.num_components: raise RuntimeError( diff --git a/tests/test_0068_integration_point_slcn_tensor.py b/tests/test_0068_integration_point_slcn_tensor.py index a51f61b1..51948410 100644 --- a/tests/test_0068_integration_point_slcn_tensor.py +++ b/tests/test_0068_integration_point_slcn_tensor.py @@ -82,7 +82,7 @@ def test_the_storage_order_is_what_the_symbol_reconstructs(vtype, dim, expected) ) var = uw.discretisation.MeshVariable(f"sv{vtype.value}{dim}", mesh, vtype=vtype, degree=1) - columns = _storage_components(vtype, tuple(var.sym.shape), mesh.dim) + columns = _storage_components(vtype, tuple(var.sym.shape)) assert columns == expected assert len(columns) == var.num_components @@ -131,7 +131,7 @@ def test_a_symmetric_tensor_history_transports_every_component(): 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), 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) @@ -191,7 +191,7 @@ def test_the_history_symbol_participates_in_expressions(vtype): 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), mesh.dim) + columns = _storage_components(vtype, tuple(var.sym.shape)) with uw.synchronised_array_update(): var.data[...] = 2.0 ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( @@ -217,3 +217,33 @@ def test_the_refusal_is_gone_but_the_rule_check_is_not(): 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" From f25149f9ef20fc89dce7e7803853ef83cc0b8dc7 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 10 Sep 2026 14:38:54 -0700 Subject: [PATCH 3/4] Review (#720): guard the setter, separate TENSOR from SYM_TENSOR, light 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 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 13 +- scripts/test.sh | 7 +- src/underworld3/systems/ddt.py | 82 ++++- tests/test_0066_integration_point_slcn.py | 303 +++++++++++++++++- ...test_0068_integration_point_slcn_tensor.py | 249 -------------- 5 files changed, 389 insertions(+), 265 deletions(-) delete mode 100644 tests/test_0068_integration_point_slcn_tensor.py diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index bb52368e..d313da6c 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -90,7 +90,7 @@ 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_0068_integration_point_slcn_tensor.py` pins it against +2-D — and `tests/test_0066_integration_point_slcn.py` pins it against what the variable's own `.sym` reconstructs. ## Implementation @@ -193,7 +193,16 @@ For the same history carried on **particles** rather than at the rule, use viscoelastic stress history (see below). The choice between them is where the state lives, not what shape it can take. -Tests: `tests/test_0068_integration_point_slcn_tensor.py`. +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 diff --git a/scripts/test.sh b/scripts/test.sh index f58ca36f..aca52192 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -94,7 +94,12 @@ 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_0070 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. Verified + # passing (83 tests) before wiring in. + $PYTEST tests/test_005[1-9]*py tests/test_006*py tests/test_0070*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 4f370fe1..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. @@ -4244,7 +4251,7 @@ def _storage_components(vtype, shape): 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_0068_integration_point_slcn_tensor.py`` asserts that round trip, so + ``test_0066_integration_point_slcn.py`` asserts that round trip, so a change of convention fails there rather than silently transposing a stress. @@ -4344,14 +4351,7 @@ def __init__( self._psi_meshVar = None self._psi_fn = psi_fn if isinstance(psi_fn, sympy.Matrix) else sympy.Matrix([[psi_fn]]) - expected = _psi_shape_for(vtype, mesh.cdim) - if expected is not None and tuple(self._psi_fn.shape) != expected: - raise ValueError( - f"IntegrationPointSemiLagrangian: psi_fn has shape " - f"{tuple(self._psi_fn.shape)} but vtype={vtype} on a cdim=" - f"{mesh.cdim} mesh needs {expected}. Pass the vtype that " - "matches the field, or reshape psi_fn." - ) + self._check_psi_shape(self._psi_fn) self._init_history_tracking(order) self._check_rule_oversampling(degree) @@ -4484,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 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) diff --git a/tests/test_0068_integration_point_slcn_tensor.py b/tests/test_0068_integration_point_slcn_tensor.py deleted file mode 100644 index 51948410..00000000 --- a/tests/test_0068_integration_point_slcn_tensor.py +++ /dev/null @@ -1,249 +0,0 @@ -"""Vector and tensor histories at the integration points. - -``IntegrationPointSemiLagrangian`` used to refuse anything but a scalar. It -now carries a vector or a tensor, which is what a Navier-Stokes momentum -history or a viscoelastic stress history needs. - -The properties checked here are the same ones the scalar case rests on — each -slot holds the snapshot evaluated exactly at the traced departure point — plus -the one the shape introduces: a field with N independent components is stored -in N columns, not in one per matrix entry, and the packing has to survive the -round trip. A symmetric tensor in 2-D is 2x2 symbolically and 3 columns in -storage; getting that wrong transposes a stress silently. -""" - -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] - -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" From 586224f6b0d2752bc133e0f5df3ff31098fa9d84 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 10 Sep 2026 14:52:33 -0700 Subject: [PATCH 4/4] Take the 006x-007x test band whole, so a sibling test is not dark again 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 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- scripts/test.sh | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/scripts/test.sh b/scripts/test.sh index aca52192..6dc40415 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -94,12 +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 - # test_006[2-9] and test_0070 matched NO batch glob and so never ran in - # CI: the whole integration-point suite (0064-0067), swarm repopulation, + # 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. Verified - # passing (83 tests) before wiring in. - $PYTEST tests/test_005[1-9]*py tests/test_006*py tests/test_0070*py || status=1 + # 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