Is there a way to get the partial derivative for a gradient vector?

[Swani]

Hello. Dear community

I was wondering if it’s possible to get the partial derivative(Cartesian coordinate) for a gradient vector in legacy FEniCS.
(this might be a bit silly, but I’m posting it for the sake of curiosity, please forgive me)

I was trying to implement the following formula in legacy FEniCS.

u_{elem}(x,y)=x^2+y^2 \to \nabla{u_{elem}(x,y)}=\left[{2x\atop2y}\right] \to {\partial{\nabla{u_{elem}(x,y)}}\over\partial{x}}=\left[{2\atop0}\right], {\partial{\nabla{u_{elem}(x,y)}}\over\partial{y}}=\left[{0\atop2}\right]

The code looks like:

from dolfin import *

dom = RectangleMesh(Point(0.0, 0.0), Point(2.0, 2.0), 2, 2, diagonal="right")
Vh = VectorFunctionSpace(dom, "CG", 1)
CG1 = FunctionSpace(dom, "CG", 1)
u_nodal = interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1)
print(f"size of u_nodal vector: {u_nodal.vector().get_local().size}")
print(f"u_nodal vector: {u_nodal.vector().get_local()}\n")

grd_u_nodal = project(grad(u_nodal), Vh)
print(f"size of grad(u_nodal) vector: {grd_u_nodal.vector().get_local().size}")
print(f"grad(u_nodal) vector: {grd_u_nodal.vector().get_local()}\n")

x0, x1 = MeshCoordinates(dom)
print(f"dgrad(u_nodal)dx: {assemble(diff(grad(u_nodal), variable(x0))*dx)}")

The output and error:

size of u_nodal vector: 9
u_nodal vector: [4. 1. 5. 0. 2. 8. 1. 5. 4.]

size of grad(u_nodal) vector: 18
grad(u_nodal) vector: [0.85714286 3.14285714 0.57142857 2.28571429 1.71428571 3.42857143
 0.28571429 0.28571429 2.         2.         3.71428571 3.71428571
 2.28571429 0.57142857 3.42857143 1.71428571 3.14285714 0.85714286]

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

I’m thinking that the error is probably caused by a difference between the array shapes, is there any way to fix this?
(And the reason for the decimals in the gradient array is numerical difference?)

I hope this clears up some of my questions, however minor.

[francesco-ballarin]

Anything wrong with u_nodal.dx(0).dx(0)?

Of course, you’d better be using a higher degree than 1 in your Vh, otherwise the second derivative will surely be zero.

[dokken]

You cannot assemble vector-expressions by themselfs. In Legacy Dolfin (and DOLFINx).
You can print grd_u_nodal because you have used project, that creates an inner product for the VectorFunctionSpace.

Please note that you can use

import dolfin

mesh = dolfin.UnitSquareMesh(10, 10)

x, y = dolfin.SpatialCoordinate(mesh)

u = x**2 + y**2

grad_u = dolfin.grad(u)

grad_u_ex = dolfin.as_vector((2 * x, 2 * y))

error = dolfin.assemble(
    dolfin.inner(grad_u - grad_u_ex, grad_u - grad_u_ex) * dolfin.dx
)
print(error)

Hessian_u = dolfin.grad(grad_u)

Hessian_u_ex = 2 * dolfin.Identity(2)

error = dolfin.assemble(
    dolfin.inner(Hessian_u - Hessian_u_ex, Hessian_u - Hessian_u_ex) * dolfin.dx
)
print(error)

which gives you exact accuracy, for instance if u_elem = x3 + y3 then

import dolfin

mesh = dolfin.UnitSquareMesh(10, 10)

x, y = dolfin.SpatialCoordinate(mesh)

u = x**3 + y**3

grad_u = dolfin.grad(u)

grad_u_ex = dolfin.as_vector((3 * x**2, 3 * y**2))

error = dolfin.assemble(
    dolfin.inner(grad_u - grad_u_ex, grad_u - grad_u_ex) * dolfin.dx
)
print(error)

Hessian_u = dolfin.grad(grad_u)

Hessian_u_ex = dolfin.as_tensor(((6 * x, 0), (0, 6 * y)))
error = dolfin.assemble(
    dolfin.inner(Hessian_u - Hessian_u_ex, Hessian_u - Hessian_u_ex) * dolfin.dx
)
print(error)

gives you

0.0
0.0

[Swani]

Many responses!
Thank you for always being available to reply my post. You’re a great help in solving my curiosity!

Dear @francesco-ballarin,
From comments, what does this mean?

@francesco-ballarin

u_nodal.dx(0).dx(0)

The output looked like that.
(I assume it was outputting a UFL form, am I right?)

(grad((grad(f_1988))[0]))[0]

Dear @dokken,
In reply, you used grad(u) and grad_u_ex to show that there was no difference between the two parameters. It was an informative(as was the Hessian code result also).

Furthermore, while checking your reply, I have some additional questions.
a. In the code, you showed the Hessian, shouldn’t the Hessian( \nabla^{2}u_{elem}(x,y) ) be viewed differently from the partial derivative of the gradient( {\partial\nabla{u_{elem}(x,y)}\over\partial{x} \ or \ \partial{y}} )?
(If my knowledge is wrong, please correct me.)

