Issue with 1D wave equation example #327
Description
Issue
Looking at the 1D wave equation example, I am not sure if the presented solution is correct. I discussed this on a julia discourse thread. To briefly summarize the solution in that NeuralPDE.jl example isn't what I would expect to see, and it doesn't match the results from a matlab wave-solver when I try to replicate the results.
Then when I went to convert the NeuralPDE example to a "purely" ModelingToolkit version I got a solution that looked like it was from a diffusion problem. I did have a similar problem when I wrote a custom wave-solver using a spectral method, which ended up being due to improperly defined initial conditions (i.e. du(0,x)/dx = 0 instead of the derivative of the u(0,x)). So I exchanged Dt(u(0,x)) ~ 0. for Dx(u(0,x)) ~ 1-2x in the bcs without it having any effect.
Matlab Code
clearvars;
% =========================================================================
% SIMULATION
% =========================================================================
% create the computational grid
Lx =1;
dx = 0.1; % grid point spacing in the x direction [m]
Nx = round(Lx/dx); % number of grid points in the x (row) direction
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
x = 0:dx:Lx-dx;
p0 = x.*(1-x);
source.p0=p0';
% define a sensor
sensor.mask = zeros(1,Nx);
% run the simulation
args = {'PMLInside', false, 'RecordMovie', true};
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor,args{:});
Matlab Results
Keep in mind the solution below has a larger spatial dim than the NeuralPDE example, due to the fact that this wave solver is using a spectral method and has to have a perfectly-matched-layer to avoid "wrap around" effects.
Converted NeuralPDE 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)
Converted NeuralPDE Results
