Updated code for acausal built-in components + added new test model TestHeatTransfer2 where code size is independent from number of temperature nodes.

Author: MartinOtter

Project.toml

@@ -1,7 +1,7 @@
 authors = ["Hilding Elmqvist <[email protected]>", "Martin Otter <[email protected]>"]
 name = "Modia"
 uuid = "cb905087-75eb-5f27-8515-1ce0ec8e839e"
-version = "0.9.4"
+version = "0.10.0-dev"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"

docs/src/index.md

@@ -42,6 +42,17 @@ functionalities of these packages.
 
 ## Release Notes
 
+### Version 0.10.0
+
+- 
+
+**Non-backwards** compatible changes (should usually not influende user models)
+
+- `_buildFunction = <functionName>` changed to `_buildFunction = Par(functionName = <functionName>)` and
+  additional argument `unitless` added to `<functionName>`.
+- `_instantiateFunction = Par(..)` changed to `_initSegmentFunction = Par(..)`
+
+
 ### Version 0.9.4
 
 - Precompile statements included (compilation of Modia package takes more time, but startup of Modia model simulations is faster).

models/HeatTransfer.jl

@@ -70,4 +70,40 @@ InsulatedRod = Model(
     ]
 )
 
+
+include("HeatTransfer/InsulatedRod2.jl")
+
+"""
+    insulatedRod = InsulatedRod2(; L, A, rho=7500.0u"kg/m^3", lambda=74.0u"W/(m*K)", c=450.0u"J/(kg*K), T0=293.15u"K", nT=1") 
     
+Generate a Model(..) instance of an insulated rod with length `L` and cross sectional area `A` that
+models 1D heat transfer from port_a to port_b. The rod is discretized with `nT` internal temperature nodes that are initialized with `T0`.
+This is an acausal built-in component where the code size does not depend on `nT` and
+`nT` can still be changed after model code is generated and compiled.
+
+For more details see Appendix B1 of [DOI: 10.3390/electronics12030500](https://doi.org/10.3390/electronics12030500).
+
+# Optional Arguments:
+- `rho`: density of rod material
+- `lambda`: thermal conductivity of rod material
+- `c`: specific heat capacity of rod material
+- `T0`: initial temperature of internal temperature nodes
+- `nT`: number of temperature nodes (`nT > 0`).
+
+# Internal variables that can be inquired/plotted:
+- `T`: Vector of temperatures at the internal temperature nodes
+- `der(T)': Vector of derivatives of `T`.
+"""
+InsulatedRod2 = Model(; 
+    _buildFunction       = Par(functionName = :(buildInsulatedRod2!)),              # Called once in @instantiateModel(..) before getDerivatives!(..) is generated 
+    _initSegmentFunction = Par(functionName = :(initInsulatedRod2ForNewSegment!)),  # Called once before initialization of a new simulation segment
+    L      = 1.0u"m",                # Length of rod
+    A      = 0.0004u"m^2",           # Rod area
+    rho    = 7500.0u"kg/m^3",        # Density of rod material
+    lambda = 74.0u"W/(m*K)",         # Thermal conductivity of rod material
+    c      = 450.0u"J/(kg*K)",       # Specific heat capacity of rod material
+    T0     = 293.15u"K",             # Initial temperature of internal temperature nodes
+    nT     = 1,                      # Number of temperature nodes (nT > 0).
+    port_a  = HeatPort,              # heat port on left side
+    port_b  = HeatPort               # heat port on right side
+)

models/HeatTransfer/InsulatedRod2.jl

