on Wave Equation with PML

Last post, I wrote some doodle work on solving the wave equation under FEM framework. And we did not consider the boundary condition, since the equation assumes valid for all plane. Now if truncated, PML is a way to save it, just with tiny modification on our previous work.

Formulation is not as simple as frequency domain using a trivial transformation to introduce absorption term. We need to use ifft to transform our Helmholtz equation back to time domain.

Think about equation:

$$u_{tt}-\nabla \cdot (\mu \nabla u)=f$$

using Laplace transform

$$\tilde{u}={\int }_0^{\infty }\exp (st)u(x,t)\mathrm{d}t$$

Thus we got Helmholtz, $s$ is complex.

$$s^2\tilde{u}=\nabla \cdot (\mu \nabla \tilde{u})+\tilde{f}$$

making coordinate transform as we do with Helmholtz:

$$x\to \tilde{x}=x+\frac{1}{s}{\int }^R\sigma (x)\mathrm{d}x$$

with conductivity $\sigma$, which is positive outside domain and vanish inside.

Extend $\tilde{u}$ onto transformed space,

$$s^2\tilde{u}=\tilde{\nabla }\cdot (\mu \tilde{\nabla }\tilde{u})+\tilde{f}$$

The transforming back,

$$\begin{array}{r}(s^2+s({\sigma }_x+{\sigma }_y)+{\sigma }_x{\sigma }_y)\tilde{u}=\nabla \cdot (\mu \nabla \tilde{u})+\frac{\partial }{\partial x}(\mu \frac{{\sigma }_y-{\sigma }_x}{s+{\sigma }_x}{\tilde{u}}_x)+\frac{\partial }{\partial y}(\mu \frac{{\sigma }_x-{\sigma }_y}{s+{\sigma }_y}{\tilde{u}}_y)+f\end{array}$$

$s$ means $\frac{\partial }{\partial t}$, therefore $u_{tt}+({\sigma }_x+{\sigma }_y)u_t+{\sigma }_x{\sigma }_{y}u=\nabla \cdot (\mu \nabla u)+\nabla \cdot \varphi +f$

where ${\varphi }_t=A\varphi +\mu B\nabla u$ and $A=diag(-{\sigma }_x,-{\sigma }_y)$, $B=diag({\sigma }_y-{\sigma }_x,{\sigma }_x-{\sigma }_y)$. Initial $\varphi$ as 0.

Using finite difference is acceptable. Central scheme with time, and work carefully with gradient on spatial grid.

Or use FEM on spatial and RK4 for time integration. It should be fast on constructing the matrices if time-independent. Then it would be easy to calculate

$${\zeta }^{\prime }=f(\zeta )$$

each step involves 4 calls on the function $f$, and each call involves 10 sparse-matrix-dense-vector multiplication, timing takes 0.05 seconds each for a matrix size of 18000 or so. That is about 12 seconds for one time step!! Unless we can find a way to reduce the calculations.

Hooray! By profiling the calculations, multiplication is just too fast and mldivide is taking most of time. But we know one thing, everything is going to be applied the same operation mldivide at last, do it at one pass with reduce the overhead at each round, the real calculation time would be only 1e-3 or so, but overhead will make it 1e-1 somehow.

Another thing is calculating the vector-mass matrix like

$${\int }_{\Omega }{\varphi }_j(A\nabla {\varphi }_i)\mathrm{d}x$$

It seems taking twice the timing than just stiffness or mass. However, we know one thing, as a sparse matrix COO-format : $I,J$ are fixed during all calculations like this, then avoiding getting extra work will make things move faster.

Update: examples are included in femex/example/PML on github. Performance is reasonable.

Problematic part: The previous overhead argument is kind of wrong in Matlab, horzcat creates a pretty large copy of the input, that will allocates lots of memory in advance.

Another thing is the ode45, it stores all solutions. That is a large cost, and porting them into plot function also takes a lot of time.

Workaround: try to write a function in C using GSL to make this faster, or write some function with C to avoid storing the solution. Or simply add a callback.

📊 Numerical Verification

The implementation using a finite difference scheme confirms the effectiveness of the PML. A moving compactly supported source is used to generate waves, which are absorbed by a cubic PML profile.

🔄 Method Comparison

To explore the best approach for time-integrating the PML system, we compared five different methods:

1. Leapfrog: 2nd-order, highly efficient, and standard for wave equations.
2. RK4: Standard 4th-order explicit (Baseline for precision, but dissipative).
3. IMEX: Implicit-Explicit 1st-order, treating the stiff damping term implicitly.
4. ETD: Exponential Time Differencing, treating the damping term exactly via analytical exponentials.
5. RKC2: A 2-stage explicit Runge-Kutta-Chebyshev scheme with extended stability.

The table below summarizes the execution times for a $150\times 150$ grid over 600 steps ($\Delta t=0.005$):

Method	Type	Time (s)	Relative Speed
IMEX	Implicit Damping	0.50s	0.92x (Fastest)
ETD	Exact Damping	0.52s	0.96x
Leapfrog	Symplectic-like	0.54s	1.00x (Baseline)
RKC2	Extended Explicit	1.30s	2.40x slower
RK4	Standard Explicit	2.58s	4.75x slower

• Comparison Code: [[codes/2015 Spring/wave_pml_all_methods.py]]
• Insights:
  • RK4 is the slowest due to its 4 derivative evaluations per step.
  • IMEX and ETD perform exceptionally well, matching Leapfrog’s speed while offering robust stability in the strongly damped PML regions.
  • Leapfrog remains a fantastic general-purpose choice, balancing simplicity, speed, and physical fidelity without requiring specialized damping treatments.

🔗 See Also

• on PML
• on PML for 3D Wave Equation