Hello everybody
I try to solve the Helmholtz equation with an absorbent subdomain.
$$\Delta p + k^2 z p = 0$$
with
$$p = p_r + ip_i$$
and
$$z = z_r + iz_i$$
The boundary conditions are periodic on the y and z surfaces of the BoxMesh. On the left (x=0) the boundary condition is fully reflective and on the right (x=Lx) I used a combination of an incoming wave and the Sommerfeld radiation condition:
$$\vec{n}\cdot\vec{\nabla}p = -ikp + 2p_0ik\exp{ikx}$$
The code I wrote runs without any problems, but the result is wrong as soon as z_i is not equal to zero in any part of the mesh.
I'm not sure if I made a mistake while splitting the equation into real and imaginary parts, or if I just have an error in the code.
Please help me
import math as m
import numpy as np
from dolfin import *
################################################################################
# set measurements and create mesh
thickness = 0.1
length = 1.1
Lx = thickness + length
Ly = 0.01
Lz = 0.01
Nx = 600
Ny = 1
Nz = 1
mesh = BoxMesh(0., 0., 0., Lx, Ly, Lz, Nx, Ny, Nz)
c = 343.4
rho = 1.204
xi = 15000.
kappa = 1.5
sigma = 0.99
# simulation name
simname = 'test'
################################################################################
# set periodic boundary conditions and create constrained function space
class PeriodicDomain(SubDomain):
# set periodic boundary conditions
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two
# slave edges
return bool((near(x[1], 0) or near(x[2], 0.)) and
(not ((near(x[1], Ly) and near(x[2], 0.)) or
(near(x[1], 0) and near(x[2], Lz)))) and on_boundary)
def map(self, x, y):
if near(x[1], Ly) and near(x[2], Lz):
y[0] = x[0]
y[1] = x[1] - Ly
y[2] = x[2] - Lz
elif near(x[1], Ly):
y[0] = x[0]
y[1] = x[1] - Ly
y[2] = x[2]
elif near(x[2], Lz):
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - Lz
else:
y[0] = -1000
y[1] = -1000
y[2] = -1000
periodicdomain = PeriodicDomain()
# create by periodic domain constrained function space
W = FunctionSpace(mesh, 'Lagrange', 2, constrained_domain=periodicdomain)
W2 = W * W
################################################################################
# define subdomains
class Absorber(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= thickness
class Air(SubDomain):
def inside(self, x, on_boundary):
return x[0] >= thickness
subdomains = MeshFunction('size_t', mesh, 3)
subdomain0 = Absorber()
subdomain0.mark(subdomains, 0)
subdomain1 = Air()
subdomain1.mark(subdomains, 1)
V0 = FunctionSpace(mesh, 'DG', 0)
zr = Function(V0)
zi = Function(V0)
zr_values = [kappa, 1.]
zi_values = [-xi*sigma / (rho), 0.]
for cell_no in range(len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
zr.vector()[cell_no] = zr_values[subdomain_no]
for cell_no in range(len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
zi.vector()[cell_no] = zi_values[subdomain_no]
################################################################################
# define ABC Boundary Condition
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
class ABC(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], Lx)# or near(x[0], 0.))
abc = ABC()
abc.mark(boundaries, 1)
################################################################################
# define expressions for variational problem and run for all frequencies
frequencies = np.arange(50, 5624, 20)
for frequency in frequencies:
omega = 2.*m.pi*frequency
agemo = Constant(1./omega)
k = omega/c
p0r = 1.
p0i = 1.
f1 = Expression('2*k_*(p0r_*sin(k_*(x[0]-L_)) + p0i_*cos(k_*(x[0]-L_)))', k_=k, p0r_=p0r, p0i_=p0i, L_=Lx)
f2 = Expression('2*k_*(p0r_*cos(k_*(x[0]-L_)) - p0i_*sin(k_*(x[0]-L_)))', k_=k, p0r_=p0r, p0i_=p0i, L_=Lx)
ksquared = Constant(k**2.)
k = Constant(k)
(pr, pi) = TrialFunction(W2)
(w1, w2) = TestFunction(W2)
a = (pr*zr*w1*ksquared - zi*agemo*pi*w1*ksquared +
zi*agemo*pr*w2*ksquared + zr*pi*w2*ksquared -
inner(nabla_grad(pr), nabla_grad(w1)) -
inner(nabla_grad(pi), nabla_grad(w2)))*dx \
+ (k*pi*w1 - k*pr*w2)*ds(1)
L = f1*w1*ds(1) - f2*w2*ds(1)
A = assemble(a, exterior_facet_domains=boundaries)
b = assemble(L, exterior_facet_domains=boundaries)
# compute solution
p = Function(W2)
P = p.vector()
solve(A, P, b)
pr, pi = p.split()
# write solution to file in VTK format
result_path = 'result/'
file = File(result_path + simname + '_' + str(frequency) + '_pr.pvd')
file << pr
file = File(result_path + simname + '_' + str(frequency) + '_pi.pvd')
file << pi
