Error: Expecting the left-hand side to be a bilinear form (not rank 1)

General
1

Hello,
I am new to Fenics. I am trying to solve a simple 2d problem(code below).

\nabla^2 v_\alpha =\xi v_\alpha + g_\alpha (x,y) (alpha=x or y)

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

a=Constant(100)
xi=Constant(100)


mesh = UnitSquareMesh(32, 32)
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

v1, v2 = TestFunctions(V)

u = Function(V)
u1, u2 = split(u)

g=Expression('exp(-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.05)',degree=2)

F=-dot(grad(u1),grad(v1))-dot(grad(u2),grad(v2))-xi*(u1*v1+u2*v2)-g*v2

a=-dot(grad(u1),grad(v1))*dx-dot(grad(u2),grad(v2))*dx
L=xi*(u1*v1+u2*v2)*dx+g*v2*dx

be = Constant(0)
def boundary_L(x, on_boundary):
    return on_boundary and near(x[0], 0, 1e-14)

def boundary_R(x, on_boundary):
    return on_boundary and near(x[0], 1, 1e-14)

def boundary_B(x, on_boundary):
    return on_boundary and near(x[1], 0, 1e-14)

def boundary_T(x, on_boundary):
    return on_boundary and near(x[1], 1, 1e-14)

def boundary(x,on_boundary):
    return on_boundary

bcs=[DirichletBC(V.sub(0), be, boundary_L),DirichletBC(V.sub(0), be, boundary_R),DirichletBC(V.sub(1), be, boundary_B),DirichletBC(V.sub(1), be, boundary_T)]

solve(a==L,u,bcs)

However, I am getting the error: Unable to define linear variational problem a(u, v) == L(v) for all v.
With reason: Expecting the left-hand side to be a bilinear form (not rank 1).

What am I doing wrong?

2

Since you are trying to solve a linear problem, then u=TrialFunction(V)

one would then define uh = Function(V) prior to solve and call:
solve(a==L, uh, bcs)

3

Hello,

Thanks a lot for your answer. I have one more question. Are couple terms like,

u1*u2.dx(0)*v1.dx(0)

allowed in fenics or would it cause numerical instabilities? The only coupled example I could find is the chemical reaction example in the fenics book but the coupling there was not involving the gradient.

4

This all depends on your PDE. As FEniCS is a framework for solving PDEs, the stability of the variational formulation is a question of mathematical stability analysis

1 Like
5

Thanks. Sorry for not being specific.
For example, when I run the following code, I get the error Too many values to unpack.

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

a=Constant(100)
xi=Constant(100)
K1 = Constant(1)
tau = Constant(0.01)
nu = Constant(1)
q02 = Constant(0.25)
dt=0.01
k=Constant(dt)


mesh = UnitSquareMesh(32, 32)
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1,P1])
V = FunctionSpace(mesh, element)

v1, v2, s1, s2 = TestFunctions(V)

u = TrialFunction(V)
u1, u2, q1, q2 = split(u)
u_n=TrialFunction(V)
u1n,u2n,q1n,q2n=split(u_n)

g=Expression('exp(-((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1)',degree=2)

be = Constant(0)
def boundary_L(x, on_boundary):
    return on_boundary and near(x[0], 0, 1e-14)

def boundary_R(x, on_boundary):
    return on_boundary and near(x[0], 1, 1e-14)

def boundary_B(x, on_boundary):
    return on_boundary and near(x[1], 0, 1e-14)

def boundary_T(x, on_boundary):
    return on_boundary and near(x[1], 1, 1e-14)

def boundary(x,on_boundary):
    return on_boundary

bcs=[DirichletBC(V.sub(0), be, boundary_L),DirichletBC(V.sub(0), be, boundary_R),DirichletBC(V.sub(1), be, boundary_B),\
     DirichletBC(V.sub(1), be, boundary_T), DirichletBC(V.sub(2), be, boundary), DirichletBC(V.sub(3), be, boundary)]


F=-dot(grad(u1),grad(v1))*dx-dot(grad(u2),grad(v2))*dx-a*q1*u1.dx(0)*v1.dx(0)*dx-a*(q1*u2.dx(1)).dx(0)*v1*dx-xi*(u1*v1+u2*v2)*dx-g*v2*dx

a=lhs(F)
L=rhs(F)

uh=Function(V)
solve(a==L,uh,bcs)

Does this count as a non-linear problem? Or is there any way to run this? I would really appreciate the answer. I have not written the full equation for the sake of simplicity. But this is the third and fourth term in F is causing the problem.

6

Yesm, this is a non-linear problem. To make this work you should replace

with u = Function(V), and keep the residual form F (do not split into a, L) and solve
with
solve(F==0, u, bcs)

1 Like
7

Hello, Thank you for your answer. It worked. But now I have a time dynamic on the solution and I want to set the initial guess for the newton iteration to be the solution of the last step. I tried to follow the answer from a similar question, How to use the solution from the last step as initial guess of a nonlinear variation problem? - FEniCS Q&A (fenicsproject.org). But ran into errors. Specifically, I tried the following (continuing from the previous question in this thread),

u = TrialFunction(V)
u1, u2, q1, q2 = split(u)
u_n=Function(V)
u1n,u2n,q1n,q2n=split(u_n)

F=-dot(grad(u1),grad(v1))*dx-dot(grad(u2),grad(v2))*dx-a*(q1*(u1.dx(0)-u2.dx(1))+q2*(u1.dx(1)+u2.dx(0)))*(v1.dx(0)+v2.dx(1))*dx-xi*(u1*v1+u2*v2)*dx-g*v2*dx\
    +((q1-q1n)/k)*s1*dx+((q2-q2n)/k)*s2*dx+K1*dot(grad(q1),grad(s1))*dx+K1*dot(grad(q2),grad(s2))*dx\
    -((q02-q1*q1-q2*q2)/tau)*(q1*s1+q2*s2)*dx-nu*((u1.dx(0)-u2.dx(1))*s1+(u1.dx(1)+u2.dx(0))*s2)*dx

a=lhs(F)
L=rhs(F)
A = assemble(a)
b = assemble(L)
solver = KrylovSolver("cg", "ilu")
solver.set_operator(A)
sol = Function(V)

But ran into the error in the step A = assemble(a). It says too many values to unpack. Can you please let me know how to achieve this? To be clear, I want to set u_n as the initial guess for the solver.

Thank you.

8

Please make a minimal reproducible example. As your current code is just a snippet, which is really hard for anyone to decipher (or get an error message from).