I’ve been trying to generate a system of equations in vector format using MTK and MOL, and I keep getting errors related to broadcast-created equations missing LHS/RHS. I’ve adapted a tutorial problem to show the issue.
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, DataInterpolations
# Parameters, variables, and derivatives
n_comp = 2
@parameters p[1:n_comp]
@variables t, x, u(..)[1:n_comp]
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
params = [p => [1.5, 2.0]]
# 1D PDE and boundary conditions
eq = Dt.(u(t, x)) .~ p .* Dxx.(u(t, x))
eqs_collected = vcat(collect.(eq)...) #tried with and without vcat(collect.()...)
bcs = [u(0, x) .~ p .* cos(x),
u(t, 0) .~ sin(t),
u(t, 1) .~ exp(-t) * cos(1),
Dx.(u(t,0)) .~ 0.0]
bcs_collected = vcat(collect.(bcs)...)
# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
x ∈ Interval(0.0, 1.0)]
# PDE system
@named pdesys = PDESystem(eqs_collected, bcs_collected, domains, [t, x], [u(t, x)], params)
# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization) #error actually occurs here when running
# Solve ODE problem
using OrdinaryDiffEq
sol = solve(prob, Tsit5(), saveat=0.2)
The actual error I get is
ERROR: type Term has no field lhs.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] throw_no_field(#unused#::Val{:Term}, s::Symbol)
@ Unityper \.julia\packages\Unityper\vtKYP\src\compactify.jl:481
[3] getproperty(x::SymbolicUtils.BasicSymbolic{Any}, s::Symbol)
@ SymbolicUtils \.julia\packages\Unityper\vtKYP\src\compactify.jl:290
[4] cardinalize_eqs!(pdesys::PDESystem)
@ MethodOfLines \.julia\packages\MethodOfLines\jB2jr\src\system_parsing\pde_system_transformation.jl:39
[5] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
@ MethodOfLines \.julia\packages\MethodOfLines\jB2jr\src\MOL_discretization.jl:29
[6] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ MethodOfLines \.julia\packages\MethodOfLines\jB2jr\src\MOL_discretization.jl:132
[7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
@ MethodOfLines \.julia\packages\MethodOfLines\jB2jr\src\MOL_discretization.jl:131
[8] top-level scope
@ Untitled-1:34
I found this thread(Cannot use a set of equations in ODESystem with variables that are vectors · Issue #1196 · SciML/ModelingToolkit.jl · GitHub) , which suggested using
eqs_collected = vcat(collect.(eqs)...)
But unfortunately, that didn’t fix the issue for the equations.
The boundary conditions form a very nice vector of equations even when they contain a differential which I thought was the issue for a while, but the equations themselves stay messy and actually cause the error.
eqs_collected looks like:
2-element Vector{SymbolicUtils.BasicSymbolic{Any}}:
(broadcast(collect, broadcast(~, broadcast(Differential(t), u(t, x)), broadcast(*, p, broadcast(ComposedFunction{Any, Any}(Differential(x), Differential(x)), u(t, x))))))[1]
(broadcast(collect, broadcast(~, broadcast(Differential(t), u(t, x)), broadcast(*, p, broadcast(ComposedFunction{Any, Any}(Differential(x), Differential(x)), u(t, x))))))[2]
while bcs_collected looks like:
8-element Vector{Equation}:
(u(0, x))[1] ~ cos(x)*p[1]
(u(0, x))[2] ~ cos(x)*p[2]
(u(t, 0))[1] ~ sin(t)
(u(t, 0))[2] ~ sin(t)
(u(t, 1))[1] ~ 0.5403023058681398exp(-t)
(u(t, 1))[2] ~ 0.5403023058681398exp(-t)
Differential(x)((u(t, 0))[1]) ~ 0.0
Differential(x)((u(t, 0))[2]) ~ 0.0
Am I making a dumb mistake somewhere, or is there a better way to do this?
edit: fixed some typos