Boundary conditions or other problems


I’m trying to solve a 2d Helmholtz scattering equation for its resonance values for n=2 and other values:
[tex]\Delta u_1+k^2(1+n)u_1=0 \in D[/tex]
[tex]\Delta u_2+k^2u_2=0 \in \mathcal{R}^2 \D[/tex]
with Sommerfeld radiation condition and equality of the u’s on the boundary of D and their partials with respect to the normal.

This turns into a quadratic eigenvalue problem and I use the following code to try and solve for k.

My code is:

#Import needed packages and programs
import netgen.gui
from netgen.geom2d import SplineGeometry
from ngsolve import *
import matplotlib.pyplot as plt
import netgen.geom2d as geom2d
from netgen.geom2d import unit_square
import numpy as np
import scipy.sparse as scs
from scipy.sparse.linalg import eigs, eigsh

#Create geometries, first circle is the scatterer, second begins the PML third truncates the infinite domain
#third circle (largest) gets a boundary condition due to the sommerfeld radiation condition. 
geo = SplineGeometry()
geo.AddCircle( (0.0, 0.0), r=1.0, leftdomain=1, rightdomain=2, bc = "scatterer")
geo.AddCircle ( (0.0, 0.0), r=1.5, leftdomain=2, rightdomain=3)
geo.AddCircle ( (0.0, 0.0), r=2.0, leftdomain=3, rightdomain=0,bc = "outerbnd")

#Set material 1 to be inner do that only those nodes inside the scatterer recieve the $n*u*v*dx("inner")$ term
geo.SetMaterial(1, "inner")
geo.SetMaterial(3, "PML")

#Generate the mesh 
ngmesh = geo.GenerateMesh(maxh=0.04)
mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle
mesh.SetPML(pml.Radial(rad=1.5,alpha=1j,origin=(0,0)), "PML")

#Set the extra value of the scatterer

#I think this tells the program what level of accuracy and where our test/trial functions live
fes = H1(mesh, order=4, complex=True, dirichlet="dir")

#create our trial function u and our test function v
u = fes.TrialFunction()
v = fes.TestFunction()

#create matrix K that is the non lambda term (sign for quad-eig)
k = BilinearForm (fes, symmetric=True)
k += grad(u)*grad(v)*dx

#create matrix C that is our lambda term (sign for quad-eig)
c = BilinearForm (fes, symmetric=True)
c += 1.j*u*v*ds('outerbnd')

#create matrix M which is our squared term and has the nn element (sign for quad-eig)
m = BilinearForm (fes, symmetric=True)
m += u*v*dx
m += nn*u*v*dx('inner')

#convert bilinearform into a matrix form

#tells the program to get ready
u = GridFunction(fes, multidim=20, name='resonances')

#moves the matrices from ngsolve into scipy.sparse because ngsolve hates stacking things
mcsr = scs.csr_matrix(m.mat.CSR(), shape=(m.mat.height, m.mat.width),dtype='complex')
ccsr = scs.csr_matrix(c.mat.CSR(), shape=(c.mat.height, c.mat.width),dtype='complex')
kcsr = scs.csr_matrix(k.mat.CSR(), shape=(k.mat.height, k.mat.width),dtype='complex')

#finds the size of M and creates the needed identity matrix of that size
leng = m.mat.height
I = scs.csr_matrix(scs.eye(leng), dtype='complex')

#Builds our first block matrix in csc form
B1 = scs.csr_matrix(scs.bmat([[mcsr, None], [None, I]], dtype='complex'))

#Builds our second block matrix in csc form 
B2 = scs.csr_matrix(scs.bmat([[ccsr, kcsr],[I, None]], dtype='complex'))

omegaval, omegavec = eigs(B2, k=11, M=B1, sigma=2.2-0.2*1.j)

for ome in omegaval:

This results in the following output results:

The problem is I can solve for the resonances using Hankel and Bessel functions and a Newton solver to get that a resonance values should be at 2.284220485273998-0.3840888660574083j. I have verified this will a friend and a Muller’s method solver.

First question is does the bilinear form account for the equality along the boundary (in both u and its derivative wrt the normal) or is there something special I need to do in the code

Second, does anyone have an idea as to how to fix the discrepancy?