Boundary conditions or other problems

More
3 years 1 month ago #3896 by ajf367
Hi

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:
Code:
#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 nn=2; #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 k.Assemble() c.Assemble() m.Assemble() #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: print(ome)
This results in the following output results:
(2.1782707705367237-0.5074461983082241j)
(2.2227185408797476-0.5094990190123052j)
(2.2675293712400597-0.4965898981642067j)
(2.306187573060749-0.49327372936216596j)
(2.3375387477427036-0.4800143486380862j)
(2.3939839608351714-0.44250305302787263j)
(2.250056254487997+0.10549392553008125j)
(2.2934973319036587+0.09638866757241127j)
(2.3307875643942353+0.08609849468627778j)
(2.386139969059608+0.060068995693344174j)
(2.372968202422313+0.07669589505376706j)

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?
Time to create page: 0.101 seconds