Issue
I was looking at using ModelingToolkit.jl and was playing around by trying to solve the 2D wave equation.
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.
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
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{:});
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!
.
