Compute derivatives with UFL when using Mixed Elements

variational formulation
1

Hi all,

I am trying to compute a Jacobian matrix using ufl.derivative.
I don’t have any problems when I use simple function spaces, but I don’t get the correct result when I am dealing with Mixed Elements.
Here is a minimal example:

import ufl
from dolfinx import fem, mesh, io, plot
from mpi4py import MPI

# Define mesh
domain =  mesh.create_rectangle(MPI.COMM_WORLD, [[0, 0],[1, 1]], [10, 10])
n = ufl.FacetNormal(domain)

# Function Spaces
el1 = ufl.FiniteElement("CG", domain.ufl_cell(), 2)
el2 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el2, el1])
W = fem.FunctionSpace(domain, mel)
num_subs = W.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
    space_i, map_i = W.sub(i).collapse()
    spaces.append(space_i)
    maps.append(map_i)
V = spaces[0]
Q = spaces[1]

# Trial and Test Functions
u, p = ufl.TrialFunctions(W)
v, q = ufl.TestFunctions(W)

# Define function
u_n = fem.Function(V)
p_n = fem.Function(Q)
F = ufl.inner(u_n,q)*ufl.dx + ufl.inner(p_n,v)*ufl.dx

# Take derivative of F w.r.t. u_n
dF_u = ufl.derivative(F, u_n, u)

# Assemble matrix (the derivative is a Jacobian matrix)
A = fem.petsc.assemble_matrix(fem.form(dF_u))
A.assemble()
C = A.convert('dense')
print(C.getDenseArray())

You can see that it returns a matrix full of zeros, but it should have some non-zero coefficients.

Here is a similar minimal example, but using simple elements:

import ufl
from dolfinx import fem, mesh, io, plot
from mpi4py import MPI

# Define mesh
domain =  mesh.create_rectangle(MPI.COMM_WORLD, [[0, 0],[1, 1]], [10, 10])
n = ufl.FacetNormal(domain)

# Function Spaces
el = ufl.FiniteElement("CG", domain.ufl_cell(), 2)
V = fem.FunctionSpace(domain, el)
Q = fem.FunctionSpace(domain, el)

# Trial and Test Functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
q = ufl.TestFunction(V)

# Define function
u_n = fem.Function(V)
p_n = fem.Function(Q)
F = ufl.inner(u_n,q)*ufl.dx + ufl.inner(p_n,v)*ufl.dx

# Take derivative of F w.r.t. u_n
dF_u = ufl.derivative(F, u_n, u)

# Assemble matrix (the derivative is a Jacobian matrix)
A = fem.petsc.assemble_matrix(fem.form(dF_u))
A.assemble()
C = A.convert('dense')
print(C.getDenseArray())

You can see that, in this case, the Jacobian matrix contains non-zero coefficients.

Does anybody know what is the problem with the mixed-space formulation?

Thanks in advance!

2

I think it does have nonzero elements. The diagonal blocks should be zeros by the definition of your weak form. The off-diagonal blocks should have nonzero entries. You can simply do a check on the dense array

...
C_dense = C.getDenseArray()
np.testing.assert_allclose(C_dense, 0.)

On another note, it looks like you want to define a bilinear form based on a mixed-formulation. If that is the case, what you would want to do is:

import ufl
from dolfinx import fem, mesh, io, plot
from mpi4py import MPI

# Define mesh
domain =  mesh.create_rectangle(MPI.COMM_WORLD, [[0, 0],[1, 1]], [10, 10])
n = ufl.FacetNormal(domain)

# Function Spaces
el1 = ufl.FiniteElement("CG", domain.ufl_cell(), 2)
el2 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el2, el1])
W = fem.FunctionSpace(domain, mel)

# Trial and Test Functions
w = ufl.TrialFunction(W)
u, p = ufl.split(w)
v, q = ufl.TestFunctions(W)

# Define function
w_n = fem.Function(W)
u_n, p_n = ufl.split(w_n)

F = ufl.inner(u_n,q)*ufl.dx + ufl.inner(p_n,v)*ufl.dx

# Take derivative of F w.r.t. u_n
dF_u = ufl.derivative(F, w_n, w)

# Assemble matrix (the derivative is a Jacobian matrix)
A = fem.petsc.assemble_matrix(fem.form(dF_u))
A.assemble()
3 Likes
3

Thank you! It worked.