11# ── Phase 3: ODE simulation with DifferentialEquations / MTK ──────────────────
22
3- import DifferentialEquations: solve, Rodas5P, ReturnCode
3+ import DifferentialEquations
44import Logging
55import ModelingToolkit
66import Printf: @sprintf
77
8+ """ Module-level default simulation settings. Modify via `configure_simulate!`."""
9+ const _SIM_SETTINGS = SimulateSettings (solver = DifferentialEquations. Rodas5P ())
10+
11+ """
12+ configure_simulate!(; solver, saveat_n) → SimulateSettings
13+
14+ Update the module-level simulation settings in-place and return them.
15+
16+ # Keyword arguments
17+ - `solver` — any SciML ODE/DAE algorithm instance (e.g. `Rodas5P`, `FBDF()`).
18+ - `saveat_n` — number of uniform time points for purely algebraic systems.
19+
20+ # Example
21+
22+ ```julia
23+ using OrdinaryDiffEqBDF
24+ configure_simulate!(solver = FBDF())
25+ ```
26+ """
27+ function configure_simulate! (;
28+ solver :: Union{Any,Nothing} = nothing ,
29+ saveat_n :: Union{Int,Nothing} = nothing ,
30+ )
31+ isnothing (solver) || (_SIM_SETTINGS. solver = solver)
32+ isnothing (saveat_n) || (_SIM_SETTINGS. saveat_n = saveat_n)
33+ return _SIM_SETTINGS
34+ end
35+
36+ """
37+ simulate_settings() → SimulateSettings
38+
39+ Return the current module-level simulation settings.
840"""
9- run_simulate(ode_prob, model_dir, model; cmp_signals, csv_max_size_mb) → (success, time, error, sol)
41+ simulate_settings () = _SIM_SETTINGS
1042
11- Solve `ode_prob` with Rodas5P (stiff solver). On success, also writes the
43+ """
44+ run_simulate(ode_prob, model_dir, model; settings, cmp_signals, csv_max_size_mb) → (success, time, error, sol)
45+
46+ Solve `ode_prob` using the algorithm in `settings.solver`. On success, also writes the
1247solution as a CSV file `<Short>_sim.csv` in `model_dir`.
1348Writes a `<model>_sim.log` file in `model_dir`.
1449Returns `nothing` as the fourth element on failure.
@@ -20,36 +55,74 @@ of signals will be compared.
2055CSV files larger than `csv_max_size_mb` MiB are replaced with a
2156`<Short>_sim.csv.toobig` marker so that the report can note the omission.
2257"""
23- function run_simulate (ode_prob, model_dir:: String ,
58+ function run_simulate (ode_prob,
59+ model_dir:: String ,
2460 model:: String ;
25- cmp_signals :: Vector{String} = String[],
26- csv_max_size_mb:: Int = CSV_MAX_SIZE_MB):: Tuple{Bool,Float64,String,Any}
27- sim_success = false
28- sim_time = 0.0
29- sim_error = " "
30- sol = nothing
61+ settings :: SimulateSettings = _SIM_SETTINGS,
62+ cmp_signals :: Vector{String} = String[],
63+ csv_max_size_mb:: Int = CSV_MAX_SIZE_MB):: Tuple{Bool,Float64,String,Any}
64+
65+ sim_success = false
66+ sim_time = 0.0
67+ sim_error = " "
68+ sol = nothing
69+ solver_settings_string = " "
3170
3271 log_file = open (joinpath (model_dir, " $(model) _sim.log" " w"
3372 println (log_file, " Model: $model "
3473 logger = Logging. SimpleLogger (log_file, Logging. Debug)
3574 t0 = time ()
75+
76+ solver = settings. solver
3677 try
37- # Rodas5P handles stiff DAE-like systems well.
3878 # Redirect all library log output (including Symbolics/MTK warnings)
3979 # to the log file so they don't clutter stdout.
4080 sol = Logging. with_logger (logger) do
4181 # Overwrite saveat, always use dense output.
4282 # For stateless models (no unknowns) the adaptive solver takes no
4383 # internal steps and sol.t would be empty with saveat=[].
4484 # Supply explicit time points so observed variables can be evaluated.
45- sys = ode_prob. f. sys
46- saveat = isempty (ModelingToolkit. unknowns (sys)) ?
47- collect (range (ode_prob. tspan[1 ], ode_prob. tspan[end ]; length = 500 )) :
48- Float64[]
49- solve (ode_prob, Rodas5P (); saveat = saveat, dense = true )
85+ sys = ode_prob. f. sys
86+ n_unknowns = length (ModelingToolkit. unknowns (sys))
87+
88+ kwargs = if n_unknowns == 0
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 )
93+ else
94+ (saveat = Float64[], dense = true )
95+ end
96+
97+ # Log solver settings — init returns NullODEIntegrator (no .opts)
98+ # when the problem has no unknowns (u::Nothing), so only inspect
99+ # opts when a real integrator is returned.
100+ # Use our own `saveat` vector for the log: integ.opts.saveat is a
101+ # BinaryHeap which does not support iterate/minimum/maximum.
102+ integ = DifferentialEquations. init (ode_prob, solver; kwargs... )
103+ saveat_s = kwargs. saveat
104+ solver_settings_string = if hasproperty (integ, :opts )
105+ sv_str = isempty (saveat_s) ? " []" : " $(length (saveat_s)) points in [$(first (saveat_s)) , $(last (saveat_s)) ]"
106+ """
107+ Solver $(parentmodule (typeof (solver))) .$(nameof (typeof (solver)))
108+ saveat: $sv_str
109+ abstol: $(@sprintf (" %.2e" . opts. abstol))
110+ reltol: $(@sprintf (" %.2e" . opts. reltol))
111+ adaptive: $(integ. opts. adaptive)
112+ dense: $(integ. opts. dense)
113+ """
114+ else
115+ sv_str = isempty (saveat_s) ? " []" : " $(length (saveat_s)) points in [$(first (saveat_s)) , $(last (saveat_s)) ]"
116+ " Solver (NullODEIntegrator — no unknowns)
117+ saveat: $sv_str
118+ dense: true"
119+ end
120+
121+ # Solve
122+ DifferentialEquations. solve (ode_prob, solver; kwargs... )
50123 end
51124 sim_time = time () - t0
52- if sol. retcode == ReturnCode. Success
125+ if sol. retcode == DifferentialEquations . ReturnCode. Success
53126 sys = sol. prob. f. sys
54127 n_vars = length (ModelingToolkit. unknowns (sys))
55128 n_obs = length (ModelingToolkit. observed (sys))
@@ -67,7 +140,8 @@ function run_simulate(ode_prob, model_dir::String,
67140 sim_time = time () - t0
68141 sim_error = sprint (showerror, e, catch_backtrace ())
69142 end
70- println (log_file, " Time: $(round (sim_time; digits= 3 )) s"
143+ println (log_file, solver_settings_string)
144+ println (log_file, " Time: $(round (sim_time; digits= 3 )) s"
71145 println (log_file, " Success: $sim_success "
72146 isempty (sim_error) || println (log_file, " \n --- Error ---\n $sim_error "
73147 close (log_file)
0 commit comments