Quadrature strategy does not converge to solution

[marcofrancis]

Hi, I want to solve a fairly simple PDE:

# Parameters
τ_min= 0.
τ_max = 1.0
τ_span = (τ_min,τ_max)

ω_min = -2.
ω_max = 2.

#Grid if needed
dω = 0.05
dt_NN = 0.01

μY = 0.03
σ = 0.03
λ = 0.1
ωb = 0

γ = 1.5
k = 0.1

#Initial Condition
h(ω) = exp.(-(γ-1)*5*ω.^2)


@parameters τ ω
@variables f(..)

Dω = Differential(ω)
Dωω = Differential(ω)^2
Dτ = Differential(τ)

μω = λ*(ω-ωb)
σω = σ

# PDE
eq  = 1/100*Dτ(f(τ,ω)) ~ 0.5*σω^2*Dωω(f(τ,ω)) - μω*Dω(f(τ,ω)) - k*f(τ,ω) + h(ω) 

# Initial and boundary conditions
bcs = [f(τ_min,ω) ~ h(ω),
       Dω(f(τ,ω_max)) ~ 0,
       Dω(f(τ,ω_min)) ~ 0]

# Space and time domains
domains = [τ ∈ Interval(τ_min,τ_max),
           ω ∈ Interval(ω_min,ω_max)]

τs,ωs = [infimum(d.domain):dω:supremum(d.domain) for d in domains]


@named pde_system = PDESystem(eq,bcs,domains,[τ,ω],[f(τ, ω)])

# NN parameters
dim =2
hls = dim+50

If I use a grid strategy everything works as expected and I get quite close to the actual solution. If I use a quadrature strategy:

    
    chain = Lux.Chain(Dense(dim,hls,Lux.σ),
                      Dense(hls,hls,Lux.σ),
                      Dense(hls,1))

    strategy = QuadratureTraining(; quadrature_alg = CubatureJLp(), 
    reltol = 1e-5, abstol = 1e-8,
    maxiters = 1_000)

    discretization = PhysicsInformedNN(chain,
                                   strategy)

    prob = discretize(pde_system,discretization)

    callback = function (p,l)
        println("Current loss is: $l")
        return false
    end

 
    res = Optimization.solve(prob,Adam(0.01);callback = callback,maxiters=5000)
    prob = remake(prob,u0=res.u)

The solution (NN line) is completely off: Initial Condition: Value at terminal time:

[KirillZubov]

QuadratureStrategy is still unstable. Try NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh()) I'm not sure that it is can be more stable it needs to dig into the problem. GridTraining is most stable at this moment in time.

[KirillZubov]

I ran the script with Grid(predict) and quadrature strategy(predict2) it is quite the same

strategy = NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh(),
                                        reltol = 1e-3, abstol = 1e-3,
                                        maxiters = 50)

[marcofrancis]

I ran the script with Grid(predict) and quadrature strategy(predict2) it is quite the same

strategy = NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh(),
                                        reltol = 1e-3, abstol = 1e-3,
                                        maxiters = 50)

Could you share the script you used? So, that I can understand what's going wrong with mine, thanks!

[KirillZubov]

the same that you, only I changed the strategy to

strategy = NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh(),
                                        reltol = 1e-3, abstol = 1e-3,
                                        maxiters = 50)

[KirillZubov]

using Flux, NeuralPDE, Test
using Optimization, OptimizationOptimJL, OptimizationOptimisers
using Integrals, IntegralsCubature
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum
using DomainSets
import Lux

τ_min = 0.0
τ_max = 1.0
τ_span = (τ_min, τ_max)

ω_min = -2.0
ω_max = 2.0

#Grid if needed
dω = 0.05
dt_NN = 0.01

μY = 0.03
σ_ = 0.03
λ = 0.1
ωb = 0

γ = 1.5
k = 0.1

#Initial Condition
h(ω) = exp.(-(γ - 1) * 5 * ω .^ 2)

@parameters τ ω
@variables f(..)

Dω = Differential(ω)
Dωω = Differential(ω)^2
Dτ = Differential(τ)

μω = λ * (ω - ωb)
σω = σ_

# PDE
eq = 1 / 100 * Dτ(f(τ, ω)) ~ 0.5 * σω^2 * Dωω(f(τ, ω)) - μω * Dω(f(τ, ω)) - k * f(τ, ω) +
                             h(ω)