@@ -0,0 +1,138 @@
+# License for this file: MIT (expat)
+# Copyright 2023, DLR Institute of System Dynamics and Control
+#
+# This file is included in Modia/models/HeatTranser.jl
+
+
+# Structure holding the internal memory of the insulated rod
+mutable struct InsulatedRodStruct{FloatType}
+  # Parameters
+  path::String               # path name of instance
+  Ge::FloatType              # = lambda*A/dx
+  Ge2::FloatType             # = 2*Ge
+  k::FloatType               # = Ge / (c*rho*A*dx)
+  T_init::Vector{FloatType}  # initial states
+
+  # Internal states and state deriatives
+  T::Vector{FloatType}       # states
+  der_T::Vector{FloatType}   # state derivatives
+
+  # Start index of state and state derivative vectors
+  T_startIndex::Int
+
+  InsulatedRodStruct{FloatType}(path::String) where {FloatType} = new(path,false)
+end
+
+
+# Called once before initialization of first simulation segment
+function initInsulatedRod2!(obj::InsulatedRodStruct{FloatType}; L, A, rho, lambda, c, T0, nT)::Nothing where {FloatType}
+    #L, A, rho=7500.0u"kg/m^3", lambda=74.0u"W/(m*K)", c=450.0u"J/(kg*K)", T0=293.15u"K", nT=1)::Nothing where {FloatType}
+
+    # Convert to SI units, strip units and check that values are positives
+    path = obj.path
+    #=
+    @show L
+    @Modia.strippedPositive!(path, L)
+    @Modia.strippedPositive!(path, A)
+    @Modia.strippedPositive!(path, rho)
+    @Modia.strippedPositive!(path, lambda)
+    @Modia.strippedPositive!(path, c)
+    @Modia.strippedPositive!(path, T0)
+    @Modia.strippedPositive!(path, nT)
+    @show L
+    @show A
+    @show rho
+    @show lambda
+    =#
+    L      = stripUnit(L)
+    A      = stripUnit(A)
+    rho    = stripUnit(rho)
+    lambda = stripUnit(lambda)
+    c      = stripUnit(c)
+    T0     = stripUnit(T0)
+
+    # Compute derived parameters
+    dx         = L/nT
+    obj.Ge     = lambda*A/dx
+    obj.Ge2    = 2*obj.Ge
+    obj.k      = obj.Ge / ( c*rho*A*dx )
+    obj.T_init = fill(T0, nT)
+
+    # Allocate and initialize internal states
+    obj.T      = copy(obj.T_init)
+    obj.der_T  = zeros(nT)
+    return nothing
+end
+
+
+# Called from @instantiateModel(..) before getDerivatives!(..) is generated
+function buildInsulatedRod2!(model::AbstractDict, FloatType, TimeType, unitless::Bool, buildDict, path)
+    pathAsString = Modia.modelPathAsString(path)
+    extraModelCode = Model(  # extra code to be merged with model
+        insRod  = Var(hideResult=true),  # InsulatedRodStruct instance (an equation to return this variable is generated by _buildFunction)
+        success = Var(hideResult=true),  # Dummy return argument
+        equations = :[
+            # copy states into insRod.T and return insRod
+            insRod = openInsulatedRod!(instantiatedModel, $pathAsString)
+      
+            # equations at the boundaries
+            port_a.Q_flow = getGe2(insRod, Val($unitless))*(port_a.T - getT1(  insRod, Val($unitless)))
+            port_b.Q_flow = getGe2(insRod, Val($unitless))*(port_b.T - getTend(insRod, Val($unitless)))
+
+            # compute insRod.der_T and copy it into state derivatives
+            success = computeInsulatedRodDerivatives!(instantiatedModel, insRod, port_a.T, port_b.T)   # instantiatedModel is provided in the generated code
+        ])
+
+  # Store build info in buildDict
+  buildDict[pathAsString] = InsulatedRodStruct{FloatType}(pathAsString)
+  return extraModelCode
+end
+
+
+# Called once before initialization of a new simulation segment
+function initInsulatedRod2ForNewSegment!(instantiatedModel::SimulationModel{FloatType,TimeType},
+                                         parameters::AbstractDict, path::String; log=false)::Nothing where {FloatType,TimeType}
+    obj::InsulatedRodStruct{FloatType} = instantiatedModel.buildDict[path]
+
+    if Modia.isFirstInitialOfAllSegments(instantiatedModel)
+        initInsulatedRod2!(obj; parameters...)
+    end
+
+    # Define a vector of new state and state derivative variables
+    obj.T_startIndex = Modia.new_x_segmented_variable!(instantiatedModel, path*".T", path*".der(T)", obj.T, "K")
+    return nothing
+end
+
+
+# Open an initialized InsulatedRod2 model and return a reference to it
+function openInsulatedRod!(instantiatedModel::SimulationModel{FloatType,TimeType}, path::String)::InsulatedRodStruct{FloatType} where {FloatType,TimeType}
+    obj::InsulatedRodStruct{FloatType} = instantiatedModel.buildDict[path]
+    Modia.get_Vector_x_segmented_value!(instantiatedModel, obj.T_startIndex, obj.T)
+    return obj
+end
+
+
+# Functions to inquire values from InsulatedRodStruct
+getT1(  insRod::InsulatedRodStruct, Unitless::Val{true})  = insRod.T[1] 
+getT1(  insRod::InsulatedRodStruct, Unitless::Val{false}) = insRod.T[1]*u"K"
+
+getTend(insRod::InsulatedRodStruct, Unitless::Val{true})  = insRod.T[end] 
+getTend(insRod::InsulatedRodStruct, Unitless::Val{false}) = insRod.T[end]*u"K"
+
+getGe2( insRod::InsulatedRodStruct, Unitless::Val{true})  = insRod.Ge2
+getGe2( insRod::InsulatedRodStruct, Unitless::Val{false}) = insRod.Ge2*u"W/K"
+
+T_grad1(T,Ta,i) = if i == 1        ; (Ta - T[1])*2 else T[i-1] - T[i] end
+T_grad2(T,Tb,i) = if i == length(T); (T[i] - Tb)*2 else T[i] - T[i+1] end
+
+function computeInsulatedRodDerivatives!(instantiatedModel, obj::InsulatedRodStruct{FloatType}, Ta, Tb)::Bool where {FloatType}
+    Ta::FloatType = stripUnit(Ta)
+    Tb::FloatType = stripUnit(Tb)
+    T  = obj.T
+    k  = obj.k
+    for i in 1:length(T)
+        obj.der_T[i] = k*(T_grad1(T,Ta,i) - T_grad2(T,Tb,i))
+    end
+    Modia.add_der_x_segmented_value!(instantiatedModel, obj.T_startIndex, obj.der_T)
+    return true
+end

