From cca5b507ceb1a1417ce2618e64a731bb4df36eb9 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 12 Oct 2026 18:03:29 +0000 Subject: [PATCH] Give each matrix its own input vector in test_discrete_gradient The MatrixCSR branch of test_discrete_gradient built one input vector on the column index map of DOLFINx's discrete gradient and passed it to the mult of both matrices. MatrixCSR.mult indexes the vector with the matrix's own column numbering, so the vector has to be on that matrix's column map. On DOLFINx nightly the two maps differ: discrete_gradient uses V's map as is, while create_matrix, which interpolation_matrix uses, adds ghosts in another order. scifem's matrix then read past the end of the vector. Those columns have zero coefficients, so the product was usually right, and NaN when the memory read held one, which is why the failures came and went, and only in parallel. Build `u` on each matrix's column map instead. Co-Authored-By: Claude Opus 5.5 --- tests/test_interpolation.py | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 01ec1a9..3f144e0 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -134,16 +134,23 @@ def test_discrete_gradient(degree, use_petsc, cell_type): G_ref = dolfinx.fem.discrete_gradient(V, W) - # Built in matrices has to use a special input vector, with additional ghosts. - try: - _x = dolfinx.la.vector(G_ref.index_map(1), G_ref.block_size[1]) - except AttributeError: - # Bug in DOLFINx 0.9.0 - _x = dolfinx.la.vector(G_ref.index_map(1), G_ref.bs[1]) - num_owned_dofs = V.dofmap.index_map.size_local * V.dofmap.index_map_bs - _x.array[:num_owned_dofs] = u.x.array[:num_owned_dofs] - _x.scatter_forward() + + def input_vector(matrix): + """`u` on the column index map of `matrix`, which `MatrixCSR.mult` indexes into. + + The column map can have more ghosts than `V`'s, and two matrices on `V` need not + share one: on DOLFINx nightly, `discrete_gradient` uses `V`'s map while + `create_matrix` adds ghosts. A vector on another matrix's map is read out of bounds. + """ + try: + x = dolfinx.la.vector(matrix.index_map(1), matrix.block_size[1]) + except AttributeError: + # Bug in DOLFINx 0.9.0 + x = dolfinx.la.vector(matrix.index_map(1), matrix.bs[1]) + x.array[:num_owned_dofs] = u.x.array[:num_owned_dofs] + x.scatter_forward() + return x if use_petsc: A = scifem.interpolation.petsc_interpolation_matrix(expr, W) @@ -153,7 +160,7 @@ def test_discrete_gradient(degree, use_petsc, cell_type): if not hasattr(dolfinx.la.MatrixCSR, "mult"): pytest.skip("MatrixCSR has no mult method") A = scifem.interpolation.interpolation_matrix(expr, W) - A.mult(_x, w.x) + A.mult(input_vector(A), w.x) w.x.scatter_forward() w_ref = dolfinx.fem.Function(W) @@ -166,7 +173,7 @@ def test_discrete_gradient(degree, use_petsc, cell_type): G_ref.assemble() G_ref.mult(u.x.petsc_vec, w_ref.x.petsc_vec) else: - G_ref.mult(_x, w_ref.x) + G_ref.mult(input_vector(G_ref), w_ref.x) w_ref.x.scatter_forward() np.testing.assert_allclose(w.x.array, w_ref.x.array, rtol=1e-11, atol=1e-12)