SciML / MethodOfLines.jl

Automatic Finite Difference PDE solving with Julia SciML
https://docs.sciml.ai/MethodOfLines/stable/
MIT License
161 stars 31 forks source link

Defining an array of variables causes error with ModelingToolkit.discretize() #43

Closed ghost closed 2 years ago

ghost commented 2 years ago

The code below generates an error message at ModelingToolkit.discretize()

begin
    # Dependencies
    using DiffEqOperators, MethodOfLines, OrdinaryDiffEq, ModelingToolkit, DomainSets, NonlinearSolve

    N = 6 # number of dependent variables

    # Variables, parameters, and derivatives
    @parameters x
    @variables u[1:N](..)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    # Domain edges
    x_min= 0.
    x_max = 1.

    # Discretization parameters
    dx = 0.1
    order = 2

    # Equations
    eqs  = Vector{ModelingToolkit.Equation}(undef, N)
    for i = 1:N
        eqs[i] = Dxx(u[i](x)) ~ u[i](x)
    end

    # Initial and boundary conditions
    bcs = Vector{ModelingToolkit.Equation}(undef, 2*N)
    for i = 1:N
        bcs[i] = Dx(u[i](x_min)) ~ 0.
    end 

    for i = 1:N
        bcs[i+N] = u[i](x_max) ~ rand()
    end

    # Space and time domains
    domains = [x ∈ Interval(x_min, x_max)]

    # PDE system
    @named pdesys = PDESystem(eqs, bcs, domains, [x], [u[i](x) for i = 1:N])

    # Method of lines discretization
    discretization = MOLFiniteDifference([x=>dx], nothing, approx_order=order)
    prob = ModelingToolkit.discretize(pdesys,discretization)

    # # Solution of the ODE system
    sol = NonlinearSolve.solve(prob,NewtonRaphson())
end

This is the full error message

MethodError: no method matching nameof(::SymbolicUtils.Term{SymbolicUtils.FnType{Tuple, Real}, Nothing})

Closest candidates are:

nameof(!Matched::SymbolicUtils.Sym) at C:\Users\User\.julia\packages\SymbolicUtils\v2ZkM\src\types.jl:144

nameof(!Matched::ModelingToolkit.AbstractSystem) at C:\Users\User\.julia\packages\ModelingToolkit\Uaky4\src\systems\abstractsystem.jl:140

nameof(!Matched::DataType) at C:\Users\User\AppData\Local\Programs\Julia-1.7.0\share\julia\base\reflection.jl:223

...

(::MethodOfLines.var"#12#27"{Nothing})(::SymbolicUtils.Term{Real, Base.ImmutableDict{DataType, Any}})@discretize_vars.jl:43
iterate@generator.jl:47[inlined]
_collect(::Vector{Any}, ::Base.Generator{Vector{Any}, MethodOfLines.var"#12#27"{Nothing}}, ::Base.EltypeUnknown, ::Base.HasShape{1})@array.jl:744
collect_similar(::Vector{Any}, ::Base.Generator{Vector{Any}, MethodOfLines.var"#12#27"{Nothing}})@array.jl:653
map(::Function, ::Vector{Any})@abstractarray.jl:2849
MethodOfLines.DiscreteSpace(::Vector{Symbolics.VarDomainPairing}, ::Vector{Any}, ::Vector{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}}, ::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})@discretize_vars.jl:41
symbolic_discretize(::ModelingToolkit.PDESystem, ::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})@MOL_discretization.jl:37
discretize(::ModelingToolkit.PDESystem, ::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})@MOL_discretization.jl:120
top-level scope@Local: 45[inlined]
xtalax commented 2 years ago

It seems that variable arrays are currently unsupported, I'll look for a workaround and try to patch in the next update.

xtalax commented 2 years ago

Found a fix, should be merged soon