DOLFIN

compute Directional derivative

Asked by Melanie Jahny

I'd like to compute the directional derivative of a function in direction u:
F(u) = \nabla_u (p^T * h(u))
with h(u) = (y, u \cdot \nabla \phi )_{L_2}

and u = \sum_k u_k \psi_k
y = \sum_i y_i \phi_i

That means h(u) reads:
h(u)_ij = \sum_k u_k \sum_T \int_T \phi_j \nabla \phi_i \cdot psi_k dx * y_j

DF(u) in direction of u_k should be a scalar for each u_k:
p^T * \sum_T \int_T \phi \nabla \phi \cdot psi_k dx * y

That means, DF(u) containing in each entry the directional derivative should be a vector
of length u.

My question is how to implement the single phi_k in the derivative.

Here is the code i've tried:

x_max = 10.0 # maximum x-coordinate
y_max = 10.0 # maximum y-coordinate

# divergence of mesh
x_div = 10
y_div = 10

mesh = Rectangle(0,0, x_max, y_max, x_div, y_div)
mesh.order()

V = FunctionSpace(mesh, "Lagrange",2)
U = FunctionSpace(mesh, "Raviart-Thomas", 1)

y = Function(V)
phi = TrialFunction(V)

u = Function(U)
p = Function(V)

h = y*dot(grad(phi),u) *dx

# I've tried to use a function on Function Space W to extract the \psi - Functions,
but it always results in a vector of dimension of y instead of u.

W = FunctionSpace(mesh, "DG", 0)
w = Function(W)

DF_k = inner(p,y)*dot(grad(phi),w*u)*dx
vec = assemble(DF_k) # VECTOR OF DIMENSION OF Y INSTEAD OF U ...

I don't find my mistake, so please could you help me?! Thanks a lot in advance.

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Martin Sandve Alnæs
Solved:
Last query:
Last reply:
Revision history for this message
Martin Sandve Alnæs (martinal) said : 		#1

Change the space phi lives in to U?

Martin

Den 20. sep. 2011 kl. 18:05 skrev Melanie Jahny <email address hidden>:

> New question #171784 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/171784
>
> I'd like to compute the directional derivative of a function in direction u:
> F(u) = \nabla_u (p^T * h(u))
> with h(u) = (y, u \cdot \nabla \phi )_{L_2}
>
> and u = \sum_k u_k \psi_k
> y = \sum_i y_i \phi_i
>
> That means h(u) reads:
> h(u)_ij = \sum_k u_k \sum_T \int_T \phi_j \nabla \phi_i \cdot psi_k dx * y_j
>
> DF(u) in direction of u_k should be a scalar for each u_k:
> p^T * \sum_T \int_T \phi \nabla \phi \cdot psi_k dx * y
>
> That means, DF(u) containing in each entry the directional derivative should be a vector
> of length u.
>
> My question is how to implement the single phi_k in the derivative.
>
> Here is the code i've tried:
>
> x_max = 10.0 # maximum x-coordinate
> y_max = 10.0 # maximum y-coordinate
>
>
> # divergence of mesh
> x_div = 10
> y_div = 10
>
> mesh = Rectangle(0,0, x_max, y_max, x_div, y_div)
> mesh.order()
>
> V = FunctionSpace(mesh, "Lagrange",2)
> U = FunctionSpace(mesh, "Raviart-Thomas", 1)
>
> y = Function(V)
> phi = TrialFunction(V)
>
> u = Function(U)
> p = Function(V)
>
> h = y*dot(grad(phi),u) *dx
>
> # I've tried to use a function on Function Space W to extract the \psi - Functions,
> but it always results in a vector of dimension of y instead of u.
>
> W = FunctionSpace(mesh, "DG", 0)
> w = Function(W)
>
>
> DF_k = inner(p,y)*dot(grad(phi),w*u)*dx
> vec = assemble(DF_k) # VECTOR OF DIMENSION OF Y INSTEAD OF U ...
>
> I don't find my mistake, so please could you help me?! Thanks a lot in advance.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Revision history for this message
Melanie Jahny (melanie-jahny) said : 		#2

Thanks Martin for your answer, but I cannot change the spaces of y, phi, u and p, it's fixed for
my problem formulation.
(I tried instead, but this causes an error when I change the space of phi; Then grad(phi) becomes a 2x2 matrix).
Is it the right way to extract the single \psi - Functions with a function on W (W = FunctionSpace(mesh, "DG", 0) )? I'm not sure about that...

Revision history for this message
Best Martin Sandve Alnæs (martinal) said : 		#3

Hi, sorry I mixed up phi in psi when comparing your maths to the code.

But without considering your problem, only the technical part,
tracking dimensions from variational form to assembled vector is quite
easy:

v1 = TestFunction(V1)
v2 = TestFunction(V2)
L1 = ...*v1*dx # Assembles to a vector of dim(V1)
L2 = ...*v2*dx # Assembles to a vector of dim(V2)

If you make \phi a TrialFunction like in your code, then you will
assemble DF_k "for all \phi=\phi_k, k=1,...". If you give u as a
Function, then u is a given function with specific values.

Do you have \phi as a Function? Or do you actually intend to assemble
the matrix of all directional derivatives? In that case, you can do

  du = TestFunction(U)
  J = derivative(DF_k, u, du)
  A = assemble(J) # Matrix, d(DF_k)/du

Not sure if this is transposed of whay you want, so maybe you need to
make du the TrialFunction and phi the TestFunction instead.

Martin

On 20 September 2011 19:10, Melanie Jahny
<email address hidden> wrote:
> Question #171784 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171784
>
>    Status: Answered => Open
>
> Melanie Jahny is still having a problem:
> Thanks Martin for your answer, but I cannot change the spaces of y, phi, u and p, it's fixed for
> my problem formulation.
> (I tried instead, but this causes an error when I change the space of phi;  Then grad(phi) becomes a 2x2 matrix).
> Is it the right way to extract the single \psi - Functions with a function on W (W = FunctionSpace(mesh, "DG", 0) )? I'm not sure about that...
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Revision history for this message
Melanie Jahny (melanie-jahny) said : 		#4

Thanks a lot, Martin, for your answer. I have always been confused about the
difference between Function and TrialFunction.
And thank you for the advice of using the "derivative" - Function.
I think I solved it now, using it and defining phi as TestFunction and psi
as TrialFunction:

y = Function(V)
phi = TestFunction(V)
u = Function(U)
psi = TrialFunction(U)
h = y*dot(grad(phi),u) *dx
Dh = assemble(derivative(h, u, psi))