adjustments for new FMIImport/FMICore versions (#187)

Project.toml

@@ -1,14 +1,15 @@
 name = "FMI"
 uuid = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
 authors = ["TT <tobias.thummerer@informatik.uni-augsburg.de>", "LM <lars.mikelsons@informatik.uni-augsburg.de>", "JK <josef.kircher@student.uni-augsburg.de>"]
-version = "0.12.3"
+version = "0.12.4"
 
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FMIExport = "31b88311-cab6-44ed-ba9c-fe5a9abbd67a"
 FMIImport = "9fcbc62e-52a0-44e9-a616-1359a0008194"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 ThreadPools = "b189fb0b-2eb5-4ed4-bc0c-d34c51242431"
@@ -17,7 +18,7 @@ ThreadPools = "b189fb0b-2eb5-4ed4-bc0c-d34c51242431"
 DiffEqCallbacks = "2.26.0"
 DifferentialEquations = "7.7.0"
 FMIExport = "0.2.0"
-FMIImport = "0.15.6"
+FMIImport = "0.15.7"
 ProgressMeter = "1.7.0"
 Requires = "1.3.0"
 ThreadPools = "2.1.1"

docs/src/faq.md

@@ -3,22 +3,24 @@
 
 This list some common - often numerical - errors, that can be fixed by better understanding the ODE-Problem inside your FMU.
 
+----------
 
-
-## Solving non-linear system
+## Solving non-linear system failes
 ### Description
 Error message or warning, that solving of a non-linear system failed, close to the simulation start time.
 
 ### Example
-- `Solving non-linear system 101 failed at time=3e-05.`
+```
+Solving non-linear system 101 failed at time=3e-05.
+```
 
 ### Reason
 This could be, because the first step of the integration is accepted by the solver's error estimation, but shouldn't. This is usually, if the first step is picked to large by the solver's start step size heuristics.
 
 ### Fix
 - Try a small start value for the integration with keyword `dt`.
 
-
+----------
 
 ## Access denied
 ### Description
@@ -35,4 +37,29 @@ Access denied
 This is because your OS doesn't allow to interact with the binaries shipped with *FMI.jl*. 
 
 ### Fix
-This can easily be solved by fixing the binary's permission options or is automatically fixed if Julia runs with admin privileges.
+This can easily be solved by fixing the binary's permission options or is automatically fixed if Julia runs with admin privileges.
+
+----------
+
+## Double Callback Crossing
+### Description
+Error message, that solving failed because of `double callback crossing`.
+
+### Example 
+```
+ERROR:
+Double callback crossing floating pointer reducer errored. Report this issue.
+```
+
+### Reason
+This is because the event instant (time point) of an FMU event indicator can't be found precisely.
+
+### Fix
+This can be solved by allowing for more interpolation points during searching of the zero crossing:
+```
+fmu.executionConfig.rootSearchInterpolationPoints = 1000 # default value is 10
+```
+This will have negative performance impact on systems with extreme amount of events (thousands per second). 
+For systems with only a few events there won't be a noticable slow down. 
+
+----------

src/FMI2/sim.jl

@@ -15,7 +15,9 @@ import FMIImport: prepareSolveFMU, finishSolveFMU, handleEvents
 import FMIImport: undual
 
 using FMIImport.ChainRulesCore
-import FMIImport.ForwardDiff
+
+import FMIImport.ReverseDiff 
+import LinearAlgebra: eigvals
 
 import ProgressMeter
 import ThreadPools
@@ -168,6 +170,64 @@ function saveValues(c::FMU2Component, recordValues, x, t, integrator, inputFunct
     return (fmiGet(c, recordValues)...,)
 end
 
+function saveEventIndicators(c::FMU2Component, recordEventIndicators, x, t, integrator, inputFunction, inputValues)
+
+    @assert c.state == fmi2ComponentStateContinuousTimeMode "saveEventIndicators(...): Must be in continuous time mode."
+
+    c.solution.evals_saveeventindicators += 1
+
+    #x_old = fmi2GetContinuousStates(c)
+    #t_old = c.t
+
+    if !c.fmu.isZeroState
+        fmi2SetContinuousStates(c, x)
+    end
+    fmi2SetTime(c, t) 
+    if inputFunction != nothing
+        fmi2SetReal(c, inputValues, inputFunction(c, x, t)) 
+    end
+
+    #fmi2SetContinuousStates(c, x_old)
+    #fmi2SetTime(c, t_old)
+
+    out = zeros(fmi2Real, c.fmu.modelDescription.numberOfEventIndicators)
+    fmi2GetEventIndicators!(c, out)
+
+    return (out[recordEventIndicators]...,)
+end
+
+function saveEigenvalues(c::FMU2Component, x, t, integrator, inputFunction, inputValues)
+
+    @assert c.state == fmi2ComponentStateContinuousTimeMode "saveEigenvalues(...): Must be in continuous time mode."
+
+    c.solution.evals_saveeigenvalues += 1
+
+    #x_old = fmi2GetContinuousStates(c)
+    #t_old = c.t
+
+    if !c.fmu.isZeroState
+        fmi2SetContinuousStates(c, x)
+    end
+    fmi2SetTime(c, t) 
+    if inputFunction != nothing
+        fmi2SetReal(c, inputValues, inputFunction(c, x, t)) 
+    end
+
+    #fmi2SetContinuousStates(c, x_old)
+    #fmi2SetTime(c, t_old)
+
+    A = ReverseDiff.jacobian(_x -> FMI.fx(c, _x, [], t), x)
+    eigs = eigvals(A)
+
+    vals = []
+    for e in eigs 
+        push!(vals, real(e))
+        push!(vals, imag(e))
+    end
+
+    return (vals...,)
+end
+
 function fx(c::FMU2Component, 
     dx::AbstractArray{<:Real},
     x::AbstractArray{<:Real}, 
@@ -246,6 +306,8 @@ Keywords:
     - solver: Any Julia-supported ODE-solver (default is the DifferentialEquations.jl default solver, currently `AutoTsit5(Rosenbrock23())`)
     - customFx: [deprecated] Ability to give a custom state derivative function ẋ=f(x,t)
     - recordValues: Array of variables (strings or variableIdentifiers) to record. Results are returned as `DiffEqCallbacks.SavedValues`
+    - recordEventIndicators: Array or Range of event indicators (identified by integers) to record
+    - recordEigenvalues: Boolean value, if eigenvalues shall be computed and recorded
     - saveat: Time points to save values at (interpolated)
     - setup: Boolean, if FMU should be setup (default=true)
     - reset: Union{Bool, :auto}, if FMU should be reset before simulation (default reset=:auto)
@@ -264,6 +326,8 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp
     solver = nothing,
     customFx = nothing,
     recordValues::fmi2ValueReferenceFormat = nothing,
+    recordEventIndicators::Union{AbstractArray{<:Integer, 1}, UnitRange{<:Integer}, Nothing} = nothing,
+    recordEigenvalues::Bool=false,
     saveat = nothing,
     x0::Union{AbstractArray{<:Real}, Nothing} = nothing,
     setup::Union{Bool, Nothing} = nothing,
@@ -275,7 +339,8 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp
     inputFunction = nothing,
     parameters::Union{Dict{<:Any, <:Any}, Nothing} = nothing,
     dtmax::Union{Real, Nothing} = nothing,
-    callbacks = [],
+    callbacksBefore = [],
+    callbacksAfter = [],
     showProgress::Bool = true,
     kwargs...)
 
@@ -287,11 +352,11 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp
     if inputFunction != nothing
         if hasmethod(inputFunction, Tuple{fmi2Real})
             _inputFunction = (c, u, t) -> inputFunction(t)
-        elseif  hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, fmi2Real})
+        elseif hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, fmi2Real})
             _inputFunction = (c, u, t) -> inputFunction(c, t)
