- Thank you received: 0
Helmholtz equation with robin condition
7 years 1 month ago #639
by gcdiwan
Helmholtz equation with robin condition was created by gcdiwan
Dear NGSolve developers,
I am a newbie to NGSolve, apologies if this sounds like a silly question.
I am trying to solve Helmholtz equation (\Delta u + \omega^2 u = 0 in \Omega with robin bc: \nabla u \cdot \mathbf{n} + i \omega u = g on \Gamma). I set my exact solution to be a plane propagating wave travelling in +ve x direction, i.e., u = \exp (i \omega x). This is imposed via robin bc, so the weak form reads:
find u \in H^1 (\Omega), s.t.
\int_\Omega \nabla u \cdot \nabla v - \omega^2 u v d\mathbf{x} + i\omega \int_{\Gamma} u v = \int_{\Gamma} g v
for all v in H^1(\Omega)
The NGSolve example code I am following is here . That example code has a non-zero source term and a zero right hand side for the robin bc. While I have opposite of that, i.e., a zero source term and a non-zero term in right hand side in robin bc.
Here is my python script:
I must be doing something stupid when I assemble my rhs linear form as f.vec are all zeros. Perhaps the way I use normal vector is not correct or I don't understand how to assemble the rhs when source term is not present. What am doing wrong here? Many thanks.
I am a newbie to NGSolve, apologies if this sounds like a silly question.
I am trying to solve Helmholtz equation (\Delta u + \omega^2 u = 0 in \Omega with robin bc: \nabla u \cdot \mathbf{n} + i \omega u = g on \Gamma). I set my exact solution to be a plane propagating wave travelling in +ve x direction, i.e., u = \exp (i \omega x). This is imposed via robin bc, so the weak form reads:
find u \in H^1 (\Omega), s.t.
\int_\Omega \nabla u \cdot \nabla v - \omega^2 u v d\mathbf{x} + i\omega \int_{\Gamma} u v = \int_{\Gamma} g v
for all v in H^1(\Omega)
The NGSolve example code I am following is here . That example code has a non-zero source term and a zero right hand side for the robin bc. While I have opposite of that, i.e., a zero source term and a non-zero term in right hand side in robin bc.
Here is my python script:
Code:
# following: https://ngsolve.org/docu/nightly/i-tutorials/unit-1.7-helmholtz/helmholtz.html
from ngsolve import *
from netgen.geom2d import unit_square
#
ngsglobals.msg_level = 1
# generate a triangular mesh of mesh-size 0.2
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
# set wavenumber
waveno = 1.0
# H1-conforming finite element space
fes = H1(mesh, order=1, complex=True)
# define trial- and test-functions
u = fes.TrialFunction()
v = fes.TestFunction()
# Forms
a = BilinearForm(fes)
a += SymbolicBFI(grad(u)*grad(v)-waveno**2*u*v)
a += SymbolicBFI(waveno*1j*u*v, definedon=mesh.Boundaries("outer"))
a.Assemble()
# is this correct way to get normals
normals = specialcf.normal(mesh.dim)
# rhs
f = LinearForm(fes)
f += SymbolicLFI((normals[0]+1.0)*(exp(1j*waveno*x)*1j*waveno) * v, definedon=mesh.Boundaries("outer"))
f.Assemble()
#solve
gfu = GridFunction(fes, name="u")
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu)
I must be doing something stupid when I assemble my rhs linear form as f.vec are all zeros. Perhaps the way I use normal vector is not correct or I don't understand how to assemble the rhs when source term is not present. What am doing wrong here? Many thanks.
- christopher
-
- Offline
- Administrator
-
Less
More
- Thank you received: 101
7 years 3 weeks ago - 7 years 3 weeks ago #667
by christopher
Replied by christopher on topic Helmholtz equation with robin condition
Hi, from a first glance I can see that you use the wrong name for the boundary conditions. The unit square has bcnames "top", "bottom", "right" and "left" and not "outer". This may be the reason why your terms do not do anything. Try replacing "outer" with "top|left|right|bottom".
You can get the available boundary names of a mesh with
Best
Christopher
You can get the available boundary names of a mesh with
Code:
print(mesh.GetBoundaries())
Christopher
Last edit: 7 years 3 weeks ago by christopher.
Time to create page: 0.116 seconds