# Issue solving 2D wave equation with ModelingToolkit.jl

**URL:** <https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq\
**Created:** [June 15, 2021, 6:05pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963 "2021-06-15T18:05:00Z")\
**Posts on this page:** 20\
**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:** [June 15, 2021, 6:05pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/1 "2021-06-15T18:05:00Z")

</div>

# Issue

I was looking at using ModelingToolkit.jl and was playing around by trying to solve the [2D wave equation](https://en.wikipedia.org/wiki/Wave_equation).

\frac{d^{2}u}{dt^{2}} = c^{2}\frac{d^{2}u}{dx^{2}}

However, the solution it outputs isn’t what I would expect. I had set the initial spatial conditions to zero, with a time varying boundary condition of u(t,0) = sin(t) \land \frac{du(t,0)}{dx} = cos(t). So I would’ve expected the wave to propagate from left to right, but it seems to become attenuated as it travels.

![1Dwave](https://global.discourse-cdn.com/julialang/original/3X/f/f/ffbf3cc3a5956ac1d16d9dbef9a549dc4a9fade7.gif)

## Code

```julia-auto
using Plots, DifferentialEquations, ModelingToolkit, DiffEqOperators

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dxx = Differential(x)^2
Dtt = Differential(t)^2
Dt = Differential(t)
Dx = Differential(x)

#2D PDE
C=1
eq = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))

# Initial and boundary conditions
bcs = [u(t,0) ~ sin(t),# for all t > 0
       Dx(u(t,0)) ~ cos(t) ]

# Space and time domains
domains = [t ∈ (0.0,50.0),
           x ∈ (0.0,10.0)]

# PDE system
pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x=>dx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
sol = solve(prob)

# Plot results
anim = @animate for i ∈ 1:length(sol.t)
plot(sol.u[i], label = "wave", ylims =[-1, 1])
end every 5

gif(anim, "1Dwave.gif", fps = 10)

```

# Additional Info

Seeing this I started to take a look at the [NeuralODE.jl wave example](https://neuralpde.sciml.ai/stable/pinn/wave/), and unless I am misreading it I wouldn’t expect the solution it is showing either. I tried to replicate the example using an [acoustic wave solver for matlab](http://www.k-wave.org/) which gave me different results.

The NeuralODE.jl example looks like the initial waveform just oscillates between ~0.2 and ~-0.2 without ever propagating outwards. Whereas the matlab wave solver shows the initial waveform moving to the boundaries where it is ultimately absorbed my the perfectly-matched-layer (since that toolbox is using a spectral method to solve the wave equation).

The question I suppose would be am I missing something? e.g. I defined the wave equation incorrectly. Or is ModelingToolkit.jl/ NeuralODE.jl giving an incorrect solution for the 2D wave equation?

## NeuralODE.jl Example Results

 ![](https://global.discourse-cdn.com/julialang/original/3X/c/6/c6438e7d0ff71b87482a99659096bea5a6b1f67d.png)

## Matlab Results

![11-Jun-2021-01-05-28-kspaceFirst](https://global.discourse-cdn.com/julialang/original/3X/0/2/025e99b32f49b0d76aa910db4dbb826eb8170c0f.gif)

## Matlab Code

```matlab
clearvars;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Lx =1; % number of grid points in the x (row) direction
dx = 0.1; % grid point spacing in the x direction [m]    
Nx = round(Lx/dx);

kgrid = kWaveGrid(Nx, dx);

% define the properties of the propagation medium
medium.sound_speed = 1; % [m/s] 

kgrid.makeTime(medium.sound_speed, [], 5);

% create initial pressure distribution using makeDisc
x = 0:dx:Lx-dx;
p0 = x.*(1-x);
source.p0=p0';
% define a centered cir
sensor.mask = zeros(1,Nx);

% run the simulation
args = {'PMLInside', false, 'RecordMovie', true};
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor,args{:});

```

---

<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:** [June 15, 2021, 6:37pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/2 "2021-06-15T18:37:15Z")

</div>

The 2D wave equation is a perfect candidate for GPU-based stencil calculations 🙂  
You may be interested in having a look at this [2D wave equation example](https://github.com/omlins/ParallelStencil.jl#acoustic-wave-2-d-app) solved using [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl#--parallelstenciljl-) on GPUs. It’s in the split V-P formulation but changing it to your needs would be straight-forward (only one kernel would be needed). Do not hesitate to ping me or @samo if you have questions.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 15, 2021, 7:11pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/3 "2021-06-15T19:11:44Z")

</div>

@apattyn, could you please indicate where the following syntax is defined for a second derivative: `Dtt = Differential(t)^2` ?  
Thanks.

(Edit): sorry my bad, in the [manual](https://mtk.sciml.ai/stable/tutorials/higher_order/) it is indicated the `D(D(t))` is equivalent to syntax `Differential(t)^2`

---

<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:** [June 15, 2021, 7:14pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/4 "2021-06-15T19:14:42Z")

</div>

> [@apattyn](#):
>
> ```julia
> eq = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))
> 
> ```

@rafael.guerra \frac{d^{2}u}{dt^{2}} is used above and defined in

> [@apattyn](#):
>
> ```julia
> @parameters t x
> @variables u(..)
> Dxx = Differential(x)^2
> Dtt = Differential(t)^2
> 
> ```

---

<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:** [June 15, 2021, 7:19pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/5 "2021-06-15T19:19:20Z")

</div>

Thanks for that link 🙂 That does seem to be an interesting library, but for now I was interested in learning ModelingToolbox.jl since it works nicely with the broader ecosystem of SciML. But I’ll keep that in mind for future projects!

---

<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:** [June 15, 2021, 7:48pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/6 "2021-06-15T19:48:19Z")

</div>

I also converted the [NeuralODE.jl example](https://neuralpde.sciml.ai/stable/pinn/wave/) into a purely ModelingToolkit.jl script and it looks like a solution to a diffusion problem rather than a wave problem.

This makes me believe its a problem with ModelingToolkit.jl or possibly a boundary condition. I say the latter may be an issue since for a numerical solution to the 2D wave equation I wrote using a spectral method, I saw a similar issue when the initial condition for \frac{du(0,x)}{dx} was improperly defined.

![1Dwave](https://global.discourse-cdn.com/julialang/original/3X/2/f/2f4d411c9871ba53bb63ce4b3ceb5c7d9499cb27.gif)

# Converted NeuralODE.jl Example

```julia
using Plots, DifferentialEquations, ModelingToolkit, DiffEqOperators

@parameters t, x
@variables u(..)
Dxx = Differential(x)^2
Dtt = Differential(t)^2
Dt = Differential(t)
Dx = Differential(x)

#2D PDE
C=1
eq = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))

# Initial and boundary conditions
bcs = [u(t,0) ~ 0.,# for all t > 0
       u(t,1) ~ 0.,# for all t > 0
       u(0,x) ~ x*(1. - x), #for all 0 < x < 1
       Dt(u(0,x)) ~ 0.0, #for all 0 < x < 1
]

# Space and time domains
domains = [t ∈ (0.0,1.0),
           x ∈ (0.0,1.0)]
# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x=>dx],t)

# PDE system
pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
sol = solve(prob)

# Plot results
anim = @animate for i ∈ 1:length(sol.t)
plot(sol.u[i], label = "wave", ylims =[-0.25, 0.25])
end every 5

gif(anim, "1Dwave.gif", fps = 10)

```

---

<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:** [June 15, 2021, 7:58pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/7 "2021-06-15T19:58:36Z")

</div>

Changing `Dt(u(0,x)) ~ 0.0` to `Dx(u(0,x)) ~ 1-2x` in the converted NeuralODE.jl example above doesn’t seem to have any affect on the solution.

---

<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:** [June 15, 2021, 8:25pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/8 "2021-06-15T20:25:41Z")

</div>

Taking a look at the docs you linked to they use `ode_order_lowering(...)` and I am wondering if that is required? It doesn’t seem to work for PDEs (which makes sense given the name):

```julia
pdesys = ode_order_lowering(pdesys)
ERROR: MethodError: no method matching ode_order_lowering(::PDESystem)
...
...

```

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 15, 2021, 8:52pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/9 "2021-06-15T20:52:52Z")

</div>

Check Chris video [here](https://youtu.be/okGybBmihOE) on solving PDEs.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 15, 2021, 9:43pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/10 "2021-06-15T21:43:43Z")

</div>

> [@apattyn](#):
>
> The NeuralODE.jl example looks like the initial waveform just oscillates between ~0.2 and ~-0.2 without ever propagating outwards. Whereas the matlab wave solver shows the initial waveform moving to the boundaries where it is ultimately absorbed my the perfectly-matched-layer (since that toolbox is using a spectral method to solve the wave equation).
> 
> The question I suppose would be am I missing something? e.g. I defined the wave equation incorrectly. Or is ModelingToolkit.jl/ NeuralODE.jl giving an incorrect solution for the 2D wave equation?

It looks like it’s the right solution to 2 digits. Plot the slides as an animation if you’re not sure.

> [@apattyn](#):
>
> ```julia
> MOLFiniteDifference([x=>dx],t)
> 
> ```

The documentation on this says it’s not complete.

> [Note: This uses a currently unreleased interface that is still a work in progress. Use at your own risk!](http://diffeqoperators.sciml.ai/dev/symbolic_tutorials/mol_heat/#Note:-This-uses-a-currently-unreleased-interface-that-is-still-a-work-in-progress.-Use-at-your-own-risk!)

Specifically, upwinding has not merged yet:

> <https://github.com/SciML/DiffEqOperators.jl/pull/390>

So yes, it’s as incorrect as the documentation page says it will be right now. We will remove that note when it’s released, but for right now, I mean, the docs said it’s not fully working yet 😅.

(Specifically, as of today, arbitrary even derivatives work fine, it’s the odd ones that need that to merge)

But anyways, @luraess mentioned, if you have a linear PDE like this then the most efficient thing to do would be to semi-discretize it with ParallelStencil.jl and then solve the ODEs using whichever method (likely one of the SSP integrators like [https://diffeq.sciml.ai/stable/solvers/ode\_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)](https://diffeq.sciml.ai/stable/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws))). The coming DiffEqOperators interface is mostly focused on nonlinear and it’s not complete right now.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 15, 2021, 9:44pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/11 "2021-06-15T21:44:13Z")

</div>

> [@apattyn](#):
>
> Taking a look at the docs you linked to they use `ode_order_lowering(...)` and I am wondering if that is required? It doesn’t seem to work for PDEs (which makes sense given the name):

Wait, where do you see this?

---

<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:** [June 15, 2021, 10:23pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/12 "2021-06-15T22:23:05Z")

</div>

> [@ChrisRackauckas](#):
>
> It looks like it’s the right solution to 2 digits. Plot the slides as an animation if you’re not sure

Hey Chris, thanks for looking at this. Looking at the NeuralODE.jl example the x-axis is time and the y-axis is the spatial dim. Based on that it doesn’t match what I was seeing in the Matlab wave-solver. We get a negative solution @ t=1 whereas in the matlab wave-solver it is always positive and propagates outwards rather than up and down.

Edit: To be specific the NeuralODE.jl solution presented isn’t matching what I see using a matlab wave-solver and its not what I would expect.

---

<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:** [June 15, 2021, 10:24pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/13 "2021-06-15T22:24:03Z")

</div>

I was referring to [https://mtk.sciml.ai/stable/tutorials/higher\_order/](https://mtk.sciml.ai/stable/tutorials/higher_order/)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/b/0b5292734f88e0c5a8d4daff4a3fdc8c58406cbf.png)

---

<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:** [June 16, 2021, 12:06am UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/14 "2021-06-16T00:06:28Z")

</div>

Based on the recommendation from Chris I’ll take a look at ParallelStencil.jl. A couple of questions before I make a new thread though. Does it work with AMD GPUs for acceleration? Then are your slack channels bridged to gitter/matrix?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 16, 2021, 6:42am UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/15 "2021-06-16T06:42:10Z")

</div>

Just one side remark: the title and posts all refer to 2D, but the equations and code displayed, with a single space variable, treat the 1D wave equation case.

---

<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:** [June 16, 2021, 7:05am UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/16 "2021-06-16T07:05:51Z")

</div>

> Does it work with AMD GPUs for acceleration?

These are coming steps but for now we have full support for Base.Threads and CUDA.jl.

> Then are your slack channels bridged to gitter/matrix?

No ParallelStencil specific Slack channel currently exists. However, we recommend posting in `#distrbuted` or `#gpu` or DM Sam or myself (and I don’t know whether those Slack channels are bridged).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 16, 2021, 8:51am UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/17 "2021-06-16T08:51:36Z")

</div>

Can you open an issue on NeuralPDE.jl with this example to look deeper into? I don’t get how it matches the analytical solution plot there and says its a low error, but also is incorrect.

---

<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:** [June 16, 2021, 5:22pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/18 "2021-06-16T17:22:00Z")

</div>

Will do. I was confused by that as well, but considering other wave-solvers are giving different results there may be an issue where the analytical solution given isn’t for the described problem.

Edit: [https://github.com/SciML/NeuralPDE.jl/issues/327](https://github.com/SciML/NeuralPDE.jl/issues/327)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [June 17, 2021, 9:46am UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/19 "2021-06-17T09:46:33Z")

</div>

Tried to adapt some “simple” Fortran code from a course I took to produce the result below, where sinusoidal boundary condition is in the center. Is it what you would expect to get?

![FiniteDifference_1D_wave_equation_sin_bc](https://global.discourse-cdn.com/julialang/original/3X/1/0/104ff15595549766e370d6245a9617931aac1edc.gif)

_NB: will post the code here after further review and comparison to analytical solution._

---

<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:** [June 17, 2021, 4:30pm UTC](https://discourse.julialang.org/t/issue-solving-2d-wave-equation-with-modelingtoolkit-jl/62963/20 "2021-06-17T16:30:25Z")

</div>

Yeah, I was able to get something similar using a Fourier-based method I had coded up a little while ago.

Since ModelingToolkit.jl is WIP I was going to just try and produce something similar with ParallelStencil.jl. I was interested in ModelingToolkit since it makes BCs, etc straight forward to handle and that it works with the SciML ecosystem.

![wave_numerical](https://global.discourse-cdn.com/julialang/original/3X/6/5/6503e2049a9017440b8c1af48ff332ed8e95208b.gif)