11# ── Phase 3: ODE simulation with DifferentialEquations / MTK ──────────────────
22
33import DifferentialEquations
4- import LinearAlgebra
5- import OrdinaryDiffEqBDF
64import Logging
75import ModelingToolkit
86import Printf: @sprintf
@@ -16,7 +14,7 @@ const _SIM_SETTINGS = SimulateSettings(solver = DifferentialEquations.Rodas5Pr()
1614Update the module-level simulation settings in-place and return them.
1715
1816# Keyword arguments
19- - `solver` — any SciML ODE/DAE algorithm instance (e.g. `FBDF() `, `Rodas5Pr ()`).
17+ - `solver` — any SciML ODE/DAE algorithm instance (e.g. `Rodas5Pr `, `FBDF ()`).
2018- `saveat_n` — number of uniform time points for purely algebraic systems.
2119
2220# Example
@@ -63,6 +61,7 @@ function run_simulate(ode_prob,
6361 settings :: SimulateSettings = _SIM_SETTINGS,
6462 cmp_signals :: Vector{String} = String[],
6563 csv_max_size_mb:: Int = CSV_MAX_SIZE_MB):: Tuple{Bool,Float64,String,Any}
64+
6665 sim_success = false
6766 sim_time = 0.0
6867 sim_error = " "
@@ -84,31 +83,13 @@ function run_simulate(ode_prob,
8483 # internal steps and sol.t would be empty with saveat=[].
8584 # Supply explicit time points so observed variables can be evaluated.
8685 sys = ode_prob. f. sys
87- M = ode_prob. f. mass_matrix
88- unknowns = ModelingToolkit. unknowns (sys)
89- n_unknowns = length (unknowns)
90- n_diff = if M isa LinearAlgebra. UniformScaling
91- n_unknowns
92- else
93- count (! iszero, LinearAlgebra. diag (M))
94- end
86+ n_unknowns = length (ModelingToolkit. unknowns (sys))
9587
9688 kwargs = if n_unknowns == 0
97- # No unknowns at all (e.g. BusUsage): the solver takes no
98- # internal steps with saveat=[], leaving sol.t empty.
99- # Use a fixed grid + adaptive=false so observed variables
100- # can be evaluated.
101- t0_s, t1_s = ode_prob. tspan
102- saveat_s = collect (range (t0_s, t1_s; length = settings. saveat_n))
103- dt_s = saveat_s[2 ] - saveat_s[1 ]
104- (saveat = saveat_s, adaptive = false , dt = dt_s, dense = false )
105- elseif n_diff == 0
106- # Algebraic unknowns only (e.g. CharacteristicIdealDiodes):
107- # the solver must take adaptive steps to track discontinuities.
108- # Keep saveat=[] + dense=true so the solver drives its own
109- # step selection; dense output is unreliable but the solution
110- # values at each step are correct.
111- (saveat = Float64[], dense = true )
89+ # No unknowns at all (e.g. BusUsage):
90+ # Supply explicit time points so observed variables can be evaluated.
91+ saveat_s = collect (range (ode_prob. tspan[1 ], ode_prob. tspan[end ]; length = settings. saveat_n))
92+ (saveat = saveat_s, dense = true )
11293 else
11394 (saveat = Float64[], dense = true )
11495 end
@@ -119,9 +100,9 @@ function run_simulate(ode_prob,
119100 # Use our own `saveat` vector for the log: integ.opts.saveat is a
120101 # BinaryHeap which does not support iterate/minimum/maximum.
121102 integ = DifferentialEquations. init (ode_prob, solver; kwargs... )
122- saveat = kwargs. saveat
103+ saveat_s = kwargs. saveat
123104 solver_settings_string = if hasproperty (integ, :opts )
124- sv_str = isempty (saveat ) ? " []" : " $(length (saveat )) points in [$(first (saveat )) , $(last (saveat )) ]"
105+ sv_str = isempty (saveat_s ) ? " []" : " $(length (saveat_s )) points in [$(first (saveat_s )) , $(last (saveat_s )) ]"
125106 """
126107 Solver $(parentmodule (typeof (solver))) .$(nameof (typeof (solver)))
127108 saveat: $sv_str
@@ -131,14 +112,14 @@ function run_simulate(ode_prob,
131112 dense: $(integ. opts. dense)
132113 """
133114 else
134- sv_str = isempty (saveat ) ? " []" : " $(length (saveat )) points in [$(first (saveat )) , $(last (saveat )) ]"
115+ sv_str = isempty (saveat_s ) ? " []" : " $(length (saveat_s )) points in [$(first (saveat_s )) , $(last (saveat_s )) ]"
135116 " Solver (NullODEIntegrator — no unknowns)
136117 saveat: $sv_str
137118 dense: true"
138119 end
139120
140121 # Solve
141- DifferentialEquations. solve (ode_prob, OrdinaryDiffEqBDF . FBDF () ; kwargs... )
122+ DifferentialEquations. solve (ode_prob, solver ; kwargs... )
142123 end
143124 sim_time = time () - t0
144125 if sol. retcode == DifferentialEquations. ReturnCode. Success
0 commit comments