Incorrect PETSc.Mat block size when using enriched_element

dolfinx
1

General info (even though I do not think it matters much):

  • MacOS 26.6.1
  • Python 3.14.6 (conda)
  • fenics-dolfinx 0.11.0 (conda)
  • petsc4py 3.25.4 (conda)

In my application, I am using the MINI element (P1 Lagrange enriched with bubble) for the (Navier-)Stokes equations. As I would like to use gamg or hypre as a multigrid preconditioner (which, at the moment, is working mediocre at best), I encountered a comment about gamg in the PETSc documentation (PCGAMG — PETSc 3.25.4 documentation) that describes the requirement “to indicate the number of degrees of freedom per grid point”. I therefore investigated and found that the submatrix (0,0) from the nested PETSc.Mat matrix only had a block_size of 1. I have broken my problem down further (the nested approach does not actually seem to matter) to pinpoint the issue (adapted from Stokes equations using Taylor-Hood elements — DOLFINx Python 0.12.0.dev0 ). A MWE can be found in the following (I know it is not a fully working solver (e.g., boundary conditions are missing), but I only included the absolutely necessary parts to reproduce the issue).

Am I doing something wrong, misunderstanding something, is this information actually not required, or is there actually an issue?

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

from dolfinx.mesh import CellType, create_rectangle
from basix.ufl import element, enriched_element
from dolfinx import default_real_type
from dolfinx.fem import functionspace, Constant
from ufl import TestFunction, TrialFunction, inner, dot, grad, dx
from dolfinx.fem.petsc import LinearProblem

# Create mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], (32, 32), CellType.triangle
)
gdim = msh.geometry.dim

# some definitions for the linear problem later
petsc_options_prefix = "block_size_test_"
petsc_options = {
    "ksp_type": "gmres",
    "ksp_rtol": 1e-9,
    "pc_type": "gamg",
}

# source term, zero vector for both approaches
f = Constant(msh, [PETSc.ScalarType(0.)] * gdim)

def use_P2():
    # P2 element
    P2 = element("Lagrange", msh.basix_cell(), degree=2, shape=(gdim,), dtype=default_real_type)
    V_P2 = functionspace(msh, P2)

    # velocity trial and test function
    u_P2 = TrialFunction(V_P2)
    v_P2 = TestFunction(V_P2)

    # variational problem
    a_ufl_P2 = inner(grad(u_P2), grad(v_P2)) * dx
    L_ufl_P2 = dot(f, v_P2) * dx

    # define solver
    problem_P2 = LinearProblem(
        a_ufl_P2,
        L_ufl_P2,
        # P=a_p,
        petsc_options_prefix = petsc_options_prefix + "P2_",
        petsc_options = petsc_options,
    )

    # print block size with P2 element
    print("Block size with P2 element:", problem_P2.A.getBlockSize())

    return

def use_P1plus():
    # enriched P1 element
    P1base = element("Lagrange", msh.basix_cell(), degree=1, shape=(gdim,), dtype=default_real_type)
    P1enrich = element("bubble", msh.basix_cell(), degree=3, shape=(gdim,), dtype=default_real_type)
    P1plus = enriched_element([P1base, P1enrich])
    V_P1plus = functionspace(msh, P1plus)

    # velocity trial and test function
    u_P1plus = TrialFunction(V_P1plus)
    v_P1plus = TestFunction(V_P1plus)

    # variational problem
    a_ufl_P1plus = inner(grad(u_P1plus), grad(v_P1plus)) * dx
    L_ufl_P1plus = dot(f, v_P1plus) * dx

    # define solver
    problem_P1plus = LinearProblem(
        a_ufl_P1plus,
        L_ufl_P1plus,
        # P=a_p,
        petsc_options_prefix=petsc_options_prefix + "P1plus_",
        petsc_options=petsc_options,
    )

    # print block size with enriched P1 element
    print("Block size with enriched P1 element:", problem_P1plus.A.getBlockSize())

    return

# run with P2 element
use_P2()

# run with enriched P1 element
use_P1plus()