I used to study Computational Fluid Dynamics (CFD), and one of my favorite techniques was Perfectly Matched Layers (PML) for absorbing boundaries because it just doesn’t seem like such a good solution should exist.
The basic problem setup is that we seek to simulate a wave equation and avoid waves reflecting off the boundaries and back into the simulation domain. Simply truncating the domain causes waves to reflect artificially, contaminating the solution.
Perfectly Matched Layers (PML) solve this problem by adding absorbing regions at the domain boundaries that eliminate reflections. You can imagine this like some sort of foam or tacky, viscous rubber damping sound waves. However, rubber behaves very differently from air and, so although it will dampen sound waves, it also causes its own spurious reflections.
The PML technique avoids these spurious reflections by “perfectly” matching the impedance of the absorbing layer to the impedance of the medium. In very hand-wavy terms, it accomplishes this by extending the solution into the complex domain and then applying damping to the imaginary component. The result is that the waves decay exponentially as they move. Moreover, this decay is frequency-independent and occurs both in the outgoing and incoming directions, so that any waves that aren’t fully damped and reflect off the boundaries are futher damped as they move back into the domain.
This notebook demonstrates PML for the scalar wave equation, starting with a 1D example before extending to 2D.
The scalar wave equation in general is given by
For one dimension, we have
We could apply PML to this equation, but it works better to start out with what’s probably a less familiar but strictly equivalent form containing only first-order derivatives,
Here, is the state variable (e.g. pressure) while functions like a velocity.
The key insight of PML is to work in the Fourier time domain and apply a complex coordinate stretching to the spatial derivatives. In practical terms, this means replacing spatial derivatives with a modified operator.
We begin by transforming to the Fourier time domain, replacing with
The PML technique replaces spatial derivatives with a complex coordinate mapping
where is the absorption strength. This stretching causes waves to decay exponentially in regions where . Notice that for (no absorption), this reduces to .
Applying this mapping gives
Multiplying both sides by gives
Transforming back to the time domain gives the final 1D PML equations
In the interior where , these reduce to the original wave equation. In the PML regions, causes exponential absorption.
The plot below shows a 1D wave propagating with PML absorbing boundaries on both ends. The shaded regions indicate where . Notice how waves exit the domain smoothly without reflecting.
The 2D case follows the same principles but with a key difference: we need separate absorption coefficients and for each spatial direction. When we apply the complex coordinate mapping independently to both and derivatives, cross-terms emerge that don’t appear in 1D. These cross-terms require an auxiliary field to handle the coupling between the two spatial directions within the PML regions.
The derivation below shows how these auxiliary variables arise naturally from the Fourier domain analysis. In regions where only one direction has absorption (e.g., but ), the equations simplify. The full complexity only appears in corner regions where both directions have non-zero absorption.
We start with the scalar wave equation in first-order form,
Following the same approach as in 1D, we transform to the Fourier time domain, replacing with
The PML technique replaces spatial derivatives with complex coordinate mappings
Applying the mapping to the velocity equations and multiplying through by the denominators yields
Transforming back to the time domain immediately gives
The equation for is more involved. Applying the PML mapping gives
Expanding and collecting terms yields
Except for the last term, we can transform this back into the time domain to obtain an equation analogous to the 1D case. However, the term in parentheses contains factors of which correspond to time integration. In other words, the equation mostly translates back into the time domain just fine, but we end up with one additional quantity we need to track and integrate along with the other state variables. We define the term in parentheses as the auxiliary quantity
Multiplying through by and transforming back to the time domain gives
The update equation for then becomes
Together with the velocity equations, we have the complete 2D PML system
These equations require tracking three fields: the scalar field , the auxiliary field , and the vector field . In practice, this naturally lends itself to a two-buffer implementation: one storing and the other storing . The two buffers are offset by half a timestep to avoid some of the common instabilities encountered when simulating wave equations.
The following WebGL implementation simulates these equations on the GPU using fragment shaders. Use the controls below to pause/resume the simulation, adjust contrast, and visualize the PML absorption regions (red for , green for ).
Although it’s numerically stable, there are some high-frequency oscillations, indicating that we probably didn’t make the best choice of numerical scheme for this wave equation. There’s plenty more to explore in designing numerical schemes for PDEs, so we’ll ignore that for now.
Note also that we add a slider for the PML width and exponent. As a baseline, we ease the strength of from 0 to 1 over the width of the layer, then tack on an exponent on top of that so that we can ease it quadratically or even just set it (nearly) constant as the exponent approaches zero. In a perfect world, the layer is perfectly matched and we could just let it be constant. However, we’ve chosen a second order scheme with significant truncation error, which means that the impedance of the absorbing layer is not actually perfectly matched and still introduces spurious reflections. Easing the strength of the layer more gently reduces these reflections, though it also requires a thicker absorbing layer.
While the resulting equations are somewhat involved, PML provides a remarkably effective solution to the absorbing boundary problem. I’m a little uncertain how to deal with irregular domains with more complicated boundaries, but that’s an exercise for the future!
This derivation extends the technique from Steven G. Johnson’s PML notes by applying the PML mapping to both x and y boundaries simultaneously.