PDESystem Internal System Compatibility

[jbiffl]

The definition for a PDESystem includes a field for internal systems, (PDESystem · ModelingToolkit.jl) but it looks like discretization of the PDESystem does not look at these systems, causing errors when discretization is attempted. From what I can tell, the equations, variables, and parameters from these internal systems are not pulled into the main PDESystem in the way that they are for ODESystems on system definition from eqs, ivs, dvs, etc.
ODESystem(eqs, t, vars, pars; systems = [sys])
vs
PDESystem(eqs, bcs, domains, [ivs], [dvs], pars; systems = [sys])

Currently, the only reference to PDESystem.systems that I can find is in Give PDESystem a systems field by xtalax · Pull Request #1588 · SciML/ModelingToolkit.jl · GitHub which seems related to visualization.

[ChrisRackauckas]

@xtalax we should make sure it's supported or throw

[ts-tkd]

Any update on this? I'm trying to implement @ChrisRackauckas 's ADS Paper via MTK & the new (great work btw!) ModelingToolkitNeuralNets.NeuralNetworkBlock and I suspect I hit the same issue (code snippet below). I guess I'll have to discretize by hand for now ? Or did you find a workaround @jbiffl ?

Also: Happy to draft a PR if someone can provide me me with an overview of what needs to be changed (maybe @xtalax ? )

Example:

using OrdinaryDiffEq
using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitNeuralNets
using MethodOfLines
using DomainSets
using Lux

@parameters t, x
@variables c(..), q(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

@named nn_in = RealInput(nin = 1)
@named nn_out = RealOutput(nout = 2)

#PDE, domains etc 
eq  = [Dt(c(t, x)) ~ Dx(c(t,x)) + Dxx(c(t, x)) - Dt(q(t,x)),
        Dt(q(t,x))  ~  nn_in.u[1],
        c ~ nn_out.u[1],
        q ~ nn_out.u[2]
    ]
bcs = [c(tstart, x) ~ 0.0,
        c(t, xstart) ~ 1.0 +   Dx(c(t, xstart)),
        Dx(c(t, xend)) ~ 0.0,
        
        q(tstart,x) ~ 3.0,
        Dx(q(t, xstart)) ~ 0.0,
        Dx(q(t, xend)) ~ 0.0
        ]

domains = [t ∈ Interval(tstart, tend),
           x ∈ Interval(xstart, xend)]
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [c(t, x),q(t,x)], systems = [nn_in, nn_out])

#discretize, convert to ODE sys
dx = 0.3
order = 2
discretization = MOLFiniteDifference([x => dx], t)

ode_prob = discretize(pdesys,discretization)  ### fails

# from here on connect NN Block, optimize etc..

[ChrisRackauckas]

I'm not sure how this issue is related to doing UDEs with MethodOfLines.jl (and this reply seems to be on the wrong repo?). It's something we need to work on but I wouldn't expect MOL to handle it right now.