Attribute _blocks lost when using SNES

LucasPa · March 30, 2026, 9:48am

Hello,

Below is a MWE I created to show an issue with PETSc.SNES on a mixed formulation. The vector F_vec below is created with F_vec = dolfinx.fem.petsc.create_vector(F) and therefore receives the attribute _blocks. However this attribute is somehow lost. An error is then raised in function_fbecause dolfinx.fem.petsc.assemble_vector requires the attribute.

import dolfinx
import scifem
import ufl
from mpi4py import MPI
from petsc4py import PETSc

msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 4, 4)

Vspace  = dolfinx.fem.functionspace(msh, ("DG", 2))
f = dolfinx.fem.Function(Vspace)
f_test = ufl.TestFunction(Vspace)

Rspace = scifem.create_real_functionspace(msh)
r = dolfinx.fem.Function(Rspace)
r_test = ufl.TestFunction(Rspace)

form = ufl.inner(f, f_test) * ufl.dx
form_r = ufl.inner(r, r_test) * ufl.dx
J = dolfinx.fem.form(ufl.derivative(form, f))

form_block = [form, form_r]

Jff = ufl.derivative(form, f)
Jfr = ufl.derivative(form, r)
Jrf = ufl.derivative(form_r, f)
Jrr = ufl.derivative(form_r, r)
J_block = [[Jff, Jfr], [Jrf, Jrr]]

F = dolfinx.fem.form(form_block)
J = dolfinx.fem.form(J_block)

F_vec = dolfinx.fem.petsc.create_vector(F)
J_mat = dolfinx.fem.petsc.assemble_matrix(J)

def function_f(solver, x, F_vec) -> None:
    """Assemble the residual."""
    print('Attribute blocks ',F_vec.getAttr('_blocks'))
    dolfinx.fem.petsc.assemble_vector(F_vec, F)
    F_vec.set(1.0)
    return
def function_j(solver, x, J, P):
    """Assemble the Jacobian matrix."""
    J.zeroEntries()        # ensure it's zero
    J.assemble()
    J.shift(1.0)           # A = A + 1 * I
    return
solution = dolfinx.fem.petsc.create_vector(F)
snes = PETSc.SNES().create(msh.comm)
snes.setFunction(function_f, F_vec)
snes.setJacobian(function_j, J=J_mat, P=None)
snes.setObjective(lambda x,y : 1.0)
snes.solve(None, solution)

The output is:

Attribute blocks  ((0, 192, 193), (193, 193, 193))
Attribute blocks  None
Traceback (most recent call last):
  File "petsc4py/PETSc/PETSc.pyx", line 348, in petsc4py.PETSc.PetscPythonErrorHandler
  File "petsc4py/PETSc/PETSc.pyx", line 348, in petsc4py.PETSc.PetscPythonErrorHandler
  File "petsc4py/PETSc/PETSc.pyx", line 348, in petsc4py.PETSc.PetscPythonErrorHandler
  [Previous line repeated 2 more times]
  File "petsc4py/PETSc/petscsnes.pxi", line 334, in petsc4py.PETSc.SNES_Function
  File "/stck/fpascal/Python/Code/FENICS/PbBlocks/issue.py", line 38, in function_f
    dolfinx.fem.petsc.assemble_vector(F_vec, F)
  File "/opt/tools/python/3.10.8-gnu831/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/opt/Softwares/FE_libraries/dolfinx/Dist/real/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 339, in _
    offset0, offset1 = b.getAttr("_blocks")
TypeError: cannot unpack non-iterable NoneType object

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/stck/fpascal/Python/Code/FENICS/PbBlocks/issue.py", line 52, in <module>
    snes.solve(None, solution)
  File "petsc4py/PETSc/SNES.pyx", line 1724, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code -1
