📂Numerical Analysis

Mur Absorbing Boundary Condition

Overview

Methods for numerically computing wave propagation include the finite difference method (FDM) and the $k$-space method. These methods assume that waves propagate to infinity, but in an actual simulation, waves propagate within a finite grid. As a result, the problem of waves reflecting at the boundaries arises. To resolve this, one can either set the grid to be much larger than the region one actually wants to simulate, or impose an appropriate boundary condition so that reflection of the waves does not occur. Such a boundary condition is called an absorbing boundary condition (ABC).

Formula

Consider the following wave equation.

$$ \dfrac{\partial ^{2}}{\partial t^{2}} u(x,t) = c^{2} \dfrac{\partial ^{2}}{\partial x^{2}} u(x,t) \quad \text{for } (x, t) \in \mathbb{R} \times [0, \infty) $$

Suppose we simulate this on the finite grid space $[0, 1] \times [0, T]$. Dividing the $x$-axis space into $M$ parts gives $\Delta x = 1/M$, and dividing time into $N$ parts gives $\Delta t = T/N$. Let the discretization of the solution $u$ of the wave equation be $U_{i}^{n} = u(i \ast \Delta x, n \ast \Delta t)$. Here $i = 0, 1, \cdots, M$ and $n = 0, 1, \cdots, N$. The absorbing boundary conditions are as follows.

$$ U_{0}^{n+1} = U_{1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{1}^{n+1} - U_{0}^{n} \right) $$

$$ U_{M}^{n+1} = U_{M-1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{M-1}^{n+1} - U_{M}^{n} \right) $$

Explanation

The results of applying the absorbing boundary conditions in one and two dimensions are as follows.

Derivation¹

Let us rearrange the wave equation as follows.

$$ \begin{align*} && \dfrac{\partial ^{2}}{\partial t^{2}} u(x,t) &= c^{2} \dfrac{\partial ^{2}}{\partial x^{2}} u(x,t) \\ \implies && \dfrac{\partial ^{2}}{\partial t^{2}} u(x,t) - c^{2} \dfrac{\partial ^{2}}{\partial x^{2}} u(x,t) &= 0 \\ \implies && \left(\dfrac{\partial }{\partial t} - c \dfrac{\partial }{\partial x} \right) \left(\dfrac{\partial }{\partial t} + c \dfrac{\partial }{\partial x} \right) u(x,t) &= 0 \\ \end{align*} $$

Here, the solution $u(x,t)$ of the wave equation satisfies one of the following two equations.

$$ \left(\dfrac{\partial }{\partial t} - c \dfrac{\partial }{\partial x} \right) u(x,t) = 0 \tag{1} $$

$$ \left(\dfrac{\partial }{\partial t} + c \dfrac{\partial }{\partial x} \right) u(x,t) = 0 \tag{2} $$

The $u_{+}(x, t) = f(x + ct)$ satisfying $(1)$ is a wave moving in the $-x$ direction. Conversely, the $u_{-}(x, t) = g(x - ct)$ satisfying $(2)$ is a wave moving in the $+x$ direction. At the boundary $x = 0$, only waves going to the left should exist for $U$, so $U$ must satisfy $(1)$. Computing the derivative of $u$ with respect to $t$ at $x = 0$,

$$ \dfrac{\partial u}{\partial t}|_{x = 0} \approx \dfrac{U_{0}^{n+1} - U_{0}^{n}}{\Delta t} \approx \dfrac{\dfrac{U_{0}^{n+1} + U_{1}^{n+1}}{2} - \dfrac{U_{0}^{n} + U_{1}^{n}}{2}}{\Delta t} $$

In the second approximation, setting the function value as the average of two adjacent points is for a better approximation². Likewise, the derivative with respect to $x$ is obtained as follows.

$$ c\dfrac{\partial u}{\partial x}|_{x = 0} \approx c \dfrac{U_{1}^{n} - U_{0}^{n}}{\Delta x} \approx c\dfrac{\dfrac{U_{1}^{n+1} + U_{1}^{n}}{2} - \dfrac{U_{0}^{n+1} + U_{0}^{n}}{2}}{\Delta x} $$

