From 0b5436428533dc2ba54f08722b360326a466b7b8 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 10 Sep 2026 18:13:04 -0700 Subject: [PATCH] Symmetric-tensor and deferred writes through .array were silently lost Two defects on the write path the style charter makes mandatory for new code (section 7: "new code uses the array property"). Both are silent, both are reachable from the interface docs/developer/subsystems/data-access.md recommends, and both predate this branch. **1. A symmetric tensor's off-diagonal write vanished.** (i, j) and (j, i) share ONE stored column, and the pack-back loop visits every (i, j), so the pair's second visit overwrote the first: var.array[:, 0, 1] = 30.0 -> stored 0.0 silently discarded var.array[:, 1, 0] = 30.0 -> stored 30.0 worked Lower beat upper only because it came last in the loop. In 3-D all three upper off-diagonals vanished. Reading .array back showed the write had not happened, so a stress or strain-rate history assembled component by component lost every shear term. The half that changed is now mirrored onto the half that did not, and setting the two corners to DIFFERENT values in one assignment is refused rather than resolved by loop order. The swarm variant failed loudly instead of silently -- its pack was a flat reshape that ignored symmetric storage, so (N, 2, 2) became (N, 4) and the assignment could not broadcast into (N, 3). A symmetric tensor on a swarm could not be written through .array at all. Now routed through _data_layout, like the mesh path. **2. Deferred writes on a swarm kept only the last.** The swarm view read its current values straight from the PETSc field. synchronised_array_update() defers the pack, so inside that context the field still held the pre-context values and every write started from them: with uw.synchronised_array_update(): velocity.array[:, 0, 0] = 1.0 velocity.array[:, 0, 1] = 2.0 # -> [0., 2.] the first write was overwritten Mesh variables were immune: their view reads and writes the canonical array. The swarm view now does too. **Why the suite missed both.** It was not thin -- 219 of 399 test files touch .array and 167 writes use the three-index form. But 148 of those are [:, N, N] with literal integers; across the whole suite there are exactly TWO off-diagonal writes and both pick [:, 1, 0], the corner that survived. The one file that writes all four corners is on a VarType.TENSOR, which stores them independently and cannot expose it. Likewise eight blocks write one variable several times inside a synchronised context -- every one on a MESH variable. Both defects sat in the single cell of the shape x carrier x index matrix that nobody had filled. The new tests are parametrised over that product rather than written one case at a time: 6 cases from 2 lines, 3 of which were failing. Full level_1 and tier_a: 1187 passed. New file 19/19 at np=1 and np=2. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../discretisation_mesh_variables.py | 42 +++ src/underworld3/swarm.py | 93 ++++++- ...test_0509_symmetric_tensor_array_writes.py | 244 ++++++++++++++++++ 3 files changed, 370 insertions(+), 9 deletions(-) create mode 100644 tests/test_0509_symmetric_tensor_array_writes.py diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index cbea9a68..a2622df2 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -2577,6 +2577,15 @@ def __setitem__(self, key, value): # Step 3: Assign the (now non-dimensional) value modified_data[key] = value + # A symmetric tensor is (dim, dim) here and stores only its + # independent components, so both entries of an off-diagonal + # pair map to ONE column. The pack loop below writes every + # (i, j), which means the pair's second visit overwrites the + # first: setting array[:, 0, 1] alone was silently discarded + # (the stale [1, 0] won), while array[:, 1, 0] alone worked. + # Mirror whichever half the caller actually changed. + self._mirror_symmetric_pairs(unpacked, modified_data) + # Route the write through the canonical array (see # SimpleMeshArrayView.__setitem__: a direct pack is a # per-write collective). _data_layout maps the structured @@ -2592,6 +2601,39 @@ def __setitem__(self, key, value): flat_data[:, self.parent._data_layout(i, j)] = modified_data[:, i, j] self.parent.data[...] = flat_data + def _mirror_symmetric_pairs(self, before, after): + """Carry an off-diagonal write across to its mirror entry. + + Only for a symmetric variable, whose (i, j) and (j, i) share a + stored column. Writing one half and leaving the other stale is + how the write got lost, so the half that changed is copied onto + the half that did not. Changing BOTH halves to different values + asks for something the storage cannot hold, and is refused + rather than resolved by the loop order. + """ + import underworld3 as uw + + if self.parent.vtype != uw.VarType.SYM_TENSOR: + return + rows, cols = self.parent.shape + for i in range(rows): + for j in range(i + 1, cols): + upper_moved = not numpy.array_equal(after[:, i, j], before[:, i, j]) + lower_moved = not numpy.array_equal(after[:, j, i], before[:, j, i]) + if upper_moved and not lower_moved: + after[:, j, i] = after[:, i, j] + elif lower_moved and not upper_moved: + after[:, i, j] = after[:, j, i] + elif upper_moved and lower_moved and not numpy.array_equal( + after[:, i, j], after[:, j, i] + ): + raise ValueError( + f"'{self.parent.name}' is a symmetric tensor: " + f"components [{i}, {j}] and [{j}, {i}] share one " + "stored value and cannot be set to different " + "values in a single assignment." + ) + @property def shape(self): return self._get_array_data().shape diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 7a0ff74c..76fefc2e 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -879,9 +879,22 @@ def __setitem__(self, key, value): f" 3. For non-dimensional values, use: {self.parent.name}.data[...] = value\n" ) - # Get current NON-DIMENSIONAL array data from PETSc - # Note: We use unpack directly here, not _get_array_data() which dimensionalizes - array_data = self.parent.unpack_uw_data_from_petsc(squeeze=False) + # Current NON-DIMENSIONAL values, read from the CANONICAL + # array rather than from the PETSc field. + # + # Reading PETSc here made every write inside + # uw.synchronised_array_update() overwrite the one before it: + # that context defers the PETSc pack, so the field still held + # the pre-context values and each __setitem__ started from + # them. Two component writes to one swarm variable kept only + # the last, silently -- and writing a vector or tensor one + # component at a time inside that context is exactly what the + # data-access guide recommends. Mesh variables were never + # affected: their view writes through the canonical array. + # (Not _get_array_data(), which dimensionalises.) + array_data = self.parent._unpack_data_to_array_format( + np.asarray(self.parent.data) + ) # Create a copy to modify (avoid modifying view directly) modified_data = array_data.copy() @@ -926,6 +939,10 @@ def __setitem__(self, key, value): # Update the specific elements modified_data[key] = value + # A symmetric tensor stores one value per off-diagonal PAIR, + # so a write to one half has to be carried to the other; see + # SwarmVariable._mirror_symmetric_pairs. + self.parent._mirror_symmetric_pairs(array_data, modified_data) # Pack back to canonical data format packed_data = self.parent._pack_array_to_data_format(modified_data) self.parent.data[:] = packed_data @@ -1009,19 +1026,77 @@ def delay_callback(self, description="array operation"): return TensorSwarmArrayView(self) + def _mirror_symmetric_pairs(self, before, after): + """Carry an off-diagonal write across to its mirror entry. + + Only for a symmetric variable, whose ``(i, j)`` and ``(j, i)`` share + one stored value. Writing one half and leaving the other stale loses + the write; changing both to different values asks for something the + storage cannot hold and is refused. + """ + if self.vtype != uw.VarType.SYM_TENSOR or after.ndim < 3: + return + rows, cols = after.shape[1], after.shape[2] + for i in range(rows): + for j in range(i + 1, cols): + upper_moved = not np.array_equal(after[:, i, j], before[:, i, j]) + lower_moved = not np.array_equal(after[:, j, i], before[:, j, i]) + if upper_moved and not lower_moved: + after[:, j, i] = after[:, i, j] + elif lower_moved and not upper_moved: + after[:, i, j] = after[:, j, i] + elif upper_moved and lower_moved and not np.array_equal( + after[:, i, j], after[:, j, i] + ): + raise ValueError( + f"'{self.clean_name}' is a symmetric tensor: components " + f"[{i}, {j}] and [{j}, {i}] share one stored value and " + "cannot be set to different values in a single assignment." + ) + + def _unpack_data_to_array_format(self, flat_data): + """Canonical ``(N, components)`` -> array ``(N, a, b)``. + + The inverse of :meth:`_pack_array_to_data_format`, and the read half + of a read-modify-write on the canonical array. Follows + ``_data_layout``, so a symmetric tensor's shared off-diagonal column + appears at both ``(i, j)`` and ``(j, i)``. + """ + flat_data = np.asarray(flat_data) + shape = self.shape + unpacked = np.empty((flat_data.shape[0], *shape), dtype=flat_data.dtype) + for i in range(shape[0]): + for j in range(shape[1]): + unpacked[:, i, j] = flat_data[:, self._data_layout(i, j)] + return unpacked + def _pack_array_to_data_format(self, array_data): - """Convert array format (N,a,b) back to canonical data format (N,components)""" - # Use existing pack logic but return numpy array instead of writing to PETSc - # This is a pure conversion method - no PETSc access + """Convert array format (N,a,b) back to canonical data format (N,components) + + A flat reshape is wrong for a symmetric tensor: the ``(N, dim, dim)`` + view has ``dim*dim`` entries and storage holds only the independent + ones, so reshaping produced ``(N, 4)`` for a ``(N, 3)`` variable and + the assignment failed to broadcast. ``_data_layout`` is the mapping + that ``pack_uw_data_to_petsc`` uses; follow it. + """ # Empty-partition guard: an N=0 array has total size 0, so numpy cannot # infer the -1 component dimension ("cannot reshape array of size 0 into # shape (0,newaxis)"). This bites a rank that owns no local particles # during a parallel read_timestep. Compute the component count from the # trailing dims explicitly. if array_data.size == 0: - ncomp = int(np.prod(array_data.shape[1:])) if array_data.ndim > 1 else 1 - return array_data.reshape(array_data.shape[0], ncomp) - return array_data.reshape(array_data.shape[0], -1) + return array_data.reshape(array_data.shape[0], self.num_components) + + if array_data.ndim < 3 or array_data.shape[1] * array_data.shape[2] == self.num_components: + return array_data.reshape(array_data.shape[0], -1) + + packed = np.empty( + (array_data.shape[0], self.num_components), dtype=array_data.dtype + ) + for i in range(array_data.shape[1]): + for j in range(array_data.shape[2]): + packed[:, self._data_layout(i, j)] = array_data[:, i, j] + return packed # Legacy methods preserved for backward compatibility (now do nothing) def use_legacy_array(self): diff --git a/tests/test_0509_symmetric_tensor_array_writes.py b/tests/test_0509_symmetric_tensor_array_writes.py new file mode 100644 index 00000000..cc1d4bbd --- /dev/null +++ b/tests/test_0509_symmetric_tensor_array_writes.py @@ -0,0 +1,244 @@ +"""Writing one off-diagonal of a symmetric tensor through ``.array``. + +A symmetric tensor is ``(dim, dim)`` in the ``.array`` view and stores only +its independent components, so ``[i, j]`` and ``[j, i]`` share ONE stored +value. The pack-back loop visits every ``(i, j)`` and writes that shared +column twice, so the pair's second visit overwrote the first: + + var.array[:, 0, 1] = 30.0 -> stored 0.0 silently discarded + var.array[:, 1, 0] = 30.0 -> stored 30.0 worked + +Lower beat upper purely because it came last in the loop. In 3-D all three +upper off-diagonals vanished. Nothing raised, and reading ``.array`` back +showed the write had not happened -- so a stress or strain-rate history +assembled component by component through the documented interface lost every +shear term. + +The charter (section 7) makes ``.array`` the interface for new code, which is +what makes this worth a guard rather than a note. +""" + +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +def _box(dim=2): + if dim == 2: + return uw.meshing.UnstructuredSimplexBox(cellSize=0.4, qdegree=3) + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), + cellSize=0.5, qdegree=3, + ) + + +@pytest.mark.parametrize("corner", [(0, 1), (1, 0)]) +def test_either_off_diagonal_write_lands(corner): + """Both halves of the pair must work, and both must mirror.""" + mesh = _box() + i, j = corner + var = uw.discretisation.MeshVariable( + f"sym{i}{j}", mesh, vtype=uw.VarType.SYM_TENSOR, degree=1) + + with uw.synchronised_array_update(): + var.array[:, i, j] = 30.0 + + read_back = np.asarray(var.array) + assert np.allclose(read_back[:, i, j], 30.0), read_back[0] + assert np.allclose(read_back[:, j, i], 30.0), "the write did not mirror" + assert np.allclose(np.asarray(var.data)[:, 2], 30.0), "it never reached storage" + + +def test_an_off_diagonal_write_lands_without_the_context_manager(): + mesh = _box() + var = uw.discretisation.MeshVariable( + "symplain", mesh, vtype=uw.VarType.SYM_TENSOR, degree=1) + var.array[:, 0, 1] = 30.0 + assert np.allclose(np.asarray(var.array)[:, 0, 1], 30.0) + + +def test_every_upper_off_diagonal_lands_in_three_dimensions(): + """2-D has one pair and 3-D has three; all of them were lost.""" + mesh = _box(dim=3) + var = uw.discretisation.MeshVariable( + "sym3d", mesh, vtype=uw.VarType.SYM_TENSOR, degree=1) + + with uw.synchronised_array_update(): + var.array[:, 0, 0] = 1.0 + var.array[:, 1, 1] = 2.0 + var.array[:, 2, 2] = 3.0 + var.array[:, 0, 1] = 4.0 + var.array[:, 0, 2] = 5.0 + var.array[:, 1, 2] = 6.0 + + read_back = np.asarray(var.array)[0] + expected = np.array([[1.0, 4.0, 5.0], [4.0, 2.0, 6.0], [5.0, 6.0, 3.0]]) + assert np.allclose(read_back, expected), read_back + + +def test_asking_for_an_asymmetric_value_is_refused(): + """[i, j] and [j, i] are one stored number. Setting them to different + values in one assignment cannot be honoured, and used to be resolved + silently by whichever the loop wrote last.""" + mesh = _box() + var = uw.discretisation.MeshVariable( + "symbad", mesh, vtype=uw.VarType.SYM_TENSOR, degree=1) + values = np.zeros((np.asarray(var.array).shape[0], 2, 2)) + values[:, 0, 1] = 7.0 + values[:, 1, 0] = 9.0 + + with pytest.raises(ValueError, match="symmetric tensor"): + with uw.synchronised_array_update(): + var.array[...] = values + + +def test_a_full_tensor_still_holds_both_corners(): + """The mirroring must not touch VarType.TENSOR, which stores all four.""" + mesh = _box() + var = uw.discretisation.MeshVariable( + "fullt", mesh, vtype=uw.VarType.TENSOR, degree=1) + + with uw.synchronised_array_update(): + var.array[:, 0, 1] = 7.0 + var.array[:, 1, 0] = 9.0 + + read_back = np.asarray(var.array)[0] + assert read_back[0, 1] == pytest.approx(7.0) + assert read_back[1, 0] == pytest.approx(9.0) + + +def test_a_whole_symmetric_assignment_is_unchanged(): + """The path that always worked must keep working.""" + mesh = _box() + var = uw.discretisation.MeshVariable( + "symwhole", mesh, vtype=uw.VarType.SYM_TENSOR, degree=1) + values = np.zeros((np.asarray(var.array).shape[0], 2, 2)) + values[:, 0, 0], values[:, 1, 1] = 1.0, 2.0 + values[:, 0, 1] = values[:, 1, 0] = 3.0 + + with uw.synchronised_array_update(): + var.array[...] = values + + assert np.allclose(np.asarray(var.data)[0], [1.0, 2.0, 3.0]) + + +# --------------------------------------------------------------------------- +# The same guarantee on the other two carriers. The swarm variant failed +# differently -- loudly rather than silently -- because its pack was a flat +# reshape that ignored symmetric storage entirely: +# +# ValueError: could not broadcast input array from shape (156,4) +# into shape (156,3) +# +# so a symmetric tensor on a swarm could not be written through .array at all. +# --------------------------------------------------------------------------- + + +def test_a_swarm_symmetric_tensor_takes_an_off_diagonal_write(): + mesh = _box() + swarm = uw.swarm.Swarm(mesh) + stress = uw.swarm.SwarmVariable("tau", swarm, vtype=uw.VarType.SYM_TENSOR) + swarm.populate(fill_param=2) + + with uw.synchronised_array_update(): + stress.array[:, 0, 1] = 30.0 + + read_back = np.asarray(stress.array) + assert np.asarray(stress.data).shape[1] == 3, "symmetric storage is 3 wide in 2-D" + assert np.allclose(read_back[:, 0, 1], 30.0) + assert np.allclose(read_back[:, 1, 0], 30.0) + + +def test_an_integration_point_symmetric_tensor_takes_an_off_diagonal_write(): + mesh = _box() + stress = uw.discretisation.IntegrationPointVariable( + "tauq", mesh, vtype=uw.VarType.SYM_TENSOR) + + with uw.synchronised_array_update(): + stress.array[:, 0, 1] = 30.0 + + read_back = np.asarray(stress.array) + assert np.allclose(read_back[:, 0, 1], 30.0) + assert np.allclose(read_back[:, 1, 0], 30.0) + + +def test_a_swarm_vector_round_trips_unchanged(): + """The pack change must not disturb the shapes that already worked.""" + mesh = _box() + swarm = uw.swarm.Swarm(mesh) + velocity = uw.swarm.SwarmVariable("vel", swarm, vtype=uw.VarType.VECTOR) + swarm.populate(fill_param=2) + + with uw.synchronised_array_update(): + velocity.array[:, 0, 0] = 1.0 + velocity.array[:, 0, 1] = 2.0 + + assert np.allclose(np.asarray(velocity.data)[0], [1.0, 2.0]) + + +# --------------------------------------------------------------------------- +# Deferred writes. +# +# `synchronised_array_update()` defers the PETSc pack, and the swarm view read +# its "current" values straight from the PETSc field -- so inside that context +# every write started from the pre-context state and overwrote the one before +# it. Only the last survived, silently. Mesh variables were immune, because +# their view reads and writes the canonical array. +# +# Writing a vector or tensor one component at a time inside that context is +# exactly what docs/developer/subsystems/data-access.md recommends. +# --------------------------------------------------------------------------- + + +def _carrier(mesh, vtype, tag, carrier): + """The same variable on a mesh, a swarm, or the integration points. + + Returns the swarm alongside it: a SwarmVariable holds only a weak + reference to its swarm, so a local one is collected when the helper + returns and the variable then refuses to be read. + """ + if carrier == "mesh": + return uw.discretisation.MeshVariable( + f"m{tag}", mesh, vtype=vtype, degree=1), None + if carrier == "integration_point": + return uw.discretisation.IntegrationPointVariable( + f"q{tag}", mesh, vtype=vtype), None + swarm = uw.swarm.Swarm(mesh) + variable = uw.swarm.SwarmVariable(f"s{tag}", swarm, vtype=vtype) + swarm.populate(fill_param=2) + return variable, swarm + + +@pytest.mark.parametrize("carrier", ["mesh", "swarm", "integration_point"]) +def test_component_writes_in_one_context_all_survive(carrier): + """Every component written inside a single deferred context must land.""" + mesh = _box() + var, _swarm = _carrier(mesh, uw.VarType.VECTOR, f"v{carrier[:2]}", carrier) + + with uw.synchronised_array_update(): + var.array[:, 0, 0] = 1.0 + var.array[:, 0, 1] = 2.0 + + assert np.allclose(np.asarray(var.data)[0], [1.0, 2.0]), ( + f"{carrier}: an earlier write in the context was overwritten") + + +@pytest.mark.parametrize("carrier", ["mesh", "swarm", "integration_point"]) +@pytest.mark.parametrize("corner", [(0, 1), (1, 0)]) +def test_a_symmetric_off_diagonal_write_lands_on_every_carrier(carrier, corner): + """The product of shape, carrier and index — the cell of this matrix that + nobody had filled is where both defects lived.""" + mesh = _box() + i, j = corner + var, _swarm = _carrier( + mesh, uw.VarType.SYM_TENSOR, f"t{carrier[:2]}{i}{j}", carrier) + + with uw.synchronised_array_update(): + var.array[:, i, j] = 30.0 + + read_back = np.asarray(var.array) + assert np.allclose(read_back[:, i, j], 30.0), f"{carrier} {corner}: lost" + assert np.allclose(read_back[:, j, i], 30.0), f"{carrier} {corner}: not mirrored"