src/CodeGeneration.jl

@@ -336,7 +336,7 @@ mutable struct SimulationModel{FloatType,TimeType}
 
     # Available after propagateEvaluateAndInstantiate!(..) called
     instantiateFunctions::Vector{Tuple{Union{Expr,Symbol},OrderedDict{Symbol,Any},String}}
-                                                # All definitions `_instantiateFunction = Par(functionName = XXX)` in the model to call
+                                                # All definitions `_initSegmentFunction = Par(functionName = XXX)` in the model to call
                                                 # `XXX(instantiatedModel, submodel, submodelPath)` in the order occurring  during evaluation
                                                 # of the parameters where, instantiatedFunctions[i] = (XXX, submodel, submodelPath)
     nsegments::Int                               # Current simulation segment

src/EvaluateParameters.jl

@@ -210,7 +210,7 @@ function propagateEvaluateAndInstantiate2!(m::SimulationModel{FloatType,TimeType
     end
     current = OrderedDict{Symbol,Any}()   # should be Map()
 
-    # Determine, whether "parameters" has a ":_constructor"  or "_instantiateFunction" key and handle this specially
+    # Determine, whether "parameters" has a ":_constructor"  or "_initSegmentFunction" key and handle this specially
     constructor          = nothing
     instantiateFunction  = nothing
     usePath              = false
@@ -248,13 +248,13 @@ function propagateEvaluateAndInstantiate2!(m::SimulationModel{FloatType,TimeType
             end            
         end
 
-    elseif haskey(parameters, :_instantiateFunction)
-        # For example: obj = (_instantiateFunction = Par(functionName = :(instantiateLinearStateSpace!))
-        _instantiateFunction = parameters[:_instantiateFunction]
-        if haskey(_instantiateFunction, :functionName)
-            instantiateFunction = _instantiateFunction[:functionName]
+    elseif haskey(parameters, :_initSegmentFunction)
+        # For example: obj = (_initSegmentFunction = Par(functionName = :(instantiateLinearStateSpace!))
+        _initSegmentFunction = parameters[:_initSegmentFunction]
+        if haskey(_initSegmentFunction, :functionName)
+            instantiateFunction = _initSegmentFunction[:functionName]
         else
-            @warn "Model $path has key :_instantiateFunction but its value has no key :functionName"
+            @warn "Model $path has key :_initSegmentFunction but its value has no key :functionName"
         end
 
     elseif haskey(parameters, :value)
@@ -269,7 +269,7 @@ function propagateEvaluateAndInstantiate2!(m::SimulationModel{FloatType,TimeType
         if log
             println(" 2:    ... key = $k, value = $v")
         end
-        if k == :_constructor || k == :_buildFunction || k == :_buildOption || k == :_instantiateFunction || 
+        if k == :_constructor || k == :_buildFunction || k == :_buildOption || k == :_initSegmentFunction || 
            k == :_path || k == :_instantiatedModel || (k == :_class && !isnothing(constructor))
             if log
                 println(" 3:    ... key = $k")

src/EventHandler.jl

@@ -73,7 +73,7 @@ mutable struct EventHandler{FloatType,TimeType}
     nz::Int                 # Number of event indicators
     nzInvariant::Int          # Number of event indicators defined in visible model equations
                             # More event indicators can be defined by objects that are not visible in the generated code, i.e. nz >= nzInvariant
-                            # These event indicators are defined in propagateEvaluateAndInstantiate!(..) via _instantiateFunction(..)
+                            # These event indicators are defined in propagateEvaluateAndInstantiate!(..) via _initSegmentFunction(..)
     z::Vector{FloatType}    # Vector of event indicators (zero crossings). If one of z[i] passes
                             # zero, that is beforeEvent(z[i])*z[i] < 0, an event is triggered
                             # (allocated during instanciation according to nz).

src/Modia.jl

@@ -9,8 +9,8 @@ Main module of Modia.
 module Modia
 
 const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
-const Version = "0.9.4"
-const Date = "2023-05-21"
+const Version = "0.10.0"
+const Date = "2023-05-22"
 const modelsPath = joinpath(Modia.path, "models")
 
 print(" \n\nWelcome to ")
@@ -170,6 +170,35 @@ The function is defined as: `stripUnit(v) = ustrip.(upreferred.(v))`.
 """
 stripUnit(v) = ustrip.(upreferred.(v))
 
+
+"""
+    @strippedPositive!(path, name)
+
+Convert `name` to its preferred units (default are the SI units), strip units and check that value is positive.
+In case of error, use `string(name)` and `path` in the error message:
+
+# Example
+```
+using Unitful
+L1 =  2.0u"m"
+@strippedPositive!("insulatedRod", L1)
+# L1 = 2.0
+
+L2 = -2.0u"m"
+@strippedPositive!("insulatedRod", L2)
+# error message:
+# Error from
+#   insulatedRod = ...(..., L2 = -2.0u"m",...): L2 > 0 required.)
+```
+"""
+macro strippedPositive!(path, name)
+    nameAsString = string(name)
+    expr = :( $name = strippedPositive!($path, $(esc(name)), $nameAsString) )   
+    return expr 
+end      
+strippedPositive!(path::String, value, name) = stripUnit(value) > 0 ? stripUnit(value) : error("\nError from\n   $path = ...(..., $name = $value, ...): $name > 0 required")
+
+
 """
     str = modelPathAsString(modelPath::Union{Expr,Symbol,Nothing})