-        elseif  hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, AbstractArray{fmi2Real,1}})
+        elseif hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, AbstractArray{fmi2Real,1}})
             _inputFunction = (c, u, t) -> inputFunction(c, u)
-        elseif  hasmethod(inputFunction, Tuple{AbstractArray{fmi2Real,1}, fmi2Real})
+        elseif hasmethod(inputFunction, Tuple{AbstractArray{fmi2Real,1}, fmi2Real})
             _inputFunction = (c, u, t) -> inputFunction(u, t)
         else 
             _inputFunction = inputFunction
@@ -303,6 +368,7 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp
     inputValueReferences = prepareValueReference(fmu, inputValueReferences)
 
     savingValues = (length(recordValues) > 0)
+    savingEventIndicators = !isnothing(recordEventIndicators) 
     hasInputs = (length(inputValueReferences) > 0)
     hasParameters = (parameters !== nothing)
     hasStartState = (x0 !== nothing)
@@ -311,7 +377,7 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp
 
     cbs = []
 
-    for cb in callbacks
+    for cb in callbacksBefore
         push!(cbs, cb)
     end
 
@@ -416,11 +482,12 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp
         push!(cbs, stepCb)
     end
 
-    if savingValues
+    if savingValues 
         dtypes = collect(fmi2DataTypeForValueReference(c.fmu.modelDescription, vr) for vr in recordValues)
