Deal with Error: too many values to unpack while solving mixed element space

Linear Algebra
1

Hello,
I am trying to use Lagrange multiplier to impose \mathbf{u.n}=0 on \partial \Omega while solving the poisson equation \Delta\mathbf{u}=\mathbf{f}. First I wanted to know wether the Lagrange multipliers provide a soft condition. Because having normal \mathbf{u} non zero on the boundaries ain’t acceptable.
Since my mesh is sometimes round I can’t use Dirichlet boundary conditions.
After multiple tries and research my code is the following:

#With a given mesh
#Declare the spaces
Vx = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
Vy = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
Lm= FiniteElement("Lagrange", mesh.ufl_cell(), 1)

# Make a mixed space
TH = MixedElement([Vx,Vy,Lm])
Esp= FunctionSpace(mesh, TH)


n = FacetNormal(mesh)

(ux,uy, l) = TrialFunctions(Esp)
(vx,vy, m) = TestFunctions(Esp)

a = (dot(grad(ux), grad(vx))+dot(grad(uy), grad(vy)))*dx\
    -(dot(grad(ux),n)*vx+dot(grad(uy),n)*vy)*ds\
    +(ux*n[0]+uy*n[1])*m*l*ds #The part I'm not sure about...
    
L = (f[0]*vx+f[1]*vy)*dx

solu = Function(Esp)
#Work fine tile here
solve(a==L, solu)

Then I get the error:

ValueError: too many values to unpack (expected 2)

The solution I found is to use a penalisation coefficient P=10^{15} before the term \int\mathbf{u.n}\,\mathrm{d}s which I multiply with the symetric expression for the test function to get the bilinear form:

P\int\mathbf{u.n}\,\mathrm{d}s\int\mathbf{v.n}\,\mathrm{d}s

VectorFunctionSpace(mesh, 'P', 2)
n = FacetNormal(mesh)

u = TrialFunction(V)
v = TestFunction(V)

a = (dot(grad(u[0]), grad(v[0]))+dot(grad(u[1]), grad(v[1])))*dx\
    -(dot(grad(u[0]),n)*v[0]+dot(grad(u[1]),n)*v[1])*ds\
    +P*dot(u,n)*dot(v,n)*ds
L = dot(f,v)*dx

u = Function(V)
solve(a==L, u)

This solution seems to work, yet it might be due to the simplicity of poisson equation.
I need to have proper methods before moving onto more complex equations (the Toner-Tu equation for active matter if you know…)

Any kind of comment or help would be appreciated.
Thanks a lot in advance!

PS: Sorry for the grammar mistakes I am just a french…

2

The problem is that +(ux*n[0]+uy*n[1])*m*l*ds is nonlinear in the trial function (ux,uy,l) (because it contains products of ux/uy with l), while the linear solver expects to have a be bilinear. To see the correct (linear) formulation, start from the Lagrangian

\int_{\partial\Omega}\lambda(\mathbf{u}\cdot\mathbf{n} - g)\,ds

(where the boundary data g happens to be zero in your case), and take a Gateaux derivative with respect to the trial function (\mathbf{u},\lambda), in the direction of the test function (\mathbf{v},\mu). This will give you two terms:

\int\lambda\mathbf{v}\cdot\mathbf{n}\,ds + \int\mu(\mathbf{u}\cdot\mathbf{n}-g)\,ds\text{ .}

The terms

\int\lambda\mathbf{v}\cdot\mathbf{n}\,ds + \int\mu\mathbf{u}\cdot\mathbf{n}\,ds

are bilinear in the trial and test functions, and thus part of the “left-hand side” form a, while

\int\mu(-g)\,ds

is only linear in the test function, and would need to be moved to the right-hand side form l (picking up a minus sign from changing sides).

1 Like
3

Thank you so much, I didn’t know the Gateaux derivative method.
This has totally answered my question!

4

The code is working but not the physics…
Since my code is:

 V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Lm= FiniteElement("Real", mesh.ufl_cell(), 0)

Esp= FunctionSpace(mesh, V*Lm)


n = FacetNormal(mesh)

(u, lamb) = TrialFunctions(Esp)
(v, mu) = TestFunctions(Esp)

a = (dot(grad(u[0]), grad(v[0]))+dot(grad(u[1]), grad(v[1])))*dx\
    -(dot(grad(u[0]),n)*v[0]+dot(grad(u[1]),n)*v[1])*ds\
    +lamb*dot(v,n)*ds+mu*dot(u,n)*ds
L = dot(f,v)*dx

solu = Function(Esp)
solve(a==L, solu)
u,lamb=solu.split(deepcopy=True)
print(u.vector()[:])
print(lamb.vector()[:])

What I get is not what I wanted…
The condition \mathbf{u.n}=0 ain’t respected on the boundaries…

When I try to change the finite element family to Lagrange with degree 1 what I get is only nan…
Am I still doing something wrong or this problem is not meant to be solve like that?

5

You need to remove the second line of a, i.e., the “consistency term”. If you view the full problem as an energy minimization, the original energy functional would be

J(u) = \frac{1}{2}\int_\Omega(\nabla u:\nabla u)\,d\Omega - \int_\Omega fu\,d\Omega

to which you would add a constraint Lagrangian to use the Lagrange multiplier approach. Taking the Gateaux derivative of J(u) in the direction v gives you the integrated-by-parts Poisson operator without the boundary term.

Another way to look at it is that the multiplier \lambda is the boundary data for a weakly-enforced Neumann/flux boundary condition, while the constraint equation weighted against \mu closes the system and ensures that the flux boundary data \lambda is selected such that the Dirichlet boundary condition u\cdot n = 0 holds.

EDIT: I should add that this does still leave a tacit homogeneous Neumann BC on the tangential flux. Often that is what is physically-relevant when you are constraining only the normal component of a vector-valued unknown, but it can take some care. For example, the divergence of viscous stress in an incompressible flow is frequently simplified to a Laplacian of the velocity, but integrating the Laplacian term by parts means that the natural flux no longer has the interpretation of traction.

6

Ok thanks I see what you mean I will try to read more about all of this, and I will try your solution. And I have no other constraint over the tangential flux, so I hope this will work.

Thanks a lot for your help!

7

I would be careful with assuming that there is “no other constraint over the tangential flux”. You could truly impose no BC on the tangential component of flux by keeping the consistency term with the substitution

\mathbf{v}~\to~(\mathbf{I} - \mathbf{n}\otimes\mathbf{n})\mathbf{v}\text{ ,}

i.e., replacing the test function with its tangential component. However, this would not clearly be well-posed; try checking the bilinear form for coercivity over the subspace satisfying the constraint.

Often I see students confuse homogeneous flux BCs with “no BCs”, since they are usually imposed by omitting a term from the variational form and colloquially called “do-nothing BCs”. For instance, a free-slip boundary condition in continuum mechanics, where the normal component of velocity/displacement is constrained to zero, also includes the condition that the tangential traction is zero, even though it doesn’t explicitly appear in a typical finite element formulation.

8

OK, I may not understood all what you told me so I will speak with my teacher about this issue.
What I am studying is active matter and a fluid description of it (Toner Tu equation). The microscopic system is composed of self propelled particles, these particles roll so I think that there is no other constraint. It may appear not physically consistent since there are not all the term are not yet taken into account, but I wanted to deal with simple equations before moving onto the more complex one.