Hi all,
I’m not sure if bug or not, but setting the fully_determined keyword argument to true on mtkcompile results in a compilation error, whereas setting it to false results in what seems to be a valid system, with an equal number of equations and unknowns, which can then be solved.
Here is code to reproduce.
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using DataInterpolationsND
itp_dims = (LinearInterpolationDimension([0.1, 10.0]), LinearInterpolationDimension([300.0, 310.0]))
pars = @parameters begin
V_internal = 10.0, []
(liq_PT_2_U_itp::NDInterpolation)(..) = NDInterpolation([112.55 154.35; 111.74 153.25], itp_dims), []
(liq_PT_2_D_itp::NDInterpolation)(..) = NDInterpolation([996.56 993.38; 1001.0 997.7], itp_dims), []
(gas_PT_2_D_itp::NDInterpolation)(..) = NDInterpolation([939.84 971.00; 940.13 971.30], itp_dims), []
(gas_PT_2_U_itp::NDInterpolation)(..) = NDInterpolation([0.16039 0.15522; 1.5971 1.5458], itp_dims), []
end
p_init = 0.3
T_liq_init = 300
T_gas_init = 300
m_liq_init = 1000
d_liq_init = liq_PT_2_D_itp(p_init, T_liq_init)
V_liq_init = m_liq_init / d_liq_init
V_ull_init = V_internal - V_liq_init
U_liq_init = m_liq_init * liq_PT_2_U_itp(p_init, T_liq_init)
d_gas_init = gas_PT_2_D_itp(p_init, T_gas_init)
m_gas_init = d_gas_init * V_ull_init
U_gas_init = m_gas_init * gas_PT_2_U_itp(p_init, T_gas_init)
vars = @variables begin
V_ull(t), [description = "ullage volume"]
V_liq(t), [description = "liquid volume"]
p(t) = p_init, [description = "pressure"]
m_gas(t) = m_gas_init, [description = "fluid mass in volume"]
U_gas(t) = U_gas_init, [description = "extensive internal energy of gas"]
T_gas(t), [description = "temperature of gas"]
d_gas(t), [description = "density of gas"]
m_liq(t) = m_liq_init, [description = "fluid mass in volume"]
U_liq(t) = U_liq_init, [description = "extensive internal energy of liquid"]
T_liq(t), [description = "temperature of liquid"]
d_liq(t), [description = "density of liquid"]
V_liq_dot(t), [description = "liquid volume rate of change"]
end
eqs = [
V_internal ~ V_ull + V_liq,
d_gas ~ m_gas / V_ull,
d_gas ~ gas_PT_2_D_itp(p, T_gas),
U_gas ~ m_gas * gas_PT_2_U_itp(p, T_gas),
D(m_gas) ~ 0.0,
D(U_gas) ~ 1.0e3 * p * V_liq_dot,
d_liq ~ liq_PT_2_D_itp(p, T_liq),
U_liq ~ m_liq * liq_PT_2_U_itp(p, T_liq),
D(m_liq) ~ 10.0,
D(U_liq) ~ 10.0 * 150.0 - 1.0e3 * p * V_liq_dot,
# chain rule D(m_liq/d_liq) to get volume rate of change
V_liq_dot ~ (D(m_liq) * d_liq - m_liq * D(d_liq)) / d_liq^2,
V_liq ~ m_liq / d_liq,
]
@named _sys = System(eqs, t, vars, pars)
# this works, resulting in system with 10 equations and 10 unknowns.
sys = mtkcompile(_sys; fully_determined=false)
# this errors out, resulting in taking derivative wrt T_gas and p_gas instead of time:
# Encountered operator `Differential(p(t), 1)` which has different independent variable than the one used in the system `t`.
# sys = mtkcompile(_sys; fully_determined=true)
guesses = [T_gas => 300, T_liq => 300, d_liq => 1000, V_ull => 9.0]
prob = ODEProblem(sys, [], (0, 10); guesses=guesses, use_scc=false)
# problem solves
sol = solve(prob; abstol=1e-6)
sol.retcode
the stack trace of error
ERROR: Encountered operator `Differential(p(t), 1)` which has different independent variable than the one used in the system `t`.
Context:
-U_gasˍt(t) + (T_gasˍt(t)*Differential(T_gas(t), 1)(gas_PT_2_U_itp(p(t), T_gas(t))) + pˍt(t)*Differential(p(t), 1)(gas_PT_2_U_itp(p(t), T_gas(t))))*m_gas(t)
Stacktrace:
[1] validate_operator(op::Differential, args::SymbolicUtils.SmallVec{…}, iv::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}; context::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…})
@ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:722
[2] validate_operator
@ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:721 [inlined]
[3] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, expr::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, iv::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, ::Type{…}; depth::Int64)
@ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:833
[4] collect_vars!
@ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:798 [inlined]
[5] #collect_vars!#37
@ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:892 [inlined]
[6] collect_vars!
@ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:887 [inlined]
[7] collect_scoped_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, sys::System, iv::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, ::Type{…}; depth::Int64)
@ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:688
[8] discover_globalscoped(sys::System)
@ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\abstractsystem.jl:643
[9] complete(sys::System; split::Bool, flatten::Bool, add_initial_parameters::Bool)
@ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\abstractsystem.jl:680
[10] complete
@ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\abstractsystem.jl:677 [inlined]
[11] mtkcompile(sys::System; additional_passes::Tuple{}, inputs::Vector{…}, outputs::Vector{…}, disturbance_inputs::Vector{…}, split::Bool, kwargs::@Kwargs{…})
@ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\systems.jl:104
[12] top-level scope
@ c:\Users\mark.garnett\EngineModeling\ignore\tank_debug_standalone.jl:65
installed packages
[a93c6f00] DataFrames v1.8.2
[4f1ef021] DataInterpolationsND v0.1.2
[8bb1440f] DelimitedFiles v1.9.1
⌅ [0c46a032] DifferentialEquations v7.17.0
⌅ [5b8099bc] DomainSets v0.7.18
⌃ [e9467ef8] GLMakie v0.13.10
⌃ [c27321d9] Glob v1.4.0
[a98d9a8b] Interpolations v0.16.2
⌃ [7ed4a6bd] LinearSolve v3.76.0
⌃ [961ee093] ModelingToolkit v11.26.0
[16a59e39] ModelingToolkitStandardLibrary v2.28.0
[15e1cf62] NPZ v0.4.3
⌃ [8913a72c] NonlinearSolve v4.19.0
⌅ [1dea7af3] OrdinaryDiffEq v6.111.0
⌅ [127b3ac7] OrdinaryDiffEqNonlinearSolve v1.28.0
[91a5bcdd] Plots v1.41.6
[efcf1570] Setfield v1.1.2
[1986cc42] Unitful v1.28.0
[de0858da] Printf v1.11.0