Hello, I am new to Julia, so apologies if I am missing something basic, but I wasn’t getting any replies in New To Julia, so I’m not sure that was the right place for this post. I am trying to solve a system of PDEs using MethodOfLines and ModelingToolkit, but I can’t seem to get past the discretization. The actual module I’m writing is quite long and makes use of a commercial library for my equation of state, so I created a MWE that I’ll be referencing for the rest of this post. Here is the MWE:
module MWE
using ModelingToolkit, DifferentialEquations, OrdinaryDiffEq, DomainSets, MethodOfLines
function h(f,g)
f*g + g
end
@register_symbolic h(f,g)
@parameters t z a b
@variables f(..) g(..)
Dt = Differential(t)
Dz = Differential(z)
eq = [
Dt(h(f(t,z),g(t,z))) ~ a*Dz(f(t,z)) + b*Dz(g(t,z)),
Dt(f(t,z)) ~ Dz(g(t,z))
]
fₒ(t) = 13 + 15*t
gₒ(t) = 10 + 30*t
bc = [
f(t,0) ~ fₒ(t),
g(t,0) ~ gₒ(t),
f(0,z) ~ 10,
g(0,z) ~ 20
]
ps = [
a => 5,
b => 10
]
domains = [
t ∈ Interval(0, 10),
z ∈ Interval(0, 30)
]
@named sys = PDESystem(eq,bc,domains,[t,z],[f(t,z),g(t,z)],ps)
dz = 0.1
discretization = MOLFiniteDifference([z=>dz],t)
prob = discretize(sys,discretization)
sol = solve(prob,Tsit5(),saveat=0.1)
end
While running, I get these warnings:
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Int64},).
└ @ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Int64},).
└ @ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Int64},).
└ @ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Int64},).
└ @ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\validation.jl:154
┌ Warning: Incompatible term found. Adding auxiliary equation var"⟦Main.MWE.h(f(t, z), g(t, z))⟧"(t, z) ~ Main.MWE.h(f(t, z), g(t, z)) to the system.
└ @ MethodOfLines C:\Users\erich\.julia\packages\MethodOfLines\byu6o\src\system_parsing\pde_system_transformation.jl:227
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Int64},).
└ @ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Int64},).
└ @ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\validation.jl:154
The full error code after these warnings is extremely long as it seems to display all of the equations for each discretized section, but after those, there is this error:
ERROR: LoadError: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[(var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[1], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[2], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[3], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[4], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[5], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[6], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[7], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[8], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[9], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[10],
... ##THE PATTERN CONTINUES THROUGH, I CUT OUT THE MIDDLE##
...(var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[291], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[292], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[293], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[294], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[295], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[296], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[297], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[298], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[299], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[300], (var"var\"⟦Main.MWE.h(f(t, z), g(t, z))⟧\""(t))[301]] are missing from the variable map.
Stacktrace:
[1] throw_missingvars(vars::Vector{SymbolicUtils.BasicSymbolic{Real}})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\variables.jl:122
[2] _varmap_to_vars(varmap::Dict{Any, Any}, varlist::Vector{SymbolicUtils.BasicSymbolic{Real}}; defaults::Dict{Any, Any}, check::Bool, toterm::typeof(ModelingToolkit.default_toterm))
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\variables.jl:116
[3] varmap_to_vars(varmap::Vector{Pair}, varlist::Vector{SymbolicUtils.BasicSymbolic{Real}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\variables.jl:85
[4] get_u0_p(sys::ModelingToolkit.ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; use_union::Bool, tofloat::Bool, symbolic_u0::Bool)
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:786
[5] get_u0_p
@ C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:768 [inlined]
[6] process_DEProblem(constructor::Type, sys::ModelingToolkit.ODESystem, u0map::Vector{…}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), kwargs::@Kwargs{…})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:811
[7] (SciMLBase.ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ModelingToolkit.ODESystem, u0map::Vector{Pair}, tspan::Tuple{Int64, Int64}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:936
[8] ODEProblem
@ C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:929 [inlined]
[9] (SciMLBase.ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ModelingToolkit.ODESystem, u0map::Vector{Pair}, tspan::Tuple{Int64, Int64})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:929
[10] (SciMLBase.ODEProblem{true})(::ModelingToolkit.ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:916
[11] (SciMLBase.ODEProblem{true})(::ModelingToolkit.ODESystem, ::Vector{Pair}, ::Vararg{Any})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:915
[12] SciMLBase.ODEProblem(::ModelingToolkit.ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:912
[13] SciMLBase.ODEProblem(::ModelingToolkit.ODESystem, ::Vector{Pair}, ::Vararg{Any})
@ ModelingToolkit C:\Users\erich\.julia\packages\ModelingToolkit\Gpzyo\src\systems\diffeqs\abstractodesystem.jl:911
[14] discretize(pdesys::ModelingToolkit.PDESystem, discretization::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::@Kwargs{})
@ PDEBase C:\Users\erich\.julia\packages\PDEBase\Aqj4G\src\discretization_state.jl:74
[15] discretize(pdesys::ModelingToolkit.PDESystem, discretization::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
@ PDEBase C:\Users\erich\.julia\packages\PDEBase\Aqj4G\src\discretization_state.jl:55
[16] top-level scope
@ C:\Users\erich\Documents\GitHub\AmmoniaSynthesis\AmmoniaPlant\src\MWE.jl:46
[17] include(fname::String)
@ Base.MainInclude .\client.jl:489
[18] top-level scope
@ REPL[16]:1
in expression starting at C:\Users\erich\Documents\GitHub\AmmoniaSynthesis\AmmoniaPlant\src\MWE.jl:1
Some type information was truncated. Use `show(err)` to see complete types.
It seems the issue is related to having a derivative of a user-defined function. In my original code, I currently do that a few times in the equations, and the missing vars issue seems to happen in each instance, so I believe that is the origin of the problem, but I’m not sure how to fix it. The user-defined functions in my original code are quite complicated, so I’d rather like to avoid deriving those functions by hand and ending up with messy and unreadable code, not to mention some of them include the commercial library I mentioned, which does take multiple of my variables as input. Does anyone have any suggestions on how to fix this issue?
Edit: Forgot to add my package list:
[6e4b80f9] BenchmarkTools v1.5.0
[479239e8] Catalyst v13.5.1
[a93c6f00] DataFrames v1.6.1
[82cc6244] DataInterpolations v4.7.0
[0c46a032] DifferentialEquations v7.13.0
⌅ [5b8099bc] DomainSets v0.6.7
[06fc5a27] DynamicQuantities v0.13.0
[94925ecb] MethodOfLines v0.10.7
⌅ [961ee093] ModelingToolkit v8.75.0
[1dea7af3] OrdinaryDiffEq v6.74.0
[438e738f] PyCall v1.96.4
[56ddb016] Logging
Of note, it seems some errors related to the variable map have been fixed in newer versions of ModelingToolkit, but I can’t update past v8.75 due to MethodOfLines. That being said, I’m not sure whether or not that would fix my specific issue. If anyone has any suggestions on how I could test that, please let me know, but I can’t think of a way to recreate the scenario without MethodOfLines.
Also, I have tried solving this with NeuralPDE, but I had the same error and same ModelingToolkit update restriction. I’m not aware of any other packages that are designed to work well with ModelingToolkit, but I would try this with a different package if possible.