Defining an array of symbolic variables causes error with NeuralPDE.discretize()

General Usage
,
1

I am trying to solve a PDE with variables [t, z1, …, zn] using NeuralPDE.jl. Here’s my code:

using Test, Flux, Optim, DiffEqFlux, Optimization
using Random, NeuralPDE, DifferentialEquations
using Statistics, Distributions, LinearAlgebra
import ModelingToolkit: Interval 
import DomainSets: UnitInterval
Random.seed!(100)

order = 3

@parameters t z[1:order]
z = collect(z)
@variables u(..)
Dt = Differential(t)
α = 1.2
β = 1.1

function dW(t, z)
    val = 0 
    for i in 1:size(z)[1]
        val += z[i]*cos((i-1/2)π*t)
    end
    val = √2*val
end
eq  = Dt(u(t,z))  ~ α*u(t,z) + β*u(t,z)*dW(t,z)
bcs = [u(0, z) ~ 1.0]

# Space and time domains
domains = Vector{Symbolics.VarDomainPairing}(undef, order + 1)
domains[1] = t ∈ Interval(0.0, 1.0)
for i in 1:order
    domains[i+1] = z[i] ∈ Interval(0.0, 1.0)
end

# number of dimensions
dim = order + 1
chain = Flux.Chain(Dense(dim,16,Flux.σ),Dense(16,16,Flux.σ),Dense(16,1))
# Initial parameters of Neural network
initθ = Float64.(DiffEqFlux.initial_params(chain))

# Discretization
dx = 0.05
discretization = PhysicsInformedNN(chain,GridTraining(dx),init_params =initθ)

@named pde_system = PDESystem(eq,bcs,domains,[t,z],[u(t,z)])
prob = discretize(pde_system,discretization)

When I run this, I get the following error:

ERROR: MethodError: no method matching nameof(::Vector{Num})
Closest candidates are:
  nameof(::Sym) at ~/.julia/packages/SymbolicUtils/vnuIf/src/types.jl:144
  nameof(::ModelingToolkit.AbstractSystem) at ~/.julia/packages/ModelingToolkit/iHLWM/src/systems/abstractsystem.jl:139
  nameof(::DataType) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/reflection.jl:223
  ...
