format markdown

sdesai1287 · Feb 12, 2023 · f09a1ca · f09a1ca
1 parent f5642a4
commit f09a1ca
Show file tree

Hide file tree

Showing 25 changed files with 747 additions and 657 deletions.
diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
@@ -1 +1,2 @@
-style = "sciml"
+style = "sciml"
+format_markdown = true
diff --git a/LICENSE.md b/LICENSE.md
@@ -19,4 +19,3 @@ The NeuralPDE.jl package is licensed under the MIT "Expat" License:
 > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 > SOFTWARE.
-> 
diff --git a/README.md b/README.md
@@ -7,7 +7,7 @@
 [![Build Status](https://github.com/SciML/NeuralPDE.jl/workflows/CI/badge.svg)](https://github.com/SciML/NeuralPDE.jl/actions?query=workflow%3ACI)
 [![Build status](https://badge.buildkite.com/fa31256f4b8a4f95fe5ab90c3bf4ef56055a2afe675435c182.svg?branch=master)](https://buildkite.com/julialang/neuralpde-dot-jl)
 
-[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
+[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
 [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
 
 NeuralPDE.jl is a solver package which consists of neural network solvers for
@@ -29,19 +29,19 @@ the documentation, which contains the unreleased features.
 
 ## Features
 
-- Physics-Informed Neural Networks for ODE, SDE, RODE, and PDE solving
-- Ability to define extra loss functions to mix xDE solving with data fitting (scientific machine learning)
-- Automated construction of Physics-Informed loss functions from a high level symbolic interface
-- Sophisticated techniques like quadrature training strategies, adaptive loss functions, and neural adapters
-  to accelerate training
-- Integrated logging suite for handling connections to TensorBoard
-- Handling of (partial) integro-differential equations and various stochastic equations
-- Specialized forms for solving `ODEProblem`s with neural networks
-- Compatability with [Flux.jl](https://docs.sciml.ai/Flux.jl/stable/) and [Lux.jl](https://docs.sciml.ai/Lux/stable/)
-  for all of the GPU-powered machine learning layers available from those libraries.
-- Compatability with [NeuralOperators.jl](https://docs.sciml.ai/NeuralOperators/stable/) for
-  mixing DeepONets and other neural operators (Fourier Neural Operators, Graph Neural Operators,
-  etc.) with physics-informed loss functions
+  - Physics-Informed Neural Networks for ODE, SDE, RODE, and PDE solving
+  - Ability to define extra loss functions to mix xDE solving with data fitting (scientific machine learning)
+  - Automated construction of Physics-Informed loss functions from a high level symbolic interface
+  - Sophisticated techniques like quadrature training strategies, adaptive loss functions, and neural adapters
+    to accelerate training
+  - Integrated logging suite for handling connections to TensorBoard
+  - Handling of (partial) integro-differential equations and various stochastic equations
+  - Specialized forms for solving `ODEProblem`s with neural networks
+  - Compatability with [Flux.jl](https://docs.sciml.ai/Flux.jl/stable/) and [Lux.jl](https://docs.sciml.ai/Lux/stable/)
+    for all of the GPU-powered machine learning layers available from those libraries.
+  - Compatability with [NeuralOperators.jl](https://docs.sciml.ai/NeuralOperators/stable/) for
+    mixing DeepONets and other neural operators (Fourier Neural Operators, Graph Neural Operators,
+    etc.) with physics-informed loss functions
 
 ## Example: Solving 2D Poisson Equation via Physics-Informed Neural Networks
 
@@ -55,52 +55,54 @@ Dxx = Differential(x)^2
 Dyy = Differential(y)^2
 
 # 2D PDE
-eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)
+eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ -sin(pi * x) * sin(pi * y)
 
 # Boundary conditions
-bcs = [u(0,y) ~ 0.0, u(1,y) ~ 0,
-       u(x,0) ~ 0.0, u(x,1) ~ 0]
+bcs = [u(0, y) ~ 0.0, u(1, y) ~ 0,
+    u(x, 0) ~ 0.0, u(x, 1) ~ 0]
 # Space and time domains
-domains = [x ∈ Interval(0.0,1.0),
-           y ∈ Interval(0.0,1.0)]
+domains = [x ∈ Interval(0.0, 1.0),
+    y ∈ Interval(0.0, 1.0)]
 # Discretization
 dx = 0.1
 
 # Neural network
 dim = 2 # number of dimensions
-chain = Lux.Chain(Dense(dim,16,Lux.σ),Dense(16,16,Flux.σ),Dense(16,1))
+chain = Lux.Chain(Dense(dim, 16, Lux.σ), Dense(16, 16, Flux.σ), Dense(16, 1))
 
 discretization = PhysicsInformedNN(chain, QuadratureTraining())
 
-@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
-prob = discretize(pde_system,discretization)
+@named pde_system = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])
+prob = discretize(pde_system, discretization)
 
-callback = function (p,l)
+callback = function (p, l)
     println("Current loss is: $l")
     return false
 end
 
-res = Optimization.solve(prob, ADAM(0.1); callback = callback, maxiters=4000)
-prob = remake(prob,u0=res.minimizer)
-res = Optimization.solve(prob, ADAM(0.01); callback = callback, maxiters=2000)
+res = Optimization.solve(prob, ADAM(0.1); callback = callback, maxiters = 4000)
+prob = remake(prob, u0 = res.minimizer)
+res = Optimization.solve(prob, ADAM(0.01); callback = callback, maxiters = 2000)
 phi = discretization.phi
 ```
 
 And some analysis:
 
 ```julia
-xs,ys = [infimum(d.domain):dx/10:supremum(d.domain) for d in domains]
-analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)
+xs, ys = [infimum(d.domain):(dx / 10):supremum(d.domain) for d in domains]
+analytic_sol_func(x, y) = (sin(pi * x) * sin(pi * y)) / (2pi^2)
 
-u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
-u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
+u_predict = reshape([first(phi([x, y], res.minimizer)) for x in xs for y in ys],
+                    (length(xs), length(ys)))
+u_real = reshape([analytic_sol_func(x, y) for x in xs for y in ys],
+                 (length(xs), length(ys)))
 diff_u = abs.(u_predict .- u_real)
 
 using Plots
-p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
-p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
-p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
-plot(p1,p2,p3)
+p1 = plot(xs, ys, u_real, linetype = :contourf, title = "analytic");
+p2 = plot(xs, ys, u_predict, linetype = :contourf, title = "predict");
+p3 = plot(xs, ys, diff_u, linetype = :contourf, title = "error");
+plot(p1, p2, p3)
 ```
 
 ![image](https://user-images.githubusercontent.com/12683885/90962648-2db35980-e4ba-11ea-8e58-f4f07c77bcb9.png)

diff --git a/docs/src/developer/debugging.md b/docs/src/developer/debugging.md
@@ -15,41 +15,40 @@ Dxx = Differential(x)^2
 Dtt = Differential(t)^2
 Dt = Differential(t)
 #2D PDE
-C=1
-eq  = Dtt(u(x,t)) ~ C^2*Dxx(u(x,t))
+C = 1
+eq = Dtt(u(x, t)) ~ C^2 * Dxx(u(x, t))
 
 # Initial and boundary conditions
-bcs = [u(0,t) ~ 0.,
-       u(1,t) ~ 0.,
-       u(x,0) ~ x*(1. - x),
-       Dt(u(x,0)) ~ 0. ]
+bcs = [u(0, t) ~ 0.0,
+    u(1, t) ~ 0.0,
+    u(x, 0) ~ x * (1.0 - x),
+    Dt(u(x, 0)) ~ 0.0]
 
 # Space and time domains
-domains = [x ∈ Interval(0.0,1.0),
-           t ∈ Interval(0.0,1.0)]
+domains = [x ∈ Interval(0.0, 1.0),
+    t ∈ Interval(0.0, 1.0)]
 
 # Neural network
-chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
+chain = FastChain(FastDense(2, 16, Flux.σ), FastDense(16, 16, Flux.σ), FastDense(16, 1))
 init_params = DiffEqFlux.initial_params(chain)
 
 eltypeθ = eltype(init_params)
 phi = NeuralPDE.get_phi(chain)
 derivative = NeuralPDE.get_numeric_derivative()
 
-u_ = (cord, θ, phi)->sum(phi(cord, θ))
+u_ = (cord, θ, phi) -> sum(phi(cord, θ))
 
-phi([1,2], init_params)
+phi([1, 2], init_params)
 
 phi_ = (p) -> phi(p, init_params)[1]
-dphi = Zygote.gradient(phi_,[1.,2.])
+dphi = Zygote.gradient(phi_, [1.0, 2.0])
 
-dphi1 = derivative(phi,u_,[1.,2.],[[ 0.0049215667, 0.0]],1,init_params)
-dphi2 = derivative(phi,u_,[1.,2.],[[0.0,  0.0049215667]],1,init_params)
-isapprox(dphi[1][1], dphi1, atol=1e-8)
-isapprox(dphi[1][2], dphi2, atol=1e-8)
+dphi1 = derivative(phi, u_, [1.0, 2.0], [[0.0049215667, 0.0]], 1, init_params)
+dphi2 = derivative(phi, u_, [1.0, 2.0], [[0.0, 0.0049215667]], 1, init_params)
+isapprox(dphi[1][1], dphi1, atol = 1e-8)
+isapprox(dphi[1][2], dphi2, atol = 1e-8)
 
-
-indvars = [x,t]
+indvars = [x, t]
 depvars = [u(x, t)]
 dict_depvars_input = Dict(:u => [:x, :t])
 dim = length(domains)
@@ -58,8 +57,12 @@ multioutput = chain isa AbstractArray
 strategy = NeuralPDE.GridTraining(dx)
 integral = NeuralPDE.get_numeric_integral(strategy, indvars, multioutput, chain, derivative)
 
-_pde_loss_function = NeuralPDE.build_loss_function(eq,indvars,depvars,phi,derivative,integral,multioutput,init_params,strategy)
+_pde_loss_function = NeuralPDE.build_loss_function(eq, indvars, depvars, phi, derivative,
+                                                   integral, multioutput, init_params,
+                                                   strategy)
+```
 
+```
 julia> expr_pde_loss_function = NeuralPDE.build_symbolic_loss_function(eq,indvars,depvars,dict_depvars_input,phi,derivative,integral,multioutput,init_params,strategy)
 
 :((cord, var"##θ#529", phi, derivative, integral, u)->begin
@@ -76,11 +79,17 @@ julia> bc_indvars = NeuralPDE.get_variables(bcs,indvars,depvars)
  [:t]
  [:x]
  [:x]
+```
 
-_bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars,
-                                                     phi,derivative,integral,multioutput,init_params,strategy,
-                                                     bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]
+```julia
+_bc_loss_functions = [NeuralPDE.build_loss_function(bc, indvars, depvars,
+                                                    phi, derivative, integral, multioutput,
+                                                    init_params, strategy,
+                                                    bc_indvars = bc_indvar)
+                      for (bc, bc_indvar) in zip(bcs, bc_indvars)]
+```
 
+```
 julia> expr_bc_loss_functions = [NeuralPDE.build_symbolic_loss_function(bc,indvars,depvars,dict_depvars_input,
                                                                         phi,derivative,integral,multioutput,init_params,strategy,
                                                                         bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]
@@ -113,10 +122,15 @@ julia> expr_bc_loss_functions = [NeuralPDE.build_symbolic_loss_function(bc,indva
               end
           end
       end)
+```
 
-train_sets = NeuralPDE.generate_training_sets(domains,dx,[eq],bcs,eltypeθ,indvars,depvars)
-pde_train_set,bcs_train_set = train_sets
+```julia
+train_sets = NeuralPDE.generate_training_sets(domains, dx, [eq], bcs, eltypeθ, indvars,
+                                              depvars)
+pde_train_set, bcs_train_set = train_sets
+```
 
+```
 julia> pde_train_set
 1-element Array{Array{Float32,2},1}:
  [0.1 0.2 … 0.8 0.9; 0.1 0.1 … 1.0 1.0]
@@ -128,10 +142,14 @@ julia> bcs_train_set
  [1.0 1.0 … 1.0 1.0; 0.0 0.1 … 0.9 1.0]
  [0.0 0.1 … 0.9 1.0; 0.0 0.0 … 0.0 0.0]
  [0.0 0.1 … 0.9 1.0; 0.0 0.0 … 0.0 0.0]
+```
 
+```julia
+pde_bounds, bcs_bounds = NeuralPDE.get_bounds(domains, [eq], bcs, eltypeθ, indvars, depvars,
+                                              NeuralPDE.StochasticTraining(100))
+```
 
-pde_bounds, bcs_bounds = NeuralPDE.get_bounds(domains,[eq],bcs,eltypeθ,indvars,depvars,NeuralPDE.StochasticTraining(100))
-
+```
 julia> pde_bounds
 1-element Vector{Vector{Any}}:
  [Float32[0.01, 0.99], Float32[0.01, 0.99]]
@@ -142,13 +160,14 @@ julia> bcs_bounds
  [1, Float32[0.0, 1.0]]
  [Float32[0.0, 1.0], 0]
  [Float32[0.0, 1.0], 0]
+```
 
-discretization = NeuralPDE.PhysicsInformedNN(chain,strategy)
-
-@named pde_system = PDESystem(eq,bcs,domains,indvars,depvars)
-prob = NeuralPDE.discretize(pde_system,discretization)
+```julia
+discretization = NeuralPDE.PhysicsInformedNN(chain, strategy)
 
-expr_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)
-expr_pde_loss_function , expr_bc_loss_functions = expr_prob
+@named pde_system = PDESystem(eq, bcs, domains, indvars, depvars)
+prob = NeuralPDE.discretize(pde_system, discretization)
 
+expr_prob = NeuralPDE.symbolic_discretize(pde_system, discretization)
+expr_pde_loss_function, expr_bc_loss_functions = expr_prob
 ```
diff --git a/docs/src/examples/3rd.md b/docs/src/examples/3rd.md
@@ -25,29 +25,29 @@ import ModelingToolkit: Interval, infimum, supremum
 Dxxx = Differential(x)^3
 Dx = Differential(x)
 # ODE
-eq = Dxxx(u(x)) ~ cos(pi*x)
+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]
+bcs = [u(0.0) ~ 0.0,
+    u(1.0) ~ cos(pi),
+    Dx(u(1.0)) ~ 1.0]
 
 # Space and time domains
-domains = [x ∈ Interval(0.0,1.0)]
+domains = [x ∈ Interval(0.0, 1.0)]
 
 # Neural network
-chain = Lux.Chain(Dense(1,8,Lux.σ),Dense(8,1))
+chain = Lux.Chain(Dense(1, 8, Lux.σ), Dense(8, 1))
 
 discretization = PhysicsInformedNN(chain, QuasiRandomTraining(20))
-@named pde_system = PDESystem(eq,bcs,domains,[x],[u(x)])
-prob = discretize(pde_system,discretization)
+@named pde_system = PDESystem(eq, bcs, domains, [x], [u(x)])
+prob = discretize(pde_system, discretization)
 
-callback = function (p,l)
+callback = function (p, l)
     println("Current loss is: $l")
     return false
 end
 
-res = Optimization.solve(prob, Adam(0.01); callback = callback, maxiters=2000)
+res = Optimization.solve(prob, Adam(0.01); callback = callback, maxiters = 2000)
 phi = discretization.phi
 ```
 
@@ -56,16 +56,16 @@ We can plot the predicted solution of the ODE and its analytical solution.
 ```@example 3rdDerivative
 using Plots
 
-analytic_sol_func(x) = (π*x*(-x+(π^2)*(2*x-3)+1)-sin(π*x))/(π^3)
+analytic_sol_func(x) = (π * x * (-x + (π^2) * (2 * x - 3) + 1) - sin(π * x)) / (π^3)
 
 dx = 0.05
-xs = [infimum(d.domain):dx/10:supremum(d.domain) for d in domains][1]
-u_real  = [analytic_sol_func(x) for x in xs]
-u_predict  = [first(phi(x,res.u)) for x in xs]
+xs = [infimum(d.domain):(dx / 10):supremum(d.domain) for d in domains][1]
+u_real = [analytic_sol_func(x) for x in xs]
+u_predict = [first(phi(x, res.u)) for x in xs]
 
 x_plot = collect(xs)
-plot(x_plot ,u_real,title = "real")
-plot!(x_plot ,u_predict,title = "predict")
+plot(x_plot, u_real, title = "real")
+plot!(x_plot, u_predict, title = "predict")
 ```
 
 ![hodeplot](https://user-images.githubusercontent.com/12683885/90276340-69bc3e00-de6c-11ea-89a7-7d291123a38b.png)
diff --git a/docs/src/examples/heterogeneous.md b/docs/src/examples/heterogeneous.md
@@ -19,29 +19,31 @@ Dx = Differential(x)
 Dy = Differential(y)
 
 # 2D PDE
-eq  = p(x) + q(y) + Dx(r(x, y)) + Dy(s(y, x)) ~ 0
+eq = p(x) + q(y) + Dx(r(x, y)) + Dy(s(y, x)) ~ 0
 
 # Initial and boundary conditions
-bcs = [p(1) ~ 0.f0, q(-1) ~ 0.0f0,
-       r(x, -1) ~ 0.f0, r(1, y) ~ 0.0f0,
-       s(y, 1) ~ 0.0f0, s(-1, x) ~ 0.0f0]
+bcs = [p(1) ~ 0.0f0, q(-1) ~ 0.0f0,
+    r(x, -1) ~ 0.0f0, r(1, y) ~ 0.0f0,
+    s(y, 1) ~ 0.0f0, s(-1, x) ~ 0.0f0]
 
 # Space and time domains
 domains = [x ∈ Interval(0.0, 1.0),
-           y ∈ Interval(0.0, 1.0)]
+    y ∈ Interval(0.0, 1.0)]
 
 numhid = 3
-chains = [[Lux.Chain(Dense(1, numhid, Lux.σ), Dense(numhid, numhid, Lux.σ), Dense(numhid, 1)) for i in 1:2];
-                        [Lux.Chain(Dense(2, numhid, Lux.σ), Dense(numhid, numhid, Lux.σ), Dense(numhid, 1)) for i in 1:2]]
+chains = [[Lux.Chain(Dense(1, numhid, Lux.σ), Dense(numhid, numhid, Lux.σ),
+                     Dense(numhid, 1)) for i in 1:2]
+          [Lux.Chain(Dense(2, numhid, Lux.σ), Dense(numhid, numhid, Lux.σ),
+                     Dense(numhid, 1)) for i in 1:2]]
 discretization = NeuralPDE.PhysicsInformedNN(chains, QuadratureTraining())
 
-@named pde_system = PDESystem(eq, bcs, domains, [x,y], [p(x), q(y), r(x, y), s(y, x)])
+@named pde_system = PDESystem(eq, bcs, domains, [x, y], [p(x), q(y), r(x, y), s(y, x)])
 prob = SciMLBase.discretize(pde_system, discretization)
 
-callback = function (p,l)
+callback = function (p, l)
     println("Current loss is: $l")
     return false
 end
 
-res = Optimization.solve(prob, BFGS(); callback = callback, maxiters=100)
+res = Optimization.solve(prob, BFGS(); callback = callback, maxiters = 100)
 ```