Mixed problem, no solution

Please provide the version of DOLFIN you are using and the output of what
DOLFIN tells you when you actually run the code. If you also provide a minimal
runnable code snippet that would also be very useful!

Johan

On Monday January 10 2011 10:00:04 B. Emek Abali wrote:
> New question #140910 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/140910
>
> setting a mixed space (velocity field v, pressure field p) as
>
> V = VectorFunctionSpace(mesh, 'CG', 1)
> Q = FunctionSpace(mesh, 'DG', 0)
> ME = MixedFunctionSpace([V,Q])
>
> and by describing
>
> null = Constant(0.0)
> class TopDrag(Expression):
> def eval(self, vec, x):
> vec[0] = 10.0
> vec[1] = 0.0
> vec[2] = 0.0
> def dim(self):
> return 3
>
> vOnTop = TopDrag()
>
> wanted to reach:
>
> a given velocity on a subdomain (not here but defined in a proper way with
> compile_subdomains() ) DirichletBC(ME.sub(0), vOnTop, top)
>
> zero velocity in one direction
> DirichletBC(ME.sub(0).sub(1), null, front)
>
> and a given pressure on a subdomain
> DirichletBC(ME.sub(1), pIn, left)
>
> which I couldn't :)
>
> Any sharp eyes where I am not thinking as FEniCS wants?

using the Ubuntu repository updated,

Dependencies:
10.09-2~ppa1~lucid1 - libdolfin0-dev (2 0.9.9) python-dolfin (2 0.9.9) dolfin-doc (2 0.9.9) dolfin-bin (2 0.9.9) python-ferari (2 0.2.0-1) python-ffc (2 0.9.4) python-fiat (2 0.9.2-1) python-instant (2 0.9.8-1) ufc (2 1.4.2) python-ufc (2 1.4.2) ufc-doc (2 1.4.2) python-viper (2 0.4.6-1) python-ufl (2 0.5.4) python-ufl-doc (2 0.5.4) syfi-dev (2 0.6.1) python-syfi0 (2 0.6.1) syfi-doc (2 0.6.1) syfi-bin (2 0.6.1)
9.12.1 - dolfin-dev (2 0.9.5) python-pydolfin0 (2 0.9.5) dolfin-doc (0 (null)) dolfin-bin (0 (null)) python-ferari (2 0.0.2) python-ffc (2 0.7.1) python-fiat (2 0.3.5) python-instant (2 0.9.6) ufc (2 1.2.0) python-ufc (2 1.2.0) ufc-doc (0 (null)) python-viper (2 0.4.5) python-ufl (2 0.4.1) syfi-dev (2 0.6.0) python-syfi0 (2 0.6.0) syfi-doc (0 (null)) syfi-bin (0 (null))

it is an iterative nonlinear problem, by cutting all the leaves away, this code leads still to NaN solutions

from dolfin import *
import numpy

def FormF(v,w,p,q):
 F = w[i]*v[i].dx(j)*v[j]*dx
 return F

def FormG(v,w,dv,p,q,dp):
 G = w[i]*dv[i].dx(j)*v[j]*dx
 return G

mesh = Box(0.0, -0.5, -0.5, 2.0, 0.5, 0.5, 10, 5, 5)
V = VectorFunctionSpace(mesh, 'CG', 1)
Q = FunctionSpace(mesh, 'CG', 1)
ME = MixedFunctionSpace([V,Q])

subdomains = MeshFunction("uint", mesh, mesh.topology().dim()-1)
right, left, back, front, top, bottom = compile_subdomains(["(fabs(x[0] - 2.0) < DOLFIN_EPS) && on_boundary", "(fabs(x[0] + 0.0) < DOLFIN_EPS) && on_boundary", "(fabs(x[1] - 0.5) < DOLFIN_EPS) && on_boundary", "(fabs(x[1] + 0.5) < DOLFIN_EPS) && on_boundary", "(fabs(x[2] - 0.5) < DOLFIN_EPS) && on_boundary", "(fabs(x[2] + 0.5) < DOLFIN_EPS) && on_boundary"])

null = Constant(0.0)
class TopDrag(Expression):
 def eval(self, vec, x):
  vec[0] = 1.0
  vec[1] = 0.0
  vec[2] = 0.0
 def dim(self):
  return 3

class BottomHold(Expression):
 def eval(self, vec, x):
  vec[0] = 0.0
  vec[1] = 0.0
  vec[2] = 0.0
 def dim(self):
  return 3

vOnTop = TopDrag()
vOnBottom = BottomHold()
pOut = 100.0 #right side pressure
pIn = 100.0 #left side pressure

bc = [DirichletBC(ME.sub(0), vOnBottom, bottom), DirichletBC(ME.sub(0), vOnTop, top), DirichletBC(ME.sub(0).sub(1), null, front), DirichletBC(ME.sub(0).sub(1), null, back), DirichletBC(ME.sub(1), pIn, left), DirichletBC(ME.sub(1), pOut, right)]