# Initial and boundary conditions
bcs = [f(τ_min, ω) ~ h(ω),
    Dω(f(τ, ω_max)) ~ 0,
    Dω(f(τ, ω_min)) ~ 0]

# Space and time domains
domains = [τ ∈ Interval(τ_min, τ_max),
    ω ∈ Interval(ω_min, ω_max)]

τs, ωs = [infimum(d.domain):dω:supremum(d.domain) for d in domains]

@named pde_system = PDESystem(eq, bcs, domains, [τ, ω], [f(τ, ω)])

# NN parameters
dim = 2
hls = dim + 50

chain = Lux.Chain(Lux.Dense(dim, hls, Lux.σ),
                  Lux.Dense(hls, hls, Lux.σ),
                  Lux.Dense(hls, 1))

chain = Lux.Chain(Lux.Dense(2, 12, Lux.tanh), Lux.Dense(12, 12, Lux.tanh), Lux.Dense(12, 1))

# strategy = NeuralPDE.QuadratureTraining(; quadrature_alg = CubatureJLp(),
#                               reltol = 1e-5, abstol = 1e-8,
#                               maxiters = 1000)

strategy = NeuralPDE.QuadratureTraining(quadrature_alg = CubatureJLh(),
                                        reltol = 1e-3, abstol = 1e-3,
                                        maxiters = 50)
# strategy = NeuralPDE.GridTraining(0.1)
discretization = NeuralPDE.PhysicsInformedNN(chain, strategy)

prob = NeuralPDE.discretize(pde_system, discretization)

callback = function (p, l)
    println("Current loss is: $l")
    return false
end

res = Optimization.solve(prob, OptimizationOptimisers.Adam(0.01); callback = callback,
                         maxiters = 5000)

grid_strategy = NeuralPDE.GridTraining(0.01)
discretization = NeuralPDE.PhysicsInformedNN(chain, grid_strategy)

prob = NeuralPDE.discretize(pde_system, discretization)

callback = function (p, l)
    println("Current loss is: $l")
    return false
end

res2 = Optimization.solve(prob, OptimizationOptimisers.Adam(0.01); callback = callback,
                          maxiters = 5000)

phi = discretization.phi
ωs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][2]
u_real = [h() for t in ωs]
u_predict = [first(phi([0, t], res.minimizer)) for t in ωs]
u_predict2 = [first(phi([0, t], res2.minimizer)) for t in ωs]

using Plots
t_plot = collect(ωs)
plot(t_plot, u_real)
plot!(t_plot, u_predict)
plot!(t_plot, u_predict2)


xs, ys = [infimum(d.domain):dx:supremum(d.domain)
          for (dx, d) in zip([0.01, 0.04], domains)]
u_predict = reshape([first(Array(phi([x, y], res.minimizer))) for x in xs for y in ys],
                    (length(xs), length(ys)))
u_predict2 = reshape([first(Array(phi([x, y], res2.minimizer))) for x in xs for y in ys],
                     (length(xs), length(ys)))

p1 = plot(xs, ys, u_predict, linetype = :contourf, title = "predict");
p2 = plot(xs, ys, u_predict2, linetype = :contourf, title = "predict2");

plot(p1, p2)

[marcofrancis]

Mmm, somehow if I run the code you pasted (only the second chain is trained, right?) I still get this for the initial data: Here's the list of packages in my environment: [052768ef] CUDA v3.12.0 [de52edbc] Integrals v3.1.2 [c31f79ba] IntegralsCubature v0.2.0 [033835bb] JLD2 v0.4.22 ⌃ [b2108857] Lux v0.4.18 ⌃ [961ee093] ModelingToolkit v8.19.0 [315f7962] NeuralPDE v5.2.0 [7f7a1694] Optimization v3.8.2 [36348300] OptimizationOptimJL v0.1.2 [42dfb2eb] OptimizationOptimisers v0.1.0 ⌃ [91a5bcdd] Plots v1.31.7 [8a4e6c94] QuasiMonteCarlo v0.2.9 [9a3f8284] Random

[KirillZubov]

The first is just too big, PINNs in general are still quite sensitive and unstable. with too big FFN is usually underfitting.