Problem with 2D Poisson-Eqn. PDE

[RainerHeintzmann]

Trying to reproduce the Poisson equation example from the JuliaCon talk, I hit a dead end, see also here.
I locally added compat entries for DiffEqOperators.jl so that it accepts the current version of MTK, but for 2-dimensional PDEs always the same problem seems to be appearing:

julia> depvar_ops = map(x->operation(x.val),pde_system.depvars)
ERROR: type CallWithMetadata has no field val
Stacktrace:
 [1] getproperty(x::Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}, f::Symbol)
   @ Base .\Base.jl:33
 [2] (::var"#7#8")(x::Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}})
   @ Main .\REPL[19]:1
 [3] iterate
   @ .\generator.jl:47 [inlined]
 [4] _collect
   @ .\array.jl:691 [inlined]
 [5] collect_similar(cont::Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}, itr::Base.Generator{Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}, var"#7#8"})
   @ Base .\array.jl:606
 [6] map(f::Function, A::Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}})
   @ Base .\abstractarray.jl:2294
 [7] top-level scope
   @ REPL[19]:1

Here is the Poinsson example that causes this problem:

using ModelingToolkit, DiffEqOperators, DifferentialEquations, DomainSets
import ModelingToolkit: Interval, infimum, supremum

@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = (Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y))

bcs = [u(0,y) ~ 0.f0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
       u(x,0) ~ 0.f0, u(x,1) ~ -sin(pi*x)*sin(pi*1) ]

domains = [x ∈ Interval(0.0,1.0),
           y ∈ Interval(0.0,1.0)]  # \in

pde_system = PDESystem(eq, bcs, domains, [x,y], [u], name=:Poisson_System)

dx = 0.05; dy = 0.05;
discretization = MOLFiniteDifference([x=>dx, y=>dy], nothing, centered_order=2)

using ModelingToolkit: operation
depvar_ops = map(x->operation(x.val),pde_system.depvars)   # copied from DiffEqOperators.jl, which shows the problem.

prob = discretize(pde_system, discretization)  # Here the above problem appears
sol = solve(prob)

For a 1D problem, like the diffusion example, the code runs fine.
Any ideas how to fix this for this 2D Poisson equation example?

[ChrisRackauckas]

What version are you on? We're in version chaos right now 😅

[RainerHeintzmann]

I tried the latest released versions that are supposed to work together, MTK: 5.23 and DiffEqOperators.jl: 4.30.0, but then I changed to the latest MTK 6.3.1 by adding the entry to the compat entry Project.tolm of DiffEqOperators.jl.
Which version(s) did you use for your tutorial demo?

[RainerHeintzmann]

Since the approach with DiffEqOperators.jl was unsuccessful, I tried the Physics-Informed NN approach. That one could actually do the discretization but failed then during the solution step. The extra code:

using NeuralPDE, Flux, DiffEqFlux, GalacticOptim
dim = 2
chain = FastChain(FastDense(dim,16,Flux.σ),
                  FastDense(16,1^6,Flux.σ), FastDense(16,1))

# discretization dx = 0.05
discretization = PhysicsInformedNN(chain, GridTraining(dx))
prob = discretize(pde_system, discretization)
# Optimizer 
opt = Optim.BFGS()
sol = GalacticOptim.solve(prob, opt, maxiters=10)

and the corresponding error message:

julia> sol = GalacticOptim.solve(prob, opt, maxiters=10)
ERROR: DimensionMismatch("A has dimensions (1,16) but B has dimensions (1,361)")
Stacktrace:
  [1] gemm_wrapper!(C::Matrix{Float32}, tA::Char, tB::Char, A::SubArray{Float32, 2, Vector{Float32}, Tuple{Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}, true}, B::Matrix{Float32}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\matmul.jl:643
  [2] mul!
    @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\matmul.jl:169 [inlined]
  [3] mul!
    @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\matmul.jl:275 [inlined]
  [4] *(A::SubArray{Float32, 2, Vector{Float32}, Tuple{Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}}}, true}, B::Matrix{Float32})
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\matmul.jl:160
  [5] adjoint
    @ ~\.julia\packages\DiffEqFlux\N7blG\src\fast_layers.jl:63 [inlined]
  [6] adjoint(__context__::Zygote.Context, f::FastDense{typeof(identity), DiffEqFlux.var"#initial_params#90"{Vector{Float32}}}, x::Matrix{Float32}, p::Vector{Float32})
    @ DiffEqFlux .\none:0

Here is a list of the installed packages:

(Modeling) pkg> st
      Status `C:\Users\pi96doc\Documents\Programming\Julia\Development\Modeling\Project.toml`
  [aae7a2af] DiffEqFlux v1.43.0
  [9fdde737] DiffEqOperators v4.30.0
  [0c46a032] DifferentialEquations v6.18.0
  [5b8099bc] DomainSets v0.5.2
  [587475ba] Flux v0.12.6
  [a75be94c] GalacticOptim v2.0.3
  [961ee093] ModelingToolkit v5.26.0
  [315f7962] NeuralPDE v3.15.0

[lamorton]

We renamed PDESystem.depvars to PDESystem.dvs somewhere between versions 5.23 and 6.3.

depvar_ops = map(x->operation(x.val),pde_system.depvars)

Try changing to this & see if it works:

depvar_ops = map(x->operation(x.val),pde_system.dvs

[RainerHeintzmann]

Thanks for the hint! This did not help, same problem. Both, pde_system.depvars and pde_system.dvs exist, yet the maybe there is a problem with the type of this dependency? This is what julia tells me about it:

julia> pde_system.depvars
1-element Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}:
 u⋆

In comparison, the Diffusion example from the landing page, which does seem to work, seems to have a somewhat more defined type:

julia> pdesys.depvars
1-element Vector{Num}:
 u(t, x)

Maybe the u* is causing the issue?

[ChrisRackauckas]

chain = FastChain(FastDense(dim,16,Flux.σ),
                  FastDense(16,1^6,Flux.σ), FastDense(16,1))

that has a typo.

chain = FastChain(FastDense(dim,16,Flux.σ),
                  FastDense(16,16,Flux.σ), FastDense(16,1))

[RainerHeintzmann]

Thanks a lot! You are absolutely right about the typo. Now, the code seems to run smoothly for the NN case.
It took me a while to find out, how to interpret the result, but finally found it in the documentation :-)
Somehow the PINN solution seems to be in a different basis, and the function disretization.phi serves to translate the NN result back to the cartesian coordinate system? Is there some documentation on what phi really does?
This video, does not seem to mention this.
In the talk on ModelingToolkit.jl Chris mentions another JuliaCon2021 video on PINN, but I was unable to find this.
P.S.: I would love to use the package on solving the inhomogeneous Helmholtz equation, which is why I am testing it.
Does it also support Boundary geometries like BEAST.jl?

[ChrisRackauckas]

I don't talk about the PDE solution interface in the video because we don't have one right now. At least, it's not unified in a way that I like, so for now there's some documentation on how to interpret the PDE solution from the PINNs and the FDM etc., but I hope we can change that with some abstract PDESolution that then has common commands across the discretizers.

[ChrisRackauckas]

This got fixed.