- Thank you received: 0
pml on waveguides
4 years 11 months ago - 4 years 11 months ago #2228
by gcdiwan
pml on waveguides was created by gcdiwan
[attachment=undefined]schema.PNG[/attachment] [attachment=undefined]schema.PNG[/attachment]
Dear NGSolve developers,
I have trouble with pml when used for a waveguide problem. Consider a 2D geometry: \Omega = [0, Range] X [0, Depth] and the Helmholtz equation:
\Delta u + k^2 u = -f in \Omega
with
\frac{\partial u}{\partial n} = 0 on \Gamma_N (top and bottom boundaries).
To the left and right, we have \Gamma_R where we use the PML, so that \partial \Omega = \Gamma_N \cup \Gamma_R
I wish to prescribe f as a point source at (x0,y0) located in \Omega so I followed the pulse implementation given in the NGSolve’s the i-tutorials page ( ngsolve.org/docu/nightly/i-tutorials/uni...holtz/helmholtz.html )
For the geometry described above, I created the mesh externally in GMsh and specified the regions for PML (to left and right of the source) by using the domain keywords i used in GMsh. While I can see that the solution decays within the PML regions, I would have expected it to undergo an exponential decay due to geometrical spreading loss as it travels to left and right of the point source. Does the NGSolve pml implementation take into account the modal behaviour of waves in a confined cavity unlike fully exterior problem?
I have attached the script that I am using and a schematic of the problem to get an idea.
Many thanks for your help.
Attachment not found
I have trouble with pml when used for a waveguide problem. Consider a 2D geometry: \Omega = [0, Range] X [0, Depth] and the Helmholtz equation:
\Delta u + k^2 u = -f in \Omega
with
\frac{\partial u}{\partial n} = 0 on \Gamma_N (top and bottom boundaries).
To the left and right, we have \Gamma_R where we use the PML, so that \partial \Omega = \Gamma_N \cup \Gamma_R
I wish to prescribe f as a point source at (x0,y0) located in \Omega so I followed the pulse implementation given in the NGSolve’s the i-tutorials page ( ngsolve.org/docu/nightly/i-tutorials/uni...holtz/helmholtz.html )
For the geometry described above, I created the mesh externally in GMsh and specified the regions for PML (to left and right of the source) by using the domain keywords i used in GMsh. While I can see that the solution decays within the PML regions, I would have expected it to undergo an exponential decay due to geometrical spreading loss as it travels to left and right of the point source. Does the NGSolve pml implementation take into account the modal behaviour of waves in a confined cavity unlike fully exterior problem?
I have attached the script that I am using and a schematic of the problem to get an idea.
Many thanks for your help.
Last edit: 4 years 11 months ago by gcdiwan. Reason: added attachments
4 years 11 months ago #2229
by gcdiwan
Replied by gcdiwan on topic pml on waveguides
Not sure why the attachment is not visible in my earlier message. Here i am pasting the code below:
Code:
from ngsolve import *
from netgen.geom2d import SplineGeometry
import numpy as np
import subprocess
from netgen.read_gmsh import ReadGmsh
Rmax = 3e3
Depth = 1e2
Lpml = 1e2
# import the Gmsh file to a Netgen mesh object
mesh = ReadGmsh('../../../meshes/waveguide2d_sym.msh')
# wrap it into an NGSolve mesh
from ngsolve import *
mesh = Mesh(mesh)
Draw(mesh)
pcR = pml.HalfSpace((Rmax,0.),(1,0),2j)
pcL = pml.HalfSpace((0,0.),(-1,0),2j)
# pc2 = pml.Cartesian((Rmax,0.),(Rmax+Lpml,Depth),2j) # doesn't absorb!
mesh.SetPML(pcR,"pmlR")
mesh.SetPML(pcL,"pmlL")
fes = H1(mesh, order=2, complex=True)
u, v = fes.TnT()
# Wavenumber & source
cspeed = 1500.
freq = 20.
x0=0.
y0=36.
omega = 2.*np.pi*freq / cspeed
ScaleFactor = 1e1
a = BilinearForm(fes, symmetric=True)
a += grad(u)*grad(v)*dx - omega*omega*u*v*dx
a.Assemble()
pulse = ScaleFactor *exp(-(omega**2)*((x-x0)*(x-x0) + (y-y0)*(y-y0)))
f = LinearForm(fes)
f += -pulse * v * dx
f.Assemble()
gfu = GridFunction(fes, name="u")
gfu.vec.data = a.mat.Inverse() * f.vec
vtk = VTKOutput(ma=mesh, coefs=[gfu.real], names = ["ureal"],filename="wgfem_real",subdivision=3)
vtk.Do()
vtk = VTKOutput(ma=mesh, coefs=[gfu.imag], names = ["uimag"],filename="wgfem_imag",subdivision=3)
vtk.Do()
- christopher
-
- Offline
- Administrator
-
Less
More
- Thank you received: 101
4 years 11 months ago #2235
by christopher
Replied by christopher on topic pml on waveguides
can you attach the mesh as well? In your example the left side pml starts at 0 and the right at Rmax, the cartesian at Rmax and Rmax+Lpml
Cartesian specifies a box, and outside of the box there are pml layers
Best
Christopher
Cartesian specifies a box, and outside of the box there are pml layers
Best
Christopher
4 years 11 months ago #2236
by gcdiwan
Replied by gcdiwan on topic pml on waveguides
Hi Christopher,
Thanks for your reply.
I have attached my geo file below.
Thanks again for your help.
Thanks for your reply.
Yes, but i have commented that specific line for cartesian pml as i did not understand how to use that or why it did not show any reduction in wave amplitude in the pml region. Most likely i have got the halfspace pml wrong as well as by definition it should be applied on one side alone, right?Cartesian specifies a box, and outside of the box there are pml layers
I have attached my geo file below.
Thanks again for your help.
Attachments:
Time to create page: 0.104 seconds