Switch solver from Rodas5P to Rodas5Pr and log resolved solver settings

.github/workflows/CI.yml

@@ -49,7 +49,13 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v3
-      - uses: julia-actions/julia-buildpkg@v1
+      - name: Build package
+        run: |
+          julia --project=. -e '
+            using Pkg
+            Pkg.add(name="BaseModelica", rev="main")
+            Pkg.instantiate()
+          '
       - uses: julia-actions/julia-runtest@v1
 
   sanity:
@@ -82,7 +88,13 @@ jobs:
           version: '1.12'
           arch: x64
       - uses: julia-actions/cache@v3
-      - uses: julia-actions/julia-buildpkg@v1
+      - name: Build package
+        run: |
+          julia --project=. -e '
+            using Pkg
+            Pkg.add(name="BaseModelica", rev="main")
+            Pkg.instantiate()
+          '
 
       - name: Sanity check ChuaCircuit
         run: |

Project.toml

@@ -7,6 +7,7 @@ authors = ["AnHeuermann"]
 BaseModelica = "a17d5099-185d-4ff5-b5d3-51aa4569e56d"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OMJulia = "0f4fe800-344e-11e9-2949-fb537ad918e1"

src/BaseModelicaLibraryTesting.jl

@@ -4,7 +4,7 @@ import Pkg
 import OMJulia
 import OMJulia: sendExpression
 import BaseModelica
-import DifferentialEquations: solve, Rodas5P, ReturnCode
+import DifferentialEquations
 import ModelingToolkit
 import Dates: now
 import Printf: @sprintf
@@ -21,11 +21,12 @@ include("pipeline.jl")
 # ── Public API ─────────────────────────────────────────────────────────────────
 
 # Shared types and constants
-export ModelResult, CompareSettings, RunInfo
+export ModelResult, CompareSettings, SimulateSettings, RunInfo
 export LIBRARY, LIBRARY_VERSION, CMP_REL_TOL, CMP_ABS_TOL
 
 # Comparison configuration
 export configure_comparison!, compare_settings
+export configure_simulate!, simulate_settings
 
 # Pipeline phases
 export run_export       # Phase 1: Base Modelica export via OMC

src/pipeline.jl

@@ -57,8 +57,12 @@ end
 
 Run the four-phase pipeline for a single model and return its result.
 """
-function test_model(omc::OMJulia.OMCSession, model::String, results_root::String,
-                    ref_root::String; csv_max_size_mb::Int = CSV_MAX_SIZE_MB)::ModelResult
+function test_model(omc::OMJulia.OMCSession,
+                    model::String,
+                    results_root::String,
+                    ref_root::String;
+                    sim_settings   ::SimulateSettings = _SIM_SETTINGS,
+                    csv_max_size_mb::Int              = CSV_MAX_SIZE_MB)::ModelResult
     model_dir = joinpath(results_root, "files", model)
     mkpath(model_dir)
 
@@ -93,6 +97,7 @@ function test_model(omc::OMJulia.OMCSession, model::String, results_root::String
 
     # Phase 3 ──────────────────────────────────────────────────────────────────
     sim_ok, sim_t, sim_err, sol = run_simulate(ode_prob, model_dir, model;
+                                               settings = sim_settings,
                                                csv_max_size_mb, cmp_signals)
 
     # Phase 4 (optional) ───────────────────────────────────────────────────────
@@ -132,6 +137,7 @@ function main(;
     results_root     :: String                = "",
     ref_root         :: String                = get(ENV, "MAPLIB_REF", ""),
     bm_options       :: String                = get(ENV, "BM_OPTIONS", "scalarize,moveBindings,inlineFunctions"),
+    sim_settings     :: SimulateSettings      = _SIM_SETTINGS,
     csv_max_size_mb  :: Int                   = CSV_MAX_SIZE_MB,
 )
     t0 = time()
@@ -200,7 +206,7 @@ function main(;
 
         for (i, model) in enumerate(models)
             @info "[$i/$(length(models))] $model"
-            result = test_model(omc, model, results_root, ref_root; csv_max_size_mb)
+            result = test_model(omc, model, results_root, ref_root; sim_settings, csv_max_size_mb)
             push!(results, result)
 
             phase = if result.sim_success && result.cmp_total > 0

src/simulate.jl

@@ -1,14 +1,49 @@
 # ── Phase 3: ODE simulation with DifferentialEquations / MTK ──────────────────
 
-import DifferentialEquations: solve, Rodas5P, ReturnCode
+import DifferentialEquations
 import Logging
 import ModelingToolkit
 import Printf: @sprintf
 
+"""Module-level default simulation settings.  Modify via `configure_simulate!`."""
+const _SIM_SETTINGS = SimulateSettings(solver = DifferentialEquations.Rodas5Pr())
+
+"""
+    configure_simulate!(; solver, saveat_n) → SimulateSettings
+
+Update the module-level simulation settings in-place and return them.
+
+# Keyword arguments
+- `solver`   — any SciML ODE/DAE algorithm instance (e.g. `Rodas5Pr`, `FBDF()`).
+- `saveat_n` — number of uniform time points for purely algebraic systems.
+
+# Example
+
+```julia
+using OrdinaryDiffEqBDF
+configure_simulate!(solver = FBDF())
+```
+"""
+function configure_simulate!(;
+    solver   :: Union{Any,Nothing} = nothing,
+    saveat_n :: Union{Int,Nothing} = nothing,
+)
+    isnothing(solver)   || (_SIM_SETTINGS.solver   = solver)
+    isnothing(saveat_n) || (_SIM_SETTINGS.saveat_n = saveat_n)
+    return _SIM_SETTINGS
+end
+
+"""
+    simulate_settings() → SimulateSettings
+
+Return the current module-level simulation settings.
 """
-    run_simulate(ode_prob, model_dir, model; cmp_signals, csv_max_size_mb) → (success, time, error, sol)
+simulate_settings() = _SIM_SETTINGS
 
-Solve `ode_prob` with Rodas5P (stiff solver).  On success, also writes the
+"""
+    run_simulate(ode_prob, model_dir, model; settings, cmp_signals, csv_max_size_mb) → (success, time, error, sol)
+
+Solve `ode_prob` using the algorithm in `settings.solver`.  On success, also writes the
 solution as a CSV file `<Short>_sim.csv` in `model_dir`.
 Writes a `<model>_sim.log` file in `model_dir`.
 Returns `nothing` as the fourth element on failure.
