Poisson-Neumann Condition (different source types)

Question

from dolfin import *
mesh = UnitSquareMesh(32, 32)
Mh = FunctionSpace(mesh, "Lagrange", 2)
Vh = VectorFunctionSpace(mesh, "Lagrange", 2)

t = 0.0
THETA = 0.05
U_exat = Expression(["THETAexp(-t)pow(x[0],2)pow((1-x[0]),2)(2x[1]pow((1-x[1]),2)-2pow(x[1],2)(1-x[1]))","-THETAexp(-t)pow(x[1],2)pow((1-x[1]),2)(2x[0]pow((1-x[0]),2)-2pow(x[0],2)(1-x[0]))"],THETA=THETA,t=0.0)
f = Expression("pow(x[0],2)pow((1-x[0]),2)(2x[1]pow((1-x[1]),2)-2pow(x[1],2)(1-x[1]))")
U_exat_pro = project(U_exat, Vh)
g = Expression("0")

p = TrialFunction(Mh)
q = TestFunction(Mh)

F = inner(grad(p), grad(q))dx + inner(f,q)dx + gqds + 1e-12pq*dx
a = lhs(F)
L = rhs(F)

rho = Function(Mh)
solve(a == L, rho)
plot(rho,title="rho",interactive=True)
plot(project(grad(rho),Vh),title="grad(rho)",interactive=True)

I find the different source type: div(u) and f, give different results.
Namely, substitute F with the form as follows,
F = inner(grad(p), grad(q))dx + inner(div(U_exat_pro),q)dx + gqds + 1e-12pq*dx
1, The source is f, the solution obeys the boundary condition
2, The source is div(u), the solution does not obey the boundary condition
Can you tell me why ????
Thanks so much.

Comments

Hi there,

can you please format your code properly, e.g. like below, so that it is easier to read? There is a corresponding button for that. Thanks!

Apart from that, looking at your code as it is, there are currently no boundary conditions implemented anywhere.

from dolfin import *

mesh = UnitSquareMesh(32, 32)
Mh = FunctionSpace(mesh, "Lagrange", 2)
Vh = VectorFunctionSpace(mesh, "Lagrange", 2)

---

from dolfin import *

mesh = UnitSquareMesh(32, 32)
Mh = FunctionSpace(mesh, "Lagrange", 2)
Vh = VectorFunctionSpace(mesh, "Lagrange", 2)

t = 0.0
THETA = 0.05
U_exat = Expression(["THETAexp(-t)pow(x[0],2)pow((1-x[0]),2)(2x[1]pow((1-x[1]),2)-2pow(x[1],2)(1-x[1]))","-THETAexp(-t)pow(x[1],2)pow((1-x[1]),2)(2x[0]pow((1-x[0]),2)-2pow(x[0],2)(1-x[0]))"],THETA=THETA,t=0.0)
f = Expression("pow(x[0],2)pow((1-x[0]),2)(2x[1]pow((1-x[1]),2)-2pow(x[1],2)(1-x[1]))")
U_exat_pro = project(U_exat, Vh)
g = Expression("0")

p = TrialFunction(Mh)
q = TestFunction(Mh)

F = inner(grad(p), grad(q))dx + inner(f,q)dx + gqds + 1e-12pq*dx
a = lhs(F)
L = rhs(F)

rho = Function(Mh)
solve(a == L, rho)
plot(rho,title="rho",interactive=True)
plot(project(grad(rho),Vh),title="grad(rho)",interactive=True)

I find the different source type: div(u) and f, give different results.
Namely, substitute F with the form as follows,
F = inner(grad(p), grad(q))dx + inner(div(U_exat_pro),q)dx + gqds + 1e-12pq*dx

1, The source is f, the solution obeys the boundary condition
2, The source is div(u), the solution does not obey the boundary condition
Can you tell me why ????
Thanks so much.

---

You haven't incorporated boundary conditions to your solver, so right now it is a neumann problem with null boundary condition. You should add the conditions as mentioned before and it will work. Wether or not your source term obeys the dirichlet boundary condition is of no significance. Apart from that, having a 1e-12 constant should make your problem prone to numerical instabilities, but I don't know the effect of this.

Best regards