Since $U$ must satisfy $(1)$ at $x = 0$,

$$ \dfrac{\dfrac{U_{0}^{n+1} + U_{1}^{n+1}}{2} - \dfrac{U_{0}^{n} + U_{1}^{n}}{2}}{\Delta t} - c\dfrac{\dfrac{U_{1}^{n+1} + U_{1}^{n}}{2} - \dfrac{U_{0}^{n+1} + U_{0}^{n}}{2}}{\Delta x} = 0 $$

$$ \implies \dfrac{ U_{0}^{n+1} + U_{1}^{n+1} - U_{0}^{n} - U_{1}^{n}}{2 \Delta t} - c\dfrac{ U_{1}^{n+1} + U_{1}^{n} - U_{0}^{n+1} - U_{0}^{n}}{2 \Delta x} = 0 $$

$$ \implies U_{0}^{n+1} + U_{1}^{n+1} - U_{0}^{n} - U_{1}^{n} - \dfrac{c\Delta t}{\Delta x} \left( U_{1}^{n+1} + U_{1}^{n} - U_{0}^{n+1} - U_{0}^{n} \right) = 0 $$

$$ \implies \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{0}^{n+1} = \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{1}^{n} - (U_{1}^{n+1} - U_{0}^{n}) + \dfrac{c\Delta t}{\Delta x} \left( U_{1}^{n+1} - U_{0}^{n} \right) $$

$$ \implies \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{0}^{n+1} = \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{1}^{n} + \left(\dfrac{c \Delta t}{\Delta x} - 1 \right) \left( U_{1}^{n+1} - U_{0}^{n} \right) $$

$$ \begin{align*} \implies U_{0}^{n+1} &= U_{1}^{n} + \dfrac{\left(\dfrac{c \Delta t}{\Delta x} - 1 \right)}{\left(\dfrac{c \Delta t}{\Delta x} + 1 \right)} \left( U_{1}^{n+1} - U_{0}^{n} \right) \\ &= U_{1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{1}^{n+1} - U_{0}^{n} \right) \end{align*} $$

Likewise, the absorbing boundary condition at $x=1$ is as follows.

$$ U_{M}^{n+1} = U_{M-1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{M-1}^{n+1} - U_{M}^{n} \right) $$

■

Code

The code that simulates the absorbing boundary condition for the one-dimensional wave equation is as follows.

using Plots

x = LinRange(-1, 1, 300)
Δx = x[2] - x[1]

t = LinRange(0, 4, 600)
Δt = t[2] - t[1]

μ = Δt / Δx
C = (μ - 1) / (μ + 1)

U = zeros(length(x), length(t))

U[:, 1] = exp.(-30 * x.^2)
U[:, 2] = exp.(-30 * x.^2)

for i ∈ 3:length(t) # without ABCs
    U[2:end-1, i] = (μ^2)*U[3:end, i-1] + 2(1-μ^2)*U[2:end-1, i-1] + (μ^2)*U[1:end-2, i-1] - U[2:end-1, i-2]
end

V = zeros(length(x), length(t))
V[:, 1] = exp.(-30 * x.^2)
V[:, 2] = exp.(-30 * x.^2)

for i ∈ 3:length(t) # with ABCs
    V[2:end-1, i] = (μ^2)*V[3:end, i-1] + 2(1-μ^2)*V[2:end-1, i-1] + (μ^2)*V[1:end-2, i-1] - V[2:end-1, i-2]
    V[1, i] = V[2, i-1] + C*(V[2, i] - V[1, i-1]) # ABCs
    V[end, i] = V[end-1, i-1] + C*(V[end-1, i] - V[end, i-1]) # ABCs
end

