BUG: Getri throws nonzero info when inverting a tensor with SLATE for DQ-2

[miguelcoolchips]

Describe the bug
I am trying to invert a block-diagonal matrix using Slate's infrastructure. It works great for simple mass matrix and for a grad-grad matrix in most cases except for DQ-2.

Steps to Reproduce

import firedrake as fd
import time
from firedrake import (
    inner,
    grad,
    dx,
)

N = 100
mesh = fd.UnitSquareMesh(N, N, quadrilateral=True)
V = fd.VectorFunctionSpace(mesh, "DQ", 2)  # Works fine with DG
v = fd.TestFunction(V)
u = fd.TrialFunction(V)
# F = inner(u, v) * dx  # Works for any order
F = inner(grad(u), grad(v)) * dx  # Works for 1, 3, 4, 5...
time_start = time.time()
fd.assemble(fd.Tensor(F).inv)
time_end = time.time()
print(f"Time taken for inverse: {time_end - time_start} seconds")

Expected behavior
No errors

Error message
Getri throws nonzero info.[1] 31916 abort

Environment:

• OS: MacOS
• Python version: 3.11.8
• Output of firedrake-status

Firedrake Configuration:
   package_manager: True
   minimal_petsc: False
   mpicc: mpicc
   mpicxx: mpicxx
   mpif90: mpif90
   mpiexec: mpiexec
   mpihome: /opt/homebrew/
   disable_ssh: False
   honour_petsc_dir: False
   with_parmetis: False
   slepc: True
   packages: []
   honour_pythonpath: False
   opencascade: False
   torch: False
   jax: False
   petsc_int_type: int32
   cache_dir: /Users/miguel/repos/firedrake-20250204/.cache
   complex: False
   remove_build_files: False
   with_blas: download
   netgen: False
Additions:
   None
Environment:
   PYTHONPATH: None
   PETSC_ARCH: None
   PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|fiat                |master                        |8839f876  |False     |
|firedrake           |master                        |96501362c |False     |
|h5py                |firedrake                     |06d2da80  |False     |
|libsupermesh        |master                        |f87cbfd   |True      |
|loopy               |main                          |27aead57  |False     |
|petsc               |firedrake                     |4354988f8b6|False     |
|pyadjoint           |master                        |da58d45   |False     |
|pytest-mpi          |main                          |f5668e4   |False     |
|slepc               |firedrake                     |cfa875895 |False     |
|ufl                 |master                        |5da0ab98  |False     |
---------------------------------------------------------------------------

[pbrubeck]

This is the expected behavior, as the laplacian with Neumann boundary conditions is singular.

[pbrubeck]

We should be throwing this error for any degree

[pbrubeck]

Why do you need DG? You could use SLATE to solve Poisson with Continuous Lagrange elements. This is standard static condensation, where the cell problems are posed on the interior with Dirichlet Boundary conditions on the facets, see this example: https://github.com/firedrakeproject/firedrake/blob/master/tests/firedrake/slate/test_cg_poisson.py

[miguelcoolchips]

Oh good point! I am trying to invert the block diagonals of a DQ discretization. I tried to stabilize the form using Nitsche boundary conditions, but I am still getting the same problem

n = fd.FacetNormal(mesh)
h = fd.CellSize(mesh)
alpha = fd.Constant(1e4 * k**2)
u_b = fd.Constant((0, 0))
F = (
    inner(grad(u), grad(v)) * dx
    - inner(dot(grad(u), n), v) * fd.ds  # Consistency term
    - inner(u, dot(grad(v), n)) * fd.ds  # Symmetry term
    + alpha / h * inner(u, v) * fd.ds  # Penalty term
)

F += (
    inner(u_b, dot(grad(v), n)) * fd.ds  # Symmetry term for BC
    - alpha / h * inner(u_b, v) * fd.ds  # Penalty term for BC
)

I think that because I am lacking the facet terms, the boundary conditions are not propagated to the element's interiors. Doesn't the HDG method do something along these lines? Do I need to impose boundary conditions on the elements somehow?

[pbrubeck]

I guess you are using IPDG to precondition a pressure Laplacian?

fd.ds is the measure on the boundary of the domain, so it only applies to the exterior facets. For interior-penalty DG, you need to include the constistency terms along the interior the facets with the measure fd.dS. The interior facet terms should always restrict the integrand to either the "+" or "-" sides of the facet (which are arbitrarily labeled).

[miguelcoolchips]

Yes, I want to use the IPDG to precondition a pressure Laplacian (I think). The SIMPLE method uses the inverse of the matrix diagonal as an approximation of the inverse of the original matrix. I want to use the element-block diagonal instead. Do you have a reference on using IPDG to precondition a pressure Laplacian? I could not find much