-        fmusol.values = SavedValues(Float64, Tuple{dtypes...})
+        fmusol.values = SavedValues(fmi2Real, Tuple{dtypes...})
         fmusol.valueReferences = copy(recordValues)
 
+        savingCB = nothing
         if saveat === nothing
             savingCB = SavingCallback((u,t,integrator) -> saveValues(c, recordValues, u, t, integrator, _inputFunction, inputValueReferences), 
                                     fmusol.values)
@@ -433,6 +500,45 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp
         push!(cbs, savingCB)
     end
 
+    if savingEventIndicators
+        dtypes = collect(fmi2Real for ei in recordEventIndicators)
+        fmusol.eventIndicators = SavedValues(fmi2Real, Tuple{dtypes...})
+        fmusol.recordEventIndicators = copy(recordEventIndicators)
+
+        savingCB = nothing
+        if saveat === nothing
+            savingCB = SavingCallback((u,t,integrator) -> saveEventIndicators(c, recordEventIndicators, u, t, integrator, _inputFunction, inputValueReferences), 
+                                    fmusol.eventIndicators)
+        else
+            savingCB = SavingCallback((u,t,integrator) -> saveEventIndicators(c, recordEventIndicators, u, t, integrator, _inputFunction, inputValueReferences), 
+                                    fmusol.eventIndicators, 
+                                    saveat=saveat)
+        end
+
+        push!(cbs, savingCB)
+    end
+
+    if recordEigenvalues
+        dtypes = collect(Float64 for _ in 1:2*length(c.fmu.modelDescription.stateValueReferences))
+        fmusol.eigenvalues = SavedValues(fmi2Real, Tuple{dtypes...})
+
+        savingCB = nothing
+        if saveat === nothing
+            savingCB = SavingCallback((u,t,integrator) -> saveEigenvalues(c, u, t, integrator, _inputFunction, inputValueReferences), 
+                                    fmusol.eigenvalues)
+        else
+            savingCB = SavingCallback((u,t,integrator) -> saveEigenvalues(c, u, t, integrator, _inputFunction, inputValueReferences), 
+                                    fmusol.eigenvalues, 
+                                    saveat=saveat)
+        end
+
+        push!(cbs, savingCB)
+    end
+
+    for cb in callbacksAfter
+        push!(cbs, cb)
+    end
+
     # if auto_dt == true
     #     @assert solver !== nothing "fmi2SimulateME(...): `auto_dt=true` but no solver specified, this is not allowed."
     #     tmpIntegrator = init(c.problem, solver)

src/extensions/Plots.jl

@@ -41,7 +41,7 @@ function fmiPlot(solution::FMUSolution; kwargs...)
     fmiPlot!(fig, solution; kwargs...)
     return fig
 end
-function fmiPlot!(fig, solution::FMUSolution;
+function fmiPlot!(fig::Plots.Plot, solution::FMUSolution;
     states::Union{Bool, Nothing}=nothing,
     values::Union{Bool, Nothing}=nothing,
     stateEvents::Union{Bool, Nothing}=nothing,
@@ -209,9 +209,9 @@ Extended the original plot-command by plotting FMUs.
 
 For further information seek `?fmiPlot`.
 """
-function Plots.plot(solution::FMUSolution, args...; kwargs...)
-    fmiPlot(solution, args...; kwargs...)
+function Plots.plot(solution::FMUSolution; kwargs...)
+    fmiPlot(solution; kwargs...)
 end
-function Plots.plot!(fig, solution::FMUSolution, args...; kwargs...)
-    fmiPlot!(fig, solution, args...; kwargs...)
+function Plots.plot!(fig::Plots.Plot, solution::FMUSolution; kwargs...)
+    fmiPlot!(fig, solution; kwargs...)
 end