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"