BUG: Greater error than expected in quadrilateral GLL elements

[Olender]

Describe the bug
We were modifying our code to work with the latest Firedrake version, and every test with SEM quads now shows greater error.

Steps to Reproduce
Steps to reproduce the behavior:
The error initially showed in an acoustic wave propagation. Below I have a reduced version with an analytical solution:

import finat
import FIAT
from firedrake import *
import numpy as np


def gauss_lobatto_legendre_line_rule(degree):
    fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule
    fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1)
    finat_ps = finat.point_set.GaussLobattoLegendrePointSet
    points = finat_ps(fiat_rule.get_points())
    weights = fiat_rule.get_weights()
    return finat.quadrature.QuadratureRule(points, weights)

def gauss_lobatto_legendre_cube_rule(dimension, degree):
    make_tensor_rule = finat.quadrature.TensorProductQuadratureRule
    result = gauss_lobatto_legendre_line_rule(degree)
    for _ in range(1, dimension):
        line_rule = gauss_lobatto_legendre_line_rule(degree)
        result = make_tensor_rule([result, line_rule])
    return result


def analytical_solution(t, V, mesh_z, mesh_x):
    analytical = Function(V)
    x = mesh_z
    y = mesh_x
    analytical.interpolate(x * (x + 1) * y * (y - 1) * t)

    return analytical


def isDiag(M):
    i, j = np.nonzero(M)
    return np.all(i == j)

degree = 4
mesh = RectangleMesh(50, 50, 1.0, 1.0, quadrilateral=True)
mesh.coordinates.dat.data[:, 0] *= -1.0
mesh_z, mesh_x = SpatialCoordinate(mesh)
element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral')
V = FunctionSpace(mesh, element)
quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree)

u = TrialFunction(V)
v = TestFunction(V)

c = Function(V, name="velocity")
c.interpolate(1 + sin(pi*-mesh_z)*sin(pi*mesh_x))
u_n = Function(V)
u_nm1 = Function(V)
dt = 0.0005
t = 0.0
final_time = 1.0
u_nm1.assign(analytical_solution((t - 2 * dt), V, mesh_z, mesh_x))
u_n.assign(analytical_solution((t - dt), V, mesh_z, mesh_x))

q_xy = Function(V)
q_xy.interpolate(-(mesh_z**2) - mesh_z - mesh_x**2 + mesh_x)
q = q_xy * Constant(2 * t)

nt = int((final_time - t) / dt) + 1

m1 = (
    1
    / (c * c)
    * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2))
    * v
    * dx(scheme=quad_rule)
)
a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule)
le = q * v * dx(scheme=quad_rule)

form = m1 + a - le

B = Cofunction(V.dual())

boundary_ids = (1, 2, 3, 4)
bcs = DirichletBC(V, 0.0, boundary_ids)
A = assemble(lhs(form), bcs=bcs)
solver_parameters = {
    "ksp_type": "preonly",
    "pc_type": "jacobi",
}
solver = LinearSolver(
    A, solver_parameters=solver_parameters
)
As = solver.A
petsc_matrix = As.petscmat
diagonal = petsc_matrix.getDiagonal()
Mdiag = diagonal.array

u_np1 = Function(V)
for step in range(nt):
    q = q_xy * Constant(2 * t)
    m1 = (
        1
        / (c * c)
        * ((u - 2.0 * u_n + u_nm1) / Constant(dt**2))
        * v
        * dx(scheme=quad_rule)
    )
    a = dot(grad(u_n), grad(v)) * dx(scheme=quad_rule)
    le = q * v * dx(scheme=quad_rule)

    form = m1 + a - le

    B = assemble(rhs(form), tensor=B)

    solver.solve(u_np1, B)

    if (step - 1) % 100 == 0:
        print(f"Time : {t}")

    u_nm1.assign(u_n)
    u_n.assign(u_np1)

    t = step * float(dt)

u_an = analytical_solution(t, V, mesh_z, mesh_x)
error = errornorm(u_n, u_an)
print(f"Error: {error}")

print("END")

In order to debug, I have also done a simple u * v * dx assemble:

import finat
import FIAT
from firedrake import *
import numpy as np


def gauss_lobatto_legendre_line_rule(degree):
    fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule
    fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1)
    finat_ps = finat.point_set.GaussLobattoLegendrePointSet
    points = finat_ps(fiat_rule.get_points())
    weights = fiat_rule.get_weights()
    return finat.quadrature.QuadratureRule(points, weights)

def gauss_lobatto_legendre_cube_rule(dimension, degree):
    make_tensor_rule = finat.quadrature.TensorProductQuadratureRule
    result = gauss_lobatto_legendre_line_rule(degree)
    for _ in range(1, dimension):
        line_rule = gauss_lobatto_legendre_line_rule(degree)
        result = make_tensor_rule([result, line_rule])
    return result