b. Since grad(u) is a UFL form, vector array extraction seems to be difficult, so I tried vector array extraction with as_vector, but with difficulty. Is there a way to extract either grad_u or grad_u_ex as a vector array? As far as I know, there doesn’t seem to be a way.
While this is an important example for noticing that the two results are not different, I think there is a difficulty in utilising the gradient(or Hessian) array.
(The reason I want to extract the values into a vector array is to compare them later with other operations).

c. This is a minor question, but is there a specific reason why you used inner product norm(similar as L^{2} norm squared) in the error checking?

[dokken]

@Swani

Furthermore, while checking your reply, I have some additional questions.
a. In the code, you showed the Hessian, shouldn’t the Hessian( \nabla^{2}u_{elem}(x,y) ∇2uelem(x,y) \nabla^{2}u_{elem}(x,y) ) be viewed differently from the partial derivative of the gradient( {\partial\nabla{u_{elem}(x,y)}\over\partial{x} \ or \ \partial{y}} ∂∇uelem(x,y)∂x or ∂y{\partial\nabla{u_{elem}(x,y)}\over\partial{x} \ or \ \partial{y}} )?
(If my knowledge is wrong, please correct me.)

You can do that by callinggrad_u.dx(0) or grad_u.dx(1).

@Swani

b. Since grad(u) is a UFL form, vector array extraction seems to be difficult, so I tried vector array extraction with as_vector, but with difficulty. Is there a way to extract either grad_u or grad_u_ex as a vector array? As far as I know, there doesn’t seem to be a way.
While this is an important example for noticing that the two results are not different, I think there is a difficulty in utilising the gradient(or Hessian) array.
(The reason I want to extract the values into a vector array is to compare them later with other operations).

Please be clear on how you want to do this comparison.
As the examples you have shown has been simple, analytic expressions, you could simply interpolate the analytical expression into the function space you want to work with. See for instance

Defining real part of complex equation using Expression() in FEniCS

Legacy dolfin doesn’t really have complex support. You can use sympy for this though: from dolfin import * mesh = UnitSquareMesh(10, 10) import sympy x, y, t = sympy.symbols("x[0], x[1] t", real=True) v_in = y * (1 - y) * sympy.exp(1j * t) v_real, _ = v_in.as_real_imag() expr = Expression(str(v_real), degree=2, t=0.0) print(assemble(expr * dx(domain=mesh))) expr.t = 2 print(assemble(expr * dx(domain=mesh)))

Convert sympy to ufl?

Using strings and eval: import sympy x = sympy.Symbol("x") cos = sympy.cos sympy_input = cos(x) + x ** 2 import fenics mesh = fenics.UnitIntervalMesh(10) cos = fenics.cos x = fenics.SpatialCoordinate(mesh)[0] ufl_output = eval(str(sympy_input)) print(fenics.assemble(ufl_output*fenics.dx)) However, is there a specfic reason for using sympy? UFL-forms can be differentiated etc, so in most use cases there is no need for using sympy. EDIT: You can also use: https://github.com/MiroK/ulfy

@Swani

c. This is a minor question, but is there a specific reason why you used inner product norm(similar as L^{2} L2L^{2} norm squared) in the error checking?

As pointed out in the previous part of the question, what I have described gives you a symbolic expression that you can evaluate when integrating. Thus I chose an L2 norm, as it Seems natural in this setting

[Swani]

Thank you for comment, @dokken! It was very helpful. I really appreciate it!
(Apologies for the late reply. )

In case it helps anyone reading this in the future, here’s my understanding based on the code you answered(bottom of my reply, “MWE”).
To wrap things up, I noticed a difference in the u array results between the code I originally attached and the code you replied to, and between projection and interpolation.

Despite being a simple mesh(I changed the mesh size 10 \times 10 to 2 \times 2 to make it easier to see).

Code:

from dolfin import *

mesh = UnitSquareMesh(2, 2)
CG1 = FunctionSpace(mesh, "CG", 1)

x, y = SpatialCoordinate(mesh)
u = x**2 + y**2
print(f"size of u vector: {project(u, CG1).vector().get_local().size}")
print(f"u vector: {project(u, CG1).vector().get_local()}")
print(f"interpolate u: {interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1).vector().get_local()}\n")

Output:

size of u vector: 9
u vector: [ 0.93571429  0.16428571  1.16428571 -0.09285714  0.42142857  1.90714286
  0.16428571  1.16428571  0.93571429]
interpolate u: [1.   0.25 1.25 0.   0.5  2.   0.25 1.25 1.  ]

There was a difference.
(When I say “gradient operation” later, I mean the same thing as in the code I referenced, but only interpolating the u part, or running it with a difference of u=x**2+y**2)

For me, intuitively, the interpolated result of the two values looks right, but when I did a “gradient operation” with interpolate u, it was different from the exact solution.
On the other hand, when I checked the u array through project, it intuitively looked like there was an error, but when I did the “gradient operation”, it matched the exact solution.

