Hi everybody
I am working on a problem of instability analysis but I am having a weird problem that I use as a very fine mesh and it only works if I make a very great reduction that works.
thanks for your incoming answers
from dolfin import *
from mshr import *
# Constants related to the geometry
xmin = -10.; xmax = 40.
ymin = -10.0; ymax = 10.
xcenter = 0.0; ycenter = 0.0; radius = 0.5
# create domain mesh
domain = Rectangle(Point(xmin,ymin), Point(xmax,ymax)) - Circle(Point(0.0,0.0),0.5)
domain.set_subdomain(1, Rectangle(Point(xmin,-3.), Point(xmax,3.)))
domain.set_subdomain(2, Rectangle(Point(-2.,-1.5), Point(7., 1.5)))
mesh = generate_mesh(domain, 50)
# mesh = refine(mesh)
# refine mesh TWICE around the cylinder
for refinements in [0,1]:
cell_markers = CellFunction("bool", mesh)
cell_markers.set_all(False)
for cell in cells(mesh):
p = cell.midpoint()
if (xmin < p[0] < xmax) and \
(-2.75 < p[1] < 2.75 ) :
cell_markers[cell] = True
mesh = refine(mesh, cell_markers)
# refine mesh tree times around the cylinder
for refinements in [0,1]:
cell_markers = CellFunction("bool", mesh)
cell_markers.set_all(False)
for cell in cells(mesh):
p = cell.midpoint()
if (-2. < p[0] < 7.) and \
(-1.25 < p[1] < 1.25 ) :
cell_markers[cell] = True
mesh = refine(mesh, cell_markers)
File("mesh.xml.gz") << mesh
# Plot the mesh
plot(mesh)
interactive()
the eigenvalue code:
from dolfin import *
# Constants related to the geometry
bmarg = 1.e-3 + DOLFIN_EPS
xmin = -10; xmax = 40
ymin = -10.0; ymax = 10
xcenter = 0.0; ycenter = 0.0; radius = 0.5
mesh = Mesh("mesh.xml.gz")
plot(mesh, title = "Mesh domain")
#boundary conditions using a mesh function
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
# Inflow boundary
class InflowBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] < xmin + bmarg
# No-slip boundary
class NoslipBoundary(SubDomain):
def inside(self, x, on_boundary):
dx = x[0] - xcenter
dy = x[1] - ycenter
r = sqrt(dxdx + dydy)
return on_boundary and r < radius + bmarg
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V*Q
# no-slip velocity b.c
noslipBoundary = NoslipBoundary()
g0 = Constant( (0., 0.) )
bc0 = DirichletBC(W.sub(0), g0, noslipBoundary)
# inlet velocity b.c
inflowBoundary = InflowBoundary()
g1 = Constant( (1., 0.) )
bc1 = DirichletBC(W.sub(0), g1, inflowBoundary)
# collect b.c.
bcs = [bc0, bc1]
# define function
u = Function(V)
p = Function(Q)
# Define variational problem for initial guess (stokes solution)
v, q = TestFunctions(W)
u, p = TrialFunctions(W)
f = Constant( (0.0, 0.0) )
F = (inner(sym(grad(u)), grad(v)) - pdiv(v) + qdiv(u) )dx + \
inner(v, f)dx
a = lhs(F); L = rhs(F)
# Assemble system
A, b = assemble_system(a, L, bcs)
up = Function(W)
solve(A, up.vector(), b, 'lu')
u_a, p_a = up.split()
plot(p_a, title = "initial guess pressure")
plot(u_a, title = "initial guess velocity")
# Save solution to file
file = File("initial_guess.pvd")
file << up
# Boundary conditions
# inlet velocity b.c
g11 = Constant( (0., 0.) )
bc11 = DirichletBC(W.sub(0), g11, inflowBoundary)
# collect b.c.
bcs1 = [bc0, bc11]
Re = 40.0
Remax = 81.0
while Re < Remax:
# define function
du = Function(V)
dp = Function(Q)
# Define variational problem for convergent solution up1
du, dp = TrialFunctions(W)
F1 = (inner(grad(du)*u_a, v) + inner(grad(u_a)*du, v) )*dx + \
(1/Re)*inner(grad(du), grad(v))*dx + \
(- dp*div(v) + q*div(du))*dx + \
(inner(grad(u_a)*u_a, v))*dx + (1/Re)*inner(grad(u_a), grad(v))*dx + \
(- p_a*div(v) + q*div(u_a))*dx
a1 = lhs(F1); L1 = rhs(F1)
# Newton iteration
err = 1.0
tol = 2.e-15
iter = 0
maxiter = 10
up1 = Function(W)
u1, p1 = up1.split()
while err > tol and iter < maxiter:
A1, b1 = assemble_system(a1, L1, bcs1)
solve(A1, up1.vector(), b1, 'lu')
err = norm(up1, 'L2')/25200
#print " iteration=%g, res=%g", iter, err
iter += 1
up.vector()[:] += up1.vector()
# save
File("base_flow.pvd") << up
File("velocity.pvd") << u_a
File("pressure.pvd") << p_a
U0, P0 = up.split() # base flow
############ normal modes problems u = U0 + du1exp(iwt), p = P0 + dp1exp(-iwt)#####################
u2, p2 = TrialFunctions(W)
a2 = (inner(grad(u2)*U0, v)+ inner(grad(U0)*u2, v) + (1/Re)*inner(grad(u2), grad(v)) - p2*div(v) + q*div(u2) )*dx
m = - inner(u2, v)*dx
A0 = PETScMatrix()
M = PETScMatrix()
assemble(a2, tensor=A0)
assemble(m, tensor=M)
for bc in bcs:
(bc.apply(A0) and bc.apply(M))
#create eigensolver
eigensolver = SLEPcEigenSolver(A0, M)
# eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectral_shift"] = 4.0
eigensolver.parameters['verbose'] = True
print 'solving'
num = 2
eigensolver.solve(num)
i = 1
for i in range(num):
#extract next eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
print("Reynolds:, Eigenvalue:", Re, i, r, c/(2*pi))
Re += 1.0
interactive()
