Is there a way to implement perfectly matched layers?

[JakobWagner]

I am currently working on an Acoustic Simulation (Helmholtz Equation) in dolfinx (With the nightly build on docker). Due to the nature of the PDE every boundary is a Neumann boundary an thus always creates reflections on the “outlet” side of the mesh.
Is there a way to implement PMLs in Dolfin-X? Is there another way of doing this?
Thanks in advance for any help

[bleyerj]

Hi,
this has been discussed here

How to treat complex damping functions for PMLs?

Hi everyone, I’m building knowledge on perfectly-matched-layers and how to define them on FEniCS, but I’m stuck now. Consider the weak form of the Helmholtz equation in the bounded domain x = [-a^*, a^*] and y = [-b^*, b^*], where a<x<a^* and b<y<b^* define the PML. \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)…

[Samuel_Groth]

Hi Jakob, a simple way to achieve this is to introduce an absorbing layer around the edge of your domain. In this region you make the wavenumber complex, which leads to damping of the waves. To reduce reflections further, you can “switch” on the imaginary part of the wavenumber gradually inside the layer, e.g., as a quadratic function. This isn’t a true PML, which requires carrying around some auxiliary variables. Therefore, you may have to make it a few wavelengths thick to reduce reflections sufficiently. If you want to read more about PMLs, I recommend these notes by Steven Johnson at MIT: https://math.mit.edu/~stevenj/18.369/pml.pdf. Let me know how you get on. I have implemented this myself so can share some code if required, but I think it’s good to give it a go yourself first. Regards, Sam.

[JakobWagner]

thank you for sharing that. I dit not see i when i was looking for it.

[JakobWagner]

Thank you - That is just what i was looking for. I’ll give it a try and keep you posted on my progress.

[JakobWagner]

Is it as easy as taking Eq. (7) and implement it like this for +x?

p = 2 # degree of sigma function in PML transformation
sigma = dolfinx.Function(V)
sigma.interpolate(lambda x: np.maximum(0, x[0] - x_start) ** p)

And just replacing

ufl.dx

with:

(1 / (1 - 1j * sigma)) * ufl.dx

Many Thanks in advance,
Jakob

It seems to work for lower frequencys

But not work at all for higher ones

[Samuel_Groth]

Hi Jakob, I think the transformation should be a bit more intricate than that. But what I was advocating in my last post was to not worry about the PML transformation and instead use an absorbing layer in which you modify your wavenumber k → k + 1j*\sigma, i.e., make it complex. Then you have to be careful about the choice of the constant in the function sigma, it’s chosen to reduce the “round-trip” reflection of a wave in the layer to a desired level. Equation (10) in “The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers” by Oskooi et al (2008) tells you how to compute this constant.

So the only modification to your code should be to add some imaginary component to your wavenumber at the right edge of your domain. If you would like to see a code of mine, check out waves-fenicsx/demo0_scattering_by_circle.py at master · samuelpgroth/waves-fenicsx · GitHub
Here I’m solving the 2D scattering problem of a plane wave incident on a circle. In it I have an absorbing border around my square domain.

I’d recommend reading the introduction sections to Oskooi et al (2008) as well. It takes a bit of time to get one’s head around the PML and absorbing layer stuff but it’s worth the effort!