Issue solving 2D wave equation with ModelingToolkit.jl

Issue

I was looking at using ModelingToolkit.jl and was playing around by trying to solve the 2D 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

Code

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, 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 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

Matlab Results

11-Jun-2021-01-05-28-kspaceFirst

Matlab Code

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{:});

1 Like

The 2D wave equation is a perfect candidate for GPU-based stencil calculations :slight_smile:
You may be interested in having a look at this 2D wave equation example solved using ParallelStencil.jl 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.

3 Likes

@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 it is indicated the D(D(t)) is equivalent to syntax Differential(t)^2

1 Like

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

Thanks for that link :slightly_smiling_face: 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!

1 Like

I also converted the NeuralODE.jl example 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

Converted NeuralODE.jl Example

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)

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.

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):

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

Check Chris video here on solving PDEs.

1 Like

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

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!

Specifically, upwinding has not merged yet:

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 :sweat_smile:.

(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 ODE Solvers · DifferentialEquations.jl). The coming DiffEqOperators interface is mostly focused on nonlinear and it’s not complete right now.

3 Likes

Wait, where do you see this?

1 Like

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.

I was referring to Automatic Transformation of Nth Order ODEs to 1st Order ODEs · ModelingToolkit.jl

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?

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.

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).

1 Like

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.

1 Like

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: Issue with 1D wave equation example · Issue #327 · SciML/NeuralPDE.jl · GitHub

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

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

1 Like

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

1 Like