anim = @animate for i ∈ 1:2:Int(length(t)/3)
    p1 = plot(x, U[:, i], xlims=(-1, 1), ylims=(-1,1), legend=false, dpi=300, title = "w/o ABCs")
    p2 = plot(x, V[:, i], xlims=(-1, 1), ylims=(-1,1), legend=false, dpi=300, title = "w/ ABCs")
    plot(p1, p2, layout=(2,1))
end

gif(anim, fps=40)

The code that simulates the absorbing boundary condition in two dimensions is as follows.

using Plots

# 파라미터 설정
nx, ny = 200, 200  # 격자 크기
nt = 400  # 시간 스텝 수
Lx, Ly = 10.0, 10.0  # 공간 영역 크기
T = 10.0  # 총 시간
c = 1.0  # 파동 속도

dx = Lx / (nx - 1)
dy = Ly / (ny - 1)
dt = T / (nt - 1)
κ = c * dt / dx

# 초기 조건 설정
u = zeros(Float64, nx, ny, nt)
v = zeros(Float64, nx, ny, nt)

for i ∈ 1:nx, j ∈ 1:ny
    if 13^2 < (i-160)^2 + (j-160)^2 < 15^2
        u[i, j, 1] = 1
        u[i, j, 2] = 1

        v[i, j, 1] = 1
        v[i, j, 2] = 1
    end
    if (i-40)^2 + (j-40)^2 < 4^2
        u[i, j, 1] = 1
        u[i, j, 2] = 1
        
        v[i, j, 1] = 1
        v[i, j, 2] = 1
    end
end

for n in 2:(nt - 1)
    u_xx = (u[3:end, 2:end-1, n] - 2 * u[2:end-1, 2:end-1, n] + u[1:end-2, 2:end-1, n]) / dx^2
    u_yy = (u[2:end-1, 3:end, n] - 2 * u[2:end-1, 2:end-1, n] + u[2:end-1, 1:end-2, n]) / dy^2
    u[2:end-1, 2:end-1, n + 1] = 2 * u[2:end-1, 2:end-1, n] - u[2:end-1, 2:end-1, n - 1] + c^2 * dt^2 * (u_xx + u_yy)

    v_xx = (v[3:end, 2:end-1, n] - 2 * v[2:end-1, 2:end-1, n] + v[1:end-2, 2:end-1, n]) / dx^2
    v_yy = (v[2:end-1, 3:end, n] - 2 * v[2:end-1, 2:end-1, n] + v[2:end-1, 1:end-2, n]) / dy^2
    v[2:end-1, 2:end-1, n + 1] = 2 * v[2:end-1, 2:end-1, n] - v[2:end-1, 2:end-1, n - 1] + c^2 * dt^2 * (v_xx + v_yy)

    u[1, 2:end-1, n+1] = u[2, 2:end-1, n] + ((κ-1)/(κ+1))*(u[2, 2:end-1, n] - u[1, 2:end-1, n-1])
    u[end, 2:end-1, n+1] = u[end-1, 2:end-1, n] + ((κ-1)/(κ+1))*(u[end-1, 2:end-1, n] - u[end, 2:end-1, n-1])
    u[2:end-1, 1, n+1] = u[2:end-1, 2, n] + ((κ-1)/(κ+1))*(u[2:end-1, 2, n] - u[2:end-1, 1, n-1])
    u[2:end-1, end, n+1] = u[2:end-1, end-1, n] + ((κ-1)/(κ+1))*(u[2:end-1, end-1, n] - u[2:end-1, end, n-1])
end

b2r = cgrad([:blue, :white, :red])

# 애니메이션 생성
anim = @animate for n in 1:nt
    p1 = heatmap(u[:, :, n], c=b2r, clim=(-1,1), legend=false, xticks=false, yticks=false, framestyle=:box, ratio=1, dpi=300, title="w/ ABCs")
    p2 = heatmap(v[:, :, n], c=b2r, clim=(-1,1), legend=false, xticks=false, yticks=false, framestyle=:box, ratio=1, dpi=300, title="w/o ABCs")
    plot(p1, p2, layout=(1,2), size=(728, 728/2))
end

gif(anim, fps=50)