Problem with arity in bilinear form with complex mode

[marazzaf]

Hi y'all,

I get a strange arity mismatch when using derivative on a energy.

Here's a snippet giving the error:

import sys
from firedrake.petsc import PETSc

L, H = 15,15
N = 50
mesh = RectangleMesh(N, N, L/2, H/2, originX=-L/2,originY=-H/2, diagonal='crossed')
h = L/N #mesh size
V = FunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, 'CG', 1)
Z = V * W
PETSc.Sys.Print('Nb dof: %i' % Z.dim())

#Elastic energy
sol = Function(Z, name='sol')
xi,lmbda = split(sol)
E = inner(xi, xi) * dx + inner(grad(xi), grad(xi)) * dx

#Weak formulation
v,mu = TestFunctions(Z)

#Weak formulation
a = inner(grad(mu), grad(xi)) * dx #constraint

#Lagrangian part
L = E - inner(grad(xi), grad(lmbda)) * dx
a += derivative(L, xi, v)

#Dirichlet BC
val_max = .3
x = SpatialCoordinate(mesh)
xi_D = -2*x[1]/H * val_max
bcs = [DirichletBC(Z.sub(0), xi_D, 3), DirichletBC(Z.sub(0), xi_D, 4), DirichletBC(Z.sub(1), Constant(0), 3), DirichletBC(Z.sub(1), Constant(0), 4)]

#Newton solver
solve(a == 0, sol, bcs=bcs, solver_parameters={'snes_monitor': None, 'snes_max_it': 25})

[colinjcotter]

I don’t think you need the third argument of derivative. On 11 Mar 2024, at 18:36, marazzaf ***@***.***> wrote:  Hi y'all, I get a strange arity mismatch when using derivative on a energy. Here's a snippet giving the error: import sys from firedrake.petsc import PETSc L, H = 15,15 N = 50 mesh = RectangleMesh(N, N, L/2, H/2, originX=-L/2,originY=-H/2, diagonal='crossed') h = L/N #mesh size V = FunctionSpace(mesh, "CG", 2) W = FunctionSpace(mesh, 'CG', 1) Z = V * W PETSc.Sys.Print('Nb dof: %i' % Z.dim()) #Elastic energy sol = Function(Z, name='sol') xi,lmbda = split(sol) E = inner(xi, xi) * dx + inner(grad(xi), grad(xi)) * dx #Weak formulation v,mu = TestFunctions(Z) #Weak formulation a = inner(grad(mu), grad(xi)) * dx #constraint #Lagrangian part L = E - inner(grad(xi), grad(lmbda)) * dx a += derivative(L, xi, v) #Dirichlet BC val_max = .3 x = SpatialCoordinate(mesh) xi_D = -2*x[1]/H * val_max bcs = [DirichletBC(Z.sub(0), xi_D, 3), DirichletBC(Z.sub(0), xi_D, 4), DirichletBC(Z.sub(1), Constant(0), 3), DirichletBC(Z.sub(1), Constant(0), 4)] #Newton solver solve(a == 0, sol, bcs=bcs, solver_parameters={'snes_monitor': None, 'snes_max_it': 25}) — Reply to this email directly, view it on GitHub<#3446>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4U437EH5Y3SHU2FCK3YXX2RDAVCNFSM6AAAAABEQ277T2VHI2DSMVQWIX3LMV43ERDJONRXK43TNFXW4OZWGM2TQOJWGM>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[marazzaf]

If I remove the 3rd argument in 'derivative', I get the following error:

File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/home/marazzato/firedrake_complex/src/firedrake/firedrake/ufl_expr.py", line 283, in derivative
raise ValueError("Shapes of u and du do not match.\n"
ValueError: Shapes of u and du do not match.
If you passed an indexed part of split(u) into derivative, you need to provide an appropriate du as well.

[colinjcotter]

I think this is because you need to take the derivative using the full unsplit variable. On 11 Mar 2024, at 21:39, marazzaf ***@***.***> wrote:  If I remove the 3rd argument in 'derivative', I get the following error: File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func File "/home/marazzato/firedrake_complex/src/firedrake/firedrake/ufl_expr.py", line 283, in derivative raise ValueError("Shapes of u and du do not match.\n" ValueError: Shapes of u and du do not match. If you passed an indexed part of split(u) into derivative, you need to provide an appropriate du as well. — Reply to this email directly, view it on GitHub<#3446 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4WCTGPO77IHQHR7A43YXYQB3AVCNFSM6AAAAABEQ277T2VHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DONJRHEYDO>. You are receiving this because you commented.Message ID: ***@***.***>

[marazzaf]

I doubt that is the issue.

The following code still gives me the arity error:

N = 10
mesh = RectangleMesh(N, N, L/2, H/2, originX=-L/2,originY=-H/2, diagonal='crossed')
h = L/N #mesh size
V = FunctionSpace(mesh, "CG", 1)
PETSc.Sys.Print('Nb dof: %i' % V.dim())

#Elastic energy
xi = Function(V, name='sol')
E = inner(xi, xi) * dx + inner(grad(xi), grad(xi)) * dx

#Weak formulation
v = TestFunction(V)
a = derivative(E, xi, v)

#Dirichlet BC
val_max = .3
x = SpatialCoordinate(mesh)
xi_D = -2*x[1]/H * val_max
bcs = [DirichletBC(V, xi_D, 3), DirichletBC(V, xi_D, 4)]

#Newton solver
solve(a == 0, xi, bcs=bcs, solver_parameters={'snes_monitor': None, 'snes_max_it': 25})

[dham]

I cannot reproduce the above error. Using a newly updated Firedrake and adding the following lines to actually define the missing symbols, this runs to completion.

from firedrake import *

L = 1
H = 1

[marazzaf]

And you are running in complex mode?

Because my firedrake is fairly recent too.

[dham]

Nothing in your post suggested that you were running in complex mode. Taking derivatives of Lagrangians in complex mode doesn't work for maths reasons. Take inner(u, u)*dx as a simple example. This is equivalent to dot(u, conj(u))*dx but conj is not complex differentiable.

[marazzaf]

You mean that the title was not explicit enough?

Also, for the energy that you suggest, the following post explains why there is no Frechet derivative: https://math.stackexchange.com/questions/4096030/is-there-a-fr%C3%A9chet-derivative-for-the-norm-on-a-complex-hilbert-space

I guess that settles the issue.