Is it possible to pass arbitrary types into MTK?

Specific Domains Modelling & Simulations
1

Based on the following docs (https://docs.sciml.ai/ModelingToolkit/stable/comparison/), which says ‘Modelica is a declarative programming language. ModelingToolkit.jl is a declarative symbolic modeling language’, I would imagine no, but then I imagine there might be a smart workaround. As an example, I am trying to do fluid modeling while obtaining fluid properties within MTK. To do so, I need to pass a parameter that is not a number (either a string or custom object). Is it possible to do so?

Please find my minimum example, where it is possible to define the fluid model outside of the mtk object, but it is not possible to pass it as a parameter. Is there a way I can make it possible to pass either a string or an object into the mtkmodel?

1. This works using CoolProp

## Coolprop working

# Helper functions
function calc_ρ(T, p)
    return PropsSI("D", "T", T, "P", p, fluid_name)
end

@register_symbolic calc_ρ(T::Real, p::Real);

# Fluid parameters
fluid_name = "Hydrogen"
T_in = 70               # K
p_in = 350e5           # Pa

# MTK model
@mtkmodel DensityCalc begin

    @parameters begin
        T
        p
    end

    @variables begin
        ρ(t)
    end

    @equations begin
        ρ ~ calc_ρ(T, p)
    end
end

# Solve
@mtkbuild sys = DensityCalc(; T=T_in, p=p_in)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, Tsit5())

# Display
println(sol[sys.ρ][end])
println(typeof(fluid_name))

67.30013202474123
String

2. This works using Clapeyron

## Clapeyron working

# Helper functions
function calc_ρ_Clap(T, p)
    return mass_density(fluid_model, p, T)
end

@register_symbolic calc_ρ_Clap(T::Real, p::Real);

# Fluid parameters
fluid_model = SingleFluid("Hydrogen")
T_in = 70               # K
p_in = 350e5           # Pa

# MTK model
@mtkmodel DensityCalc begin

    @parameters begin
        T
        p
    end

    @variables begin
        ρ(t)
    end

    @equations begin
        ρ ~ calc_ρ_Clap(T, p)
    end
end

# Solve
@mtkbuild sys = DensityCalc(; T=T_in, p=p_in)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, Tsit5())

# Display
println(sol[sys.ρ][end])
println(typeof(fluid_model))

3. Coolprop not working

## Coolprop not working

# Helper functions
function calc_ρ(T, p, fluid_name)
    return PropsSI("D", "T", T, "P", p, fluid_name)
end

@register_symbolic calc_ρ(T::Real, p::Real, fluid_name);

# Fluid parameters
fluid_name = "Hydrogen"
T_in = 70               # K
p_in = 350e5           # Pa

# MTK model
@mtkmodel DensityCalc begin

    @parameters begin
        T
        p
        fluid_name
    end

    @variables begin
        ρ(t)
    end

    @equations begin
        ρ ~ calc_ρ(T, p, fluid_name)
    end
end

# Solve
@mtkbuild sys = DensityCalc(; T=T_in, p=p_in, fluid_name=fluid_name)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, Tsit5())

# Display
println(sol[sys.ρ][end])
println(typeof(fluid_name))

TypeError: in keyword argument fluid_name, expected Union{Missing, Nothing, ModelingToolkit.NoValue, Real}, got a value of type String

Stacktrace:
[1] (::ModelingToolkit.Model{typeof(DensityCalc), Dict{Symbol, Any}})(; kw::@Kwargs{name::Symbol, T::Int64, p::Float64, fluid_name::String})
@ ModelingToolkit C:\Users\matth.julia\packages\ModelingToolkit\aau6A\src\systems\model_parsing.jl:25
[2] top-level scope
@ C:\Users\matth.julia\packages\ModelingToolkit\aau6A\src\systems\abstractsystem.jl:2268

4. Clapeyron not working

## Clapeyron not working

# Helper functions
function calc_ρ_Clap(T, p, fluid_model)
    return mass_density(fluid_model, p, T)
end

@register_symbolic calc_ρ_Clap(T::Real, p::Real);

# Fluid parameters
fluid_model = SingleFluid("Hydrogen")
T_in = 70               # K
p_in = 350e5           # Pa

# MTK model
@mtkmodel DensityCalc begin

    @parameters begin
        T
        p
        fluid_model
    end

    @variables begin
        ρ(t)
    end

    @equations begin
        ρ ~ calc_ρ_Clap(T, p, fluid_model)
    end
end

# Solve
@mtkbuild sys = DensityCalc(; T=T_in, p=p_in, fluid_model=fluid_model)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, Tsit5())

# Display
println(sol[sys.ρ][end])
println(typeof(fluid_model))

TypeError: in keyword argument fluid_model, expected Union{Missing, Nothing, ModelingToolkit.NoValue, Real}, got a value of type SingleFluid{EmpiricAncillary}