[0] SNESSolve() at /opt/Softwares/petsc/src/snes/interface/snes.c:4846
[0] SNESSolve_NEWTONLS() at /opt/Softwares/petsc/src/snes/impls/ls/ls.c:235
[0] SNESLineSearchApply() at /opt/Softwares/petsc/src/snes/linesearch/interface/linesearch.c:643
[0] SNESLineSearchApply_BT() at /opt/Softwares/petsc/src/snes/linesearch/impls/bt/linesearchbt.c:333
[0] SNESComputeFunction() at /opt/Softwares/petsc/src/snes/interface/snes.c:2490

The two first lines show that the attribute is removed at some point.

PS : The ugly fix is to add this code at the beginning of function_f:

 F_vec.setAttr('_blocks', solver.getFunction()[0].getAttr('_blocks'))

But maybe we don’t want to have to use this fix.

dokken · March 31, 2026, 6:52am

This is a subtle issue with PETSc, as SNES internally creates some new working vectors, but it doesn’t copy over attributes set to the original vector.
As seen in the assemble_residual used in DOLFINx, we have an internal function that assigns block data:

github.com/FEniCS/dolfinx

python/dolfinx/fem/petsc.py

da835de50


      
          if isinstance(residual, Sequence) and b.getType() != PETSc.Vec.Type.NEST:  # type: ignore[attr-defined]
              _assign_block_data(residual, b)
              _assign_block_data(residual, x)

which calls:

github.com/FEniCS/dolfinx

python/dolfinx/fem/petsc.py

