I am currently working on an implementation of PMLs in FEniCS to allow better simulation of scattering problems. PMLs are layers surrounding the region of interest (the inner domain) in which the scattered waves are absorbed so as to avoid reflections back into the inner domain.
Starting in 1D, i have tried creating an IntervalMesh and splitting it into an inner domain and a (surrounding) pml domain. Since the incoming field should only be present on the inner domain, the BCs are applied at the edge of the inner domain. Hence if the thickness of the pml layers (dpml) is nonzero, the inner domain boundaries will not coincide with the mesh boundaries. In this case, i get an all-zero solution. If dpml=0 i get the expected solution.
Below is a MWE. To get the all-zero solution, simply change dpml to a non-zero value.
from dolfin import *
import numpy as np
class LambdaDomain(SubDomain):
def __init__(self, condition):
self.condition = condition
SubDomain.__init__(self)
def inside(self, x, on_boundary):
return self.condition(x, on_boundary)
wl = 500
lx = wl*2
nx = 1000
dpml = 0.0 * wl
# Create mesh.
lx_full = int(lx + 2.0 * dpml)
nx_full = int((lx_full / lx) * nx)
mesh = IntervalMesh(nx_full, 0.0, lx_full)
# Mark boundaries.
facet_function = FacetFunction('size_t', mesh)
LambdaDomain(lambda x, on: near(x[0], dpml, lx / nx) and x[0] >= dpml).mark(facet_function, 1)
# ds = Measure("ds")[facet_function] # Original code
dS = Measure("dS")[facet_function]
# Mark domains.
cell_function = CellFunction('size_t', mesh)
LambdaDomain(lambda x, on: dpml <= x[0] <= lx + dpml).mark(cell_function, 1)
dx = Measure("dx")[cell_function]
# Define constants.
k0 = Constant(2 * np.pi / wl)
A0 = Constant(1.0)
# Create double valued space (real + imag).
V = FunctionSpace(mesh, 'CG', 2)
V2 = V * V
# Create test/trail functions (real + imag).
(u_r, u_i) = TrialFunction(V2)
(v_r, v_i) = TestFunction(V2)
# EQ: 1D scalar helmholtz.
a = k0 * k0 * (v_r * u_r + v_i * u_i) * dx(1) - \
(inner(nabla_grad(u_r), nabla_grad(v_r)) + inner(nabla_grad(u_i), nabla_grad(v_i))) * dx(1)
# BC: Incoming plane wave from left.
# L = - 2 * v_r * Expression('-A0*(k0*sin(k0*x[0]))', k0=k0, A0=A0) * ds(1) + \
# - 2 * v_i * Expression('-A0*(k0*cos(k0*x[0]))', k0=k0, A0=A0) * ds(1) # Original code
L = - 2 * avg(v_r) * avg(Expression('-A0*(k0*sin(k0*x[0]))', k0=k0, A0=A0)) * dS(1) + \
- 2 * avg(v_i) * avg(Expression('-A0*(k0*cos(k0*x[0]))', k0=k0, A0=A0)) * dS(1)
# Solve equations
A = assemble(a)
b = assemble(L)
p = Function(V2)
solve(A, p.vector(), b)
pr, pi = p.split(True)
# Plot stuff.
plot(facet_function)
plot(cell_function)
plot(pr)
interactive()
I guess i am applying the BCs wrong somehow. Can anyone point me in the right direction?
