Second weak formulation and solution in one file not working - min. example incl.

Question

Dear FEniCS users,

I met a strange problem in FEniCS. You can try it by yourself. At first try to run the code and then comment the blocks marked by (2).

from dolfin import *

mesh = UnitSquareMesh(15,15)
V =  FunctionSpace(mesh, "CG", 3)

# Boundary conditions
def down(x, on_boundary): return x[0] < (0.0 + DOLFIN_EPS)
g0 = Constant(0.0)
bc = DirichletBC(V, g0, down)

# Parameters
epsilon = Constant(0.1)
b = Expression(('-x[1]', 'x[0]'))
f = Constant(1.)

uh = TrialFunction(V)
vh = TestFunction(V)

# (1)
a = (epsilon*dot(grad(uh),grad(vh)) + vh*dot(b,grad(uh)))*dx
L = f*vh*dx

# Compute solution (1)
uh = Function(V)
solve(a == L, uh, bc)

# (2)
psih = TrialFunction(V)
vh = TestFunction(V)
psih = Function(V)
a = (epsilon*dot(grad(psih),grad(vh)) + vh*dot(b,grad(psih)))*dx
L = f*vh*dx

# Compute solution (2)
psih = Function(V)
solve(a == L, psih, bc)

plot(uh)
interactive()

Why is fenics throwing the following error in the first case?

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0
*** 
*** DOLFIN version: 1.3.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Answer 1

Your left hand side is rank 1 since it does not contain a TrialFunction. You have defined psih twice in block 2, so psih is a Function when a is created and not a TrialFunction as it probably should be.

Comments

Thank you Mikael, the working code is then:

from dolfin import *

mesh = UnitSquareMesh(15,15)
V =  FunctionSpace(mesh, "CG", 3)

# Boundary conditions
def down(x, on_boundary): return x[0] < (0.0 + DOLFIN_EPS)
g0 = Constant(0.0)
bc = DirichletBC(V, g0, down)

# Parameters
epsilon = Constant(0.1)
b = Expression(('-x[1]', 'x[0]'))
f = Constant(1.)

uh = TrialFunction(V)
vh = TestFunction(V)

# (1)
a = (epsilon*dot(grad(uh),grad(vh)) + vh*dot(b,grad(uh)))*dx
L = f*vh*dx

# Compute solution (1)
uh = Function(V)
solve(a == L, uh, bc)

# (2)
psih = TrialFunction(V)
vh = TestFunction(V)
a = (epsilon*dot(grad(psih),grad(vh)) + vh*dot(b,grad(psih)))*dx
L = f*vh*dx

# Compute solution (2)
psih = Function(V)
solve(a == L, psih, bc)

plot(uh)
interactive()

Everything works fine now.

---

But why I need to have the line

uh = Function(V)

in the first case? (Otherwise i get an error...)

---

Because you need a container for your solution. uh in solve must be a Function. You can call it anything you like, doesn't have to be the same name as the TrialFunction.

---

Thank you, so if I got it from a novice's point of view, it is necessary to have everytime you use solve() the following form (if I am naming TrialFunction and Function of solution the same as my convention):

uh = TrialFunction(V)
vh = TestFunction(V)
a = ...; L = ...
uh = Function(V)
solve(a == L, uh, bc)

Thank you Mikael and Miro!

Answer 2

Hi, the problem in (2) is due to psih defined as Function and not TrialFunction. Because of this a is not a bilinear form. The fix is

psih = TrialFunction(V)
vh = TestFunction(V)
# psih = Function(V)
a = (epsilon*dot(grad(psih),grad(vh)) + vh*dot(b,grad(psih)))*dx
L = f*vh*dx

Comments

Thank you MiroK! It really was a basic question, but then:
why I need to have the line

uh = Function(V)

in the first case? (Otherwise i get an error...)

---

As the error message says

Expecting second argument to be a Function.

But

uh = TrialFunction(V)
type(uh)

tells you that uh in an Argument instance.

---

Thank you Miro, I got it!