- Thank you received: 0

The forum is in read only mode.

#
Question about Integrate

*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))

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