I am trying to solve the vectorial helmholtz equation in a cylindrical inhomogenous waveguide, with a geometry and olverlayn permittivity shown in Fig. (a). The variational formulation for the problem with PEC boundary conditions and E split into a transverse and z term is shown in (b) and its implementation for Fenics is based on the work from Lezar. A eigenvalue spectrum is shown in ©, with the horiztonal line showing the maximum physically possible. I guess the higher values belong to spurious modes, fine. (d) shows the x-component of a mode. For some mesh parameters I get solutions that make sense (as given in the full code below). However, when making the mesh finer by adjusting

mesh = generate_mesh(domain, 15)

more spurious modes kick in and my previous solutions are overlayn with random new features (e). I haven’t seen this kind of behaviour in COMSOL before. Have I formulated an ill-posed problem maybe? Using gmsh to refine the mesh only around the core also did not help. A ‘minimal’ example is given below.

This has been bugging me for a quite a while, so if someone knows a way out I would appreciate it.

Steffen

```
from dolfin import *
import numpy as np
import fenics as fn
from mshr import *
class set_properties(UserExpression):
def __init__(self, subdomains, param, **kwargs):
super().__init__(**kwargs)
self.subdomains = subdomains
self.params = param
def eval_cell(self, values, x, cell):
values[0] = self.params[self.subdomains[cell.index]-1]
class OperatingWavenumberSquared(UserExpression):
def eval (self, values, x):
values[0] = self.__k_o_squared
def set_wavelength(self, wavelength):
self.__k_o_squared = (2*np.pi/wavelength)**2
def curl_t(w):
return Dx(w[1], 0) - Dx(w[0], 1)
vector_order = 2
nodal_order = 3
wavelength = 1.55
k0 = 2*pi/wavelength
r_clad = 40
r_core = 8
n_clad = 1.4378
n_core = 1.4457
domain = Circle(Point(0, 0), r_clad)
core = Circle(Point(0, 0), r_core)
domain.set_subdomain(1, core)
mesh = generate_mesh(domain, 15)
surface_markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
e_rs = np.array([n_core**2, n_clad**1])
e_r = set_properties(surface_markers, e_rs, degree=2)
k_0_squared = OperatingWavenumberSquared()
k_0_squared.set_wavelength(wavelength)
>
vector_element = FiniteElement("N1curl", mesh.ufl_cell(), vector_order)
nodal_element = FiniteElement("Lagrange", mesh.ufl_cell(), nodal_order)
V = FunctionSpace(mesh, vector_element * nodal_element) # NOTE: V.sub(0) and V.sub(1) are the subspaces
(N_i, L_i) = TestFunctions(V)
(N_j, L_j) = TrialFunctions(V)
u_r = 1.0
s_tt_ij = 1.0/u_r * dot(curl_t(N_i), curl_t(N_j))
t_tt_ij = e_r * dot(N_i, N_j)
s_zz_ij = 1.0/u_r * dot(grad(L_i), grad(L_j))
t_zz_ij = e_r * L_i * L_j
b_tt_ij = 1.0/u_r * dot(N_i, N_j)
b_tz_ij = 1.0/u_r * dot(N_i, grad(L_j))
b_zt_ij = 1.0/u_r * dot(grad(L_i), N_j)
a_ij = ( s_tt_ij - k_0_squared * t_tt_ij ) * dx
b_ij = ( s_zz_ij - k_0_squared * t_zz_ij + b_tt_ij + b_tz_ij + b_zt_ij ) * dx
electric_wall = DirichletBC(V, Constant((0.0, 0.0, 0.0)), DomainBoundary())
parameters["linear_algebra_backend"] = "PETSc"
A = PETScMatrix()
B = PETScMatrix()
assemble(a_ij, tensor=A)
assemble(b_ij, tensor=B)
electric_wall.apply(B)
electric_wall.apply(A)
neigs = 50
disp = np.zeros(neigs)
# Create VTK file for saving solution
vtkfile = File('solved_modes/solution.pvd')
solver = SLEPcEigenSolver(A, B)
solver.parameters['tolerance'] = 1e-12
solver.parameters["spectrum"] = "target magnitude"
solver.parameters["spectral_transform"] = "shift-and-invert"
solver.parameters["spectral_shift"] = -n_core**2*k0**2
solver.parameters["solver"] = "arnoldi"
solver.solve(neigs)
assert solver.get_number_converged() > 0
for j in range(min(neigs, solver.get_number_converged())):
r, c, rx, cx = solver.get_eigenpair(j)
disp[j] = np.sqrt(np.abs(r/k0**2))
mode = Function(V)
mode.vector()[:] = rx
(Et, Ez) = mode.split()
Et.rename("Et", "Transverse mode solution")
Ez.rename("Ez", "z-axis mode solution")
vtkfile << (Et, j)
vtkfile << (Ez, j)
```