Unusual Eigenvalues in for EM wave in Transverse Cylindrical waveguide

Question

Hello,

I am solving the eigenvalue problem for a cylindrical waveguide. I followed the reference code from chapter 34. However, I can't find the right eigenvalue, as the field has many spikes. I used the following parameters in the SlepcEigenSolver() class to make sure it's close to the target k_0:

#k_0 = 2*pi/1.55 is the wavevector (lamba = 1.55 microns = c0/f)
eigensolver = SLEPcEigenSolver(A,B)
eigensolver.parameters["spectral_shift"] = ((2*pi/1.55)**2)
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectrum"] = "target real" 

I am having a tough time trying to understand what is wrong. Any sort of guidance towards finding the right eigenvalues would appreciated.

Thanks a lot.

For reference, here is the code for reference:

from dolfin import *

# Test for PETSc and SLEPc
...

# Define mesh, function space

transverse_order = 2
axial_order = 1

# Importing GMSH file for single mode profile
mesh = Mesh("SMF.xml")
Nl = FunctionSpace(mesh, "N1curl", transverse_order)
Vl = FunctionSpace(mesh, "CG", axial_order)
combined_space = Nl * Vl

# Define basis and bilinear form
(N_i, L_i) = TestFunctions(combined_space)
(N_j, L_j) = TrialFunctions(combined_space)

# Define Dielectric constant for core radius = 1 micron
...
e_r = Dielectric()

# Relative permiability  (=1)
one_over_u_r = Constant(1.0)

# Define wave vector (in micron^(-1) units), wavelength of light = 1.550 microns

k_o_squared = Constant((2*pi/1.55)**2)

# Define Zero Boundary condition
class Zero_Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

zero = Constant ( (0.0, 0.0, 0.0) )

dirichlet_bc = DirichletBC(combined_space, zero, Zero_Boundary())

# 2D curl function
def curl_t(w):
    return Dx(w[1],0) - Dx(w[0],1)

#Define Matrices A and B such that A*e = lambda*B*e (eigenvalue: lambda)

s_tt = one_over_u_r*dot(curl_t(N_i), curl_t(N_j))
t_tt = e_r*dot(N_i, N_j)

s_zz = one_over_u_r*dot(grad(L_i), grad(L_j))
t_zz = e_r*L_i*L_j

b_tt = one_over_u_r*dot(N_i, N_j)
b_tz = one_over_u_r*dot(N_i, grad(L_j))
b_zt = one_over_u_r*dot(grad(L_i), N_j)

# define the matrix entries a_ij, b_ij with * dx resulting in integration over 
# the domain of the mesh when assemble() is called
a = ( s_tt - k_o_squared * t_tt ) * dx
b = ( s_zz - k_o_squared * t_zz + b_tt + b_tz + b_zt ) * dx

# Assembling the matrices
A = PETScMatrix()
B = PETScMatrix()
assemble(a, tensor=A)
assemble(b, tensor=B)

# Applying boundary condition on the Matrices A,B
dirichlet_bc.apply(A)
dirichlet_bc.apply(B)

# Create eigensolver with generalized Hermitian condition

eigensolver = SLEPcEigenSolver(A,B)
eigensolver.parameters["spectral_shift"] = ((2*pi/1.55)**2)
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectrum"] = "target real"


# Compute all eigenvalues of A x = \lambda B x
print "Computing eigenvalues. This can take a minute."
eigensolver.solve()
print "Eigenvalues computed: ",eigensolver.get_number_converged()

if eigensolver.get_number_converged() > 0:

    r, c, rx, cx = eigensolver.get_eigenpair(0)

    print "Eigenvalue : ",r

# Initialize function and assign eigenvector
    u = Function(combined_space,rx)
    transverse, axial = u.split()

# Plot eigenfunction
    plot(abs(axial))
    interactive()