# \[ParallelStencil\] Wave Equation

**URL:** <https://discourse.julialang.org/t/parallelstencil-wave-equation/70012>\
**Category:** Modelling & Simulations\
**Tags:** pde, physics\
**Created:** [October 19, 2021, 5:08am UTC](https://discourse.julialang.org/t/parallelstencil-wave-equation/70012 "2021-10-19T05:08:11Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![apattyn](https://avatars.discourse-cdn.com/v4/letter/a/ee7513/32.png) [@apattyn](https://discourse.julialang.org/u/apattyn)\
**Post date:** [October 19, 2021, 5:08am UTC](https://discourse.julialang.org/t/parallelstencil-wave-equation/70012/1 "2021-10-19T05:08:11Z")

</div>

I was looking at [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl#questions-comments-and-discussions), and trying to just use it to directly solve the second order form of the wave equation. I noticed in the [2D acoustic example](https://github.com/omlins/ParallelStencil.jl/blob/main/miniapps/acoustic2D.jl) 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?

---

<div class="post-metadata">

**Author:** ![luraess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luraess/32/16189_2.png) [@luraess](https://discourse.julialang.org/u/luraess)\
**Post date:** [October 19, 2021, 5:52am UTC](https://discourse.julialang.org/t/parallelstencil-wave-equation/70012/2 "2021-10-19T05:52:26Z")

</div>

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](https://github.com/omlins/ParallelStencil.jl/blob/main/miniapps/acoustic2D.jl),

```julia
# 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](https://github.com/omlins/ParallelStencil.jl/blob/bb6aa77d5a0dba649da10534328f3a759fd6ef0d/miniapps/ThermalConvection2D.jl#L52) for usage) to define your kernel.

Hope this helps!

---

<div class="post-metadata">

**Author:** ![apattyn](https://avatars.discourse-cdn.com/v4/letter/a/ee7513/32.png) [@apattyn](https://discourse.julialang.org/u/apattyn)\
**Post date:** [October 19, 2021, 9:07pm UTC](https://discourse.julialang.org/t/parallelstencil-wave-equation/70012/3 "2021-10-19T21:07:03Z")

</div>

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.