u = Function(ME)
v = u.split()[0]
p = u.split()[1]
omega = TestFunction(ME)
w, q = split(omega)
du = TrialFunction(ME)
dv, dp = split(du)

F = FormF(v,w,p,q)
G = FormG(v,w,dv,p,q,dp)
a = lhs(F-G)
L = rhs(F-G)
A = assemble(a)
b = assemble(L)
for condition in bc:
 condition.apply(A, b)

du = Function(ME)
solve(A, du.vector(), b)

print 'some value:', u((0.5,0.5,0.5))
print 'some value:', du((0.5,0.5,0.5))

I cannot spot anything right now.

You might want to write down the math of your form and maybe someone else
might give you a hint.

Johan

On Monday January 10 2011 10:36:53 B. Emek Abali wrote:
> Question #140910 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/140910
>
> B. Emek Abali posted a new comment:
> it is an iterative nonlinear problem, by cutting all the leaves away,
> this code leads still to NaN solutions
>
> from dolfin import *
> import numpy
>
> def FormF(v,w,p,q):
> F = w[i]*v[i].dx(j)*v[j]*dx
> return F
>
> def FormG(v,w,dv,p,q,dp):
> G = w[i]*dv[i].dx(j)*v[j]*dx
> return G
>
> mesh = Box(0.0, -0.5, -0.5, 2.0, 0.5, 0.5, 10, 5, 5)
> V = VectorFunctionSpace(mesh, 'CG', 1)
> Q = FunctionSpace(mesh, 'CG', 1)
> ME = MixedFunctionSpace([V,Q])
>
> subdomains = MeshFunction("uint", mesh, mesh.topology().dim()-1)
> right, left, back, front, top, bottom = compile_subdomains(["(fabs(x[0] -
> 2.0) < DOLFIN_EPS) && on_boundary", "(fabs(x[0] + 0.0) < DOLFIN_EPS) &&
> on_boundary", "(fabs(x[1] - 0.5) < DOLFIN_EPS) && on_boundary",
> "(fabs(x[1] + 0.5) < DOLFIN_EPS) && on_boundary", "(fabs(x[2] - 0.5) <
> DOLFIN_EPS) && on_boundary", "(fabs(x[2] + 0.5) < DOLFIN_EPS) &&
> on_boundary"])
>
> null = Constant(0.0)
> class TopDrag(Expression):
> def eval(self, vec, x):
> vec[0] = 1.0
> vec[1] = 0.0
> vec[2] = 0.0
> def dim(self):
> return 3
>
> class BottomHold(Expression):
> def eval(self, vec, x):
> vec[0] = 0.0
> vec[1] = 0.0
> vec[2] = 0.0
> def dim(self):
> return 3
>
> vOnTop = TopDrag()
> vOnBottom = BottomHold()
> pOut = 100.0 #right side pressure
> pIn = 100.0 #left side pressure
>
> bc = [DirichletBC(ME.sub(0), vOnBottom, bottom), DirichletBC(ME.sub(0),
> vOnTop, top), DirichletBC(ME.sub(0).sub(1), null, front),
> DirichletBC(ME.sub(0).sub(1), null, back), DirichletBC(ME.sub(1), pIn,
> left), DirichletBC(ME.sub(1), pOut, right)]
>
> u = Function(ME)
> v = u.split()[0]
> p = u.split()[1]
> omega = TestFunction(ME)
> w, q = split(omega)
> du = TrialFunction(ME)
> dv, dp = split(du)
>
> F = FormF(v,w,p,q)
> G = FormG(v,w,dv,p,q,dp)
> a = lhs(F-G)
> L = rhs(F-G)
> A = assemble(a)
> b = assemble(L)
> for condition in bc:
> condition.apply(A, b)
>
> du = Function(ME)
> solve(A, du.vector(), b)
>
> print 'some value:', u((0.5,0.5,0.5))
> print 'some value:', du((0.5,0.5,0.5))

Thanks Johan Hake, if I did understand correctly, You could not find anything odd by the boundary conditions definitions,

if yes, then my problem in fenics has been solved.

On Tuesday January 11 2011 07:32:21 B. Emek Abali wrote:
> Question #140910 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/140910
>
> Status: Answered => Solved
>
> B. Emek Abali confirmed that the question is solved:
> Thanks Johan Hake, if I did understand correctly, You could not find
> anything odd by the boundary conditions definitions,
>
> if yes, then my problem in fenics has been solved.

I am actually not sure. I tried to use the SubDomain::mark method to debug
what part of the boundary got marked. But could only get one correct plot.
Maybe you should strip your system down to a simple Poisson and then apply one
boundary condition at a time to check each one of them?

Johan