Physics-informed neural networks (PINNs) solver on Julia. GSoC 2020. Second evaluation.
Previous blog
Link to the first blog.
NeuralPDE.jl: Scientific Machine Learning for Partial Differential Equations.
NeuralPDE.jl is a solver package that consists of neural network solvers for partial differential equations using scientific machine learning (SciML) techniques such as physics-informed neural networks (PINNs) and deep BSDE solvers. This package utilizes deep neural networks and neural stochastic differential equations to solve high dimensional PDEs at a greatly reduced cost and greatly increased generality compared to classical methods.
Physics-Informed Neural Networks for automated PDE solving
Forward-Backwards Stochastic Differential Equation (FBSDE) methods for parabolic PDEs
Deep learning based solvers for optimal stopping time and Kolmogorov backwards equations
The goals of my project on GSoC is development and improvement PINNs solver for NeuralPDE.jl.
What is done to the second evaluation?
Higher-order derivatives (more 2nd) for the numeric derivative.
All basic types of boundary conditions like Neumann, Robin, and periodic boundary conditions.
System of PDEs
Add Readme and Docs
Add new tests.
Dependencies
"$VERSION"
using Pkg
Pkg.add("ModelingToolkit")
Pkg.add("DiffEqFlux")
Pkg.add("Flux")
Pkg.add("Plots")
Pkg.add("Test")
]add NeuralPDE#master
using ModelingToolkit
using NeuralPDE
using DiffEqFlux
using Flux
using Plots
using Test
Examples
Let's take a look at a couple of examples with new features that have been added.
ODE with high order derivative
Let's consider ODE with 3rd order derivative:
x θ
u(..)
Dxxx'''~x
Dx'~x
#ODE
eq = Dxxx(u(x,θ)) ~ cos(pi*x)
# Initial and boundary conditions
bcs = [u(0.,θ) ~ 0.0,
u(1.,θ) ~ cos(pi),
Dx(u(1.,θ)) ~ 1.0]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.1
discretization = NeuralPDE.PhysicsInformedNN(dx)
# Neural network and optimizer
opt = Flux.ADAM(0.01)
chain = FastChain(FastDense(1,8,Flux.σ),FastDense(8,1))
pde_system = PDESystem(eq,bcs,domains,[x],[u])
prob = NeuralPDE.discretize(pde_system,discretization)
alg = NeuralPDE.NNDE(chain,opt,autodiff=false)
phi,res = solve(prob,alg,verbose=true, maxiters=10000)
analytic_sol_func(x) = (π*x*(-x+(π^2)*(2*x-3)+1)-sin(π*x))/(π^3)
xs = [domain.domain.lower:dx/10:domain.domain.upper for domain in domains][1]
u_real = [analytic_sol_func(x) for x in xs]
u_predict = [first(phi(x,res.minimizer)) for x in xs]
u_predict ≈ u_real atol = 10.0
x_plot = collect(xs)
plot(x_plot ,u_real)
plot!(x_plot ,u_predict)
The linear PDE system
x, y, θ
u1(..), u2(..)
Dx'~x
Dy'~y
# System of pde
eqs = [Dx(u1(x,y,θ)) + 4*Dy(u2(x,y,θ)) ~ 0,
Dx(u2(x,y,θ)) + 9*Dy(u1(x,y,θ)) ~ 0]
# Initial and boundary conditions
bcs = [[u1(x,0,θ) ~ 2x, u2(x,0,θ) ~ 3x]]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.1
discretization = NeuralPDE.PhysicsInformedNN(dx)
# Neural network and optimizer
opt = Flux.ADAM(0.1)
chain = FastChain(FastDense(2,8,Flux.σ),FastDense(8,2))
pde_system = PDESystem(eqs,bcs,domains,[x,y],[u1,u2])
prob = NeuralPDE.discretize(pde_system,discretization)
alg = NeuralPDE.NNDE(chain,opt,autodiff=false)
phi,res = solve(prob,alg,verbose=true, maxiters=500)
analytic_sol_func(x,y) = [1/3*(6x - y), 1/2*(6x - y)]
xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
u_real = [analytic_sol_func(x,y) for x in xs for y in ys]
u_predict = [phi([x,y],res.minimizer) for x in xs for y in ys]
u_predict ≈ u_real atol = 10.0
2D wave equation with Neumann boundary condition
Let's solve this 2d wave equation:
x, t, θ
u(..)
Dxx''~x
Dtt''~t
Dt'~t
#2D PDE
C=1
eq = Dtt(u(x,t,θ)) ~ C^2*Dxx(u(x,t,θ))
# Initial and boundary conditions
bcs = [u(0,t,θ) ~ 0.,# for all t > 0
u(1,t,θ) ~ 0.,# for all t > 0
u(x,0,θ) ~ x*(1. - x), #for all 0 < x < 1
Dt(u(x,0,θ)) ~ 0. ] #for all 0 < x < 1]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
t ∈ IntervalDomain(0.0,1.0)]
# Discretization
dx = 0.1
discretization = NeuralPDE.PhysicsInformedNN(dx)
# Neural network and optimizer
opt = Flux.ADAM(0.02)
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
pde_system = PDESystem(eq,bcs,domains,[x,t],[u])
prob = NeuralPDE.discretize(pde_system,discretization)
alg = NeuralPDE.NNDE(chain,opt,autodiff=false)
phi,res = solve(prob,alg,verbose=true, maxiters=5000)
xs,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
analytic_sol_func(x,t) = sum([(8/(k^3*pi^3)) * sin(k*pi*x)*cos(C*k*pi*t) for k in 1:2:50000])
u_predict = reshape([first(phi([x,t],res.minimizer)) for x in xs for t in ts],(length(xs),length(ts)))
u_real = reshape([analytic_sol_func(x,t) for x in xs for t in ts], (length(xs),length(ts)))
u_predict ≈ u_real atol = 10.0
diff_u = abs.(u_predict .- u_real)
p1 = plot(xs, ts, u_real, linetype=:contourf,title = "analytic");
p2 =plot(xs, ts, u_predict, linetype=:contourf,title = "predict");
p3 = plot(xs, ts, diff_u,linetype=:contourf,title = "error");
plot(p1,p2,p3)
More examples you can find in docs or tests of NeuralPDE.
Future Goals
Setting boundary conditions of arbitrary geometry (in spherical and cartesian coordinates).
The ability to add additional conditions on the solution for the loss function (like the normalization condition)
Instead of using the entire set of points for training, in each optimization iteration, we could select randomly different residual points which will
significantly reduce training time.
Support automatic differentiation(AD)
Investigate physics-informed neural network (PINN) training strategies
Expand the test suite.
High dimensionality
Combined NN solution with grid numerical methods
Estimation problem, data-driven solution
At the moment, high-order derivatives, systems of PDE with a large number of equations, and high-dimensional equations still remain a challenge.