PML for 3D Wave Equation

In previous post, the wave equation in 2D, with PML for absorption boundary condition. The system is augmented to 4 unknown variables. And in 3D case, regardless of possible issues from PML, the system should involve more unknowns for additional dimension, however, it is actually more than that. It will bring an $s^{-1}$ term, representing an temporal integral.

By Laplace transform (or Fourier transform) to frequency domain, the equation is simply as ($s$ is complex)

$s^{2}u=c^2\Delta u$

By changing coordinate for each axis (e.g. $x,y,z$), $\tilde{x}=x+{\int }_0^x{\sigma }_x(\xi )d\xi$

we will arrive at a new system

$(s^2+s({\sigma }_x+{\sigma }_y+{\sigma }_z)+({\sigma }_x{\sigma }_y+{\sigma }_y{\sigma }_z+{\sigma }_z{\sigma }_x)+{\sigma }_x{\sigma }_y{\sigma }_z{s}^{-1})u=c^2\Delta u+\nabla \cdot \Phi$

where $\Phi =({\varphi }^1,{\varphi }^2,{\varphi }^3)$ are auxiliary functions, inverting $s$ to ${\partial }_t$, we put the equation back to time domain,

$u_{tt}+p_1{u}_t+p_{2}u+p_{3}U=c^2\Delta u+\nabla \cdot \Phi$

${\Phi }_t+\Sigma \Phi =c^2(p_1\mathbb{I}-\Sigma )\nabla u+\Gamma \nabla U$

$U_t=u$

where $\Sigma =\mathrm{diag}({\sigma }_x,{\sigma }_y,{\sigma }_z)$ and $\Gamma =\mathrm{diag}({\sigma }_y{\sigma }_z,{\sigma }_z{\sigma }_x,{\sigma }_x{\sigma }_y)$. $p_1=\mathrm{Tr}(\Sigma )$, $p_2=\mathrm{Tr}(\Gamma )$, $p_3=\det (\Sigma )$.

There are 6 unknowns to solve : $Y=(u,u_t,{\varphi }^1,{\varphi }^2,{\varphi }^3,U)$. The system’s initial condition is from wave equation, additional variables are initialized as zero.

$Y_t=LY$

can be solved with various numerical methods, $L$ is a second order operator in spatial variables, thus we can use FDM, FEM, pseudo-spectral.

Pseudo-spectral is slower in complexity, but it will involve less points, since it is more accurate on the derivatives. The system needs evaluation on $U_x,U_y,U_z,u_x,u_y,u_z,c^2\Delta u,\nabla \cdot \Phi$ respectively, which requires fft for 5 times on $u,U,\Phi$ and ifft for 8 times, each fft/ifft in theory needs $15Nlo{g}_{2}N$ in 3D, thus total flops are $195N{\log }_{2}N$ flops, for a grid as large as 200x200x200, the total flops will be 7E10 or so. On a single core machine at effective frequency 2.0GHz, the time will be 30s for one evaluation! For multi-core (say quad-core) platform, it will require (maybe) 10s for one run, 1000 time-steps with forward Euler will be about 3h, multi-step methods like RK2, RK3, RK4, the time will be much longer.

If precision is not important, using single precision will cut the timing in half.