I would like to solve the Stoke's equation for pressure and velocity fields which vary sinusoidally in time: $ u(x,y,t) = \tilde{u}(x,y) \exp( i \omega t) $. Following advice on this form and in the documentation I split the complex $\tilde{u}(x,y) $ into a real and imaginary part, select Taylor-Hood elements and solve. Using a direct solver (mumps) this yields the solution of Tuck ([E.O. Tuck, J Eng Math 3 29]) for an oscillating cantilever.
I tried to then solve with an iterative solver. However, it is not clear to me how to define the preconditioner and as an engineer the references in the Fenics book (chapter 35) are hard to follow. The best I could manage so far was a shotgun approach (commented in the code listing below for the amusement of any specialists on this forum). Can anyone point me in the right direction (either a form for the preconditioner or if it exists a simpler discussion of preconditioning)?
As a final comment using a popular but very expensive commercial tool this model did solve using a ilu precondition and gmres (fortunately for simple folk like me the user doesn't have to specify a preconditioner for this tool) so I guess the situation is not intrinsically hopeless. However, I ultimately want to use dolfin-adjoint so this route is not particularly attractive...
In any event here is my model:
import numpy as np
from dolfin import *
import mshr as mr
#--- variables
R_L = 50.
Ro = 10
Re=10
t = 0.02
a = 1
radius = R_L/np.sqrt(Re)
#--- define geometry (use polygon for bc later)
cantilever = mr.Polygon([Point(-a,-t/2),Point(a,-t/2),Point(a,t/2),Point(0,t/2)\
,Point(-a,t/2)])
ob = mr.Circle(Point(0,0),radius,segments=100)
domain = ob-cantilever
mesh = mr.generate_mesh(domain,20)
#--- mark boundaries
class outer_boundary(SubDomain):
def inside(self,x,on_boundary):
dist = 0
if on_boundary:
dist = np.abs(np.sqrt(x[0]**2+x[1]**2)-radius)
return on_boundary and dist<0.1
ob = outer_boundary()
class canti_boundary(SubDomain):
def inside(self,x,on_boundary):
dist = 0
if on_boundary:
dist = np.abs(np.sqrt(x[0]**2+x[1]**2)-radius)
return on_boundary and not dist <0.1
cb =canti_boundary()
bnd_mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bnd_mf.set_all(0)
ob.mark(bnd_mf,1)
cb.mark(bnd_mf,2)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2, dim=4) # Trial: pair of velocity fields, Test: Cauchy Tensor
Q = VectorFunctionSpace(mesh, "CG", 1, dim=2) # Trial: pair of pressure fields, Test: Divergence condition
W = V * Q
S = FunctionSpace(mesh,'CG',1) # for post only
VV = VectorFunctionSpace(mesh,'CG',1) # for post only
# No-slip boundary condition for velocity
noslip = Constant((0., 0., 0., 0.))
osc = Expression(("0", "1","0.", "0.0"))
freeS_p = Constant((0,0))
def is_canti_center(x,on_boundary):
return near(x[0],0) and near(x[1],t/2)
bc1 = DirichletBC(W.sub(0), noslip, bnd_mf, 1)
bc2 = DirichletBC(W.sub(0), osc, bnd_mf, 2)
bcp = DirichletBC(W.sub(1),freeS_p,is_canti_center,'pointwise')
# Collect boundary conditions
bcs = [bc1,bc2,bcp]
# Define variational problem
(un, pn) = TrialFunctions(W)
(vn, qn) = TestFunctions(W)
ur = as_vector([un[0],un[1]])
ui = as_vector([un[2],un[3]])
pr = pn[0]
pi = pn[1]
vr = as_vector([vn[0],vn[1]])
vi = as_vector([vn[2],vn[3]])
qr = qn[0]
qi = qn[1]
fr = Constant((0,0))
fi = Constant((0,0))
(j,l) = indices(2)
a = Dx(ur[l],j)*Dx(vr[l],j)*dx-Dx(vr[j],j)*pr*dx + qr*Dx(ur[j],j)*dx+\
Dx(ui[l],j)*Dx(vi[l],j)*dx-Dx(vi[j],j)*pi*dx + qi*Dx(ui[j],j)*dx-\
Constant(Re)*ui[j]*vr[j]*dx+Constant(Re)*ur[j]*vi[j]*dx
L = fr[j]*vr[j]*dx+fi[j]*vi[j]*dx
###---- Solution with Direct Solver
w = Function(W)
solve(a==L,w,bcs=bcs,solver_parameters={"linear_solver": "mumps"})
###/---- Solution with Direct Solver
###---- Solution with Iterative Solver ??
#A,bb = assemble_system(a,L,bcs)
#precond = (inner(nabla_grad(ur), nabla_grad(vr))+ inner(nabla_grad(ui), nabla_grad(vi)))*dx\
# + (dot(ur,vr)+dot(ui,vi)) * dx\
# + (pr*qr+pi*qi) * dx\
# + (dot(grad(pr),grad(qr))+dot(grad(pi),grad(qi))) * dx
#P, btmp = assemble_system(precond, L, bcs)
#w = Function(W)
#solver = KrylovSolver("gmres", "ilu")
#solver.set_operators(A, P)
#solver.parameters['monitor_convergence'] = True
#solver.parameters['maximum_iterations'] = 1000
#solver.solve(w.vector(), bb)
###/---- Solution with Iterative Solver
# # Split the mixed solution using a shallow copy
(un,pn) = w.split()
ur = as_vector([un[0],un[1]])
ui = as_vector([un[2],un[3]])
pr = pn[0]
pi = pn[1]
# Plot solution: in phase and quadrature compents
plot(ur,title='u: llel V')
plot(ui,title='u: llel A')
plot(pr,title='p: llel V')
plot(pi,title='p: llel A')
#-- post process
uxr = project(ur[0],S)
uyr = project(ur[1],S)
uxt = Function(S)
uxi = project(ui[0],S)
uyi = project(ui[1],S)
uyt = Function(S)
ut = Function(VV)
interactive()
