I apologize if this seems naive but I’m running the following code but every time I rerun the code the value of the integral changes. On occasions the value of the integral has positive real part, on others it has negative real parts, etc. I have the the out resonance of the Arnoldi unchanged through each run but the u still changes.

Could someone pinpoint what is occurring to cause this lack of consistency?

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

geo.AddCircle ( (0.0, 0.0), r=2.0, leftdomain=2, rightdomain=3)

geo.AddCircle ( (0.0, 0.0), r=3.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.04)

mesh = Mesh(ngmesh)

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

mesh.SetPML(pml.Radial(rad=2.0, alpha=1j, origin=(0,0)), “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.111)

Draw(u)

print ("lam: ", lam)

print(Integrate(u, mesh))[/code]