In my opinion, defining u as u=x**2+y**2 rather than using interpolate seems to be wrong to the eye, but I think it is correct considering the operation, is there any way to improve this?
(I was hoping to be able to use projection and still see results that were reasonably consistent with the interpolation results)

MWE:

from dolfin import *

mesh = UnitSquareMesh(2, 2)
CG1 = FunctionSpace(mesh, "CG", 1)
Vh = VectorFunctionSpace(mesh, "CG", 1)
Th = TensorFunctionSpace(mesh, "CG", 1)

x, y = SpatialCoordinate(mesh)
u = x**2 + y**2

# Hmm... different!
print(f"size of u vector: {project(u, CG1).vector().get_local().size}")
print(f"u vector: {project(u, CG1).vector().get_local()}")
print(f"interpolate u: {interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1).vector().get_local()}\n")

grad_u = grad(u)
grad_u_ex = as_vector((2*x, 2*y))

print(f"grad_u size: {project(grad_u, Vh).vector().size()}")
print(f"grad_u vector: {project(grad_u, Vh).vector()[:]}")
print(f"x of grad_u size: {project(grad_u.dx(0), Vh).vector().get_local().size}")
print(f"x of grad_u vector: {project(grad_u.dx(0), Vh).vector().get_local()}")
print(f"y of grad_u size: {project(grad_u.dx(1), Vh).vector().get_local().size}")
print(f"y of grad_u vector: {project(grad_u.dx(1), Vh).vector().get_local()}")
grad_error = assemble(inner(grad_u - grad_u_ex, grad_u - grad_u_ex) * dx)
print(f"error between grad_u and grad_u_ex: {grad_error}\n")

Hessian_u = grad(grad_u)
Hessian_u_ex = 2*Identity(2)

print(f"Hessian size: {project(Hessian_u, Th).vector().get_local().size}")
print(f"Hessian vector: {project(Hessian_u, Th).vector().get_local()}")
hes_error = assemble(inner(Hessian_u - Hessian_u_ex, Hessian_u - Hessian_u_ex) * dx)
print(f"error between Hessian_u and Hessian_u_ex: {hes_error}")

Output:

size of u vector: 9
u vector: [ 0.93571429  0.16428571  1.16428571 -0.09285714  0.42142857  1.90714286
  0.16428571  1.16428571  0.93571429]
interpolate u: [1.   0.25 1.25 0.   0.5  2.   0.25 1.25 1.  ]

grad_u size: 18
grad_u vector: [-9.58434720e-17  2.00000000e+00  2.51534904e-17  1.00000000e+00
  1.00000000e+00  2.00000000e+00  8.32667268e-17  6.44227717e-17
  1.00000000e+00  1.00000000e+00  2.00000000e+00  2.00000000e+00
  1.00000000e+00  7.13714802e-17  2.00000000e+00  1.00000000e+00
  2.00000000e+00 -8.32667268e-17]
x of grad_u size: 18
x of grad_u vector: [2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0.]
y of grad_u size: 18
y of grad_u vector: [0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2. 0. 2.]
error between grad_u and grad_u_ex: 0.0

Hessian size: 36
Hessian vector: [2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2.
 2. 0. 0. 2. 2. 0. 0. 2. 2. 0. 0. 2.]
error between Hessian_u and Hessian_u_ex: 0.0

[dokken]

@Swani

To wrap things up, I noticed a difference in the u array results between the code I originally attached and the code you replied to, and between projection and interpolation.

Despite being a simple mesh(I changed the mesh size 10 \times 10 10×1010 \times 10 to 2 \times 2 2×22 \times 2 to make it easier to see).

Code:

from dolfin import *

mesh = UnitSquareMesh(2, 2)
CG1 = FunctionSpace(mesh, "CG", 1)

x, y = SpatialCoordinate(mesh)
u = x**2 + y**2
print(f"size of u vector: {project(u, CG1).vector().get_local().size}")
print(f"u vector: {project(u, CG1).vector().get_local()}")
print(f"interpolate u: {interpolate(Expression("x[0]*x[0]+x[1]*x[1]", degree=2), CG1).vector().get_local()}\n")

Output:

size of u vector: 9
u vector: [ 0.93571429  0.16428571  1.16428571 -0.09285714  0.42142857  1.90714286
  0.16428571  1.16428571  0.93571429]
interpolate u: [1.   0.25 1.25 0.   0.5  2.   0.25 1.25 1.  ]

There was a difference.

Of course there will be a difference between interpolation and projection, see for instance: Lecture 8 — Lecture notes MATMEK-4270

@Swani

In my opinion, defining u as u=x**2+y**2 rather than using interpolate seems to be wrong to the eye, but I think it is correct considering the operation, is there any way to improve this?

Using Expression the string describing u will be evaluated at interpolation points, while using x, y = SpatialCoordinate(mesh) gives you an evaluation at the quadrature points used in the projection scheme.
This means that you can control the accuracy of approximation by choosing a suitable quadrature rule (up to the best possible approximation).

[Swani]

Thanks for your reply!
I’ll do some more studying based on what you’ve show me! Thanks!
I’ll leave a comment if I don’t understand something later!