SciML / NeuralPDE.jl

Physics-Informed Neural Networks (PINN) Solvers of (Partial) Differential Equations for Scientific Machine Learning (SciML) accelerated simulation
https://docs.sciml.ai/NeuralPDE/stable/
Other
991 stars 198 forks source link

Random noise as an initial condition #481

Open oashour opened 2 years ago

oashour commented 2 years ago

I am moving this over here from the Discourse as per Chris's request. A short summary is below.

Essentially I'm trying to use random noise as the initial condition of my simulation of the PDEs shown in that thread. Currently, I have

bcs = [c₁(0,x) ~ rand(),
       c₂(0,x) ~ rand(),
       c₁(t,0) ~ 0.,
       c₂(t,0) ~ 0.,
       c₁(t,30) ~ 0.,
       c₂(t,30) ~ 0.]

But this generates a single value as the initial condition instead of changing it every time as per Chris. I am not sure how to proceed, something was mentioned about solving over all initial conditions and using that as an extra dimension but I don't quite understand this, I'm very new to PINNs.

Thanks!

ChrisRackauckas commented 2 years ago

As I mentioned in the Discourse thread, I think the right thing to do is to train with the initial condition being a different independent variable and then sample across that dimension.

Maybe @KirillZubov can help get a demo of that together.

KirillZubov commented 2 years ago

@oashour you need @register_symbolic https://symbolics.juliasymbolics.org/dev/manual/functions/

@parameters t x
@variables c₁(..) c₂(..)
rand_f(x) = rand(size(x)...)
@register_symbolic rand_f(x)
Base.Broadcast.broadcasted(::typeof(rand_f), x) = rand_f(x)
vec_ = ones(1,10)
rand_f(vec_)
"""
julia> rand_f(vec_)
1×10 Matrix{Float64}:
 0.384208  0.0996561  0.990935  0.9898  0.479029  0.516312  0.635055  0.480348  0.789933  0.665574

"""

eq = c₁(t,x)+c₂(t,x) ~ 0
bcs = [c₁(0,x) ~ rand_f(x),
       c₂(0,x) ~ rand_f(x),
       c₁(t,0) ~ 0.,
       c₂(t,0) ~ 0.,
       c₁(t,30) ~ 0.,
       c₂(t,30) ~ 0.]
"""
output:
6-element Vector{Equation}:
 c₁(0, x) ~ rand_f(x)
 c₂(0, x) ~ rand_f(x)
 c₁(t, 0) ~ 0.0
 c₂(t, 0) ~ 0.0
 c₁(t, 30) ~ 0.0
 c₂(t, 30) ~ 0.0
"""

domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0)]

chains = [FastChain(FastDense(2,12,Flux.tanh),FastDense(12,12,Flux.tanh),FastDense(12,1)),
         FastChain(FastDense(2,12,Flux.tanh),FastDense(12,12,Flux.tanh),FastDense(12,1))]

initθ = map(c -> Float64.(c), DiffEqFlux.initial_params.(chains))
grid_strategy = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chains,
                                             grid_strategy;
                                             init_params = initθ)

@named pde_system = PDESystem(eq,bcs,domains,[t,x],[c₁(t,x),c₂(t,x)])
prob = NeuralPDE.discretize(pde_system,discretization)
sym_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)

prob.f(reduce(vcat, initθ), nothing)
res = GalacticOptim.solve(prob, ADAM(0.1); maxiters=500)

#two dim case
rand_f(x,y) = rand(size(vcat(x,y))...)
@register_symbolic rand_f(x,y)
Base.Broadcast.broadcasted(::typeof(rand_f), x,y) = rand_f(x,y)
rand_f.(ones(1,10),ones(1,10))
ChrisRackauckas commented 2 years ago

I'm not sure that's the right way to handle it. Why not handle the random factor as deterministic for the training and then sample on the results?

KirillZubov commented 2 years ago

@ChrisRackauckas do you mean something like this algorithm?

for every x in x_vector create a random vector - rand(N), N - size of vector
xi ~ rand(N) 
vector of results of initial conditions - c₁(0,xi) ~ rand(N) for each xi in x_vector
loss_fuctions  =  vector of loss_fuction(c₁(0,xi), rand(N)) for each xi in x_vector
choice with the best convergence from the vector of loss_fuctions
ChrisRackauckas commented 2 years ago

No. Think of it as a PDE with 2 more dimensions. Make the initial condition be c₁(0,x,r1,r2) ~ r1, c₂(0,x,r1,r2) ~ r2 and solve. Then every slice NN(:,:,r1,r2) is a solution with the given initial condition. This means that NN(:,:,rand(2)...) is a quick and effective way to sample solutions with random initial conditions.

KirillZubov commented 2 years ago

why do we use random numbers how initial conditions? is it an unknown initial condition and we try to pick one from a random set which is better? sort of Monte Carlo method.

Trying to figure out this problem https://discourse.julialang.org/t/diffeqoperators-and-differentialequations-solving-a-coupled-system-of-pdes/74722

ChrisRackauckas commented 2 years ago

There is no randomness in the way I described the training.

KirillZubov commented 2 years ago

that is, we initialize random vectors r1, r2 once and use them without generating new ones at each iteration?

ChrisRackauckas commented 2 years ago

No, they are not random. They are independent variables, like t or x.

KirillZubov commented 2 years ago

ok, got it