Absorbing boundary conditions for acoustics in Discontinuous Galerkin

Question

Note: I'm trying to implement a Discontinuous Galerkin method, as kind of a way to learn about these things.

As of now, I've taken the acoustic wave equation $c^2 \nabla \cdot \nabla u(x,t) - \frac{\partial^2 u(x,t)}{\partial t^2} = f(x,t)$, rewritten it in weak form with the test function $v(x,t)$, and used Green's first identity to obtain a flux integral from the laplace operator:

$$\int_{\Gamma_K} c^2 v(x,t) \nabla u(x,t) \cdot dS - \int_K c^2 \nabla u(x,t) \cdot \nabla v(x,t) dV - \int_K \frac{\partial^2 u(x,t)}{\partial t^2} v(x,t) dV = \int_K f(x,t) v(x,t) dV$$ where $K$ is an element of the mesh and $\Gamma_K$ is its edges / faces.

I'm not super sure how ABC's work, but I believe the reason for using Green's first identity to get the flux integral on $\Gamma_K$ is so that we can replace existing boundary conditions with absorbing boundary conditions on the boundary faces of $\partial \Omega$? If I understand correctly, these absorbing boundary conditions would also be a differential equation that are only imposed on the boundary $\partial \Omega$.

My question is: What are some good absorbing boundary conditions (differential equations) I can use for the acoustic wave equation in a 2-dimensional or 3-dimensional rectangular domain?

Answer 1

The problem with ABC in the form of derivative operators at the boundary is that the exact ABCs are non-local in space and often in time, which makes them hard to use in numerical simulations. You can derive approximate ABCS based on series expansions, but this can be tedious and it may be difficult to get stable schemes. This approach was made famous by a very widely cited paper by Enquist and Majda.

An easier and widely used approach are perfectly matched layers. You extend your computational domain $\Omega$ by an additional layer in which you modify your derivatives by a complex scaling factor

$\partial_x \rightarrow \left(1 + \frac{i \sigma(x)}{\omega} \right) \partial_x$

This modifies the dispersion relation of your system such that traveling wave modes in the exterior layer are damped. Here, $\sigma$ is chosen such that it is zero at the boundary of your physical domain $\partial \Omega$ and increases towards the end of your layer, thus slowly increasing the damping. A function like $\textrm{erf}(x) = \frac{1}{\sqrt{\pi}}\int_{0}^{x} \exp(-s^2) ds$ might work well for $\sigma$, but you will need to rescale it so that it starts at zero at the physical boundary. You may have to play a little with your parameters $\sigma$ and $\omega$ to get good results.

Note that PML provides an approximate ABC, unless the layer would be infinitely long, which is of course impossible in a numerical simulation. However, the amplitude of the waves usually decays exponentially, so that small layers can already damp your reflections to below the discretization error of your DG method in the interior.

Originally, PML was introduced by Berenger although this paper might not be the best way to start. There is also a well written Wikipedia entry and searching for perfectly matched layer certainly will provide you with good material - as said, the method is very widely used.

Answer 2

There exist Absorbing Boundary Conditions for the wave equation that are stable and that go up to any order of accuracy (limited only by the accuracy of discretization of your model), so that they are good competition for PMLs.

A commonly used ABC of the Enquist-Majda type mentioned above, which is of the second order, is:

$u_{tt}-u_{xt}-\frac{1}{2}u_{yy}=0$

(Derivatives are implied). This should be plugged into the boundary term in your weak form, but you may need to define $u_x$ on the boundary as an additional variable in your model to do that.

If you need higher accuracy I suggest the Hagstrom-Warburton boundary condition: this condition is based on the definition of auxiliary variables on the boundary, which also satisfy the wave equation there. Corners require special treatment.