Question about Integrate

More
2 years 5 months ago #4344 by ajf367
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))
More
2 years 5 months ago #4345 by ajf367
Replied by ajf367 on topic Question about Integrate
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))
print(Integrate(u, mesh))
More
2 years 5 months ago #4346 by joachim
Replied by joachim on topic Question about Integrate
The Arnold method starts with a random initial guess, which leads to random scaling of eigenvectors
The following user(s) said Thank You: ajf367
Time to create page: 0.101 seconds