- Thank you received: 0
Question about Integrate
2 years 5 months ago #4344
by ajf367
Question about Integrate was created 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?
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))
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))
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))
Time to create page: 0.101 seconds