degree = 2
# mesh = RectangleMesh(3, 2, 2.0, 1.0, quadrilateral=True)
mesh = RectangleMesh(1, 1, 1.0, 1.0, quadrilateral=True)
element = FiniteElement('CG', mesh.ufl_cell(), degree=degree, variant='spectral')
V = FunctionSpace(mesh, element)
quad_rule = gauss_lobatto_legendre_cube_rule(dimension=2, degree=degree)

u = TrialFunction(V)
v = TestFunction(V)

form = u*v*dx(scheme=quad_rule)
A = assemble(form)
M = A.M.values
Mdiag = M.diagonal()

x_mesh, y_mesh = SpatialCoordinate(mesh)
x_func = Function(V)
y_func = Function(V)
x_func.interpolate(x_mesh)
y_func.interpolate(y_mesh)
x = x_func.dat.data[:]
y = y_func.dat.data[:]

print("END")

Expected behavior
For the first code, I expected an errornorm close to 5e^-8, which was achieved in older versions of Firedrake. For the second code, I have used only a degree=2 element, because it is easier. According to x_func and y_func, the node 0 has position (0.5, 0.5) and the assembled diagonal matrix has Mdiag[0] value 0.02777. Solving for the matrix gave me 0.4444.

I have made a table below of x and y values, Mdiag values, and values I expected to see:

Node	x Value	y Value	Mdiag Value	Expected Value
0	0.5	0.5	0.0277	0.4444
1	0.0	0.5	0.0277	0.1111
2	0.5	1.0	0.1111	0.1111
3	1.0	0.5	0.1111	0.1111
4	0.5	0.0	0.0277	0.1111
5	0.0	0.0	0.0277	0.0277
6	0.0	1.0	0.1111	0.0277
7	1.0	1.0	0.4444	0.0277
8	1.0	0.0	0.1111	0.0277

Environment:

• OS: Ubuntu 20.04.6
• Python version: 3.8.10
• Output of firedrake-status:

DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
 __import__('pkg_resources').require('firedrake==0.13.0+6134.gcb77d32ac')
Firedrake Configuration:
   package_manager: True
   minimal_petsc: False
   mpicc: None
   mpicxx: None
   mpif90: None
   mpiexec: None
   disable_ssh: True
   honour_petsc_dir: False
   with_parmetis: False
   slepc: False
   packages: []
   honour_pythonpath: False
   opencascade: False
   tinyasm: False
   torch: False
   petsc_int_type: int32
   cache_dir: /home/olender/Firedrake/2024-06/firedrake/.cache
   complex: False
   remove_build_files: False
   with_blas: None
   netgen: False
Additions:
   None
Environment:
   PYTHONPATH: None
   PETSC_ARCH: None
   PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |f352325   |False     |
|PyOP2               |master                        |a52b9982  |False     |
|fiat                |master                        |fa86ed3   |False     |
|firedrake           |master                        |cb77d32ac |False     |
|h5py                |firedrake                     |6b512e5e  |False     |
|libsupermesh        |master                        |84becef   |False     |
|loopy               |main                          |967461ba  |False     |
|petsc               |firedrake                     |272f13d92c7|False     |
|pyadjoint           |master                        |849ee3e   |False     |
|pytest-mpi          |main                          |f2566a1   |False     |
|tsfc                |master                        |f80f29d   |False     |
|ufl                 |master                        |fbd288e6  |False     |
---------------------------------------------------------------------------

[connorjward]

@pbrubeck do you have any ideas?

[pbrubeck]

The CG dofs are not in lexicographic ordering anymore, but we still generate quadrature rules in lexicographic ordering. In 1D we order first the vertex dofs and then the interior ones.

[pbrubeck]

There is definitely an inconsistency with the DOF and quadrature orderings, I'll assign this bug to myself.

[pbrubeck]

A temporary fix is to reorder the GLL quadrature rule to match the new order the GLL DOFs.

def gauss_lobatto_legendre_line_rule(degree):
    fiat_make_rule = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule
    fiat_rule = fiat_make_rule(FIAT.ufc_simplex(1), degree + 1)
    finat_ps = finat.point_set.GaussLobattoLegendrePointSet
    perm = [0, degree, *list(range(1, degree))]
    points = finat_ps(fiat_rule.get_points()[perm])
    weights = fiat_rule.get_weights()[perm]
    return finat.quadrature.QuadratureRule(points, weights)

We wanted to make this quadrature rule more user-friendly, by supporting strings as the scheme kwarg:
https://github.com/firedrakeproject/tsfc/pull/306/files