da835de50


      
                  Note:
                      The function(s) do not share memory with the solution
                      vector ``x``.
                  """
                  return self._u
          
          
          # -- High-level interface for SNES ---------------------------------------
          
          
          def _assign_block_data(forms: Sequence[Form], vec: PETSc.Vec):  # type: ignore[name-defined]
              """Assign block data to a PETSc vector.
          
              Args:
                  forms: List of forms to extract block data from.
                  vec: PETSc vector to assign block data to.
              """
              maps = (
                  (
                      form.function_spaces[0].dofmaps(0).index_map,  # type: ignore[attr-defined]
                      form.function_spaces[0].dofmaps(0).index_map_bs,  # type: ignore[attr-defined]

LucasPa · March 31, 2026, 8:54am

Thank you for your answer. Please note that I created as well an issue on the PETSc gitlab . In particular I posted a MWE 100% PETSc where a vector attribute is not deleted (see below). What is the difference with the dolfinx implementation given above where the attribute disappears?

The MWE:

from petsc4py import PETSc
from mpi4py import MPI

ts = PETSc.TS().create()

F_vec = PETSc.Vec().create(comm=PETSc.COMM_SELF)
F_vec.setSizes(10)
F_vec.setUp()
F_vec.setAttr('attribute', 'value')
assert(F_vec.getAttr('attribute') == 'value')

mat = PETSc.Mat().create()
mat.setSizes((10,10))
mat.setFromOptions()

def evalFunction(ts, t, x, xdot, f):
    assert(f.getAttr('attribute') == 'value')
    return
def evalJacobian(ts, t, x, xdot, a, A, B):
    return

ts.setIFunction(evalFunction, F_vec)
ts.setIJacobian(evalJacobian, mat)

ts.setEquationType(PETSc.TS.EquationType.IMPLICIT)
ts.setProblemType(PETSc.TS.ProblemType.NONLINEAR)

x = PETSc.Vec().create(comm=PETSc.COMM_SELF)
x.setSizes(10)
x.setUp()
ts.setSolution(x)

ts.setTime(0.0)
ts.setTimeStep(1.0)
ts.setMaxSteps(5)

opts = PETSc.Options()
opts['ts_adapt_type'] = 'basic'
opts['ksp_type'] = 'preonly'
opts['pc_type'] = 'lu'
ts.setFromOptions()
ts.solve()

Matthew_Knepley · March 31, 2026, 12:02pm

@dokken This time the issue is not with PETSc making duplicates, because the vector returned by TSGetIFunction() should be the same as the one originally set. The only way this changes is someone is calling TSSetIFunction() in the interim. I want him to run in the debugger and get a stack trace when this happens. Could you do that?

dokken · March 31, 2026, 12:13pm

Hi Matt,
Isn’t the difference between the two examples here that one uses PETSc TS while the other one uses PETSc SNES?

LucasPa · March 31, 2026, 12:15pm

Please note that I get this issue both with SNES and TS.

dokken · March 31, 2026, 12:32pm

Do you have a pure PETSc snes problem that you could run where the attribute is copied?

And could you do as Matt suggests:

In general DOLFINx doesn’t do anything magical when adding the attributes:

github.com/FEniCS/dolfinx

python/dolfinx/la/petsc.py

594052ca7


      
          
          maps = [(index_map, bs) for index_map, bs in maps]
          off_owned = tuple(
              itertools.accumulate(maps, lambda off, m: off + m[0].size_local * m[1], initial=0)
          )
          off_ghost = tuple(
              itertools.accumulate(
                  maps, lambda off, m: off + m[0].num_ghosts * m[1], initial=off_owned[-1]
              )
          )
          vec.setAttr("_blocks", (off_owned, off_ghost))

github.com/FEniCS/dolfinx

python/dolfinx/la/petsc.py

594052ca7


      
              index_map, bs = maps[0]
              ghosts = index_map.ghosts.astype(PETSc.IntType)  # type: ignore[attr-defined]
              size = (index_map.size_local * bs, index_map.size_global * bs)
              b = PETSc.Vec().createGhost(ghosts, size=size, bsize=bs, comm=index_map.comm)  # type: ignore
              if kind == PETSc.Vec.Type.MPI:  # type: ignore[attr-defined]
                  _assign_block_data(maps, b)
              return b
          
          if kind is None or kind == PETSc.Vec.Type.MPI:  # type: ignore[attr-defined]
              b = dolfinx.cpp.fem.petsc.create_vector_block(maps)
              _assign_block_data(maps, b)
              return b
          elif kind == PETSc.Vec.Type.NEST:  # type: ignore[attr-defined]
              return dolfinx.cpp.fem.petsc.create_vector_nest(maps)
          else:
              raise NotImplementedError(
                  "Vector type must be specified for blocked/nested assembly."
                  f"Vector type '{kind}' not supported."
                  "Did you mean 'nest' or 'mpi'?"
              )

github.com/FEniCS/dolfinx

python/dolfinx/fem/petsc.py

594052ca7


      
          self._A = create_matrix(self.J, kind=kind)
          # Create PETSc structure for preconditioner if provided
          if self._preconditioner is not None:
              self._P_mat = create_matrix(self._preconditioner, kind=kind)
          else:
              self._P_mat = None
          
          # Determine the vector kind based on the matrix type
          kind = "nest" if self._A.getType() == PETSc.Mat.Type.NEST else kind  # type: ignore[attr-defined]
          assert kind is None or isinstance(kind, str)
          self._b = create_vector(_extract_function_spaces(self.F), kind=kind)  # type: ignore
          self._x = create_vector(_extract_function_spaces(self.F), kind=kind)  # type: ignore
          
          # Create the SNES solver and attach the corresponding Jacobian and
          # residual computation functions
          self._snes = PETSc.SNES().create(self.A.comm)  # type: ignore[attr-defined]
          self.solver.setJacobian(
              partial(assemble_jacobian, u, self.J, self.preconditioner, bcs), self.A, self.P_mat
          )
          self.solver.setFunction(partial(assemble_residual, u, self.F, self.J, bcs), self.b)

I know that @garth implemented the “_blocks” attribute so maybe he can supply more information/context

LucasPa · March 31, 2026, 1:05pm

I created a MWE 100% PETSc that triggers the problem:

from petsc4py import PETSc
from mpi4py import MPI

snes = PETSc.SNES().create()

F_vec = PETSc.Vec().create(comm=PETSc.COMM_SELF)
F_vec.setSizes(10)
F_vec.setUp()
F_vec.setAttr('_blocks', 'value')
assert(F_vec.getAttr('_blocks') == 'value')

mat = PETSc.Mat().create()
mat.setSizes((10,10))
mat.setFromOptions()

def evalFunction(snes, x, f):
    print(f.getAttr('_blocks'),'/',snes.getFunction()[0].getAttr('_blocks'))
    f.set(1.0)
    return
def evalJacobian(snes, x, A, B):
    A.assemble()
    A.shift(1.0)
    return

snes.setFunction(evalFunction, F_vec)
snes.setJacobian(evalJacobian, mat)


x = PETSc.Vec().create(comm=PETSc.COMM_SELF)
x.setSizes(10)
x.setUp()
snes.setSolution(x)


snes.setIterationNumber(1)

opts = PETSc.Options()
opts['ksp_type'] = 'preonly'
opts['pc_type'] = 'lu'
snes.setFromOptions()
snes.solve()

The first lines of the output (see the print in evalFunction) are:

value / value
None / value
None / value

This shows that the attribute _blocks “disappeared”.

dokken · March 31, 2026, 1:53pm

Could you setattr for x as well?

LucasPa · March 31, 2026, 2:03pm

The script below shows that _blocks is removed from x as well.

from petsc4py import PETSc
from mpi4py import MPI

snes = PETSc.SNES().create()

F_vec = PETSc.Vec().create(comm=PETSc.COMM_SELF)
F_vec.setSizes(10)
F_vec.setUp()
F_vec.setAttr('_blocks', 'value')
assert(F_vec.getAttr('_blocks') == 'value')

mat = PETSc.Mat().create()
mat.setSizes((10,10))
mat.setFromOptions()

def evalFunction(snes, x, f):
    print(f.getAttr('_blocks'),'/',snes.getFunction()[0].getAttr('_blocks'),'/',x.getAttr('_blocks'))
    f.set(1.0)
    return
def evalJacobian(snes, x, A, B):
    A.assemble()
    A.shift(1.0)
    return

snes.setFunction(evalFunction, F_vec)
snes.setJacobian(evalJacobian, mat)


x = PETSc.Vec().create(comm=PETSc.COMM_SELF)
x.setSizes(10)
x.setUp()
x.setAttr('_blocks', 'value')
snes.setSolution(x)


snes.setIterationNumber(1)

opts = PETSc.Options()
opts['ksp_type'] = 'preonly'
opts['pc_type'] = 'lu'
snes.setFromOptions()
snes.solve()

dokken · March 31, 2026, 3:18pm

This is more along the lines of what I’ve seen with the attributes in PETSc. @Matthew_Knepley any suggestions?

EDIT: I’ve tried to adapt DOLFINx to instead use the args and kwargs positional input arguments, rather than using functools.partial, which obfuscates the input arguments, at: Use PETSc CTX for SNES solver by jorgensd · Pull Request #4139 · FEniCS/dolfinx · GitHub

Matthew_Knepley · April 1, 2026, 1:48pm

The attributes are Python-only. I am now discussing it with Lisandro.

dokken · April 1, 2026, 4:49pm

That’s fine by me, I kind of like the option of sending args and kwargs to the function instead, as done in the PR to DOLFINx. I wasn’t aware that it was an option.

Topic		Replies	Views
Issues with using self created PETSc mats instead of dolfinx PETSc mats in SNES Linear Algebra	2	813	July 11, 2022
Using petsc4py.PETSc.SNES directly Linear Algebra	23	4820	June 13, 2024
Error using PETSc SNES General	5	440	April 16, 2024
From 0.9 to 0.10: Issues at custom SNES Interface for blocked problems dolfinx dolfinx	2	83	January 8, 2026
Setting SNES solver with petsc4py segfaults Linear Algebra	8	1402	July 5, 2023

Attribute _blocks lost when using SNES

Related topics