@@ -20,36 +55,74 @@ of signals will be compared.
 CSV files larger than `csv_max_size_mb` MiB are replaced with a
 `<Short>_sim.csv.toobig` marker so that the report can note the omission.
 """
-function run_simulate(ode_prob, model_dir::String,
+function run_simulate(ode_prob,
+                      model_dir::String,
                       model::String;
-                      cmp_signals    ::Vector{String} = String[],
-                      csv_max_size_mb::Int            = CSV_MAX_SIZE_MB)::Tuple{Bool,Float64,String,Any}
-    sim_success = false
-    sim_time    = 0.0
-    sim_error   = ""
-    sol         = nothing
+                      settings       ::SimulateSettings = _SIM_SETTINGS,
+                      cmp_signals    ::Vector{String}   = String[],
+                      csv_max_size_mb::Int              = CSV_MAX_SIZE_MB)::Tuple{Bool,Float64,String,Any}
+
+    sim_success           = false
+    sim_time              = 0.0
+    sim_error             = ""
+    sol                   = nothing
+    solver_settings_string = ""
 
     log_file = open(joinpath(model_dir, "$(model)_sim.log"), "w")
     println(log_file, "Model:   $model")
     logger = Logging.SimpleLogger(log_file, Logging.Debug)
     t0 = time()
+
+    solver = settings.solver
     try
-        # Rodas5P handles stiff DAE-like systems well.
         # Redirect all library log output (including Symbolics/MTK warnings)
         # to the log file so they don't clutter stdout.
         sol = Logging.with_logger(logger) do
             # Overwrite saveat, always use dense output.
             # For stateless models (no unknowns) the adaptive solver takes no
             # internal steps and sol.t would be empty with saveat=[].
             # Supply explicit time points so observed variables can be evaluated.
-            sys    = ode_prob.f.sys
-            saveat = isempty(ModelingToolkit.unknowns(sys)) ?
-                         collect(range(ode_prob.tspan[1], ode_prob.tspan[end]; length = 500)) :
-                         Float64[]
-            solve(ode_prob, Rodas5P(); saveat = saveat, dense = true)
+            sys        = ode_prob.f.sys
+            n_unknowns = length(ModelingToolkit.unknowns(sys))
+
+            kwargs = if n_unknowns == 0
+                # No unknowns at all (e.g. BusUsage):
+                # Supply explicit time points so observed variables can be evaluated.
+                saveat_s = collect(range(ode_prob.tspan[1], ode_prob.tspan[end]; length = settings.saveat_n))
+                (saveat = saveat_s, dense = true)
+            else
+                (saveat = Float64[], dense = true)
+            end
+
+            # Log solver settings — init returns NullODEIntegrator (no .opts)
+            # when the problem has no unknowns (u::Nothing), so only inspect
+            # opts when a real integrator is returned.
+            # Use our own `saveat` vector for the log: integ.opts.saveat is a
+            # BinaryHeap which does not support iterate/minimum/maximum.
+            integ = DifferentialEquations.init(ode_prob, solver; kwargs...)
+            saveat_s = kwargs.saveat
+            solver_settings_string = if hasproperty(integ, :opts)
+                sv_str = isempty(saveat_s) ? "[]" : "$(length(saveat_s)) points in [$(first(saveat_s)), $(last(saveat_s))]"
+                """
+                Solver $(parentmodule(typeof(solver))).$(nameof(typeof(solver)))
+                    saveat:   $sv_str
+                    abstol:   $(@sprintf("%.2e", integ.opts.abstol))
+                    reltol:   $(@sprintf("%.2e", integ.opts.reltol))
+                    adaptive: $(integ.opts.adaptive)
+                    dense:    $(integ.opts.dense)
+                """
+            else
+                sv_str = isempty(saveat_s) ? "[]" : "$(length(saveat_s)) points in [$(first(saveat_s)), $(last(saveat_s))]"
+                "Solver (NullODEIntegrator — no unknowns)
+                    saveat: $sv_str
+                    dense:  true"
+            end
+
+            # Solve
+            DifferentialEquations.solve(ode_prob, solver; kwargs...)
         end
         sim_time = time() - t0
-        if sol.retcode == ReturnCode.Success
+        if sol.retcode == DifferentialEquations.ReturnCode.Success
             sys    = sol.prob.f.sys
             n_vars = length(ModelingToolkit.unknowns(sys))
             n_obs  = length(ModelingToolkit.observed(sys))
@@ -67,7 +140,8 @@ function run_simulate(ode_prob, model_dir::String,
         sim_time  = time() - t0
         sim_error = sprint(showerror, e, catch_backtrace())
     end
-    println(log_file, "Time:    $(round(sim_time; digits=3)) s")
+    println(log_file, solver_settings_string)
+    println(log_file, "Time: $(round(sim_time; digits=3)) s")
     println(log_file, "Success: $sim_success")
     isempty(sim_error) || println(log_file, "\n--- Error ---\n$sim_error")
     close(log_file)

src/types.jl

@@ -37,6 +37,25 @@ Base.@kwdef mutable struct CompareSettings
     error_fn :: Symbol  = :mixed
 end
 
+# ── Simulation settings ────────────────────────────────────────────────────────
+
+"""
+    SimulateSettings
+
+Mutable configuration struct for ODE simulation.
+
+# Fields
+- `solver`   — any SciML ODE/DAE algorithm instance.  Default: `nothing`,
+               resolved to `Rodas5Pr()` when the module-level singleton is
+               constructed in `simulate.jl`.
+- `saveat_n` — number of evenly-spaced time points used for purely algebraic
+               systems (all mass-matrix rows zero).  Default: `500`.
+"""
+Base.@kwdef mutable struct SimulateSettings
+    solver   :: Any = nothing
+    saveat_n :: Int = 500
+end
+
 # ── Run metadata ───────────────────────────────────────────────────────────────
 
 """

