[ParallelStencil] Wave Equation

Specific Domains Modelling & Simulations
,
1

I was looking at ParallelStencil.jl, and trying to just use it to directly solve the second order form of the wave equation. I noticed in the 2D acoustic example that it is broken up into two first order equations, and I didn’t see any other examples where the PDE was kept in second order form. So I was wondering if this is a requirement of the library? Or if one is able to keep their PDEs in second order form?

2

Hi @apattyn, thanks for your message and interest in ParallelStencil.

The 2D acoustic wave example featured in the miniapps is indeed in the velocity-pressure formulation. However, the “order” in which you implement an equation is independent from ParallelStencil. It is the user’s responsibility to discretise the PDE accordingly to their need.

You could implement following in ParallelStencil, modifying acoustic2D.jl,

# define the second order kernel
@parallel function compute_P!(Pn, P, Po, dt, ρ, K, dx, dy)
    @inn(Pn) = 2.0*@inn(P) - @inn(Po) + 
               dt*dt*K/ρ*(@d2_xi(P)/dx/dx + @d2_yi(P)/dy/dy)
    return
end
# [...]
# Time loop
for it = 1:nt
    # [...]
    @parallel compute_P!(Pn, P, Po, dt, ρ, K, dx, dy)
    Po .= P
    P  .= Pn
    t += dt
    # visualisation ...
end
# [...]

If you want an even more performance oriented implementation, consider using @parallel_indices (see here for usage) to define your kernel.

Hope this helps!

2 Likes
3

Thanks for your help @luraess. That does clear things up! Just so I am clear though, the reason the wave equation is formatted that way in your example is due to us wanting to solve for the next time step, so using a center finite difference scheme we have:

\frac{d^2P}{dt^2} \approx \frac{P(t+dt)-2P(t)+P(t-dt)}{dt^2}

Therefore, at the next time step our solution is:

P(t+dt) \approx 2P(t) - P(t-dt) + dt^2(c^2\nabla^2P)

Where c is the sound speed of the medium and \nabla^2P would be defined by @d2_xi(P)/dx^2 + ... in the code.

1 Like