Gradient of a Form as a vector

1

I have the following form and corresponding boundary conditions:

from fenics import *
mesh = UnitSquareMesh(16,16)
VectorSpaces=[["Lagrange",2],["Lagrange",1]]
V = VectorElement(VectorSpaces[0][0],mesh.ufl_cell(),VectorSpaces[0][1])
Q = FiniteElement(VectorSpaces[1][0],mesh.ufl_cell(),VectorSpaces[1][1])
K = FunctionSpace(mesh,"DG",0)#function space of piecewise constants
Welem = MixedElement([V,Q])
W = FunctionSpace(mesh,Welem)

up = TrialFunction(W)
u,p = split(up)
vq = TestFunction(W)
v,q = split(vq)

F = (
    inner(dot(grad(u),u),v)*dx - #Nonlinear Advection Term
    inner(p,div(v))*dx + #Pressure Gradient term
    inner(grad(u),grad(v))*dx + #Diffusion Term
    inner(q,div(u))*dx # Divergence-Free Condition Term
)


#boundary conditions
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
down   = 'near(x[1], 0)'
up = 'near(x[1], 1)'
al = 'near(x[0], 0) || near(x[0], 1) || near(x[1], 0) || near(x[1], 1)'
# Define boundary conditions
bcu_down  = DirichletBC(W.sub(0), Constant((0, 0)), down)
bcu_up  = DirichletBC(W.sub(0), Constant((1, 0)), up)
bcp_inflow = DirichletBC(W.sub(1), Expression("0",degree=2), inflow)#0 pressu
bcp_outflow = DirichletBC(W.sub(1), Expression("0",degree=2), outflow)#
bdc = [bcp_inflow,bcp_outflow,bcu_up,bcu_down]

How would I obtain a dolfin.cpp.la.Vector vector of the gradient of F given a specific function up.

2

You can use the functions ufl.action (Form language — Unified Form Language (UFL) 2021.1.0 documentation) to create a form that uses up as the input for u, and reduced the problem to a vector assembly.

3

Do you mean I manually find the gradient form and then apply action to that ?

4

Sorry, I misunderstood your question, you can use derivative(F, u, up) and assemble the result to get the gradient of F in direction up. See Form language — Unified Form Language (UFL) 2021.1.0 documentation

5

So the modified code would be (qith s being the function that is being inputted into up):

from fenics import *
mesh = UnitSquareMesh(16,16)
VectorSpaces=[["Lagrange",2],["Lagrange",1]]
V = VectorElement(VectorSpaces[0][0],mesh.ufl_cell(),VectorSpaces[0][1])
Q = FiniteElement(VectorSpaces[1][0],mesh.ufl_cell(),VectorSpaces[1][1])
K = FunctionSpace(mesh,"DG",0)#function space of piecewise constants
Welem = MixedElement([V,Q])
W = FunctionSpace(mesh,Welem)

up = TrialFunction(W)
u,p = split(up)
vq = TestFunction(W)
v,q = split(vq)

F = (
    inner(dot(grad(u),u),v)*dx - #Nonlinear Advection Term
    inner(p,div(v))*dx + #Pressure Gradient term
    inner(grad(u),grad(v))*dx + #Diffusion Term
    inner(q,div(u))*dx # Divergence-Free Condition Term
)


#boundary conditions
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
down   = 'near(x[1], 0)'
up_ = 'near(x[1], 1)'
al = 'near(x[0], 0) || near(x[0], 1) || near(x[1], 0) || near(x[1], 1)'
# Define boundary conditions
bcu_down  = DirichletBC(W.sub(0), Constant((0, 0)), down)
bcu_up  = DirichletBC(W.sub(0), Constant((1, 0)), up_)
bcp_inflow = DirichletBC(W.sub(1), Expression("0",degree=2), inflow)#0 pressu
bcp_outflow = DirichletBC(W.sub(1), Expression("0",degree=2), outflow)#
bdc = [bcp_inflow,bcp_outflow,bcu_up,bcu_down]

s = Function(W)
action(F,s)
gradF = assemble(derivative(F,s,s))
print(type(gradF))

But, that returns a matrix, not a vector. Am I misunderstanding something ?

6

Action does not do an in-place operation. You need to do as follows:

s = Function(W)
gradF = assemble(derivative(action(F, s),s,s))
print(type(gradF))
7

Ah yes, that fixes the issue.

Thank you !

Does this work regardless of boundary conditions ?

8

What do you mean when you say regardless of boundary conditions?
If you want to take the derivative in a certain direction, let’s call it ds, then you should do the following:

s0 = Function(W) # Point to take the derivative around
ds = Function(W) # Direction to differentiate in
gradF = assemble(derivative(action(F, s),s, ds))
print(type(gradF))