test/chua_circuit.jl

@@ -0,0 +1,31 @@
+@testset "ChuaCircuit pipeline" begin
+    tmpdir    = mktempdir()
+    model_dir = joinpath(tmpdir, "files", TEST_MODEL_CHUA)
+    mkpath(model_dir)
+    bm_path = replace(abspath(joinpath(model_dir, "$TEST_MODEL_CHUA.bmo")), "\\" => "/")
+
+    omc = OMJulia.OMCSession(TEST_OMC)
+    try
+        OMJulia.sendExpression(omc, """setCommandLineOptions("--baseModelica --baseModelicaOptions=scalarize,moveBindings -d=evaluateAllParameters")""")
+        ok = OMJulia.sendExpression(omc, """loadModel(Modelica, {"4.1.0"})""")
+        @test ok == true
+
+        exp_ok, _, exp_err = run_export(omc, TEST_MODEL_CHUA, model_dir, bm_path)
+        @test exp_ok
+        exp_ok || @warn "Export error: $exp_err"
+
+        if exp_ok
+            par_ok, _, par_err, ode_prob = run_parse(bm_path, model_dir, TEST_MODEL_CHUA)
+            @test par_ok
+            par_ok || @warn "Parse error: $par_err"
+
+            if par_ok
+                sim_ok, _, sim_err, _ = run_simulate(ode_prob, model_dir, TEST_MODEL_CHUA)
+                @test sim_ok
+                sim_ok || @warn "Simulation error: $sim_err"
+            end
+        end
+    finally
+        OMJulia.quit(omc)
+    end
+end