Stacktrace:
[1] (::ModelingToolkit.Model{typeof(DensityCalc), Dict{Symbol, Any}})(; kw::@Kwargs{name::Symbol, T::Int64, p::Float64, fluid_model::SingleFluid{EmpiricAncillary}})
@ ModelingToolkit C:\Users\matth.julia\packages\ModelingToolkit\aau6A\src\systems\model_parsing.jl:25
[2] top-level scope
@ C:\Users\matth.julia\packages\ModelingToolkit\aau6A\src\systems\abstractsystem.jl:2268

2 
    @parameters begin
        T
        p
        fluid_name::String
    end
2 Likes
3

Just a mention about CoolProp. when calling PropsSI(prop,spec1,val1,spec2,val2,fluid), you are actually building the fluid model each time. there is a low level interface, but it does not have an ergonomic api.

4

Yes, I would like to make use of the low-level API. Do you have any idea how I could pass the ‘abstract object’ into MTK? This goes for both CoolProp Low-level and Clapeyron I guess. I would need to do something like this:

CoolProp high level: Simplest but slow:

    @parameters begin
        T
        p
        fluid_name::String
    end

Clapeyron and CP low-level:

    @parameters begin
        T
        p
        fluid_model::AbstractObject???
    end

Is passing a custom object type into MTK possible?

5

for Clapeyron. the type of SingleFluid is SingleFluid{EmpiricAncillary} (you can check it via typeof(model).

For coolprop, the low level interface works via handles. a handle is just a CInt that points to the model in memory, if i remember correctly. The low level api with your example is something similar to this:

function calc_ρ(handle,T,p)
  inputs =CoolProp.get_input_pair_index("PT_INPUTS") #it is just an integer
  prop = get_param_index("D") #it is just an integer
  #AbstractState_set_fractions(handle, [0.4, 0.6]); for multicomponent models
  CoolProp.AbstractState_update(handle, inputs, p, T) #this is the part that does the calculation
  return CoolProp.AbstractState_keyed_output(handle, prop)
end

handle = CoolProp.AbstractState_factory("HEOS", "Hydrogen")

Clapeyron (in the particular case of SingleFluid eos model)

    @parameters begin
        T
        p
        fluid_model::SingleFluid{EmpiricAncillary}
    end

Coolprop (low level api)

    @parameters begin
        T
        p
        handle::Cint
    end
1 Like
6

@longemen3000 @ChrisRackauckas
Thank you so much for the input, the full solution is as follows:

CoolProp High-Level:

# Helper functions
function calc_ρ(T, p, fluid_name)
    return PropsSI("D", "T", T, "P", p, fluid_name)
end

@register_symbolic calc_ρ(T::Real, p::Real, fluid_name::String);

# MTK model
@mtkmodel DensityCalc begin

    @parameters begin
        T
        p
        fluid_name::String
    end

    @variables begin
        ρ(t)
    end

    @equations begin
        ρ ~ calc_ρ(T, p, fluid_name)
    end
end

# Fluid parameters
fluid_name = "Hydrogen"
T_in = 70               # K
p_in = 350e5           # Pa

# Solve
@mtkbuild sys = DensityCalc(; T=T_in, p=p_in, fluid_name=fluid_name)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, Tsit5())

# Display
println(sol[sys.ρ][end])

67.30013202474123

Clapeyron:

# Helper functions
function calc_ρ_Clap(T, p, fluid_model)
    return mass_density(fluid_model, p, T)
end

# Register the function with the type annotation
@register_symbolic calc_ρ_Clap(T::Real, p::Real, fluid_model::SingleFluid);

# MTK model
@mtkmodel DensityCalc begin
    @parameters begin
        T
        p
        fluid_model::SingleFluid
    end

    @variables begin
        ρ(t)
    end

    @equations begin
        ρ ~ calc_ρ_Clap(T, p, fluid_model)
    end
end

# Fluid parameters
fluid_model = SingleFluid("Hydrogen")
T_in = 70               # K
p_in = 350e5           # Pa

# Solve
@mtkbuild sys = DensityCalc(; T=T_in, p=p_in, fluid_model=fluid_model)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, Tsit5())

# Display
println(sol[sys.ρ][end])

67.30013202474127

1 Like
7

Yes, MTK parameters can be declared to be of any Julia type.

1 Like
8

There is not by any chance a more general Julia type for fluid_model in Clapeyron? For example, if I would like to pass a SingleFluid today and a Peng-Robinson tomorrow, there is not a more general ‘Clapeyron fluid model type’ that could handle all of these?

9

All Clapeyron models are subtypes of EosModel. With that said, i don’t know if MTK supports abstract types in the parameter declaration

1 Like
10

It does but it has performance implications, because then it won’t strictly type the container, so it probably shouldn’t be done.

1 Like