add RealtimeClock

src/clock.jl

@@ -98,11 +98,14 @@
 abstract type AbstractClock <: AbstractDiscrete end
 
 """
-    Clock <: AbstractClock
-    Clock([t]; dt)
+    Clock()
+    Clock(dt)
+    Clock(t, dt)
 
 The default periodic clock with independent variables `t` and tick interval `dt`.
 If `dt` is left unspecified, it will be inferred (if possible).
+
+See also [`RealtimeClock`](@ref).
 """
 struct Clock <: AbstractClock
     "Independent variable"
@@ -122,3 +125,25 @@
 end
 
 is_concrete_time_domain(x) = x isa Union{AbstractClock, Continuous}
+
+"""
+    RealtimeClock()
+    RealtimeClock(dt)
+    RealtimeClock(t, dt)
+
+Similar to [`Clock`](@ref), but with with the additional property that the simulation is run no faster than real-time, i.e., a simulation with `tspan = (0, 10)` will take at least 10 seconds of wallclock time to run.
+
+To achieve real-time execution, the wall-clock duration of each integration step is measured, and a call to `Libc.systemsleep` is made to ensure that the next integration step is not started before `dt` seconds of wall-clock time has elapsed. This leads to a simulation that takes at least as long as the length of the time span, but may take longer if the computational load is high enough that one integration step takes more than `dt` seconds.
+"""
+struct RealtimeClock <: AbstractClock
+    clock::Clock
+    RealtimeClock(args...) = new(Clock(args...))
+end
+
+sampletime(c) = sampletime(c.clock)
+Base.hash(c::RealtimeClock, seed::UInt) = hash(c.dt, seed ⊻ 0x9d3d7a9a18874b90)
+Base.:(==)(c1::RealtimeClock, c2::RealtimeClock) = c1.clock == c2.clock
+function Base.getproperty(c::RealtimeClock, s::Symbol)
+    s ∈ fieldnames(typeof(c)) && return getfield(c, s)
+    getproperty(getfield(c, :clock), s)
+end

src/systems/clock_inference.jl

@@ -236,3 +236,28 @@
     defaults = Dict{Any, Any}(v => 0.0 for v in Iterators.flatten(inputs))
     return affects, clocks, svs, appended_parameters, defaults
 end
+
+function generate_discrete_callback(affect, clock, sv)
+    error("$clock is not a supported clock type.")
+end # Fallback
+
+function generate_discrete_callback(affect, clock::Clock, sv)
+    PeriodicCallback(DiscreteSaveAffect(affect, sv), clock.dt)
+end
+
+function generate_discrete_callback(affect, clock::RealtimeClock, sv)
+    start_time = time() # Set start time before simulation starts
+    function realtime_affect!(integrator, saved_values)
+        # This is a closure that closes over start_time
+        affect(integrator, saved_values)
+        execution_time = time() - start_time
+        sleep_time = clock.dt - execution_time
+        if sleep_time < 0
+            @debug "RealtimeClock failed to meet deadline dt = $(clock.dt), execution time = $execution_time at simulation time t = $(t)"
+        end
+        Libc.systemsleep(max(0, sleep_time))
+        start_time = time() # Update start time to current time
+        nothing
+    end
+    PeriodicCallback(DiscreteSaveAffect(realtime_affect!, sv), clock.dt)
+end

src/systems/diffeqs/abstractodesystem.jl

@@ -947,13 +947,7 @@
     cbs = process_events(sys; callback, has_difference, kwargs...)
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
         affects, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...)
-        discrete_cbs = map(affects, clocks, svs) do affect, clock, sv
-            if clock isa Clock
-                PeriodicCallback(DiscreteSaveAffect(affect, sv), clock.dt)
-            else
-                error("$clock is not a supported clock type.")
-            end
-        end
+        discrete_cbs = map(ModelingToolkit.generate_discrete_callback, affects, clocks, svs)
         if cbs === nothing
             if length(discrete_cbs) == 1
                 cbs = only(discrete_cbs)

test/clock.jl

@@ -310,3 +310,26 @@ if VERSION >= v"1.7"
 
     @test sol.u ≈ sol2.u
 end
+
+# RealtimeClock
+dt = 0.01
+@variables t x(t)=0 y(t)=0 u(t)=0 yd(t)=0 ud(t)=0 r(t)=1
+@parameters kp = 1
+D = Differential(t)
+d = ModelingToolkit.RealtimeClock(dt)
+# u(n + 1) := f(u(n))
+
+eqs = [yd ~ Sample(d)(y)
+    ud ~ kp * (r - yd)
+    r ~ 1.0
+
+# plant (time continuous part)
+    u ~ Hold(ud)
+    D(x) ~ -x + u
+    y ~ x]
+@named sys = ODESystem(eqs)
+prob = ODEProblem(structural_simplify(sys), [], (0.0, 3.0))
+solve(prob, Tsit5(), kwargshandle = KeywordArgSilent)
+execution_time = @elapsed solve(prob, Tsit5(), kwargshandle = KeywordArgSilent)
+@test_skip execution_time≈3 atol=10 * dt # Very crude test that can be evaluated locally but not on CI
+@test execution_time >= 3 - dt # This can be evaluated on CI as well since it only tests a lower bound