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 PkgPkg.add("ModelingToolkit")Pkg.add("DiffEqFlux")Pkg.add("Flux")Pkg.add("Plots")Pkg.add("Test")]add NeuralPDE#masterusing ModelingToolkitusing NeuralPDEusing DiffEqFluxusing Fluxusing Plotsusing TestExamples
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#ODEeq = Dxxx(u(x,θ)) ~ cos(pi*x)# Initial and boundary conditionsbcs = [u(0.,θ) ~ 0.0, u(1.,θ) ~ cos(pi), Dx(u(1.,θ)) ~ 1.0]# Space and time domainsdomains = [x ∈ IntervalDomain(0.0,1.0)]# Discretizationdx = 0.1discretization = NeuralPDE.PhysicsInformedNN(dx)# Neural network and optimizeropt = 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.0x_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 pdeeqs = [Dx(u1(x,y,θ)) + 4*Dy(u2(x,y,θ)) ~ 0, Dx(u2(x,y,θ)) + 9*Dy(u1(x,y,θ)) ~ 0]# Initial and boundary conditionsbcs = [[u1(x,0,θ) ~ 2x, u2(x,0,θ) ~ 3x]]# Space and time domainsdomains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]# Discretizationdx = 0.1discretization = NeuralPDE.PhysicsInformedNN(dx)# Neural network and optimizeropt = 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.02D wave equation with Neumann boundary condition
Let's solve this 2d wave equation:

x, t, θ u(..) Dxx''~x Dtt''~t Dt'~t#2D PDEC=1eq = Dtt(u(x,t,θ)) ~ C^2*Dxx(u(x,t,θ))# Initial and boundary conditionsbcs = [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 domainsdomains = [x ∈ IntervalDomain(0.0,1.0), t ∈ IntervalDomain(0.0,1.0)]# Discretizationdx = 0.1discretization = NeuralPDE.PhysicsInformedNN(dx)# Neural network and optimizeropt = 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.0diff_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.