I have a simple Stokes equation with a weak boundary condition on the pressure on one frontier. The code runs and results are satisfying when using a direct solver but when using an iterative solver (I took the example demo_stokes-iterative.py that runs on many other geometries) it does not converge anymore.
mesh = Mesh("Mesh2D.xml")
File("maillage2D.pvd")<<mesh
Define function spaces
V1 = FiniteElement("Lagrange",mesh.ufl_cell(),2)
V2 = FiniteElement("Lagrange",mesh.ufl_cell(),2)
P = FiniteElement("Lagrange",mesh.ufl_cell(),1)
M = MixedElement([[V1,V2], P])
W = FunctionSpace(mesh, M)
Np = 10
Lp = 2.5
Rp = 0.5
Lr = 8
Rr = 4
P0 = -1.
nu=1.
tol = 1e-6
Boundaries
def b1(x, on_boundary): return on_boundary and near(x[1], -Rp, tol)
def b2(x, on_boundary): return on_boundary and near(x[0], Lp, tol) and (x[1]<=-Rp)
def b3(x,on_boundary): return on_boundary and near(x[1], -Rr,tol)
def b4(x, on_boundary): return on_boundary and near(x[0], Lp+Lr, tol)
def b5(x,on_boundary): return on_boundary and near(x[1], Rr, tol)
def b6(x,on_boundary): return on_boundary and near(x[0],Lp,tol) and (x[1]>=Rp)
def b7(x, on_boundary): return on_boundary and near(x[1], Rp, tol)
def b8(x,on_boundary): return on_boundary and near (x[0], 0, tol)
def bound(x,on_boundary): return on_boundary
No-slip boundary condition for velocity
noslip = project(Constant((0.0, 0.0)), W.sub(0).collapse())
bc1 = DirichletBC(W.sub(0), noslip, b1)
bc7 = DirichletBC(W.sub(0), noslip, b7)
slip boundary condition for velocity u.n = 0
slip = project(Constant(0.), W.sub(0).sub(0).collapse())
bc2 = DirichletBC(W.sub(0).sub(0), slip, b2)
bc6 = DirichletBC(W.sub(0).sub(0),slip, b6)
collect boundary conditions
bcs = [bc1,bc2,bc6,bc7]
class fsi_boundary(SubDomain):
def inside(self, x,on_boundary):
return on_boundary and near(x[0], 0, tol)
define measure dsS
fsimarkers = FacetFunction("size_t", mesh)
fsimarkers.set_all(0)
fsibound = fsi_boundary()
fsibound.mark(fsimarkers, 8)
dsS = Measure('ds', domain=mesh, subdomain_data = fsimarkers)
Define variational problem
normal = project(Constant((1.,0.)), W.sub(0).collapse())
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
w= Function(W)
f = Constant((0.0, 0.0))
a = nuinner(sym(grad(u)), sym(grad(v)))dx - pdiv(v)dx + qdiv(u)dx
L = inner(f, v)dx + Constant(P0) dot(v,normal)*dsS(8)
Form for use in constructing preconditioner matrix
b = inner(sym(grad(u)), sym(grad(v)))dx +pq*dx
Assemble system
A, bb = assemble_system(a, L, bcs)
Assemble preconditioner system
P, btmp = assemble_system(b, L, bcs)
Create Krylov solver and AMG preconditioner
solver = KrylovSolver(krylov_method, "amg")
Associate operator (A) and preconditioner matrix (P)
solver.set_operators(A, P)
Solve
solver.solve(w.vector(), bb)
Get sub-functions
u, p = w.split(True)
