I have searched within this forum and within the project documentation, but couldn’t find a clear way to declare a complex funtion in FEniCS. I wish to declare the following incident plane wave in order to study a scattering problem
p^{inc}(\mathbf{x}) = e^{-i \kappa x_1}
But I don’t know how to deal with the complex part. Does the 2019.1.0 version support complex numbers natively or do I need to work with them by other means?
Besides, I’m not sure if this equation should be defined only on the boundary or within the weak form, for I have a(p,\overline{q}) = L ({\overline{q}}) and p = p^{inc} + p^{sc}.
Can somebody address these questions?
Many thanks!
Hi,
You can either split your formulation into real or complex parts. See some of the discussion on this post or if you want native support for complex variational forms, Dolfinx could be the way forward. Specifically, you may want to check the Helmholtz demo
Thank you @bhaveshshrimali. Well, I’m a bit overwhelmed by the requirements in order to install Dolphinx. I decided to try the other way around instead. Indeed, I’m dealing with the Helmholtz problem \nabla ^2 p + \kappa^2 p = 0
for which the weak form is \int_{\Omega}(\nabla p \cdot \nabla \overline{q} - \kappa^2 p \overline{q})d \Omega = \int_{\Gamma}(g \overline{q})ds
By looking through the provided link I wrote the following, where e^{-i \kappa x_1} = cos(\kappa x_1) - i sin(\kappa x_1)
from fenics import *
import matplotlib.pyplot as plt
nx = 200
ny = 200
k = 15
#mesh and function space
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1*P1)
tol = 1E-14
def LEFT(x, on_boundary):
return on_boundary and abs(x[0] + 1.0) < tol
bound_x_re = Expression(("cos(k*x[0])"), degree=1, k=k)
bound_x_im = Expression(("-sin(k*x[0])"), degree=1, k=k)
#boundary conditions
bcs = [DirichletBC(V.sub(0), bound_x_re, LEFT),
DirichletBC(V.sub(1), bound_x_im, LEFT)]
#variational problem
(p_re, p_im) = TrialFunction(V)
(q_re, q_im) = TestFunction(V)
g_re = Constant(0.0)
g_im = Constant(0.0)
a = inner(grad(p_re),grad(q_re))*dx + inner(grad(p_im),grad(q_im))*dx - k*k*p_re*q_re*dx - k*k*p_im*q_im*dx
L = dot(g_re, q_re)*ds + dot(g_im, q_im)*ds
#compute solution
p = Function(V)
solve(a == L, p, bcs)
p_re, p_im = p.split()
plot(p_re)
plt.show()
which gives
Could you tell me if there’s any issue in my implementation?
I believe the minus sign should go along with p_i and q_i, but if that’s not the case, can you please give some more details on the way you developed this expression?
Thanks!
The correct bilinear form for the Helmholtz operator in a lossless medium is \int\left (\nabla p_r \cdot\nabla q_r + \nabla p_i \cdot\nabla q_i - k^2(p_r q_r + p_i q_i)\right ) d\mathbf{x} .
The coupling between real and imaginary components comes at the exit boundaries (impedance boundary condition). You will have surface terms of the form k \int \left (p_r q_i - p_i q_r\right ) d\mathbf{S}
on these wave exit boundaries.
Your error is that you appear not to have included any exit boundaries. All the boundary conditions are reflecting (Dirichlet).
Hi @BillS
Thank you for your valuable comments. I think that FEniCS assumes Newmann BCs when no Boundary Condition is informed. Nevertheless, my model should indeed avoid the spurious reflexions on the outer boundaries.
The way I found to do this was to introduce PMLs, i.e \large{\int_{\Omega} (\frac{\gamma_y}{\gamma_x} \frac{\partial p}{\partial x} \frac{\partial q}{\partial x} + \frac{\gamma_x}{\gamma_y} \frac{\partial p}{\partial y} \frac{\partial q}{\partial y}} - \kappa^2 \gamma_x \gamma_y p q)dxdy = \int_{\Gamma}gq ds
by means of the damping functions \gamma_x and \gamma_y.
As for now I’m still struggling a bit to implement all this because I’m learning the math of the PML’s on the fly. My motivation is that my model should also account for wave propagation in solid media, for which I have already found a ready example using PML’s on the internet. Would you have any comments on my approach?
Ahh…PMLs are another story. These are anisotropic lossy pseudomaterials. In this case, you will have q_r, p_i and q_i, p_r coupling everywhere within the PML region. As you have found, splitting up the forms into real and imaginary parts is frustratingly messy.
If the input/output port regions permit only a single mode, I would use the local absorbing boundary condition. It’s much more compact (and sparse) than the PML and is simple to apply.
If you want to model radiation problems in the near field (or almost-near-field), PMLs may be the way to go. However, maybe higher-order local absorbing boundary conditions may provide satisfactory solutions on a radiation sphere that lies in the Fresnel (transition) zone.
Scattering problems in free space are best handled using integral-equation methods, in my opinion. You can use exact Green’s functions and need only calculate the fields at the scatterer surface, making the problem more tractable.
Modeling waves in solids can be fiendishly complicated: you have compression (dilatation) waves and two shear wave polarizations and coupling between everything at boundary interfaces. I would start easy - try the acoustic Helmholtz equation first. You may want to try Dolfinx, just to avoid the mess of algebra needed to express the bilinear forms in terms of real quantities.
@luc@BillS Please note some discussion following the FEniCS implementation of the Helmholtz equations and PMLs with RT Finite elements explained in this document : https://doi.org/10.5281/zenodo.3888145.
