acoustic_wave_animation.jl

# Equations:
# ∂P/∂t = -k(∂Vx/∂x + ∂Vy/∂y)
# ∂Vx/∂t = -1/ρ ∂P/∂x
# ∂Vy/∂t = -1/ρ ∂P/∂y

# => Update rules:
# Vx = Vx - Δt/ρ ΔP/Δx
# Vy = Vy - Δt/ρ ΔP/Δy
# P = P - Δt*k (Δ(Vx)/Δx + Δ(Vy)/Δy)

using Plots
using Printf

const PROBLEM_SIZE = 255

@views function acoustic_2D_animation()
    # physical parameters
    lx, ly = 40.0, 40.0  # domain extends
    k = 1.0 # bulk modulus
    ρ = 1.0 # density
    t = 0.0 # physical time
    # plotting
    nout = 20

    # spatial grid
    nx, ny = PROBLEM_SIZE, PROBLEM_SIZE
    dx, dy = lx / (nx - 1), ly / (ny - 1)
    xc, yc = (-lx / 2):dx:(lx / 2), (-ly / 2):dy:(ly / 2)

    # time discretization
    dt = min(dx, dy) / sqrt(k / ρ) / 4.1
    nt = 1000 # timesteps

    # initial condition
    P = [exp(-((ix - 1) * dx - 0.3 * lx)^2 - ((iy - 1) * dy - 0.5 * ly)^2) +
         exp(-((ix - 1) * dx - 0.7 * lx)^2 - ((iy - 1) * dy - 0.5 * ly)^2)
         for ix in 1:nx, iy in 1:ny]

    # preallocation
    Vx = zeros(nx - 1, ny - 2)
    Vy = zeros(nx - 2, ny - 1)

    # time loop
    anim = Animation(mktempdir(), String[])
    for it in 1:nt
        # -------- stencil kernel --------
        # Velocities
        Vx .= Vx .- dt ./ ρ .* diff(P[:, 2:(end - 1)], dims = 1) ./ dx
        Vy .= Vy .- dt ./ ρ .* diff(P[2:(end - 1), :], dims = 2) ./ dy
        # Pressure
        P[2:(end - 1), 2:(end - 1)] .= P[2:(end - 1), 2:(end - 1)] .-
                                       dt .* k .* (diff(Vx, dims = 1) ./ dx .+
                                        diff(Vy, dims = 2) ./ dy)
        # Time
        t += dt
        # --------------------------------

        # plotting
        if mod(it, nout) == 0
            heatmap(xc, yc, P', xlabel = "x", ylabel = "y", title = "Pressure, i=$it",
                    colormap = :viridis, clims = (0.0, 0.2))
            frame(anim)
        end
    end
    gif(anim, "acoustic_wave_animation.gif", fps = 15)
    return nothing
end

acoustic_2D_animation()