inner(u, div(w)) * dx incorrect for non-affine geometric mappings

[jpdean]

The following test assembles the matrix corresponding to inner(u, div(w)) * dx on a mesh of two square elements and two trapezium shaped elements, where u is from a "DQ" space and w is from an "RTCF" space. Due to the properties of the Piola transform, these matrices should be equal, however, the test fails.

I believe this issue is causing the vector unknown in a mixed Poisson problem to not converge in the L^2-norm (at all) on general quadrilateral meshes. There are, however, no issues on meshes of affine quadrilaterals. Note that this issue is not related to FEniCS/dolfinx#1455.

import pytest
import dolfinx
import ufl
import numpy as np
from mpi4py import MPI


def assemble_matrix(k, offset):
    mesh = create_mesh(offset)

    V = dolfinx.FunctionSpace(mesh, ("DQ", k))
    W = dolfinx.FunctionSpace(mesh, ("RTCF", k + 1))

    u = ufl.TrialFunction(V)
    w = ufl.TestFunction(W)

    a = ufl.inner(u, ufl.div(w)) * ufl.dx
    A = dolfinx.fem.assemble_matrix(a)
    A.assemble()
    return A[:, :]


def create_mesh(offset):
    """Creates a mesh of two square elements if offset = 0, or two
    trapezium elements if |offset| > 0.
    """
    x = np.array([[0, 0],
                  [1, 0],
                  [0, 0.5 + offset],
                  [1, 0.5 - offset],
                  [0, 1],
                  [1, 1]])
    cells = np.array([[0, 1, 2, 3],
                      [2, 3, 4, 5]])

    ufl_mesh = ufl.Mesh(ufl.VectorElement("Lagrange", "quadrilateral", 1))
    mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, x, ufl_mesh)
    return mesh


@pytest.mark.parametrize("k", [0, 1, 2])
def test_general_quads(k):
    # Assemble matrix on a mesh of square elements
    A_square = assemble_matrix(k, 0)
    # Assemble matrix on a mesh of trapezium elements
    A_trap = assemble_matrix(k, 0.25)

    # Due to the properties of the Piola transform, A_square and A_trap
    # should be equal
    assert np.allclose(A_square, A_trap, atol=1e-6)

[jpdean]

This issue appears to be related to the J_deriv_* terms that are present in the generated code. Whilst they have the correct values, I am unsure why they are needed here. Setting them to zero causes the above test to pass and also fixes the issue with the convergence rate of the vector unknown in a mixed Poisson problem.

[jpdean]

The J_deriv_* arise from the product rule. However, if my hand calculations are correct, the terms involving them should not actually contribute anything to the matrix. This would explain why setting them to zero gives the correct entries. I am currently unsure why they are contributing non-zero values.

The following is the simplest test that I can think of that demonstrates the problem:

import dolfinx
import ufl
import numpy as np
from mpi4py import MPI


def assemble_vector(k, offset):
    mesh = create_mesh(offset)

    V = dolfinx.FunctionSpace(mesh, ("RTCF", k + 1))
    v = ufl.TestFunction(V)

    l = ufl.inner(dolfinx.Constant(mesh, 1), ufl.div(v)) * ufl.dx
    L = dolfinx.fem.assemble_vector(l)

    return L[:]


def create_mesh(offset):
    """Creates a mesh of a single square element if offset = 0, or a
    trapezium element if |offset| > 0.
    """
    x = np.array([[0, 0],
                  [1, 0],
                  [0, 0.5 + offset],
                  [1, 0.5 - offset]])
    cells = np.array([[0, 1, 2, 3]])

    ufl_mesh = ufl.Mesh(ufl.VectorElement("Lagrange", "quadrilateral", 1))
    mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, x, ufl_mesh)
    return mesh


@pytest.mark.parametrize("k", [0, 1, 2])
def test_general_quads_vec(k):
    L_square = assemble_vector(k, 0)
    L_trap = assemble_vector(k, 0.25)

    assert np.allclose(L_square, L_trap, atol=1e-6)

[jpdean]

This issue also seems to affect curved triangles

[jpdean]

I think I've figured out the problem. The ufl for inner(u, div(w)) * dx eventually gets expanded out to include terms involving Grad(Jacobian(...)). However, for some reason, in the generated code, this gets interpreted as ReferenceGrad(Jacobian(...)). In other words, J_deriv_* in the generated code are the derivatives of the Jacobian with respect to the reference coordinate system, however, they are being used as though they are the derivatives with respect to the physical coordinate system. To test this hypothesis, I modified the generated code to compute the derivatives of the Jacobian with respect to the global coordinates and used these values in place of the J_deriv_* terms, and it seems to fix the issue (at least on the example I tried). I will investigate further.

[michalhabera]

@jpdean Can you move this issue to ffcx?

When you say it gets expanded in UFL, did you check after all preprocessings steps (apply pullbacks and derivatives)?

[jpdean]

Sure! I think it could be a UFL issue though? When I check preprocessed_form in https://github.com/FEniCS/ufl/blob/main/ufl/algorithms/compute_form_data.py#L411, it contains Grad(Jacobian(...)) terms, but by this stage they should have been converted to ReferenceGrad(Jacobian(...)), right?

After talking to @mscroggs, it looks like ffcx should just be checking there aren't any global derivatives and throwing and error if so: see this

ffcx/ffcx/ir/elementtables.py

Line 210 in df0284b

# Grad(Jacobian(...)) should be a local derivative

and compare with this https://github.com/FEniCS/ffcx/blame/82f9c336677bbd516da2cf0524540d9a888b5a90/ffcx/ir/elementtables.py#L325-L326.

I'll create a branch adding this check back in and see what else breaks. @michalhabera are you happy for me to move this issue to UFL rather than ffcx?

[jpdean]

I'll move this issue to FFCx and create a new issue in UFL demonstrating the problem.

[jpdean]

I've managed to fix this issue by removing the line https://github.com/FEniCS/ffcx/blob/main/ffcx/analysis.py#L163. PR #374 does this and adds back the check to make sure there are no global derivatives of the Jacobian present. @mscroggs, is there a reason why preserve_geometry_types=(ufl.classes.Jacobian, ) is needed?