Thanks for the suggestion. I forgot to mention that my eigenvalue problem is complex. I was too scared to follow your advice, since I cannot even pronounce sesquilinear, let alone understanding what it means. In the end my colleague who wrote the code just followed the recommendation from this post.
Below is a (not so minimal) example. Is it straight forward to solve it with SLEPc? The usual problems that we solve are more complicated, so I do not want to make any violent changes to our well tested code.
from dolfin import * # overwrites local variables with name conflict!
import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
from scipy.sparse.linalg import eigs
Re = Constant(896.425)
eta = 1.135
d = eta-1.0
Nr = 20
velocity_order = 2
kz = Constant(1/2/d)
k2pi = kz*2*np.pi
num_eig = 10
class InnerCylinder(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],1.0)
class OuterCylinder(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],eta)
mesh = IntervalMesh(Nr,1.0,eta)
x = SpatialCoordinate(mesh)
r = x[0]
# Function space
V0 = FunctionSpace(mesh, "CG", velocity_order )
Q0 = FunctionSpace(mesh, "CG", velocity_order-1)
Vp = VectorElement("CG", mesh.ufl_cell(), velocity_order , dim=3)
Qp = FiniteElement("CG", mesh.ufl_cell(), velocity_order-1)
Wp = FunctionSpace(mesh, MixedElement(Vp,Qp))
up, pp = TrialFunctions(Wp)
vp, qp = TestFunctions(Wp)
ur,uphi,uz = split(up) # here the order of velocity components is more-less arbitrary
vr,vphi,vz = split(vp) # but for 2D and 3D meshes it must match the order of coordinates in FEniCS !
# Analytical basic state
c1,c2 = -1/(eta**2-1), eta**2/(eta**2-1)
Uphi_exact = Expression('c1*x[0]+c2/x[0]',c1=c1,c2=c2,degree=2)
P_exact = Expression('pow(A,2)*pow(x[0],2)/2+2*A*B*std::log(x[0])-pow(B,2)/2/pow(x[0],2)-(pow(A,2)-pow(B,2))/2',
A=c1,B=c2,degree=1)
Uphi,P = Function(V0), Function(Q0)
Uphi.interpolate(Uphi_exact)
P.interpolate(P_exact)
# Boundary conditions
bcs = [DirichletBC(Wp.sub(0), Constant((0.0,0.0,0.0)), OuterCylinder())
,DirichletBC(Wp.sub(0), Constant((0.0,0.0,0.0)), InnerCylinder())
]
# Variational forms
Ji = -vz*k2pi*pp *r*dx # z-momentum: pressure gradient
Ji+= -uz*k2pi*qp *r*dx # continuity, z-direction
J = ( 2*Uphi*uphi*vr/r +pp*(vr.dx(0)+vr/r) ) *r*dx # r: convection + pressure
J += -vphi*ur *( Uphi.dx(0) +Uphi/r ) *r*dx # phi: convection
J += 1/Re *( -vr.dx(0)*ur.dx(0) -k2pi**2*ur*vr -vr*ur/r**2 ) *r*dx # r
J += 1/Re *( -vphi.dx(0)*uphi.dx(0) -k2pi**2*uphi*vphi -vphi*uphi/r**2 ) *r*dx #phi
J += 1/Re *( -vz.dx(0)*uz.dx(0) -k2pi**2*uz*vz ) *r*dx # z: Laplacian
J += -qp *( ur.dx(0) +ur/r ) *r*dx # continuity
m = dot(up,vp) *r*dx
# Create matrices in FEniCS
A, Ai, M = PETScMatrix(), PETScMatrix(), PETScMatrix() # initialize empty matrices in FEniCS
assemble(J , tensor=A ) # populate them with coefficients from the weak forms
assemble(Ji, tensor=Ai)
assemble(m , tensor=M )
# Impose Dirichlet BCs
bcinds = []
for bc in bcs:
bc.apply(A )
bc.apply(Ai)
bc.apply(M )
bcdict = bc.get_boundary_values()
bcinds.extend(bcdict.keys())
# Export matrices from FEniCS into sparse matrices for SciPy
Ar = sp.csr_matrix( A.mat().getValuesCSR()[::-1])
Ai = sp.csr_matrix(Ai.mat().getValuesCSR()[::-1])
M = sp.csr_matrix( M.mat().getValuesCSR()[::-1])
# Create shift matrix -> take care of the dirichlet boundary conditions
# otherwise one gets ficticious eigenvalues equal to 1
shift = 1.2345e6*np.ones(len(bcinds)) # large arbitrary multiplier
S = sp.csr_matrix((shift, (bcinds, bcinds)), shape=Ar.shape)
# Complex matrix
A = Ar + 1j*Ai + S
# EGV solver
v, V = eigs(A, num_eig, M, sigma=1.0)
# order eigenvalues from largest to smallest real part
idxs = np.argsort(-np.real(v))
v,V = v[idxs], V[:,idxs[0]]
print(v[0])