I can solve coupled advection problems in Stokes flow, if I use fixed point iterations, but I would like to use the builtin non-linear solver of FEniCS.
Why do I get a NaN residual in the first Newton iteration? This happens even, if I use a Constant characteristic velocity and a one-way coupling:
from dolfin import *
mesh = RectangleMesh(Point(0.,0.),Point(2.,1.),100,20)
def left(x, on_boundary):
return x[0] < 1e3*DOLFIN_EPS
def right(x, on_boundary):
return 2.-1e3*DOLFIN_EPS < x[0]
def top_bottom(x, on_boundary):
return 1.-1e3*DOLFIN_EPS < x[1] or x[1] < 1e3*DOLFIN_EPS
V = FunctionSpace(mesh,'CG',1)
W = VectorFunctionSpace(mesh,'CG',2)
Q = MixedFunctionSpace([V,W,V])
U = Function(Q)
p,v,c = split(U)
q,w,c_tst = TestFunctions(Q)
h = CellSize(mesh)
vchar = Constant(1.)
F = dot(v,grad(c))*(c_tst+dot(grad(c_tst),v)*h/vchar)*dx
F += (inner(sym(grad(v)),sym(grad(w))) + p*div(w) + q*div(v))*dx
bc1 = DirichletBC(Q.sub(1),Expression(("x[1]*(1.-x[1])","0.")),left)
bc2 = DirichletBC(Q.sub(0),Constant(0.),right)
bc3 = DirichletBC(Q.sub(1).sub(1),Constant(0.),right)
bc4 = DirichletBC(Q.sub(1),Constant((0.,0.)),top_bottom)
bc5 = DirichletBC(Q.sub(2),Expression("x[1] < 0.5"),left)
bcs = [bc1,bc2,bc3,bc4,bc5]
solve(F==0, U, bcs)
It works just fine, if I solve it decoupled !?
Q = MixedFunctionSpace([V,W])
U = Function(Q)
p,v = split(U)
q,w = TestFunctions(Q)
vchar = Constant(1.)
F = (inner(sym(grad(v)),sym(grad(w))) + p*div(w) + q*div(v))*dx
bc1 = DirichletBC(Q.sub(1),Expression(("x[1]*(1.-x[1])","0.")),left)
bc2 = DirichletBC(Q.sub(0),Constant(0.),right)
bc3 = DirichletBC(Q.sub(1).sub(1),Constant(0.),right)
bc4 = DirichletBC(Q.sub(1),Constant((0.,0.)),top_bottom)
bcs = [bc1,bc2,bc3,bc4]
solve(F==0, U, bcs)
c = Function(V)
c_tst = TestFunction(V)
F = dot(v,grad(c))*(c_tst+dot(grad(c_tst),v)*h/vchar)*dx
bcs = DirichletBC(V,Expression("x[1] < 0.5"),left)
solve(F==0, c, bcs)