Stacktrace:
  [1] _getname(x::Vector{Num}, #unused#::Dict{Any, Any})
    @ Symbolics ~/.julia/packages/Symbolics/sDAUx/src/variable.jl:351
  [2] getname(x::Vector{Num}, val::Dict{Any, Any}) (repeats 2 times)
    @ Symbolics ~/.julia/packages/Symbolics/sDAUx/src/variable.jl:364
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:670 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:643 [inlined]
  [5] getindex
    @ ./broadcast.jl:597 [inlined]
  [6] copyto_nonleaf!(dest::Vector{Symbol}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Symbolics.getname), Tuple{Base.Broadcast.Extruded{Vector{Any}, Tuple{Bool}, Tuple{Int64}}}}, iter::Base.OneTo{Int64}, state::Int64, count::Int64)
    @ Base.Broadcast ./broadcast.jl:1055
  [7] copy
    @ ./broadcast.jl:907 [inlined]
  [8] materialize
    @ ./broadcast.jl:860 [inlined]
  [9] get_vars(indvars_::Vector{Any}, depvars_::Vector{Num})
    @ NeuralPDE ~/Documents/Julia-fork/NeuralPDE.jl/src/pinns_pde_solve.jl:750
 [10] discretize_inner_functions(pde_system::PDESystem, discretization::PhysicsInformedNN{GridTraining{Float64}, Vector{Float64}, NeuralPDE.Phi{Optimisers.Restructure{Chain{Tuple{Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}}}, NamedTuple{(:layers,), Tuple{Tuple{NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}}}}}}, typeof(NeuralPDE.numeric_derivative), Bool, Nothing, NonAdaptiveLoss{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
    @ NeuralPDE ~/Documents/Julia-fork/NeuralPDE.jl/src/pinns_pde_solve.jl:1257
 [11] discretize(pde_system::PDESystem, discretization::PhysicsInformedNN{GridTraining{Float64}, Vector{Float64}, NeuralPDE.Phi{Optimisers.Restructure{Chain{Tuple{Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}}}, NamedTuple{(:layers,), Tuple{Tuple{NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}}}}}}, typeof(NeuralPDE.numeric_derivative), Bool, Nothing, NonAdaptiveLoss{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
    @ NeuralPDE ~/Documents/Julia-fork/NeuralPDE.jl/src/pinns_pde_solve.jl:1512
 [12] top-level scope
    @ ~/Documents/Julia-fork/NeuralPDE.jl/test/NNSDE_wp_tests.jl:56

The stack trace tells me that the error is coming from Symbolics.jl, but I’m not sure if the issue is with how I am initializing my variables or if something else is going on here.

2

z is a vector of Nums (vector of symbolic values), not a Num itself. But the constructor that the example shows is always just an array or list of symbolic values. So you want to do

@named pde_system = PDESystem(eq,bcs,domains,[t,z...],[u(t,z...)])

which is the same as writing out

@named pde_system = PDESystem(eq,bcs,domains,[t,z1,z2,...,zn],[u(t,z1,z2,...,zn)])
3

I suspect you need to splat z everywhere you call u:

using Flux, DiffEqFlux
using Random, NeuralPDE, DifferentialEquations
import ModelingToolkit: Interval
import DomainSets: UnitInterval
Random.seed!(100)

order = 3
@parameters z[1:order]
# z = collect(z)
@variables t u(..)
Dt = Differential(t)
α = 1.2
β = 1.1

function dW(t, z)
    val = 0
    for (i, zi) in enumerate(z)
        val += zi*cos((i-1/2)π*t)
    end
    val = √2*val
end

eq  = Dt(u(t,z...))  ~ α*u(t,z...) + β*u(t,z...)*dW(t,z)

bcs = [u(0, z...) ~ 1.0]

# Space and time domains
domains = Vector{Symbolics.VarDomainPairing}(undef, order + 1)
domains[1] = t ∈ Interval(0.0, 1.0)
for i in 1:order
    domains[i+1] = z[i] ∈ Interval(0.0, 1.0)
end

# number of dimensions
dim = order + 1
chain = Flux.Chain(Dense(dim,16,Flux.σ),Dense(16,16,Flux.σ),Dense(16,1))
# Initial parameters of Neural network
initθ = Float64.(DiffEqFlux.initial_params(chain))

# Discretization
dx = 0.05
discretization = PhysicsInformedNN(chain,GridTraining(dx),init_params =initθ)

@named pde_system = PDESystem(eq,bcs,domains,[t,z...],[u(t,z...)])
prob = discretize(pde_system,discretization)
But that still gives an error: No method matching nameof(::Term)

julia> prob = discretize(pde_system,discretization)
ERROR: MethodError: no method matching nameof(::Term{Real, Base.ImmutableDict{DataType, Any}})
Closest candidates are:
nameof(::Sym) at ~/.julia/packages/SymbolicUtils/vnuIf/src/types.jl:144
nameof(::ModelingToolkit.AbstractSystem) at ~/.julia/packages/ModelingToolkit/tMgaW/src/systems/abstractsystem.jl:139
nameof(::DataType) at /usr/local/stow/julia-1.7.3/share/julia/base/reflection.jl:223
…
Stacktrace:
[1] (::NeuralPDE.var"#40#41")(argument::Term{Real, Base.ImmutableDict{DataType, Any}})
@ NeuralPDE ./none:0
[2] iterate
@ ./generator.jl:47 [inlined]
[3] collect_to!(dest::Vector{Symbol}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, NeuralPDE.var"#40#41"}, offs::Int64, st::Int64)
@ Base ./array.jl:782
[4] collect_to_with_first!
@ ./array.jl:760 [inlined]
[5] collect(itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, NeuralPDE.var"#40#41"})
@ Base ./array.jl:734
[6] get_vars(indvars_::Vector{Num}, depvars_::Vector{Num})
@ NeuralPDE ~/.julia/packages/NeuralPDE/iNhvg/src/symbolic_utilities.jl:353
[7] symbolic_discretize(pde_system::PDESystem, discretization::PhysicsInformedNN{GridTraining{Float64}, Vector{Float64}, NeuralPDE.Phi{Optimisers.Restructure{Chain{Tuple{Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}}}, NamedTuple{(:layers,), Tuple{Tuple{NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}}}}}, Nothing}, typeof(NeuralPDE.numeric_derivative), Bool, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
@ NeuralPDE ~/.julia/packages/NeuralPDE/iNhvg/src/discretize.jl:420
[8] discretize(pde_system::PDESystem, discretization::PhysicsInformedNN{GridTraining{Float64}, Vector{Float64}, NeuralPDE.Phi{Optimisers.Restructure{Chain{Tuple{Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}}}, NamedTuple{(:layers,), Tuple{Tuple{NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}}}}}, Nothing}, typeof(NeuralPDE.numeric_derivative), Bool, Nothing, Nothing, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
@ NeuralPDE ~/.julia/packages/NeuralPDE/iNhvg/src/discretize.jl:669
[9] top-level scope
@ REPL[29]:1
[10] top-level scope
@ ~/.julia/packages/CUDA/tTK8Y/src/initialization.jl:52

4

@shashi let’s chat about this one when I get back.

5

Is there anything I can do to help with this?

6

it requires splatting z at the moment… We should eventually support it, but still only when z is not collected… Do open an issue about this so that it does not get lost.!

7

@hpieper14 open an issue on NeuralPDE.jl with this example. I think you’d just have to go step by step through the parser and make it expand z internally. Though the nameof issue might be solved by @shashi’s recent change?

8

No I think the recent changes are unrelated. It really needs to be handled at the modeling DSL level, not Symbolics level.