Hello,

I’m running some resonance experiments with a base problem and a perturbed problem. I calculate the resonance in the base case (first block of code at the bottom) and then run the perturbed case in block 2 and 4 changing the choices of epsilon. My issue is I know from a 1d example that I need epsilon in the range of 1/20-1/40 to obtain the numerical results that coincide with the analytical results that I already have. The idea is that the epsilon cases are perturbed resonance values that converge to the base resonance as epsilon goes to 0 along the path delta.

The issue is I somehow need to strengthen the code to handle the smaller epsilons but I am not sure how. I have tried expanding the PML final layer further out but that does not work. I have also tried to shrink the layers and shrink the max mesh size but have issue obtaining the base resonance (which is not perfect but close to the true values, verified through alternate means). Any help would be appreciated. I tried attaching a jupyter notebook but it would not attach. As such, I have the code posted below.

Thank you in advance.

Block 1

[code]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

#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.AddRectangle( (0.0, 0.0), (1.0, 1.0), leftdomain=1, rightdomain=2, bc = “scatterer”)

geo.AddCircle ( (0.5, 0.5), r=12.0, leftdomain=2, rightdomain=3)

geo.AddCircle ( (0.5, 0.5), r=35.0, leftdomain=3, rightdomain=0,bc = “dir”)

#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.1)

mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle

mesh.SetPML(pml.Radial(rad=12.0, alpha=1j, origin=(0.5,0.5)), “PML”)

#Set the extra value of the scatterer

nn=10;

#I think this tells the program what level of accuracy and where our test/trial functions live

fes = H1(mesh, order=3, complex=True, dirichlet=“dir”)

#create our trial function u and our test function v they are H1 functions

u = fes.TrialFunction()

v = fes.TestFunction()

#create matrix K that is the non lambda term (sign for quad-eig)

a = BilinearForm (fes, symmetric=True)

a += grad(u)*grad(v)*dx

#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

a.Assemble()

m.Assemble()

u = GridFunction(fes, multidim=1, name=‘resonances’)

with TaskManager():

lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),

u.vecs, shift=0.1502419606698783-0.19383557950824307*1.j)

Draw(u)

print ("lam: ", lam)

tplam=lam[0]

print(tplam)[/code]

Block 2:

[code]siz=5;

tpepsvec = np.zeros(siz);

tplamvec = np.zeros(siz, dtype = complex);

for i in range(siz):

eps = 1/(15+i+0.05);

tpepsvec[i]=eps;

#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.AddRectangle( (0.0, 0.0), (1.0, 1.0), leftdomain=1, rightdomain=2, bc = “scatterer”)

geo.AddCircle ( (0.5, 0.5), r=12.0, leftdomain=2, rightdomain=3)

geo.AddCircle ( (0.5, 0.5), r=35.0, leftdomain=3, rightdomain=0,bc = “dir”)

#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.1)

mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle

mesh.SetPML(pml.Radial(rad=12.0, alpha=1j, origin=(0.5,0.5)), “PML”)

#Set the extra value of the scatterer

#Set the extra value of the scatterer

f = 10+2*np.pi*sin(2*np.pi*x/eps)+2*np.pi*cos(2*np.pi*y/eps)

#I think this tells the program what level of accuracy and where our test/trial functions live

fes = H1(mesh, order=3, complex=True, dirichlet=“dir”)

#create our trial function u and our test function v they are H1 functions

u = fes.TrialFunction()

v = fes.TestFunction()

#create matrix K that is the non lambda term (sign for quad-eig)

a = BilinearForm (fes, symmetric=True)

a += grad(u)*grad(v)*dx

#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 += f*u*v*dx(‘inner’)

#convert bilinearform into a matrix form

a.Assemble()

m.Assemble()

```
u = GridFunction(fes, multidim=5, name='resonances')
with TaskManager():
lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),
u.vecs, shift=tplam)
for l in lam:
if abs(l-(tplam))<0.05:
tplamvec[i]=lam[0]
```

print(tplamvec) [/i][/i][/code]

[i][i]Block 3:

[code]N=5-1;

sq = np.zeros(N+1);

for j in range(N+1):

sq[j]=np.abs(tplamvec[j]-tplam);

print(sq)

for i in range(N+1):

ax = plt.gca()

ax.scatter(tpepsvec[i] ,sq[i] , c=‘blue’);

ax.set_yscale(‘log’);

ax.set_xscale(‘log’);

m, b = np.polyfit(np.log(tpepsvec), np.log(sq), 1);

plt.plot(tpepsvec, np.exp(b)*tpepsvec**m)

plt.title(‘\delta = 0.05, on loglog scale’)

plt.xlabel(‘Epsilon’)

plt.ylabel('Resonance values ')

plt.show()

print(‘The slope is:’, m)[/i][/i][/code]

The slope output is [/i][/i]0.29766692963062147[i][i][i][i]Block 4:

[code]siz=5;

tepsvec = np.zeros(siz);

tlamvec = np.zeros(siz, dtype = complex);

for i in range(siz):

eps = 1/(2+i+0.05);

tepsvec[i]=eps;

#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.AddRectangle( (0.0, 0.0), (1.0, 1.0), leftdomain=1, rightdomain=2, bc = “scatterer”)

geo.AddCircle ( (0.5, 0.5), r=12.0, leftdomain=2, rightdomain=3)

geo.AddCircle ( (0.5, 0.5), r=35.0, leftdomain=3, rightdomain=0,bc = “dir”)

#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.1)

mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle

mesh.SetPML(pml.Radial(rad=12.0, alpha=1j, origin=(0.5,0.5)), “PML”)

#Set the extra value of the scatterer

#Set the extra value of the scatterer

f = 10+2*np.pi*sin(2*np.pi*x/eps)+2*np.pi*cos(2*np.pi*y/eps)

#I think this tells the program what level of accuracy and where our test/trial functions live

fes = H1(mesh, order=3, complex=True, dirichlet=“dir”)

#create our trial function u and our test function v they are H1 functions

u = fes.TrialFunction()

v = fes.TestFunction()

#create matrix K that is the non lambda term (sign for quad-eig)

a = BilinearForm (fes, symmetric=True)

a += grad(u)*grad(v)*dx

#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 += f*u*v*dx(‘inner’)

#convert bilinearform into a matrix form

a.Assemble()

m.Assemble()

```
u = GridFunction(fes, multidim=5, name='resonances')
with TaskManager():
lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),
u.vecs, shift=tplam)
for l in lam:
if abs(l-(tplam))<0.05:
tlamvec[i]=lam[0]
```

print(tlamvec)

N=5-1;

sq = np.zeros(N+1);

for j in range(N+1):

sq[j]=np.abs(tlamvec[j]-tplam);

print(sq)

for i in range(N+1):

ax = plt.gca()

ax.scatter(tepsvec[i] ,sq[i] , c=‘blue’);

ax.set_yscale(‘log’);

ax.set_xscale(‘log’);

m, b = np.polyfit(np.log(tepsvec), np.log(sq), 1);

plt.plot(tepsvec, np.exp(b)*tepsvec**m)

plt.title(‘\delta = 0.05, on loglog scale’)

plt.xlabel(‘Epsilon’)

plt.ylabel('Resonance values ')

plt.show()

print(‘The slope is:’, m)

This has slope: [/i][/i][/i][/i][/code][/i][/i][/i][/i]0.9212055445929848