PML boundary conditions

Question

I set up two one-way wave equations for constant velocity $c$ in one-dimension. When I implement them I get a highly unstable (divergent) solution. I wonder if someone could give me a suggestion about what could be wrong with my reasoning. Thanks.

The idea behind the PML method Perfectly Matched Layer is to attenuate the wavefield by doing analytic continuation of the real field into the complex in a way that the field in a buffer zone is damped exponentially. This is done by doing a change of variable.

\begin{eqnarray} x' = x + \mathrm{i} \frac{f(x)}{\omega}. \tag{eq. 1} \end{eqnarray}

Here $x, f(x), \omega \in \mathbb{R}$ and $x'$ is a complex number with $\mathrm{i}=\sqrt{-1}$. The factor $1/\omega$ is set up to avoid grid dispersion.

We assume a wave equation in the complex plane and time with variables $(x', t)$ This is:

\begin{eqnarray*} u_{x'x'} - \frac{1}{c^2} u_{tt}=0. \end{eqnarray*}

After applyting the Fourier transform in time we find the Helmholtz equation:

\begin{eqnarray*} \frac{\partial^2 U(x',\omega)}{\partial x'^2} + \frac{\omega^2}{c^2} U(x',\omega)=0 \end{eqnarray*} which we can factor as:

\begin{eqnarray*} \left ( \frac{\partial }{\partial x'} + \mathrm{i} \frac{\omega}{c} \right ) \left (\frac{\partial }{\partial x'} - \mathrm{i} \frac{\omega}{c} \right ) U(x', \omega)=0. \end{eqnarray*} and find

\begin{eqnarray} \left ( \frac{\partial }{\partial x'} + \mathrm{i} \frac{\omega}{c} \right ) U(x', \omega) = 0 \quad , \quad \left ( \frac{\partial }{\partial x'} - \mathrm{i} \frac{\omega}{c} \right ) U(x',\omega) = 0. \tag{eq.2} \end{eqnarray} These are advection equations . The first represents a wave traveling forward (say, to the right) with wavespeed $c$ and the second a wave traveling backward (say, to the left) with wave speed $c$.

The change of variable $x'=x+\mathrm{i} \frac{f(x)}{\omega}$ produces, by using the chain rule,

\begin{eqnarray*} \frac{\partial U}{\partial x'} &=& \frac{\partial U}{\partial x} \frac{d x}{d x'} = \frac{1}{\left ( 1 + \frac{\mathrm{i} f_x}{\omega} \right )} \frac{\partial U}{\partial x} \end{eqnarray*}

now from the first of equations eq.2

\begin{eqnarray*} -\mathrm{i} \frac{\omega}{c} U = \frac{1}{ \left ( 1 + \mathrm{i} \frac{f_x}{\omega} \right )} U_{x}. \end{eqnarray*} or

\begin{eqnarray*} -\mathrm{i} \frac{\omega}{c} U + \frac{f_x}{c} U = U_x. \end{eqnarray*} Clearly, if $f=0$ (or constant) the equation in $x$ coincides $x'$ and the wave in $x'$ coincides with the classical wave equation for real domain $x$. We can write the equation in time (by using the Fourier inverse transform) as

\begin{eqnarray} \frac{1}{c} u_t + \frac{f_x}{c} u = u_x \tag{eq.3} \end{eqnarray}

We now verify that this equation fits our purposes. Let us assume a solution of the form

\begin{eqnarray} u(x,t)=A \mathrm{e}^{\mathrm{i}(k x-\omega t)} \mathrm{e}^{-k f(x)/\omega} \tag{eq.4} \end{eqnarray} If we substitute this solution in equation (eq.3) we find, for the left hand side

\begin{eqnarray*} \left ( -\frac{\mathrm{i} \omega }{c} + \frac{f_x}{c} \right ) u(x,t) \end{eqnarray*} and, for the right hand side

\begin{eqnarray*} \left ( \mathrm{i k} - \frac{k f'(x)}{\omega} \right ) u(x,t). \end{eqnarray*} Since $k=\omega/c$ we see that the function $u(x,t)$ defined in (eq.4) satisfies equation (eq.3).

If, instead of chosing the first of equations (eq.3), we pick the second we find, following a similar procedure

\begin{eqnarray} -\frac{1}{c} u_t + \frac{f_x}{c} u = u_x. \tag{eq.5} \end{eqnarray} Similarly, we find that

\begin{eqnarray} u(x,t)=A \mathrm{e}^{\mathrm{i}(k x+\omega t)} \mathrm{e}^{-k f(x)/\omega} \tag{eq.6} \end{eqnarray} satisfies (eq.5).

Solution (eq.4) represents a wave moving to the right with a speed $c$, wave number $k$, frequency $\omega$ and attenuation factor $f(x)$. Likewise, solution (eq.6) represents a wave moving to the left with the same parameters above. The superposition of all this "plane-wave" solutions over all frequencies and wave numbers constrained with initial and boundary conditions is the solution of the initial/boundary value problem.

In both cases (eq.4) and (eq.6), if $f(x)=0$ we get back to the initial problem with no attenuation over the real variable $x$.

Let us define $f_x=\sigma(x)$. That is $f(x)=\int_{x_0}^x \sigma(s) ds$ then, in (eq.3) we use forward differences in time and central differences in space to get:

\begin{eqnarray*} \frac{ w_{i j+1} - w_{ij}}{\Delta t} = \frac{c}{2 \Delta x} (w_{i+1 j} - w_{i-1 j}) - \sigma_i w_{ij}, \end{eqnarray*} using a forward difference in time and central in space and the CFL number $\gamma= c \Delta t/\Delta x$. We find

\begin{eqnarray*} w_{i j+1} = \frac{\gamma}{2}( w_{i+1 j} - w_{i-1 j}) - w_{ij}(\sigma_i \Delta t -1) \end{eqnarray*}

We use equation (eq.5) for the left side of the working interval (model). The corresponding finite difference system is

\begin{eqnarray*} w_{i j+1} = -\frac{\gamma}{2}( w_{i+1 j} - w_{i-1 j}) + w_{ij}(\sigma_i \Delta t-1). \end{eqnarray*}

I use these finite difference equations in the edges of the model (last and first 20% or 10%....nothing works).

----------------update------------------

I changed the finite difference scheme to forward in both, space and time. Now it is stable, but the the modeling acts as if the sides where shifted and the reflections are taking place at $[a+d,b-d]$ where $[a,b]$ is the spatial domain of $u$ and $d$ is the buffer zone for the damping. The field $u$ is zeroed in the intervals $[a,a+d]$ and $[b-d,b]$.

Answer 1

I found that the solution of the PML problem by the method proposed here should not work. Next is a detailed analysis of the problem.

I start with the von Neumann stability analysis for the FTCS system corresponding to (eq 3.) on the original question.

\begin{eqnarray*} \frac{ w_{i j+1} - w_{ij}}{\Delta t} = \frac{c}{2 \Delta x} (w_{i+1 j} - w_{i-1 j}) - \sigma_i w_{ij}, \end{eqnarray*} Call $\gamma=c \Delta t/\Delta x$ (the CFL number) and so we have \begin{eqnarray*} w_{i j+1} = \frac{\gamma}{2} ( w_{i+1 j} - w_{i-1 j}) - w_{ij}(\sigma_i \Delta -1). \end{eqnarray*}

For the von Neumann analysis we assume $w_{ij}=\mathrm{e}^{\mathrm{i} \omega t_j} \mathrm{e}^{\mathrm{i} k x_k}$ and call $G=w_{i j+1}/w_{ij}=\mathrm{e}^{\mathrm{i} \omega \Delta t }$ the complex gain factor. After replacing and dividing by $w_{ij}$ we find

\begin{eqnarray*} G = \mathrm{e}^{\mathrm{i} \omega \Delta t} = \gamma \mathrm{i} \sin k \Delta x - \sigma_i \Delta t +1. \end{eqnarray*} and, since $\gamma, \sigma_i k, \Delta x, \Delta t \in \mathbb{R}$ we find that

\begin{eqnarray*} |G| = \sqrt{(1-\sigma_i \Delta t)^2 + \gamma^2 \sin^2 k \Delta x} \end{eqnarray*} If we want $|G|<1$ then we need

\begin{eqnarray*} \gamma^2 \sin^2 k \Delta x \le 1 - (1 -\sigma_i \Delta t)^2 \end{eqnarray*} so that, we need

\begin{eqnarray*} \gamma \le \sqrt{ \frac{1-(1-\sigma_i \Delta t)^2}{|\sin k \Delta x|}} \le \sqrt{ 1 - (1 - \sigma_i \Delta t)^2} \end{eqnarray*} If $\sigma_i=0$ (case $x=x_0$, starting point of buffer zone where $f(x_0)=0$) then $\gamma \le 0$, and this impedes the election of $\gamma$ which guarantee stability. This idea then does not work.

Let us now try forward differences in both domain FTFS: \begin{eqnarray*} \frac{ w_{i j+1} - w_{ij}}{\Delta t} = \frac{c}{\Delta x} (w_{i+1 j} - w_{i j}) - \sigma_i w_{ij}, \end{eqnarray*} That is,

\begin{eqnarray*} w_{i j+1} = \gamma w_{i+1 j} - (\sigma_i \Delta t+\gamma -1) w_{ij}. \end{eqnarray*}

We perform von Neumann stability analysis:

\begin{eqnarray*} G = \mathrm{e}^{\mathrm{i} \omega \Delta t} &=& \gamma \mathrm{e}^{\mathrm{i} k \Delta x} -(\sigma_i \Delta t + \gamma -1) \\ &=& \gamma \cos k \Delta x -(\sigma_i \Delta t + \gamma -1) + \mathrm{i} \gamma \sin k \Delta x \end{eqnarray*} La absolute value of $G$ is given by

\begin{eqnarray*} |G| &=& \gamma^2 - 2 \gamma \cos k \Delta x(\sigma_i \Delta t + \gamma -1) + (\sigma_i \Delta t + \gamma -1)^2 \\ &=& \gamma^2 (2 - 2 \cos k \Delta x ) + 2 \gamma ( \sigma_i \Delta t -1 - \cos k \Delta x + \sigma_i \Delta t \cos k \Delta x) -2 \sigma_i \Delta t + 1 \\ &=& \gamma^2 (2 - 2 \cos k \Delta x ) + 2 \gamma [ \sigma_i \Delta t -1 + \cos k \Delta x( \sigma_i \Delta t -1)] -2 \sigma_i \Delta t + 1 \\ &=& \gamma^2 (2 - 2 \cos k \Delta x ) + 2 \gamma [ (\sigma_i \Delta t -1)(1 + \cos k \Delta x)] -2 \sigma_i \Delta t + 1 \\ \end{eqnarray*} If we want $|G| \le 1$ we need

\begin{eqnarray*} \gamma^2 (2 - 2 \cos k \Delta x ) + 2 \gamma [ (\sigma_i \Delta t -1)(1 + \cos k \Delta x)] -2 \sigma_i \Delta t \le 0 \end{eqnarray*} The left hand-side is a parabola convex up, since $2-2 \cos k \Delta x > 0$. The discriminant $D$ of this equation is given by

\begin{eqnarray*} D = 4 (\sigma_i \Delta t-1)^2 (1+ \cos k \Delta x)^2 +8 \sigma_i \Delta t (2 - 2 \cos k \Delta x) > 0 \end{eqnarray*} so that there exist a range of values of $\gamma$ where this parabola is under the $x$ axis. This range is given by the interval:

\begin{eqnarray*} [ -2(\sigma_i \Delta t -1)(1 + \cos k \Delta x) - \sqrt{D} \, , \, -2(\sigma_i \Delta t -1)(1 + \cos k \Delta x) + \sqrt{D}] \end{eqnarray*} $D$ can be as close as we want to $2 ( \sigma_i \Delta t -1)(1+\cos k \Delta x)$ by making $k \Delta x$ close to $0$. Hence we cannot guaranteed $\gamma > 0$ for stability. We see two classical methods for which stability fails.

We now turn our attention to implicit methods such as the Crank-Nicholson method.

Let us get back to equation (eq.3) and write:

\begin{eqnarray*} w_{i j+1} = \frac{\gamma}{4} ( w_{i+1 j+1} - w_{i-1 j+1} ) + \frac{\gamma}{4} ( w_{i+1 j} - w_{i-1 j}) + w_{i j-1} (1 -\sigma_i). \end{eqnarray*} We do a von Neumann analysis.

Assuming: $G=w_{i j+1}/w_{ij}= \mathrm{e}^{\mathrm{i} \omega \Delta t}$

\begin{eqnarray*} G &=& \frac{\gamma}{2} G \mathrm{i} \sin k \Delta x + \frac{\gamma}{2} \mathrm{i} \sin k \Delta x + (1-\sigma_i) \\ \left ( 1 - \mathrm{i} \frac{\gamma}{2} \sin k \Delta x \right) G &=& \left ( 1 + \mathrm{i} \frac{\gamma}{2} \sin k \Delta x - \sigma_i \right ) \end{eqnarray*} that is,

\begin{eqnarray*} | G | = \frac{(1-\sigma_i)^2 + (\gamma^2/4) \sin^2 k \Delta x} {1 + (\gamma^2/4) \sin^2 k \Delta x}. \end{eqnarray*} We find that if, $0<\sigma_i<2$ then $|G| \le 1$ and the system is unconditionally convergent and if $\sigma_i < 0$ the system is unconditionally divergent. The solution (eq.4) diverges for $f(x)>0$ (assuming $k>0, \omega > 0$), and if, $\sigma(x) >0$ then $f(x)=\int_{x_0}^x \sigma(s) ds >0 $ for $x_0<x$ and $f(x) \le 0$ for $x_0>x$. However the solution (eq.4) corresponds to a pulse traveling to the left, so that it diverges for $x<x_0$ (left).

We conclude that all methods tried to solve the PML problem based on the factorization of the two-way wave equation into two single-way way equations are not appropriate.