I’m writting a research in which I heavily used NeuralPDE.jl, and I would like to showcase solution vs analytical, or classical simulations vs provided by the package.
In one of these instances, I’m getting an answer I don’t understand. Because, I believe it’s wrong. The solver is not obeying my initial condition.
eq: ∂u/∂t + u(x,t)(∂u/∂x) = 0,
bc: u(0,t) = (1/(2√π))exp((-1/2)(x-3)^2)
f(x) = (1/(2*√π))*exp((-1/2)*(x-3)^2)
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
import ModelingToolkit: Interval, infimum, supremum
@parameters t, x
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
#2D PDE
eq = Dt(u(t,x)) + u(t,x)*Dx(u(t,x)) ~ 0
# Initial and boundary conditions
bcs = [u(0,x) ~ f(x),
u(0,x) - f(x) ~ 0]
# Space and time domains
domains = [t ∈ Interval(0.0,20.0),
x ∈ Interval(0.0,10.0)]
# Discretization
nx=100;
dx = 10/(nx)
endTime = 20
nt = 100
dt = endTime/nt
The neural network:
# Neural network
dim = 2 # number of dimensions
chain = FastChain(FastDense(dim,25,Flux.σ),FastDense(25,25,Flux.σ),FastDense(25,1))
# Initial parameters of Neural network
initθ = Float64.(DiffEqFlux.initial_params(chain))
# Discretization
discretization = PhysicsInformedNN(chain,GridTraining(dx),init_params=initθ)
@named pde_system = PDESystem(eq,bcs,domains,[x,t],[u(x, t)])
prob = discretize(pde_system,discretization)
cb = function (p,l)
println("Current loss is: $l")
return false
end
opt = Optim.BFGS()
res = GalacticOptim.solve(prob,opt; cb = cb, maxiters=5000)
phi = discretization.phi
Finally, the plot,
using Plots; pyplot();
xs,ts = [infimum(d.domain):dx:supremum(d.domain) for (d,dx) in zip(domains,[dx,dt])]
u_predict = reshape([Array(phi([x, t], res.minimizer))[1] for x in xs for t in ts], length(xs), length(ts))
plot(xs,ts,u_predict',st=:surface, title="Burguer equation, PINN", xlabel="X", ylabel="Time", zlabel="U")
which gives me:
And, the initial condition is pretty different from what I specified:
p1=plot(xs, u_predict[1,:], title="Predicted")
p2=plot(xs, map(x -> f(x), xs), title="Analytical")
plot(p1,p2)
By the Euler Explicit method,
nxx=100;
delta_xx = 10/(nxx - 1)
xx = range(0, stop=delta_xx*(nxx-1), length=nxx)
endTime = 20
nt = 1000
delta_t = endTime/nt
t = range(0, stop=endTime, length=nt)
f(x) = (1/(2*√π))*exp((-1/2)*(x-3)^2)
v_zero = f.(xx)
v=zeros((nxx,nt+1))
v[:,1]=copy(v_zero)
for n in 1:nt
v[:,n+1] = copy(v[:,n])
for i in 2:nxx
v[i,n+1] = v[i,n] - v[i,n] * delta_t/delta_xx * (v[i,n] - v[i-1,n])
end
end
using Plots; pyplot()
xxs = collect(xx)
ts = collect(t)
plot(collect(xx),collect(t),v'[1:1000,1:100],st=:surface, title="Burguer equation", xlabel="X", ylabel="Time", zlabel="V")
.