From fc3db87309d483a203e0a5ea6e0b610342803389 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 12 May 2026 15:59:44 +0530 Subject: [PATCH 01/24] refactor: use `IRInfo` and `IRStructure` for `full_equations` --- .../src/systems/abstractsystem.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl index 31d5deec18..00d57653f9 100644 --- a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl +++ b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl @@ -1998,6 +1998,20 @@ equations during `mtkcompile`. These equations matches generated numerical code. See also [`equations`](@ref) and [`ModelingToolkitBase.get_eqs`](@ref). """ function full_equations(sys::AbstractSystem; simplify = false) + subsys = get_systems(sys) + # Fast path using `IRInfo` + if isempty(subsys) + new_eqs = Equation[] + eqs = equations(sys) + sizehint!(new_eqs, length(eqs)) + info = get_ir_info(sys) + ir = get_irstructure(sys) + @assert length(info.eqs_idxs) == length(eqs) + for (eq, rhs_idx) in zip(eqs, info.eqs_idxs) + push!(new_eqs, eq.lhs ~ ir[rhs_idx]) + end + return new_eqs + end empty_substitutions(sys) && return equations(sys) subs = get_substitutions(sys) neweqs = map(equations(sys)) do eq From 338e20f595df2481badda96a484108be270139b8 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 13 May 2026 15:55:12 +0530 Subject: [PATCH 02/24] fix: do not mutate decomposition subsystems in `SCCNonlinearProblem` --- src/problems/sccnonlinearproblem.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl index 4c16426d72..15aebcdcb1 100644 --- a/src/problems/sccnonlinearproblem.jl +++ b/src/problems/sccnonlinearproblem.jl @@ -296,7 +296,10 @@ function build_caches!(sys::System, decomposition::SCCDecomposition) push!(banned_vars, split_indexed_var(u)[1]) end - _eqs = get_eqs(subsys) + # While we own the system and so mutation _should_ be safe, `IRInfo` exists. + # It stores the indices corresponding to equations in the `IRStructure`, which + # would be incorrect if we mutate. + _eqs = copy(get_eqs(subsys)) exprs_to_search = SymbolicT[] for i in eachindex(_eqs) push!(exprs_to_search, _eqs[i].rhs) @@ -306,6 +309,7 @@ function build_caches!(sys::System, decomposition::SCCDecomposition) for i in eachindex(_eqs) _eqs[i] = _eqs[i].lhs ~ subber(_eqs[i].rhs) end + subsys = decomposition.subsystems[i] = ConstructionBase.setproperties(subsys; eqs = _eqs) if decomposition.islinear[i] store_to_mutable_cache!(subsys, CachedLinearAb, nothing) From dfb8ad17ed02eac210744075c504186ef89c91c0 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 15 May 2026 12:01:51 +0530 Subject: [PATCH 03/24] test: minor test update to account for new `full_equations` --- lib/ModelingToolkitBase/test/serialization.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/lib/ModelingToolkitBase/test/serialization.jl b/lib/ModelingToolkitBase/test/serialization.jl index 9bffd12f96..16814a6830 100644 --- a/lib/ModelingToolkitBase/test/serialization.jl +++ b/lib/ModelingToolkitBase/test/serialization.jl @@ -29,8 +29,7 @@ str = String(take!(io)) sys = include_string(@__MODULE__, str) rc2 = expand_connections(rc_model) -@test isapprox(sys, rc2) -@test issetequal(equations(sys), equations(rc2)) +@test issetequal(full_equations(sys), full_equations(rc2)) @test issetequal(unknowns(sys), unknowns(rc2)) @test issetequal(parameters(sys), parameters(rc2)) @@ -54,5 +53,4 @@ sol_ = solve(prob_, ImplicitEuler()) probexpr = ODEProblem{true}(ss, [capacitor.v => 0.0], (0, 0.1); expr = Val{true}, missing_guess_value); prob_obs = eval(probexpr) sol_obs = solve(prob_obs, ImplicitEuler()) -@show all_obs @test sol_obs[all_obs] == sol[all_obs] From 751b59354c8b769863b4f25db1a95b96f6a86464 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Mon, 18 May 2026 21:10:13 +0200 Subject: [PATCH 04/24] avoid calling `unhack_system` twice for the same `system` This shaves off around 33% of runtime in `mtkcompile` for the model I am testing it on: Before: ``` 11.193316 seconds (136.79 M allocations: 6.120 GiB, 6.71% gc time) ``` After: ``` 7.570217 seconds (85.88 M allocations: 3.767 GiB, 4.58% gc time) ``` --- lib/ModelingToolkitBase/src/systems/abstractsystem.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl index 31d5deec18..729dc98f81 100644 --- a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl +++ b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl @@ -556,7 +556,7 @@ end supports_initialization(sys::AbstractSystem) = true -function add_initialization_parameters(sys::AbstractSystem; split = true) +function add_initialization_parameters(sys::AbstractSystem; split = true, _unhack_sys = nothing) @assert !has_systems(sys) || isempty(get_systems(sys)) supports_initialization(sys) || return sys is_initializesystem(sys) && return sys @@ -564,8 +564,7 @@ function add_initialization_parameters(sys::AbstractSystem; split = true) all_initialvars = Set{SymbolicT}() # time-independent systems don't initialize unknowns # but may initialize parameters using guesses for unknowns - eqs = equations(sys) - _sys = unhack_system(sys) + _sys = _unhack_sys === nothing ? unhack_system(sys) : _unhack_sys obs = observed(_sys) eqs = equations(_sys) for x in unknowns(_sys) @@ -691,10 +690,10 @@ function complete( end sys = newsys check_no_parameter_equations(sys) + _unhack_sys = unhack_system(sys) if add_initial_parameters - sys = add_initialization_parameters(sys; split)::T + sys = add_initialization_parameters(sys; split, _unhack_sys)::T end - _unhack_sys = unhack_system(sys) cb_alg_eqs = Equation[alg_equations(_unhack_sys); observed(_unhack_sys)] if has_continuous_events(sys) && is_time_dependent(sys) cevts = SymbolicContinuousCallback[] From 53dcd46f9e97011dd54b596c052691a80a8eb09c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 19 May 2026 18:16:58 +0530 Subject: [PATCH 05/24] refactor: add `wrap_as_any` flag to `concrete_getu` --- lib/ModelingToolkitBase/src/systems/problem_utils.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/ModelingToolkitBase/src/systems/problem_utils.jl b/lib/ModelingToolkitBase/src/systems/problem_utils.jl index d6c9400fd0..9e7f99a910 100644 --- a/lib/ModelingToolkitBase/src/systems/problem_utils.jl +++ b/lib/ModelingToolkitBase/src/systems/problem_utils.jl @@ -672,7 +672,7 @@ Note that the getter ONLY works for problem-like objects, since it generates an function. It does NOT work for solutions. """ Base.@nospecializeinfer function concrete_getu( - indp, syms; + indp, syms; wrap_as_any = false, eval_expression, eval_module, force_time_independent = false, kwargs... ) @nospecialize @@ -680,6 +680,9 @@ Base.@nospecializeinfer function concrete_getu( indp, syms; wrap_delays = false, eval_expression, eval_module, force_time_independent, kwargs... ) + if wrap_as_any + return ObservedWrapper{is_time_dependent(indp) && !force_time_independent, Any}(obsfn) + end return ObservedWrapper{is_time_dependent(indp) && !force_time_independent}(obsfn) end From eeb66ca41893919753823e37265399b2888584ca Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Tue, 19 May 2026 18:26:54 +0200 Subject: [PATCH 06/24] fix spelling typos identified by CI --- .typos.toml | 4 + docs/src/basics/MTKLanguage.md | 2 +- docs/src/internals/systems.md | 4 +- docs/src/tutorials/initialization.md | 2 +- .../src/systems/nonlinear/initializesystem.jl | 2 +- .../test/initial_values.jl | 2 +- src/systems/alias_elimination.jl | 2 +- src/systems/codegen.jl | 6 +- test_compiler_options_fallback.jl | 97 +++++++++++++++++++ 9 files changed, 111 insertions(+), 10 deletions(-) create mode 100644 test_compiler_options_fallback.jl diff --git a/.typos.toml b/.typos.toml index a933260fb5..95a9e01066 100644 --- a/.typos.toml +++ b/.typos.toml @@ -7,3 +7,7 @@ ser = "ser" isconnection = "isconnection" Ue = "Ue" Derivate = "Derivate" +missings = "missings" +alph = "alph" +fo = "fo" +shs = "shs" diff --git a/docs/src/basics/MTKLanguage.md b/docs/src/basics/MTKLanguage.md index d075330e6c..2089047ec7 100644 --- a/docs/src/basics/MTKLanguage.md +++ b/docs/src/basics/MTKLanguage.md @@ -254,7 +254,7 @@ end - Defining discrete events as described [here](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/#Discrete-events-support). - If this block is not defined in the model, no discrete events will be added. - The block can also handle `if`-statements (evaluated at compilation) and unpacked lists defined by [comprehension](https://docs.julialang.org/en/v1/manual/arrays/#man-comprehensions) - - Dicrete events have the form `*condition* => *affect*` where *condition* is a boolean expression + - Discrete events have the form `*condition* => *affect*` where *condition* is a boolean expression - Discrete parameters and other keyword arguments should be specified in a vector, as seen below. ```@example mtkmodel-example diff --git a/docs/src/internals/systems.md b/docs/src/internals/systems.md index 130fc55dd5..5f5455415e 100644 --- a/docs/src/internals/systems.md +++ b/docs/src/internals/systems.md @@ -88,7 +88,7 @@ This constructor is responsible for: - Populating `var_to_name`. - Validating bindings. - Parsing bindings/initial conditions/guesses from variable metadata. -- Converting symbolic event shorthand into `SymbolicContinousCallback`/`SymbolicDiscreteCallback`. +- Converting symbolic event shorthand into `SymbolicContinuousCallback`/`SymbolicDiscreteCallback`. !!! note "Parsing metadata" It is often useful to disable parsing bindings/initial conditions/guesses from variable metadata during @@ -200,7 +200,7 @@ The order and number of unknowns should correspond to that of equations. For exa D(x) ~ rhs ``` -Then `x` should be the `i`th unknown. Algebraic variables often have a weaker correspondance. +Then `x` should be the `i`th unknown. Algebraic variables often have a weaker correspondence. ## `AtomicArrayDict` diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index 8e3e4d308c..beaea968d4 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -99,7 +99,7 @@ allowed values on the RHS along with its behavior during initialization. For thi nomenclature "FP parameter" or "FP variable" is used to refer to parameters/variables which take continuous values, determined by their `SymbolicUtils.symtype` being `Real`, or a subtype of `AbstractFloat`, or an array thereof. The term "variable" in this table is used -to refer to both continously-varying quantities such as the unknowns of a DAE (created +to refer to both continuously-varying quantities such as the unknowns of a DAE (created via `@variables`) as well as discretely-varying quantities updated in events (created by `@discretes`). Note that all FP variables are solved for by initialization. Any quantity that initialization solves for is referred to as an "initialization unknown". diff --git a/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl b/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl index 2c86663f63..9f2e13ca40 100644 --- a/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl +++ b/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl @@ -438,7 +438,7 @@ function timevaring_initsys_process_op!( is a variable/parameter of the system, or change the set of initial conditions \ provided to the problem constructor. 2. An initial condition was provided for a bound parameter. Bound parameters cannot \ - be given intial conditions. Their values are determined solely by the binding \ + be given initial conditions. Their values are determined solely by the binding \ relation. In case neither of these is the case, please open an issue in ModelingToolkit.jl \ diff --git a/lib/ModelingToolkitBase/test/initial_values.jl b/lib/ModelingToolkitBase/test/initial_values.jl index dfbdfaf820..164c1b52bf 100644 --- a/lib/ModelingToolkitBase/test/initial_values.jl +++ b/lib/ModelingToolkitBase/test/initial_values.jl @@ -455,7 +455,7 @@ if !@isdefined(ModelingToolkit) @parameters p d eqs = [D(X) ~ p - d * X] @mtkcompile sys = System(eqs, t) - sim_cond = [X => 1.0, X => 2.0, p => 1.0, d => 2.0] # `X` provide twice, likely a misstake + sim_cond = [X => 1.0, X => 2.0, p => 1.0, d => 2.0] # `X` provide twice, likely a mistake @test_throws ["duplicate entries", "X(t)"] ODEProblem(sys, sim_cond, (0.0, 10.0)) end end diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl index b5e668febb..82123d5201 100644 --- a/src/systems/alias_elimination.jl +++ b/src/systems/alias_elimination.jl @@ -253,7 +253,7 @@ function alias_elimination!(state::TearingState; fully_determined = true, (; fullvars, sys) = state (; graph, var_to_diff, solvable_graph) = state.structure # Previously, underconstrained variables were zeroed out. This leads to significant - # unintuitive behavior, especially for intialization systems. The underconstrained variables + # unintuitive behavior, especially for initialization systems. The underconstrained variables # should constitute an underdetermined error, and hence the behavior is removed. This pass # continues to remove redundant equations, since it is essential for adding analysis points # to existing connections. diff --git a/src/systems/codegen.jl b/src/systems/codegen.jl index 1e74858be8..f69d959d6a 100644 --- a/src/systems/codegen.jl +++ b/src/systems/codegen.jl @@ -191,7 +191,7 @@ end Generate `f1` and `f2` for [`SemilinearODEFunction`](@ref) (internally represented as a `SplitFunction`). `A`, `B`, `C` are the matrices returned from [`calculate_semiquadratic_form`](@ref). This expects that the system has the necessary -extra parmameters added by [`add_semiquadratic_parameters`](@ref). +extra parameters added by [`add_semiquadratic_parameters`](@ref). ## Keyword Arguments @@ -367,7 +367,7 @@ Generate the jacobian of `f1` for [`SemilinearODEFunction`](@ref) (internally re `SplitFunction`). `A`, `B`, `C` are the matrices returned from [`calculate_semiquadratic_form`](@ref). `Cjac` is the jacobian of `C` with respect to the unknowns of the system, or `nothing` if `C === nothing`. This expects that the system has the -necessary extra parmameters added by [`add_semiquadratic_parameters`](@ref). +necessary extra parameters added by [`add_semiquadratic_parameters`](@ref). ## Keyword Arguments @@ -530,7 +530,7 @@ Return the sparsity pattern of the jacobian of `f1` for [`SemilinearODEFunction (internally represented as a `SplitFunction`). `A`, `B`, `C` are the matrices returned from [`calculate_semiquadratic_form`](@ref). `Cjac` is the jacobian of `C` with respect to the unknowns of the system, or `nothing` if `C === nothing`. This expects that the system has the -necessary extra parmameters added by [`add_semiquadratic_parameters`](@ref). +necessary extra parameters added by [`add_semiquadratic_parameters`](@ref). ## Keyword Arguments diff --git a/test_compiler_options_fallback.jl b/test_compiler_options_fallback.jl new file mode 100644 index 0000000000..ce65137e3d --- /dev/null +++ b/test_compiler_options_fallback.jl @@ -0,0 +1,97 @@ +# Test that module-level compiler options are inherited by functions eval'd in those modules, +# including RuntimeGeneratedFunctions. +# +# Run with: julia --project test_compiler_options_fallback.jl + +using RuntimeGeneratedFunctions +using InteractiveUtils: code_native, code_typed + +# Create modules with different compiler options +module EvalOpt0 + Base.Experimental.@compiler_options optimize=0 + using RuntimeGeneratedFunctions + RuntimeGeneratedFunctions.init(@__MODULE__) +end + +module EvalOpt1 + Base.Experimental.@compiler_options optimize=1 + using RuntimeGeneratedFunctions + RuntimeGeneratedFunctions.init(@__MODULE__) +end + +module EvalOpt0NoInfer + Base.Experimental.@compiler_options optimize=0 infer=false + using RuntimeGeneratedFunctions + RuntimeGeneratedFunctions.init(@__MODULE__) +end + +module EvalDefault + using RuntimeGeneratedFunctions + RuntimeGeneratedFunctions.init(@__MODULE__) +end + +# A non-trivial expression to compile +test_expr = :(function (x) + y = sin(x) + cos(x) + return y * y + exp(x) +end) + +# Test eval path +println("=== Testing eval path ===") +f_default = EvalDefault.eval(test_expr) +f_opt0 = EvalOpt0.eval(test_expr) +f_opt1 = EvalOpt1.eval(test_expr) +f_opt0_noinfer = EvalOpt0NoInfer.eval(test_expr) + +x = 1.5 +println("Default: f($x) = $(f_default(x))") +println("Opt0: f($x) = $(f_opt0(x))") +println("Opt1: f($x) = $(f_opt1(x))") +println("Opt0+no infer: f($x) = $(f_opt0_noinfer(x))") +println("Results match: ", f_default(x) == f_opt0(x) == f_opt1(x) == f_opt0_noinfer(x)) + +# Test RGF path +println("\n=== Testing RuntimeGeneratedFunction path ===") +rgf_default = RuntimeGeneratedFunction(EvalDefault, EvalDefault, test_expr) +rgf_opt0 = RuntimeGeneratedFunction(EvalOpt0, EvalOpt0, test_expr) +rgf_opt1 = RuntimeGeneratedFunction(EvalOpt1, EvalOpt1, test_expr) +rgf_opt0_noinfer = RuntimeGeneratedFunction(EvalOpt0NoInfer, EvalOpt0NoInfer, test_expr) + +println("RGF Default: f($x) = $(rgf_default(x))") +println("RGF Opt0: f($x) = $(rgf_opt0(x))") +println("RGF Opt1: f($x) = $(rgf_opt1(x))") +println("RGF Opt0+no infer: f($x) = $(rgf_opt0_noinfer(x))") +println("Results match: ", rgf_default(x) == rgf_opt0(x) == rgf_opt1(x) == rgf_opt0_noinfer(x)) + +# Check that the compiler options actually take effect by inspecting code_native. +# At opt level 0, the native code should contain more instructions (no inlining/folding). +println("\n=== Checking native code differences ===") +native_default = sprint(code_native, rgf_default, Tuple{Float64}) +native_opt0 = sprint(code_native, rgf_opt0, Tuple{Float64}) +native_opt1 = sprint(code_native, rgf_opt1, Tuple{Float64}) + +# Count call instructions — opt0 should have more calls (nothing inlined) +count_calls(s) = count(r"callq?\s"i, s) +println("Call instructions - Default: $(count_calls(native_default)), Opt0: $(count_calls(native_opt0)), Opt1: $(count_calls(native_opt1))") +println("Raw instruction count - Default: $(count('\n', native_default)), Opt0: $(count('\n', native_opt0)), Opt1: $(count('\n', native_opt1))") +if count_calls(native_opt0) > count_calls(native_opt1) + println("Opt0 has more calls than Opt1 -> module-level options ARE taking effect") +elseif count_calls(native_opt0) == count_calls(native_opt1) + println("Same call count — check raw instruction counts above for differences") +else + println("WARNING: Opt0 has fewer calls than Opt1 — unexpected") +end + +# Check inference works/doesn't work +println("\n=== Checking infer=false effect ===") +ci_default = only(code_typed(rgf_default, Tuple{Float64})) +ci_noinfer = only(code_typed(rgf_opt0_noinfer, Tuple{Float64})) +println("Default return type: $(ci_default[2])") +println("Opt0+no infer return type: $(ci_noinfer[2])") +if ci_noinfer[2] === Any + println("infer=false IS taking effect (return type is Any)") +else + println("infer=false may NOT be taking effect (return type is concrete)") +end + +println("\nDone.") From bbd1c3b22b9dbb189fc4896ffdc05bd12c03e943 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 19 May 2026 18:16:19 +0530 Subject: [PATCH 07/24] refactor: improve `isinitial` --- .../src/systems/abstractsystem.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl index 491cf1e633..aa7e6928f5 100644 --- a/lib/ModelingToolkitBase/src/systems/abstractsystem.jl +++ b/lib/ModelingToolkitBase/src/systems/abstractsystem.jl @@ -626,12 +626,17 @@ end """ Returns true if the parameter `p` is of the form `Initial(x)`. """ +function isinitial(p::SymbolicT) + p, _ = split_indexed_var(p) + Moshi.Match.@match p begin + BSImpl.Term(; f) => f isa Initial + _ => false + end +end function isinitial(p) - p = unwrap(p) - return iscall(p) && ( - operation(p) isa Initial || - operation(p) === getindex && isinitial(arguments(p)[1]) - ) + up = unwrap(p) + up === p && return false + return isinitial(up) end """ From d060e81440653e421642d854e12dc7bc929cfc89 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 15 May 2026 19:14:43 +0530 Subject: [PATCH 08/24] fix: avoid generating massive functions in `get_mtkparameters_reconstructor` This used to be the lion's share of the compile time in large models. The infrastructure added here might be useful in other places too. --- .../src/systems/problem_utils.jl | 299 ++++++++++++++---- 1 file changed, 241 insertions(+), 58 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/problem_utils.jl b/lib/ModelingToolkitBase/src/systems/problem_utils.jl index 9e7f99a910..b05cbea2dd 100644 --- a/lib/ModelingToolkitBase/src/systems/problem_utils.jl +++ b/lib/ModelingToolkitBase/src/systems/problem_utils.jl @@ -721,6 +721,223 @@ function (pca::PConstructorApplicator)(x::AbstractArray{<:AbstractArray}) return pca.p_constructor(pca.(x)) end +""" + $TYPEDEF + +Callable struct designed for use by `MTKParametersReconstructor`. Uses a fixed set of templates to +act as a very dynamic (and limited) observed function returning an array. See `__apply_copy_template` +for the supported templates. +""" +struct CopyParamsByTemplate{IsRoot, T, N} + """ + List of templates. + """ + template::T # TODO: This field is parametric because I thought we might want to specialize it in some cases. + """ + Size of the returned buffer. + """ + size::NTuple{N, Int} +end + +function CopyParamsByTemplate{IR}(temp::T, size::NTuple{N, Int}) where {IR, T, N} + return CopyParamsByTemplate{IR, T, N}(temp, size) +end + +function __apply_copy_template(valp, template) + p = parameter_values(valp) + u = state_values(valp) + if template isa ParameterIndex{SciMLStructures.Tunable, UnitRange{Int}} + if p isa MTKParameters + return p.tunable[template.idx] + else + return p[template.idx] + end + elseif template isa ParameterIndex{SciMLStructures.Initials, UnitRange{Int}} + return p.initials[template.idx] + elseif template isa ParameterIndex{SciMLStructures.Discrete, Tuple{Int, UnitRange{Int}}} + return p.discrete[template.idx[1]][template.idx[2]] + elseif template isa ParameterIndex{SciMLStructures.Constants, Tuple{Int, UnitRange{Int}}} + return p.constant[template.idx[1]][template.idx[2]] + elseif template isa ParameterIndex{Nonnumeric, Tuple{Int, UnitRange{Int}}} + return p.nonnumeric[template.idx[1]][template.idx[2]] + elseif template isa UnitRange{Int} + return u[template] + elseif template isa ObservedWrapper + return template(valp) + elseif template isa CopyParamsByTemplate + return template(valp) + elseif template isa IndepVarTemplate + return current_time(valp) + else + # MethodError because this is a manual dispatch chain + throw(MethodError(__apply_copy_template, (valp, template))) + end +end + +function (cp::CopyParamsByTemplate{IsRoot})(src) where {IsRoot} + if IsRoot + reshape(mapreduce(Base.Fix1(__apply_copy_template, src), vcat, cp.template), cp.size) + else + buffers = map(Base.Fix1(__apply_copy_template, src), cp.template) + if cp.template isa Tuple + buffers = collect(buffers) + end + reshape(buffers, cp.size) + end +end + +struct IndepVarTemplate end +const IV_TEMPLATE = IndepVarTemplate() + +Base.@nospecializeinfer function __specialize_templates(template::Vector{Any}, elem_types::Set{DataType}) + if length(template) <= 4 + return Tuple(template) + elseif length(elem_types) <= 4 + return Vector{Union{collect(elem_types)...}}(template) + else + return template + end +end + +function CopyParamsByTemplate(srcsys::AbstractSystem, syms::AbstractArray{SymbolicT}; kws...) + template = [] + elem_types = Set{DataType}() + iv = get_iv(srcsys) + for sym in syms + if iv isa SymbolicT && isequal(iv, sym) + push!(template, IV_TEMPLATE) + push!(elem_types, IndepVarTemplate) + continue + end + symidx = parameter_index(srcsys, sym) + if symidx === nothing + symidx = variable_index(srcsys, sym) + if symidx === nothing + if isempty(template) + push!(elem_types, Vector{SymbolicT}) + push!(template, SymbolicT[sym]) + continue + end + prev = template[end] + if prev isa Vector{SymbolicT} + push!(prev, sym) + else + push!(elem_types, Vector{SymbolicT}) + push!(template, SymbolicT[sym]) + end + continue + end + if isempty(template) + push!(elem_types, UnitRange{Int}) + push!(template, symidx:symidx) + continue + end + prev = template[end] + if prev isa UnitRange{Int} && last(prev) + 1 == symidx + template[end] = first(prev):symidx + else + push!(elem_types, UnitRange{Int}) + push!(template, symidx:symidx) + end + continue + elseif symidx isa Int + symidx = ParameterIndex(SciMLStructures.Tunable(), symidx) + end + portion = symidx.portion + _bufidx = symidx.idx + bufidx::UnitRange{Int} = if _bufidx isa AbstractVector{Int} + @assert isequal(vec(_bufidx), first(_bufidx):last(_bufidx)) + subidx = nothing + first(_bufidx):last(_bufidx) + elseif _bufidx isa Int + subidx = nothing + _bufidx:_bufidx + elseif _bufidx isa NTuple{2, Int} + subidx = _bufidx[1] + _bufidx[2]:_bufidx[2] + else + # Will error due to the typeassert on `bufidx` + nothing + end + if isempty(template) + pidx = if subidx === nothing + ParameterIndex(symidx.portion, bufidx) + else + ParameterIndex(symidx.portion, (subidx, bufidx)) + end + push!(template, pidx) + push!(elem_types, typeof(pidx)) + continue + end + prev = template[end] + if prev isa ParameterIndex && prev.portion === symidx.portion && (subidx === nothing && last(prev.idx) + 1 == first(bufidx) || subidx == prev.idx[1] && last(prev.idx[2]) + 1 == first(bufidx)) + if subidx === nothing + template[end] = ParameterIndex(prev.portion, first(prev.idx):last(bufidx)) + else + template[end] = ParameterIndex(prev.portion, (subidx, first(prev.idx[2]):last(bufidx))) + end + elseif subidx === nothing + push!(template, ParameterIndex(symidx.portion, bufidx)) + push!(elem_types, typeof(template[end])) + else + push!(template, ParameterIndex(symidx.portion, (subidx, bufidx))) + push!(elem_types, typeof(template[end])) + end + end + + for i in eachindex(template) + if template[i] isa Vector{SymbolicT} + template[i] = concrete_getu(srcsys, template[i]; wrap_as_any = true, kws...) + delete!(elem_types, Vector{SymbolicT}) + push!(elem_types, typeof(template[i])) + end + end + + return CopyParamsByTemplate{true}(__specialize_templates(template, elem_types), size(syms)) +end + +function CopyParamsByTemplate(srcsys::AbstractSystem, syms::AbstractArray; kws...) + template = [] + elem_types = Set{DataType}() + for sym in syms + push!(template, CopyParamsByTemplate(srcsys, sym; kws...)) + push!(elem_types, typeof(template[end])) + end + return CopyParamsByTemplate{false}(__specialize_templates(template, elem_types), size(syms)) +end + +struct MTKParametersReconstructor{T, I, D, C, N} + tunables_fn::T + initials_fn::I + discretes_fn::D + consts_fn::C + nonnumerics_fn::N + diffcache_buffer_idx::Int +end + +# TODO: make this infer when the nonnumerics are non-trivial +function (recon::MTKParametersReconstructor)(src, dst) + src_ps = parameter_values(src) + dst_ps = parameter_values(dst) + oldcache = dst_ps.caches + # I don't know why but this makes it infer properly + if recon.tunables_fn isa ComposedFunction + tunablevals = recon.tunables_fn.outer(recon.tunables_fn.inner(src)) + else + tunablevals = recon.tunables_fn(src) + end + initialvals = recon.initials_fn(src) + nonnumerics = recon.nonnumerics_fn(src)::typeof(dst_ps.nonnumeric) + (; diffcache_buffer_idx) = recon + if !iszero(diffcache_buffer_idx) + @set! nonnumerics[diffcache_buffer_idx] = DiffCacheAllocatorAPIWrapper{ForwardDiff.valtype(eltype(initialvals))}.(nonnumerics[diffcache_buffer_idx]) + end + return MTKParameters( + tunablevals, initialvals, recon.discretes_fn(src), + recon.consts_fn(src), nonnumerics, oldcache isa Tuple{} ? () : copy.(oldcache) + ) +end + """ $(TYPEDSIGNATURES) @@ -735,10 +952,10 @@ takes a value provider of `srcsys` and a value provider of `dstsys` and returns unwrapped. - `p_constructor`: The `p_constructor` argument to `process_SciMLProblem`. """ -function get_mtkparameters_reconstructor( +function MTKParametersReconstructor( srcsys::AbstractSystem, dstsys::AbstractSystem; initials = false, unwrap_initials = false, p_constructor = identity, - eval_expression = false, eval_module = @__MODULE__, force_time_independent = false, + force_time_independent = false, kwargs... ) _p_constructor = p_constructor @@ -757,32 +974,27 @@ function get_mtkparameters_reconstructor( tunable_getter = if isempty(tunable_syms) Returns(SVector{0, Float64}()) else - p_constructor ∘ concrete_getu( - srcsys, tunable_syms; eval_expression, eval_module, - force_time_independent, kwargs... - ) + p_constructor ∘ CopyParamsByTemplate(srcsys, tunable_syms; kwargs...) end initials_getter = if initials && !isempty(syms[2]) - initsyms = Vector{Any}(syms[2]) - allsyms = Set(variable_symbols(srcsys)) + initsyms = syms[2]::Vector{SymbolicT} + allsyms = Set{SymbolicT}(variable_symbols(srcsys)) if unwrap_initials for i in eachindex(initsyms) sym = initsyms[i] - innersym = if operation(sym) === getindex - sym, idxs... = arguments(sym) - only(arguments(sym))[idxs...] + arr, isarr = split_indexed_var(sym) + innersym = if isarr + sidx = get_stable_index(sym) + first(arguments(arr))[sidx] else - only(arguments(sym)) + first(arguments(arr)) end if innersym in allsyms initsyms[i] = innersym end end end - p_constructor ∘ concrete_getu( - srcsys, initsyms; eval_expression, eval_module, - force_time_independent, kwargs... - ) + p_constructor ∘ CopyParamsByTemplate(srcsys, initsyms; kwargs...) else Returns(SVector{0, Float64}()) end @@ -795,29 +1007,25 @@ function get_mtkparameters_reconstructor( p_constructor(map(x -> x.length, bufsizes)) end ) + # discretes need to be blocked arrays # the `getu` returns a tuple of arrays corresponding to `p.discretes` # `Base.Fix1(...)` applies `p_constructor` to each of the arrays in the tuple # `Base.Fix2(...)` does `BlockedArray.(tuple_of_arrs, blockarrsizes)` returning a # tuple of `BlockedArray`s Base.Fix2(Broadcast.BroadcastFunction(BlockedArray), blockarrsizes) ∘ - Base.Fix1(broadcast, p_constructor) ∘ + Base.Fix1(broadcast, p_constructor) ∘ Tuple ∘ # This `broadcast.(collect, ...)` avoids `ReshapedArray`/`SubArray`s from # appearing in the result. - concrete_getu( - srcsys, Tuple(broadcast.(collect, syms[3])); - eval_expression, eval_module, force_time_independent, kwargs... - ) + CopyParamsByTemplate(srcsys, broadcast.(collect, syms[3]); kwargs...) end - const_getter = if syms[4] == () + const_getter = if isempty(syms[4]) Returns(()) else - Base.Fix1(broadcast, p_constructor) ∘ concrete_getu( - srcsys, Tuple(syms[4]); - eval_expression, eval_module, force_time_independent, kwargs... - ) + Base.Fix1(broadcast, p_constructor) ∘ Tuple ∘ CopyParamsByTemplate(srcsys, syms[4]; kwargs...) end - nonnumeric_getter = if syms[5] == () + diffcache_buffer_idx = 0 + nonnumeric_getter = if isempty(syms[5]) Returns(()) else ic = get_index_cache(dstsys) @@ -828,44 +1036,19 @@ function get_mtkparameters_reconstructor( ) diffcache_params = SU.getmetadata(dstsys, DiffCacheParams, Dict{SymbolicT, Int}())::Dict{SymbolicT, Int} - diffcache_buffer_idx = 0 if !isempty(diffcache_params) representative = first(keys(diffcache_params)) diffcache_buffer_idx, _ = ic.nonnumeric_idx[representative] @set! buftypes[diffcache_buffer_idx] = identity + for (i, sym) in enumerate(syms[5][diffcache_buffer_idx]) + end end # nonnumerics retain the assigned buffer type without narrowing Base.Fix1(broadcast, _p_constructor) ∘ - Base.Fix1(Broadcast.BroadcastFunction(call), buftypes) ∘ - concrete_getu( - srcsys, Tuple(syms[5]); - eval_expression, eval_module, force_time_independent, kwargs... - ) - end - getters = ( - tunable_getter, initials_getter, discs_getter, const_getter, nonnumeric_getter, - ) - getter = let getters = getters, diffcache_buffer_idx = diffcache_buffer_idx - function _getter(valp, initprob) - oldcache = parameter_values(initprob).caches - tunablevals = getters[1](valp) - initialvals = getters[2](valp) - nonnumerics = getters[5](valp) - if !iszero(diffcache_buffer_idx) - @set! nonnumerics[diffcache_buffer_idx] = DiffCacheAllocatorAPIWrapper{ForwardDiff.valtype(eltype(initialvals))}.(nonnumerics[diffcache_buffer_idx]) - end - return promote_with_nothing( - promote_type_with_nothing(eltype(tunablevals), initialvals), - MTKParameters( - tunablevals, initialvals, getters[3](valp), - getters[4](valp), nonnumerics, oldcache isa Tuple{} ? () : - copy.(oldcache) - ) - ) - end + Base.Fix1(Broadcast.BroadcastFunction(call), buftypes) ∘ Tuple ∘ CopyParamsByTemplate(srcsys, syms[5]; kwargs...) end - return getter + return MTKParametersReconstructor(tunable_getter, initials_getter, discs_getter, const_getter, nonnumeric_getter, diffcache_buffer_idx) end function call(f, args...) @@ -893,7 +1076,7 @@ function ReconstructInitializeprob( kwargs... ) if is_split(dstsys) - pgetter = get_mtkparameters_reconstructor( + pgetter = MTKParametersReconstructor( srcsys, dstsys; p_constructor, eval_expression, eval_module, force_time_independent = is_steadystateprob, kwargs... ) @@ -973,7 +1156,7 @@ function construct_initializeprobpmap( ) @assert is_initializesystem(initsys) if is_split(sys) - return let getter = get_mtkparameters_reconstructor( + return let getter = MTKParametersReconstructor( initsys, sys; initials = true, unwrap_initials = true, p_constructor, eval_expression, eval_module, kwargs... ) From 53a488f6c40afcee27458b8fdf2c49f644013540 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 20 May 2026 00:03:27 +0530 Subject: [PATCH 09/24] refactor: improve `get_updated_symbolic_problem` Enables supporting ReverseDiff AD --- .../src/systems/nonlinear/initializesystem.jl | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl b/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl index 2c86663f63..292a3a6634 100644 --- a/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl +++ b/lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl @@ -812,14 +812,9 @@ function DiffEqBase.get_updated_symbolic_problem( end u0 = DiffEqBase.promote_u0(u0, buffer, t0) + u0 = ArrayInterface.restructure(u0, meta.get_updated_u0(prob, initdata.initializeprob)) - if ArrayInterface.ismutable(u0) - T = typeof(u0) - else - T = StaticArrays.similar_type(u0) - end - - return remake(prob; u0 = T(meta.get_updated_u0(prob, initdata.initializeprob)), p) + return remake(prob; u0, p) end """ From 1f026bdb459726c2e7405ed12b8a46f7fd7f18be Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 20 May 2026 00:02:47 +0530 Subject: [PATCH 10/24] feat: add MTKTrackerExt --- lib/ModelingToolkitBase/Project.toml | 2 ++ lib/ModelingToolkitBase/ext/MTKTrackerExt.jl | 14 ++++++++++++++ 2 files changed, 16 insertions(+) create mode 100644 lib/ModelingToolkitBase/ext/MTKTrackerExt.jl diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml index 5f6e9eb742..6781128148 100644 --- a/lib/ModelingToolkitBase/Project.toml +++ b/lib/ModelingToolkitBase/Project.toml @@ -73,6 +73,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" Pyomo = "0e8e1daf-01b5-4eba-a626-3897743a3816" +Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" [extensions] MTKBifurcationKitExt = "BifurcationKit" @@ -86,6 +87,7 @@ MTKLabelledArraysExt = "LabelledArrays" MTKLatexifyExt = "Latexify" MTKMooncakeExt = "Mooncake" MTKPyomoDynamicOptExt = "Pyomo" +MTKTrackerExt = "Tracker" [compat] ADTypes = "1.14.0" diff --git a/lib/ModelingToolkitBase/ext/MTKTrackerExt.jl b/lib/ModelingToolkitBase/ext/MTKTrackerExt.jl new file mode 100644 index 0000000000..fd97a518df --- /dev/null +++ b/lib/ModelingToolkitBase/ext/MTKTrackerExt.jl @@ -0,0 +1,14 @@ +module MTKTrackerExt + +import ModelingToolkitBase as MTKBase +import Tracker + +function MTKBase.promote_type_with_nothing(::Type{Tracker.TrackedReal{T}}, x::Tracker.TrackedArray{T}) where {T} + return Tracker.TrackedReal{T} +end + +function MTKBase.promote_with_nothing(::Type{Tracker.TrackedReal{T}}, x::Tracker.TrackedArray{T}) where {T} + return x +end + +end From 2845a607e1213573f45ccda1b5233a34be43b4c5 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Tue, 19 May 2026 18:27:59 +0200 Subject: [PATCH 11/24] fix trailing whitespace --- .github/workflows/Documentation.yml | 80 +++++++-------- .github/workflows/SpellCheck.yml | 2 +- LICENSE.md | 30 +++--- NEWS.md | 6 +- docs/src/API/System.md | 2 +- docs/src/API/model_building.md | 4 +- docs/src/API/variables.md | 4 +- docs/src/basics/Events.md | 2 +- docs/src/basics/InputOutput.md | 2 +- docs/src/basics/Linearization.md | 4 +- docs/src/basics/MTKLanguage.md | 6 +- docs/src/basics/PrecompileComponents.md | 2 +- docs/src/index.md | 32 +++--- docs/src/tutorials/SampledData.md | 6 +- docs/src/tutorials/attractors.md | 2 +- .../tutorials/change_independent_variable.md | 6 +- docs/src/tutorials/fmi.md | 6 +- docs/src/tutorials/initialization.md | 4 +- docs/src/tutorials/linear_analysis.md | 2 +- docs/src/tutorials/modelingtoolkitize.md | 4 +- docs/src/tutorials/nonlinear.md | 2 +- docs/src/tutorials/optimization.md | 2 +- docs/src/tutorials/stochastic_diffeq.md | 2 +- ext/MTKFMIExt.jl | 2 +- lib/ModelingToolkitBase/src/problems/docs.jl | 4 +- .../src/systems/callbacks.jl | 8 +- .../src/systems/codegen_utils.jl | 2 +- .../src/systems/imperative_affect.jl | 4 +- .../src/systems/problem_utils.jl | 2 +- lib/ModelingToolkitBase/src/variables.jl | 12 +-- .../test/input_output_handling.jl | 4 +- lib/ModelingToolkitBase/test/odesystem.jl | 2 +- lib/SciCompDSL/src/model_parsing.jl | 6 +- lib/SciCompDSL/test/model_parsing.jl | 8 +- src/ModelingToolkit.jl | 2 +- src/problems/sccnonlinearproblem.jl | 2 +- src/systems/alias_elimination.jl | 2 +- src/systems/analysis_points.jl | 2 +- src/systems/if_lifting.jl | 2 +- test/downstream/inversemodel.jl | 2 +- test_compiler_options_fallback.jl | 97 ------------------- 41 files changed, 139 insertions(+), 236 deletions(-) delete mode 100644 test_compiler_options_fallback.jl diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 73dde99d76..7d10317fa1 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -1,40 +1,40 @@ -name: Documentation - -on: - push: - branches: - - master - - v10 - - 'as/su-v4' - tags: '*' - pull_request: - -concurrency: - # Skip intermediate builds: always, but for the master branch. - # Cancel intermediate builds: always, but for the master branch. - group: ${{ github.workflow }}-${{ github.ref }} - cancel-in-progress: ${{ github.ref != 'refs/heads/master' && github.refs != 'refs/tags/*'}} - -jobs: - build: - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@latest - with: - version: 'lts' - - run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev - - name: Install dependencies - run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ -e 'using Pkg; Pkg.develop([PackageSpec(path=pwd()), PackageSpec(path=joinpath(pwd(), "lib", "ModelingToolkitBase")), PackageSpec(path=joinpath(pwd(), "lib", "SciCompDSL"))]); Pkg.instantiate()' - - name: Build and deploy - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key - JULIA_DEBUG: "Documenter" - run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ --code-coverage=user docs/make.jl - - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v6 - with: - files: lcov.info - token: ${{ secrets.CODECOV_TOKEN }} - fail_ci_if_error: true +name: Documentation + +on: + push: + branches: + - master + - v10 + - 'as/su-v4' + tags: '*' + pull_request: + +concurrency: + # Skip intermediate builds: always, but for the master branch. + # Cancel intermediate builds: always, but for the master branch. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ github.ref != 'refs/heads/master' && github.refs != 'refs/tags/*'}} + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@latest + with: + version: 'lts' + - run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev + - name: Install dependencies + run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ -e 'using Pkg; Pkg.develop([PackageSpec(path=pwd()), PackageSpec(path=joinpath(pwd(), "lib", "ModelingToolkitBase")), PackageSpec(path=joinpath(pwd(), "lib", "SciCompDSL"))]); Pkg.instantiate()' + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + JULIA_DEBUG: "Documenter" + run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ --code-coverage=user docs/make.jl + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v6 + with: + files: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + fail_ci_if_error: true diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 746b039a14..e3e08cd594 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v6 - name: Check spelling - uses: crate-ci/typos@v1.18.0 \ No newline at end of file + uses: crate-ci/typos@v1.18.0 \ No newline at end of file diff --git a/LICENSE.md b/LICENSE.md index f25fcc1542..cc09b22cb9 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -2,33 +2,33 @@ The ModelingToolkit.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2018-22: Yingbo Ma, Christopher Rackauckas, Julia Computing, and > contributors -> +> > Permission is hereby granted, free of charge, to any person obtaining a copy -> +> > of this software and associated documentation files (the "Software"), to deal -> +> > in the Software without restriction, including without limitation the rights -> +> > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -> +> > copies of the Software, and to permit persons to whom the Software is -> +> > furnished to do so, subject to the following conditions: -> +> > The above copyright notice and this permission notice shall be included in all -> +> > copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -> +> > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -> +> > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -> +> > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -> +> > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -> +> > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -> +> > SOFTWARE. diff --git a/NEWS.md b/NEWS.md index 4c8fde0d6e..c1f09df1a8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -488,7 +488,7 @@ now always performs this wrapping. ### Upgrade guide - The function `states` is renamed to `unknowns`. In a similar vein: - + + `unknown_states` is now `solved_unknowns`. + `get_states` is `get_unknowns`. + `get_unknown_states` is now `get_solved_unknowns`. @@ -498,7 +498,7 @@ now always performs this wrapping. - ModelingToolkit.jl now exports common definitions of `t` (time independent variable) and `D` (the first derivative with respect to `t`). Any models made using ModelingToolkit.jl should leverage these common definitions. There are three variants: - + + `t` and `D` use DynamicQuantities.jl units. This is the default for standard library components. + `t_unitful` and `D_unitful` use Unitful.jl units. @@ -522,7 +522,7 @@ now always performs this wrapping. parameters. Their value cannot be modified directly. Whenever a parameter value is changed, dependent parameter values are recalculated. For example, if `@parameters p1 p2 = 3p1` then `p2` can not be modified directly. If `p1` is changed, then `p2` will be updated accordingly. To restore the old behavior: - + + Pass the `split = false` keyword to `structural_simplify`. E.g. `ss = structural_simplify(sys; split = false)`. + Pass `split = false` to `@mtkbuild`. E.g. `@mtkbuild sys = ODESystem(...) split = false`. - Discrete-time system using `Difference` are unsupported. Instead, use the new `Clock`-based syntax. diff --git a/docs/src/API/System.md b/docs/src/API/System.md index 488c7e3555..b12c08b6ce 100644 --- a/docs/src/API/System.md +++ b/docs/src/API/System.md @@ -126,7 +126,7 @@ namespaced version of `var`. Note that this can also be used to access subsystem or analysis points. !!! note - + By default, top-level systems not marked as `complete` will apply their namespace. Systems marked as `complete` will not do this namespacing. This namespacing behavior can be toggled independently of whether the system is completed using [`toggle_namespacing`](@ref) and the diff --git a/docs/src/API/model_building.md b/docs/src/API/model_building.md index 398c3e12a8..735ff9bda5 100644 --- a/docs/src/API/model_building.md +++ b/docs/src/API/model_building.md @@ -242,7 +242,7 @@ has the ability to represent them. Compilation strategies can be implemented ind on top of [`mtkcompile`](@ref) using the `additional_passes` functionality. !!! warn - + These operators are considered experimental API. ```@docs; canonical = false @@ -266,7 +266,7 @@ While ModelingToolkit has the capability to represent state machines, it lacks t to compile and simulate them. !!! warn - + This functionality is considered experimental API ```@docs diff --git a/docs/src/API/variables.md b/docs/src/API/variables.md index 55f060732d..3101adbadc 100644 --- a/docs/src/API/variables.md +++ b/docs/src/API/variables.md @@ -251,7 +251,7 @@ ModelingToolkit.isconstant ``` !!! note - + [`@constants`](@ref) is a convenient way to create `@parameters` with `tunable = false` metadata @@ -422,7 +422,7 @@ system happen at discrete intervals on a clock. While ModelingToolkit cannot yet such systems, it has the capability to represent them. !!! warn - + These operators are considered experimental API. ```@docs diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 848ce411ba..ca16edc47a 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -76,7 +76,7 @@ event = SymbolicDiscreteCallback( This way, we enforce that `y` will remain the same, and `p` will change. !!! warning - + Symbolic affects come with the guarantee that the state after the callback will be consistent. However, when using [general functional affects](@ref func_affects) or [imperative affects](@ref imp_affects) one must be more careful. In diff --git a/docs/src/basics/InputOutput.md b/docs/src/basics/InputOutput.md index b1eb2905df..b99374eaa8 100644 --- a/docs/src/basics/InputOutput.md +++ b/docs/src/basics/InputOutput.md @@ -29,7 +29,7 @@ ModelingToolkit can generate the dynamics of a system, the function ``M\dot x = This function takes a vector of variables that are to be considered inputs, i.e., part of the vector ``u``. Alongside returning the function ``f``, [`ModelingToolkit.generate_control_function`](@ref) also returns the chosen state realization of the system after simplification. This vector specifies the order of the state variables ``x``, while the user-specified vector `u` specifies the order of the input variables ``u``. !!! note "Un-simplified system" - + This function expects `sys` to be un-simplified, i.e., `mtkcompile` or `@mtkcompile` should not be called on the system before passing it into this function. `generate_control_function` calls a special version of `mtkcompile` internally. ### Example: diff --git a/docs/src/basics/Linearization.md b/docs/src/basics/Linearization.md index 14781d2c0c..6a7a77d1f6 100644 --- a/docs/src/basics/Linearization.md +++ b/docs/src/basics/Linearization.md @@ -42,11 +42,11 @@ using ModelingToolkit: inputs, outputs ``` !!! note "Inputs must be unconnected" - + The model above has 4 variables but only three equations, there is no equation specifying the value of `r` since `r` is an input. This means that only unbalanced models can be linearized, or in other words, models that are balanced and can be simulated _cannot_ be linearized. To learn more about this, see [How to linearize a ModelingToolkit model (YouTube)](https://www.youtube.com/watch?v=-XOux-2XDGI&t=395s). Also see [ModelingToolkitStandardLibrary: Linear analysis](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/linear_analysis/) for utilities that make linearization of completed models easier. !!! note "Un-simplified system" - + Linearization expects `sys` to be un-simplified, i.e., `mtkcompile` or `@mtkcompile` should not be called on the system before linearizing. ## Operating point diff --git a/docs/src/basics/MTKLanguage.md b/docs/src/basics/MTKLanguage.md index 2089047ec7..c8481a41dd 100644 --- a/docs/src/basics/MTKLanguage.md +++ b/docs/src/basics/MTKLanguage.md @@ -173,13 +173,13 @@ julia> ModelingToolkit.getdefault(model_c1.v) One or more partial systems can be extended in a higher system with `@extend` statements. This can be done in two ways: - `@extend PartialSystem(var1 = value1)` - + + This is the recommended way of extending a base system. + The default values for the arguments of the base system can be declared in the `@extend` statement. + Note that all keyword arguments of the base system are added as the keyword arguments of the main system. - `@extend var_to_unpack1, var_to_unpack2 = partial_sys = PartialSystem(var1 = value1)` - + + In this method: explicitly list the variables that should be unpacked, provide a name for the partial system and declare the base system. + Note that only the arguments listed out in the declaration of the base system (here: `var1`) are added as the keyword arguments of the higher system. @@ -381,7 +381,7 @@ end ``` !!! note - + For more examples of usage, checkout [ModelingToolkitStandardLibrary.jl](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/) ## [More on `Model.structure`](@id model_structure) diff --git a/docs/src/basics/PrecompileComponents.md b/docs/src/basics/PrecompileComponents.md index d2c83f459c..d5380f6aa7 100644 --- a/docs/src/basics/PrecompileComponents.md +++ b/docs/src/basics/PrecompileComponents.md @@ -159,7 +159,7 @@ is the most reliable tool. Use `@descend` on the component constructor and look (runtime-dispatch) calls. In Cthulhu, press `T` to switch to Typed IR, `o` to enable optimizations, and `w` to highlight unstable code — this view is harder to read but gives the most accurate picture of what the compiler actually infers. The Julia flags `--trace-compile` -and `--trace-compile-timing` are also useful for identifying methods that were not +and `--trace-compile-timing` are also useful for identifying methods that were not precompiled or were invalidated at load time. !!! note "Julia 1.12 and precompilation" diff --git a/docs/src/index.md b/docs/src/index.md index ffece1504f..b10b7c3387 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -37,7 +37,7 @@ If you use ModelingToolkit in your work, please cite the following: ## Feature Summary !!! danger "ModelingToolkit version 11" - + ModelingToolkit version 11 just released. Please refer to the [changelog](https://github.com/SciML/ModelingToolkit.jl/blob/master/NEWS.md) for a summary of the changes. Some documentation pages may be broken while downstream packages update to the new version. @@ -90,7 +90,7 @@ a standard library of prebuilt components for the ModelingToolkit ecosystem. ## Model Import Formats - [CellMLToolkit.jl](https://docs.sciml.ai/CellMLToolkit/stable/): Import [CellML](https://www.cellml.org/) models into ModelingToolkit - + + Repository of more than a thousand pre-made models + Focus on biomedical models in areas such as: Calcium Dynamics, Cardiovascular Circulation, Cell Cycle, Cell Migration, Circadian Rhythms, @@ -100,10 +100,10 @@ a standard library of prebuilt components for the ModelingToolkit ecosystem. Protein Modules, Signal Transduction, and Synthetic Biology. - [SBMLToolkit.jl](https://docs.sciml.ai/SBMLToolkit/stable/): Import [SBML](http://sbml.org/) models into ModelingToolkit - + + Uses the robust libsbml library for parsing and transforming the SBML - [ReactionNetworkImporters.jl](https://docs.sciml.ai/ReactionNetworkImporters/stable/): Import various models into ModelingToolkit - + + Supports the BioNetGen `.net` file + Supports importing networks specified by stoichiometric matrices @@ -115,37 +115,37 @@ Below is an incomplete list of extension libraries one may want to be aware of: - [Catalyst.jl](https://docs.sciml.ai/Catalyst/stable/): Symbolic representations of chemical reactions - + + Symbolically build and represent large systems of chemical reactions + Generate code for ODEs, SDEs, continuous-time Markov Chains, and more + Simulate the models using the SciML ecosystem with O(1) Gillespie methods - [DataDrivenDiffEq.jl](https://docs.sciml.ai/DataDrivenDiffEq/stable/): Automatic identification of equations from data - + + Automated construction of ODEs and DAEs from data + Representations of Koopman operators and Dynamic Mode Decomposition (DMD) - [MomentClosure.jl](https://sciml.github.io/MomentClosure.jl/dev/): Automatic transformation of ReactionSystems into deterministic systems - + + Generates Systems for the moment closures + Allows for geometrically-distributed random reaction rates - [ReactionMechanismSimulator.jl](https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl): Simulating and analyzing large chemical reaction mechanisms - + + Ideal gas and dilute liquid phases. + Constant T and P and constant V adiabatic ideal gas reactors. + Constant T and V dilute liquid reactors. + Diffusion limited rates. Sensitivity analysis for all reactors. + Flux diagrams with molecular images (if molecular information is provided). - [NumCME.jl](https://github.com/voduchuy/NumCME.jl): High-performance simulation of chemical master equations (CME) - + + Transient solution of the CME + Dynamic state spaces + Accepts reaction systems defined using Catalyst.jl DSL. - [FiniteStateProjection.jl](https://github.com/SciML/FiniteStateProjection.jl): High-performance simulation of chemical master equations (CME) via finite state projections - + + Accepts reaction systems defined using Catalyst.jl DSL. ## Compatible Numerical Solvers @@ -158,21 +158,21 @@ the solver libraries which are the numerical targets of the ModelingToolkit system: - [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) - + + Multi-package interface of high performance numerical solvers for `System`, `SDESystem`, and `JumpSystem` - [NonlinearSolve.jl](https://docs.sciml.ai/NonlinearSolve/stable/) - + + High performance numerical solving of `NonlinearSystem` - [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) - + + Multi-package interface for numerical solving `OptimizationSystem` - [NeuralPDE.jl](https://docs.sciml.ai/NeuralPDE/stable/) - + + Physics-Informed Neural Network (PINN) training on `PDESystem` - [MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/) - + + Automated finite difference method (FDM) discretization of `PDESystem` ## Contributing @@ -183,7 +183,7 @@ system: - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions. - There are a few community forums: - + + The #diffeq-bridged and #sciml-bridged channels in the [Julia Slack](https://julialang.org/slack/) + The #diffeq-bridged and #sciml-bridged channels in the diff --git a/docs/src/tutorials/SampledData.md b/docs/src/tutorials/SampledData.md index 9f3f340f46..1f482ac1f0 100644 --- a/docs/src/tutorials/SampledData.md +++ b/docs/src/tutorials/SampledData.md @@ -3,11 +3,11 @@ A sampled-data system contains both continuous-time and discrete-time components, such as a continuous-time plant model and a discrete-time control system. ModelingToolkit supports the modeling and simulation of sampled-data systems by means of *clocks*. !!! danger "Experimental" - + The sampled-data interface is currently experimental and at any time subject to breaking changes **not** respecting semantic versioning. !!! note "Negative shifts" - + The initial release of the sampled-data interface **only supports negative shifts**. A clock can be seen as an *event source*, i.e., when the clock ticks, an event is generated. In response to the event the discrete-time logic is executed, for example, a control signal is computed. For basic modeling of sampled-data systems, the user does not have to interact with clocks explicitly, instead, the modeling is performed using the operators @@ -26,7 +26,7 @@ The [`ShiftIndex`](@ref) operator is used to refer to past and future values of ```math \begin{align} - x(k+1) &= 0.5x(k) + u(k) \\ + x(k+1) &= 0.5x(k) + u(k) \\ y(k) &= x(k) \end{align} ``` diff --git a/docs/src/tutorials/attractors.md b/docs/src/tutorials/attractors.md index 8b5fecbef9..02a7e4d1e7 100644 --- a/docs/src/tutorials/attractors.md +++ b/docs/src/tutorials/attractors.md @@ -14,7 +14,7 @@ Hence the "nonlocal-" component. More differences and pros & cons are discussed in the documentation of Attractors.jl. !!! note "Attractors and basins" - + This tutorial assumes that you have some familiarity with dynamical systems, specifically what are attractors and basins of attraction. If you don't have this yet, we recommend Chapter 1 of the textbook diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md index ffad0d34cc..8ed29ce74d 100644 --- a/docs/src/tutorials/change_independent_variable.md +++ b/docs/src/tutorials/change_independent_variable.md @@ -57,7 +57,7 @@ M2s # display this # hide The derivatives are now with respect to the new independent variable $x$, which can be accessed with `M2.x`. !!! info - + At this point `x`, `M1.x`, `M1s.x`, `M2.x`, `M2s.x` are *three* different variables. Meanwhile `y`, `M1.y`, `M1s.y`, `M2.y` and `M2s.y` are *four* different variables. It can be instructive to inspect these yourself to see their subtle differences. @@ -73,9 +73,9 @@ plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)! ``` !!! tip "Usage tips" - + Look up the documentation of [`change_independent_variable`](@ref) for tips on how to use it. - + For example, if you also need $t(x)$, you can tell it to add a differential equation for the old independent variable in terms of the new one using the [inverse function rule](https://en.wikipedia.org/wiki/Inverse_function_rule) (here $\mathrm{d}t/\mathrm{d}x = 1 / (\mathrm{d}x/\mathrm{d}t)$). If you know an analytical expression between the independent variables (here $t = x/v$), you can also pass it directly to the function to avoid the extra differential equation. ## 2. Alleviating stiffness by changing to logarithmic time diff --git a/docs/src/tutorials/fmi.md b/docs/src/tutorials/fmi.md index f2675429f9..4e22f898cb 100644 --- a/docs/src/tutorials/fmi.md +++ b/docs/src/tutorials/fmi.md @@ -7,7 +7,7 @@ Events, non-floating-point variables and array variables are not supported. Addi time derivatives of FMU states/outputs is not supported. !!! danger "Experimental" - + This functionality is currently experimental and subject to change without a breaking release of ModelingToolkit.jl. @@ -43,7 +43,7 @@ Note how hierarchical names in the FMU (e.g. `mass.m` or `spring.f`) are turned names, with `__` being the namespace separator (`mass__m` and `spring__f`). !!! note - + Eventually we plan to reconstruct a hierarchical system structure mirroring the one indicated by the variables in the FMU. This would allow accessing the above mentioned variables as `model.mass.m` and `model.spring.f` instead of `model.mass__m` and `model.spring__f` respectively. @@ -61,7 +61,7 @@ same reference, and hence refer to the same quantity. Correspondingly, there is `mass__vˍt(t) ~ mass__a(t)` in the system. !!! note - + Any variables and/or parameters that are not part of the FMU should be ignored, as ModelingToolkit creates them to manage the FMU. Unexpected usage of these variables/parameters can lead to errors. diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index beaea968d4..78879d6a5f 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -298,7 +298,7 @@ long enough you will see that `λ = 0` is required for this equation, but since `λ = 1` we end up with a set of equations that are impossible to satisfy. !!! note - + If you would prefer to have an error instead of a warning in the context of non-fully determined systems, pass the keyword argument `fully_determined = true` into the problem constructor. Additionally, any warning about not being fully determined can @@ -562,7 +562,7 @@ sol = solve(iprob) ``` !!! note - + For more information on solving NonlinearProblems and NonlinearLeastSquaresProblems, check out the [NonlinearSolve.jl tutorials!](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/getting_started/). diff --git a/docs/src/tutorials/linear_analysis.md b/docs/src/tutorials/linear_analysis.md index 5317b45fc9..c54aaa3c07 100644 --- a/docs/src/tutorials/linear_analysis.md +++ b/docs/src/tutorials/linear_analysis.md @@ -21,7 +21,7 @@ connect(comp1.output, :analysis_point_name, comp2.input, comp3.input, comp4.inpu ``` !!! warning "Causality" - + Analysis points are *causal*, i.e., they imply a directionality for the flow of information. The order of the connections in the connect statement is thus important, i.e., `connect(out, :name, in)` is different from `connect(in, :name, out)`. The directionality of an analysis point can be thought of as an arrow in a block diagram, where the name of the analysis point applies to the arrow itself. diff --git a/docs/src/tutorials/modelingtoolkitize.md b/docs/src/tutorials/modelingtoolkitize.md index 8cd8015630..c73935f209 100644 --- a/docs/src/tutorials/modelingtoolkitize.md +++ b/docs/src/tutorials/modelingtoolkitize.md @@ -17,7 +17,7 @@ so, ModelingToolkit analysis passes and transformations can be run as intermedia to improve a simulation code before it's passed to the solver. !!! note - + `modelingtoolkitize` does have some limitations, i.e. not all codes that work with the numerical solvers will work with `modelingtoolkitize`. Namely, it requires the ability to trace the equations with Symbolics.jl `Num` types. Generally, a code which is @@ -25,7 +25,7 @@ to improve a simulation code before it's passed to the solver. `modelingtoolkitize`. !!! warn - + `modelingtoolkitize` expressions cannot keep control flow structures (loops), and thus equations with long loops will be translated into large expressions, which can increase the compile time of the equations and reduce the SIMD vectorization achieved by LLVM. diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md index 8342eb788b..9845706844 100644 --- a/docs/src/tutorials/nonlinear.md +++ b/docs/src/tutorials/nonlinear.md @@ -6,7 +6,7 @@ This is, for example, useful for finding the steady state of an ODE. This steady state is reached when the nonlinear system of differential equations equals zero. !!! note - + The high level `@mtkmodel` macro used in the [getting started tutorial](@ref getting_started) is not yet compatible with `NonlinearSystem`. diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md index a932510308..5fb8ceb461 100644 --- a/docs/src/tutorials/optimization.md +++ b/docs/src/tutorials/optimization.md @@ -4,7 +4,7 @@ ModelingToolkit.jl is not only useful for generating initial value problems (`OD The package can also build optimization systems. !!! note - + The high level `@mtkmodel` macro used in the [getting started tutorial](@ref getting_started) is not yet compatible with `OptimizationSystem`. diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md index d7326c36c6..a9ed44c217 100644 --- a/docs/src/tutorials/stochastic_diffeq.md +++ b/docs/src/tutorials/stochastic_diffeq.md @@ -7,7 +7,7 @@ In particular, we show how to represent a as a `SDESystem`. !!! note - + The high level `@mtkmodel` macro used in the [getting started tutorial](@ref getting_started) is not yet compatible with `SDESystem`. diff --git a/ext/MTKFMIExt.jl b/ext/MTKFMIExt.jl index 5acbbae0c0..29f26a99a8 100644 --- a/ext/MTKFMIExt.jl +++ b/ext/MTKFMIExt.jl @@ -417,7 +417,7 @@ function parseFMIVariableName(name::AbstractString) name = @view name[5:(end - 1)] der = 1 else - + der = parse(Int, @view name[(idx + 1):(end - 1)]) name = @view name[5:(idx - 1)] end diff --git a/lib/ModelingToolkitBase/src/problems/docs.jl b/lib/ModelingToolkitBase/src/problems/docs.jl index 104bdcd730..9bfdaa7077 100644 --- a/lib/ModelingToolkitBase/src/problems/docs.jl +++ b/lib/ModelingToolkitBase/src/problems/docs.jl @@ -172,7 +172,7 @@ as in `x(0.5) ~ 1`). More general constraints that should hold over the entire such as `x(t)^2 + y(t)^2`, should be specified as one of the equations used to build the `System`. -If a `System` without `constraints` is specified, it will be treated as an initial value problem. +If a `System` without `constraints` is specified, it will be treated as an initial value problem. ```julia @parameters g t_c = 0.5 @@ -191,7 +191,7 @@ If a `System` without `constraints` is specified, it will be treated as an initi bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false) ``` -If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting +If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting `BVProblem` must be solved using BVDAE solvers, such as Ascher. """ diff --git a/lib/ModelingToolkitBase/src/systems/callbacks.jl b/lib/ModelingToolkitBase/src/systems/callbacks.jl index c13baa5104..9840218d53 100644 --- a/lib/ModelingToolkitBase/src/systems/callbacks.jl +++ b/lib/ModelingToolkitBase/src/systems/callbacks.jl @@ -315,7 +315,7 @@ end const Affect = Union{AffectSystem, ImperativeAffect} """ - SymbolicContinuousCallback(eqs::Vector{Equation}, affect = nothing, iv = nothing; + SymbolicContinuousCallback(eqs::Vector{Equation}, affect = nothing, iv = nothing; affect_neg = affect, initialize = nothing, finalize = nothing, rootfind = SciMLBase.LeftRootFind, initialize_save_discretes = true) @@ -360,7 +360,7 @@ and combined with the remaining `Equation`s. + `ctx` is a user-defined context object passed to `f!` when invoked. This value is aliased for each problem. * A [`ImperativeAffect`](@ref); refer to its documentation for details. -`reinitializealg` is used to set how the system will be reinitialized after the callback. +`reinitializealg` is used to set how the system will be reinitialized after the callback. - Symbolic affects have reinitialization built in. In this case the algorithm will default to SciMLBase.NoInit(), and should **not** be provided. - Functional and imperative affects will default to SciMLBase.CheckInit(), which will error if the system is not properly reinitialized after the callback. If your system is a DAE, pass in an algorithm like SciMLBase.BrownBasicFullInit() to properly re-initialize. @@ -533,7 +533,7 @@ end A callback that triggers at the first timestep that the conditions are satisfied. -The condition can be one of: +The condition can be one of: - Δt::Real - periodic events with period Δt - ts::Vector{Real} - events trigger at these preset times given by `ts` - eqs::Vector{SymbolicT} - events trigger when the condition evaluates to true @@ -542,7 +542,7 @@ The condition can be one of: interest of correctness. The callback will trigger and save at `tspan[1]` if the clock would tick at `tspan[1]`. -Arguments: +Arguments: - iv: The independent variable of the system. This must be specified if the independent variable appears in one of the equations explicitly, as in x ~ t + 1. """ struct SymbolicDiscreteCallback <: AbstractCallback diff --git a/lib/ModelingToolkitBase/src/systems/codegen_utils.jl b/lib/ModelingToolkitBase/src/systems/codegen_utils.jl index 95ea45da94..531e960df4 100644 --- a/lib/ModelingToolkitBase/src/systems/codegen_utils.jl +++ b/lib/ModelingToolkitBase/src/systems/codegen_utils.jl @@ -332,7 +332,7 @@ Base.@nospecializeinfer function build_function_wrapper( p_start += 1 p_end += 1 end - + ir_info = get_ir_info(sys) expr = ir_info.obs_subber(expr) diff --git a/lib/ModelingToolkitBase/src/systems/imperative_affect.jl b/lib/ModelingToolkitBase/src/systems/imperative_affect.jl index 4885fe020f..842bfbbfbb 100644 --- a/lib/ModelingToolkitBase/src/systems/imperative_affect.jl +++ b/lib/ModelingToolkitBase/src/systems/imperative_affect.jl @@ -3,7 +3,7 @@ `ImperativeAffect` is a helper for writing affect functions that will compute observed values and ensure that modified values are correctly written back into the system. The affect function `f` needs to have -the signature +the signature ``` f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple @@ -15,7 +15,7 @@ so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` cu `(; x = y)` or `(; x)` (which aliases `x` as itself) are allowed. The argument NamedTuples (for instance `(;x=y)`) will be populated with the declared values on function entry; if we require `(;x=y)` in `observed` and `y=2`, for example, -then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function `f`. +then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function `f`. The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write diff --git a/lib/ModelingToolkitBase/src/systems/problem_utils.jl b/lib/ModelingToolkitBase/src/systems/problem_utils.jl index d6c9400fd0..5497536df6 100644 --- a/lib/ModelingToolkitBase/src/systems/problem_utils.jl +++ b/lib/ModelingToolkitBase/src/systems/problem_utils.jl @@ -823,7 +823,7 @@ function get_mtkparameters_reconstructor( Vector{bufsize.type} end ) - + diffcache_params = SU.getmetadata(dstsys, DiffCacheParams, Dict{SymbolicT, Int}())::Dict{SymbolicT, Int} diffcache_buffer_idx = 0 if !isempty(diffcache_params) diff --git a/lib/ModelingToolkitBase/src/variables.jl b/lib/ModelingToolkitBase/src/variables.jl index 04a107cf11..84859d633e 100644 --- a/lib/ModelingToolkitBase/src/variables.jl +++ b/lib/ModelingToolkitBase/src/variables.jl @@ -366,7 +366,7 @@ end distribute_shift(eq::Equation) = distribute_shift(eq.lhs) ~ distribute_shift(eq.rhs) distribute_shift(var::Union{Num, Arr}) = distribute_shift(unwrap(var)) """ -Distribute a shift applied to a whole expression or equation. +Distribute a shift applied to a whole expression or equation. Shift(t, 1)(x + y) will become Shift(t, 1)(x) + Shift(t, 1)(y). Only shifts variables whose independent variable is the same t that appears in the Shift (i.e. constants, time-independent parameters, etc. do not get shifted). """ @@ -857,7 +857,7 @@ getmisc(x::SymbolicT) = Symbolics.getmetadata(x, VariableMisc, nothing) hasmisc(x) Determine whether a symbolic variable `x` has misc -metadata associated with it. +metadata associated with it. See also [`getmisc(x)`](@ref). """ @@ -897,15 +897,15 @@ An operator that evaluates time-dependent variables at a specific absolute time - `t::Union{SymbolicT, Number}`: The absolute time at which to evaluate the variable. # Description -`EvalAt` is used to evaluate time-dependent variables at a specific time point. This is particularly -useful in optimization problems where you need to specify constraints or costs at particular moments +`EvalAt` is used to evaluate time-dependent variables at a specific time point. This is particularly +useful in optimization problems where you need to specify constraints or costs at particular moments in time, or delay differential equations for setting a delay time. -The operator works by replacing the time argument of time-dependent variables with the specified +The operator works by replacing the time argument of time-dependent variables with the specified time `t`. For variables that don't depend on time, `EvalAt` returns them unchanged. # Behavior -- For time-dependent variables like `x(t)`, `EvalAt(τ)(x)` returns `x(τ)` +- For time-dependent variables like `x(t)`, `EvalAt(τ)(x)` returns `x(τ)` - For time-independent parameters, `EvalAt` returns them unchanged - For derivatives, `EvalAt` evaluates the derivative at the specified time - For arrays of variables, `EvalAt` is applied element-wise diff --git a/lib/ModelingToolkitBase/test/input_output_handling.jl b/lib/ModelingToolkitBase/test/input_output_handling.jl index ea4e1fc4e1..0e6fdb6a4e 100644 --- a/lib/ModelingToolkitBase/test/input_output_handling.jl +++ b/lib/ModelingToolkitBase/test/input_output_handling.jl @@ -517,8 +517,8 @@ if @isdefined(ModelingToolkit) f, x, p, simplified_system = ModelingToolkit.generate_control_function(model, [u;]) - x0 = ModelingToolkitBase.get_u0(simplified_system, []) - p = ModelingToolkitBase.get_p(simplified_system, []) + x0 = ModelingToolkitBase.get_u0(simplified_system, []) + p = ModelingToolkitBase.get_p(simplified_system, []) @test f[1](x0, zeros(2), p, 0) != f[1](x0, ones(2), p, 0) end diff --git a/lib/ModelingToolkitBase/test/odesystem.jl b/lib/ModelingToolkitBase/test/odesystem.jl index c66701aa0b..576e38476d 100644 --- a/lib/ModelingToolkitBase/test/odesystem.jl +++ b/lib/ModelingToolkitBase/test/odesystem.jl @@ -1670,7 +1670,7 @@ end using ModelingToolkitBase using ModelingToolkitBase: t_nounits as t, D_nounits as D @variables x(t)[1:3]=[0,0,1] - @variables u1(t)=0 u2(t)=0 + @variables u1(t)=0 u2(t)=0 y₁, y₂, y₃ = x k₁, k₂, k₃ = 1,1,1 eqs = [ diff --git a/lib/SciCompDSL/src/model_parsing.jl b/lib/SciCompDSL/src/model_parsing.jl index eef1bb5b1d..b0bb62e3e6 100644 --- a/lib/SciCompDSL/src/model_parsing.jl +++ b/lib/SciCompDSL/src/model_parsing.jl @@ -1344,9 +1344,9 @@ end function check_event_syntax(line) if line.args[1] == :(=>) return true - elseif (line.head == :(...)) || - (line.head == :generator) || - (line.head == :comprehension)|| + elseif (line.head == :(...)) || + (line.head == :generator) || + (line.head == :comprehension)|| (line.head == :if) return check_event_syntax(line.args[1]) else diff --git a/lib/SciCompDSL/test/model_parsing.jl b/lib/SciCompDSL/test/model_parsing.jl index 16fdd2dc1c..272e363554 100644 --- a/lib/SciCompDSL/test/model_parsing.jl +++ b/lib/SciCompDSL/test/model_parsing.jl @@ -480,7 +480,7 @@ end using ModelingToolkitBase: D_nounits @testset "Event handling in MTKModel" begin - + cond = true @mtkmodel A begin @variables begin @@ -532,9 +532,9 @@ using ModelingToolkitBase: D_nounits @mtkcompile model = M() - u0 = [model.x => 10, - model.y => 0, - model.z => 0, + u0 = [model.x => 10, + model.y => 0, + model.z => 0, model.c => 0, model.d => 0] diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 4812820c04..bdba846f91 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -184,5 +184,5 @@ function __init__() SU.hashcons(unwrap(ODE_C), true) SU.hashcons(SCC_EXPLICITFUN_CACHE_OUT, true) end - + end # module diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl index 15aebcdcb1..9f52f321fe 100644 --- a/src/problems/sccnonlinearproblem.jl +++ b/src/problems/sccnonlinearproblem.jl @@ -289,7 +289,7 @@ function build_caches!(sys::System, decomposition::SCCDecomposition) for i in eachindex(decomposition.subsystems) empty!(banned_vars) empty!(state) - + subsys = decomposition.subsystems[i] union!(banned_vars, unknowns(subsys)) for u in unknowns(subsys) diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl index 82123d5201..f2b94fb63a 100644 --- a/src/systems/alias_elimination.jl +++ b/src/systems/alias_elimination.jl @@ -157,7 +157,7 @@ function find_perfect_aliases!( state.always_present[v] = true continue end - + push!(vars_to_rm, v) subs[fullvars[v]] = fullvars[target] push!(state.additional_observed, fullvars[v] ~ fullvars[target]) diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl index ee51ffbd4a..769c56a6e0 100644 --- a/src/systems/analysis_points.jl +++ b/src/systems/analysis_points.jl @@ -18,7 +18,7 @@ const DOC_LOOP_OPENINGS = """ const DOC_SYS_MODIFIER = """ - `system_modifier`: A function taking the transformed system and applying any additional transformations, returning the modified system. The modified system - is passed to `linearization_function`. + is passed to `linearization_function`. """ """ diff --git a/src/systems/if_lifting.jl b/src/systems/if_lifting.jl index 3cd218213c..ae2ce59300 100644 --- a/src/systems/if_lifting.jl +++ b/src/systems/if_lifting.jl @@ -547,7 +547,7 @@ const CONDITION_SIMPLIFIER = Rewriters.Fixpoint( If lifting converts (nested) if statements into a series of continuous events + a logically equivalent if statement + parameters. Lifting proceeds through the following process: -* rewrite comparisons to be of the form eqn [op] 0; subtract the RHS from the LHS +* rewrite comparisons to be of the form eqn [op] 0; subtract the RHS from the LHS * replace comparisons with generated parameters; for each comparison eqn [op] 0, generate an event (dependent on op) that sets the parameter !!! warn diff --git a/test/downstream/inversemodel.jl b/test/downstream/inversemodel.jl index a33db82bbe..09a02561cd 100644 --- a/test/downstream/inversemodel.jl +++ b/test/downstream/inversemodel.jl @@ -34,7 +34,7 @@ rc = 0.25 # Reference concentration c = RealOutput() T = RealOutput() end - + xc = @something(xc, c0) xT = @something(xT, T0) vars = @variables begin diff --git a/test_compiler_options_fallback.jl b/test_compiler_options_fallback.jl deleted file mode 100644 index ce65137e3d..0000000000 --- a/test_compiler_options_fallback.jl +++ /dev/null @@ -1,97 +0,0 @@ -# Test that module-level compiler options are inherited by functions eval'd in those modules, -# including RuntimeGeneratedFunctions. -# -# Run with: julia --project test_compiler_options_fallback.jl - -using RuntimeGeneratedFunctions -using InteractiveUtils: code_native, code_typed - -# Create modules with different compiler options -module EvalOpt0 - Base.Experimental.@compiler_options optimize=0 - using RuntimeGeneratedFunctions - RuntimeGeneratedFunctions.init(@__MODULE__) -end - -module EvalOpt1 - Base.Experimental.@compiler_options optimize=1 - using RuntimeGeneratedFunctions - RuntimeGeneratedFunctions.init(@__MODULE__) -end - -module EvalOpt0NoInfer - Base.Experimental.@compiler_options optimize=0 infer=false - using RuntimeGeneratedFunctions - RuntimeGeneratedFunctions.init(@__MODULE__) -end - -module EvalDefault - using RuntimeGeneratedFunctions - RuntimeGeneratedFunctions.init(@__MODULE__) -end - -# A non-trivial expression to compile -test_expr = :(function (x) - y = sin(x) + cos(x) - return y * y + exp(x) -end) - -# Test eval path -println("=== Testing eval path ===") -f_default = EvalDefault.eval(test_expr) -f_opt0 = EvalOpt0.eval(test_expr) -f_opt1 = EvalOpt1.eval(test_expr) -f_opt0_noinfer = EvalOpt0NoInfer.eval(test_expr) - -x = 1.5 -println("Default: f($x) = $(f_default(x))") -println("Opt0: f($x) = $(f_opt0(x))") -println("Opt1: f($x) = $(f_opt1(x))") -println("Opt0+no infer: f($x) = $(f_opt0_noinfer(x))") -println("Results match: ", f_default(x) == f_opt0(x) == f_opt1(x) == f_opt0_noinfer(x)) - -# Test RGF path -println("\n=== Testing RuntimeGeneratedFunction path ===") -rgf_default = RuntimeGeneratedFunction(EvalDefault, EvalDefault, test_expr) -rgf_opt0 = RuntimeGeneratedFunction(EvalOpt0, EvalOpt0, test_expr) -rgf_opt1 = RuntimeGeneratedFunction(EvalOpt1, EvalOpt1, test_expr) -rgf_opt0_noinfer = RuntimeGeneratedFunction(EvalOpt0NoInfer, EvalOpt0NoInfer, test_expr) - -println("RGF Default: f($x) = $(rgf_default(x))") -println("RGF Opt0: f($x) = $(rgf_opt0(x))") -println("RGF Opt1: f($x) = $(rgf_opt1(x))") -println("RGF Opt0+no infer: f($x) = $(rgf_opt0_noinfer(x))") -println("Results match: ", rgf_default(x) == rgf_opt0(x) == rgf_opt1(x) == rgf_opt0_noinfer(x)) - -# Check that the compiler options actually take effect by inspecting code_native. -# At opt level 0, the native code should contain more instructions (no inlining/folding). -println("\n=== Checking native code differences ===") -native_default = sprint(code_native, rgf_default, Tuple{Float64}) -native_opt0 = sprint(code_native, rgf_opt0, Tuple{Float64}) -native_opt1 = sprint(code_native, rgf_opt1, Tuple{Float64}) - -# Count call instructions — opt0 should have more calls (nothing inlined) -count_calls(s) = count(r"callq?\s"i, s) -println("Call instructions - Default: $(count_calls(native_default)), Opt0: $(count_calls(native_opt0)), Opt1: $(count_calls(native_opt1))") -println("Raw instruction count - Default: $(count('\n', native_default)), Opt0: $(count('\n', native_opt0)), Opt1: $(count('\n', native_opt1))") -if count_calls(native_opt0) > count_calls(native_opt1) - println("Opt0 has more calls than Opt1 -> module-level options ARE taking effect") -elseif count_calls(native_opt0) == count_calls(native_opt1) - println("Same call count — check raw instruction counts above for differences") -else - println("WARNING: Opt0 has fewer calls than Opt1 — unexpected") -end - -# Check inference works/doesn't work -println("\n=== Checking infer=false effect ===") -ci_default = only(code_typed(rgf_default, Tuple{Float64})) -ci_noinfer = only(code_typed(rgf_opt0_noinfer, Tuple{Float64})) -println("Default return type: $(ci_default[2])") -println("Opt0+no infer return type: $(ci_noinfer[2])") -if ci_noinfer[2] === Any - println("infer=false IS taking effect (return type is Any)") -else - println("infer=false may NOT be taking effect (return type is concrete)") -end - -println("\nDone.") From a0b1c0da4dbd7a7488912b73b35cda3aa0841daa Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 20 May 2026 16:23:06 +0530 Subject: [PATCH 12/24] build: bump MTKBase patch version --- lib/ModelingToolkitBase/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ModelingToolkitBase/Project.toml b/lib/ModelingToolkitBase/Project.toml index 6781128148..a33f98b39e 100644 --- a/lib/ModelingToolkitBase/Project.toml +++ b/lib/ModelingToolkitBase/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkitBase" uuid = "7771a370-6774-4173-bd38-47e70ca0b839" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "1.36.2" +version = "1.36.3" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From f6a820b56778c78cefb71f312e6698a988af1f12 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 20 May 2026 16:23:11 +0530 Subject: [PATCH 13/24] build: bump MTK patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 88ffd63680..174fb9e80c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" -version = "11.26.3" +version = "11.26.4" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] [deps] From 36a6a95797f28d4ea6c782103a27c4a0b7df2e12 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:14:05 +0800 Subject: [PATCH 14/24] test: failing Q1 basic test for homotopy rewrite_trivial --- lib/ModelingToolkitBase/test/homotopy_init.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 lib/ModelingToolkitBase/test/homotopy_init.jl diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl new file mode 100644 index 0000000000..76a6f74bce --- /dev/null +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -0,0 +1,15 @@ +using Test +using ModelingToolkitBase +using ModelingToolkitBase: rewrite_trivial +using Symbolics + +@testset "homotopy operator — L0 trivial rewrite" begin + @testset "Q1 basic" begin + @variables x + @parameters p + expr = homotopy(x^2 - p, x - sqrt(p)) + rewritten = rewrite_trivial(expr) + @test isequal(Symbolics.unwrap(rewritten), Symbolics.unwrap(x^2 - p)) + @test !occursin("homotopy", repr(rewritten)) + end +end From 8c28ed9f265b62e921c9a793de0e803cfbede587 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:19:47 +0800 Subject: [PATCH 15/24] feat: register homotopy operator + rewrite_trivial (L0 trivial form) Implements Modelica spec 3.7.4 trivial form: homotopy(a, s) is a registered symbolic primitive that rewrites to actual=a. Hand-written TermInterface walk (not @rule) per #1385 feedback. Metadata + symtype propagated on rebuilt nodes. Runtime fallback methods for Real and AbstractArray inputs. --- .../src/ModelingToolkitBase.jl | 2 + .../src/systems/homotopy.jl | 71 +++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 lib/ModelingToolkitBase/src/systems/homotopy.jl diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl index da845565bf..3a71cb2124 100644 --- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl +++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl @@ -194,6 +194,7 @@ include("systems/ir_info.jl") include("systems/codegen_utils.jl") include("problems/docs.jl") include("systems/codegen.jl") +include("systems/homotopy.jl") include("systems/problem_utils.jl") include("problems/compatibility.jl") @@ -270,6 +271,7 @@ export liouville_transform, change_independent_variable, export respecialize export PDESystem export Differential, expand_derivatives, @derivatives +export homotopy export Equation export Term export SymScope, LocalScope, ParentScope, GlobalScope diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl new file mode 100644 index 0000000000..57e528fee1 --- /dev/null +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -0,0 +1,71 @@ +# Modelica homotopy operator — L0 trivial form (spec 3.7.4). +# `homotopy(actual, simplified)` is a symbolic primitive; at runtime and in +# the init pipeline it evaluates to `actual`. The `simplified` argument is +# carried symbolically so future PRs (L1 parameter sweep) can use it for +# continuation without changing call sites. + +using Symbolics +using SymbolicUtils + +@register_symbolic homotopy(actual, simplified) + +""" + homotopy(actual, simplified) + +Modelica homotopy operator (spec 3.7.4). At runtime and at L0 trivial rewrite +time, equivalent to `actual`. Future PRs may use `simplified` as a starting +point for parameter-sweep continuation during initialization. +""" +homotopy(actual::Real, simplified::Real) = actual +homotopy(actual::AbstractArray, simplified::AbstractArray) = actual + +""" + rewrite_trivial(expr) + +Recursively replace every `homotopy(a, s)` node in `expr` with `a`. Preserves +`metadata` and `symtype` on every rebuilt node. Implemented as a hand-written +walk over `TermInterface` — does NOT use `@rule` (per AayushSabharwal feedback +in #1385: per-node matcher overhead is unacceptable at scale). +""" +function rewrite_trivial(expr) + x = Symbolics.unwrap(expr) + return _rewrite_trivial(x) +end + +function _rewrite_trivial(x) + if !SymbolicUtils.iscall(x) + return x + end + op = SymbolicUtils.operation(x) + args = SymbolicUtils.arguments(x) + new_args = map(_rewrite_trivial, args) + if op === homotopy + length(new_args) == 2 || throw(ArgumentError( + "homotopy() expects exactly 2 arguments, got $(length(new_args))")) + return new_args[1] + end + return SymbolicUtils.maketerm( + typeof(x), op, new_args, + SymbolicUtils.metadata(x), + ) +end + +""" + has_homotopy(expr) + +Return `true` iff `expr` contains at least one `homotopy(...)` node. +""" +function has_homotopy(expr) + x = Symbolics.unwrap(expr) + return _has_homotopy(x) +end + +function _has_homotopy(x) + if !SymbolicUtils.iscall(x) + return false + end + if SymbolicUtils.operation(x) === homotopy + return true + end + return any(_has_homotopy, SymbolicUtils.arguments(x)) +end From 2b7b770ad609867abe31ca1fee7358c1c1369779 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:32:01 +0800 Subject: [PATCH 16/24] refactor: align homotopy.jl with MTKBase style + tighten test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Per code review: - Drop top-level 'using Symbolics' / 'using SymbolicUtils' inside the included file — MTKBase parent module already imports iscall, operation, arguments, maketerm, metadata bare (L21-23 of ModelingToolkitBase.jl) and re-exports Symbolics, so qualified prefixes are redundant noise versus surrounding files (utils.jl, codegen_utils.jl, if_lifting.jl, etc.). - Drop unreachable AbstractArray runtime fallback (@register_symbolic is scalar-only; users go through broadcast). - Drop arity ArgumentError guard in _rewrite_trivial (@register_symbolic enforces 2-arity at parse time; guard is dead code). - Soften docstring: metadata is preserved explicitly via maketerm arg; symtype is inferred by maketerm itself, not preserved by us. - Replace brittle !occursin("homotopy", repr(rewritten)) assertion with !has_homotopy(rewritten), which uses the public predicate and gives has_homotopy test coverage as a bonus. --- .../src/systems/homotopy.jl | 40 ++++++------------- lib/ModelingToolkitBase/test/homotopy_init.jl | 4 +- 2 files changed, 15 insertions(+), 29 deletions(-) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl index 57e528fee1..e01b720b6a 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -1,12 +1,3 @@ -# Modelica homotopy operator — L0 trivial form (spec 3.7.4). -# `homotopy(actual, simplified)` is a symbolic primitive; at runtime and in -# the init pipeline it evaluates to `actual`. The `simplified` argument is -# carried symbolically so future PRs (L1 parameter sweep) can use it for -# continuation without changing call sites. - -using Symbolics -using SymbolicUtils - @register_symbolic homotopy(actual, simplified) """ @@ -17,37 +8,32 @@ time, equivalent to `actual`. Future PRs may use `simplified` as a starting point for parameter-sweep continuation during initialization. """ homotopy(actual::Real, simplified::Real) = actual -homotopy(actual::AbstractArray, simplified::AbstractArray) = actual """ rewrite_trivial(expr) Recursively replace every `homotopy(a, s)` node in `expr` with `a`. Preserves -`metadata` and `symtype` on every rebuilt node. Implemented as a hand-written -walk over `TermInterface` — does NOT use `@rule` (per AayushSabharwal feedback -in #1385: per-node matcher overhead is unacceptable at scale). +`metadata` on every rebuilt node; `symtype` is inferred by `maketerm`. +Implemented as a hand-written walk over the TermInterface protocol — does NOT +use `@rule` (per AayushSabharwal feedback in #1385: per-node matcher overhead +is unacceptable at scale). """ function rewrite_trivial(expr) - x = Symbolics.unwrap(expr) + x = unwrap(expr) return _rewrite_trivial(x) end function _rewrite_trivial(x) - if !SymbolicUtils.iscall(x) + if !iscall(x) return x end - op = SymbolicUtils.operation(x) - args = SymbolicUtils.arguments(x) + op = operation(x) + args = arguments(x) new_args = map(_rewrite_trivial, args) if op === homotopy - length(new_args) == 2 || throw(ArgumentError( - "homotopy() expects exactly 2 arguments, got $(length(new_args))")) return new_args[1] end - return SymbolicUtils.maketerm( - typeof(x), op, new_args, - SymbolicUtils.metadata(x), - ) + return maketerm(typeof(x), op, new_args, metadata(x)) end """ @@ -56,16 +42,16 @@ end Return `true` iff `expr` contains at least one `homotopy(...)` node. """ function has_homotopy(expr) - x = Symbolics.unwrap(expr) + x = unwrap(expr) return _has_homotopy(x) end function _has_homotopy(x) - if !SymbolicUtils.iscall(x) + if !iscall(x) return false end - if SymbolicUtils.operation(x) === homotopy + if operation(x) === homotopy return true end - return any(_has_homotopy, SymbolicUtils.arguments(x)) + return any(_has_homotopy, arguments(x)) end diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 76a6f74bce..4c9a59a3a6 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -1,6 +1,6 @@ using Test using ModelingToolkitBase -using ModelingToolkitBase: rewrite_trivial +using ModelingToolkitBase: rewrite_trivial, has_homotopy using Symbolics @testset "homotopy operator — L0 trivial rewrite" begin @@ -10,6 +10,6 @@ using Symbolics expr = homotopy(x^2 - p, x - sqrt(p)) rewritten = rewrite_trivial(expr) @test isequal(Symbolics.unwrap(rewritten), Symbolics.unwrap(x^2 - p)) - @test !occursin("homotopy", repr(rewritten)) + @test !has_homotopy(rewritten) end end From 86b613ed9ea98d6273ca4bc7cac01694e5d5bef5 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:35:18 +0800 Subject: [PATCH 17/24] test: Q2 nested homotopy collapses to innermost actual --- lib/ModelingToolkitBase/test/homotopy_init.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 4c9a59a3a6..85d7d65144 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -12,4 +12,20 @@ using Symbolics @test isequal(Symbolics.unwrap(rewritten), Symbolics.unwrap(x^2 - p)) @test !has_homotopy(rewritten) end + + @testset "Q2 nested" begin + @variables x + @parameters p + inner = homotopy(x^2 - p, x - sqrt(p)) + outer = homotopy(inner, x - 1) + @test isequal(Symbolics.unwrap(rewrite_trivial(outer)), + Symbolics.unwrap(x^2 - p)) + + triple = homotopy(homotopy(homotopy(x, x + 1), x + 2), x + 3) + @test isequal(Symbolics.unwrap(rewrite_trivial(triple)), + Symbolics.unwrap(x)) + + @test !has_homotopy(rewrite_trivial(outer)) + @test !has_homotopy(rewrite_trivial(triple)) + end end From 414261c0701daa4d2b90efda3b0a0e14a176bbd9 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:37:27 +0800 Subject: [PATCH 18/24] test: Q3 homotopy preserved inside Base.ifelse branch structure --- lib/ModelingToolkitBase/test/homotopy_init.jl | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 85d7d65144..6977dbf08b 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -28,4 +28,27 @@ using Symbolics @test !has_homotopy(rewrite_trivial(outer)) @test !has_homotopy(rewrite_trivial(triple)) end + + @testset "Q3 Base.ifelse branches" begin + @variables x + @parameters p + cond = p > 0 + branch_expr = Base.ifelse( + cond, + homotopy(x^2 - p, x - sqrt(abs(p))), + homotopy(-(x^2) - p, x + sqrt(abs(p))), + ) + rewritten = rewrite_trivial(branch_expr) + + # Outer ifelse structure preserved + @test occursin("ifelse", repr(rewritten)) + # Both branches' homotopy nodes are gone + @test !has_homotopy(rewritten) + + # Folding the cond to true/false picks the right actual + true_folded = Symbolics.simplify(Symbolics.substitute(rewritten, Dict(cond => true))) + false_folded = Symbolics.simplify(Symbolics.substitute(rewritten, Dict(cond => false))) + @test isequal(Symbolics.unwrap(true_folded), Symbolics.unwrap(x^2 - p)) + @test isequal(Symbolics.unwrap(false_folded), Symbolics.unwrap(-(x^2) - p)) + end end From 20be7b0ba4b4cee72ab7c011c318c284278fa848 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:38:00 +0800 Subject: [PATCH 19/24] test: Q4 vectorized homotopy via broadcast --- lib/ModelingToolkitBase/test/homotopy_init.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 6977dbf08b..92be761255 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -51,4 +51,18 @@ using Symbolics @test isequal(Symbolics.unwrap(true_folded), Symbolics.unwrap(x^2 - p)) @test isequal(Symbolics.unwrap(false_folded), Symbolics.unwrap(-(x^2) - p)) end + + @testset "Q4 broadcast" begin + @variables x[1:3] p[1:3] + actuals = [x[i]^2 - p[i] for i in 1:3] + simplifieds = [x[i] - p[i] for i in 1:3] + vec_expr = homotopy.(actuals, simplifieds) + rewritten = rewrite_trivial.(vec_expr) + @test length(rewritten) == 3 + for i in 1:3 + @test isequal(Symbolics.unwrap(rewritten[i]), + Symbolics.unwrap(actuals[i])) + @test !has_homotopy(rewritten[i]) + end + end end From 4bcc771fe81dd77ad89b9d9b5772fee66a685a0a Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 14:50:29 +0800 Subject: [PATCH 20/24] test: failing integration test for homotopy in InitializationProblem Adds has_homotopy_in_equations and rewrite_trivial_in_equations helpers to homotopy.jl (internal use; not exported). Adds RED integration test asserting that an InitializationProblem built from a NonlinearSystem containing homotopy(x^2 - p, x - 1) = 0 solves the actual branch (x^2 = p) rather than the simplified branch (x = 1). Currently fails because the init pipeline does not yet apply rewrite_trivial. --- .../src/systems/homotopy.jl | 27 +++++++++++++++++++ lib/ModelingToolkitBase/test/homotopy_init.jl | 24 +++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/lib/ModelingToolkitBase/src/systems/homotopy.jl b/lib/ModelingToolkitBase/src/systems/homotopy.jl index e01b720b6a..5c0f0abfd0 100644 --- a/lib/ModelingToolkitBase/src/systems/homotopy.jl +++ b/lib/ModelingToolkitBase/src/systems/homotopy.jl @@ -55,3 +55,30 @@ function _has_homotopy(x) end return any(_has_homotopy, arguments(x)) end + +""" + has_homotopy_in_equations(eqs) + +Return `true` iff any equation in `eqs` (lhs or rhs) contains a `homotopy(...)` +node. `eqs` is an iterable of `Equation`. +""" +function has_homotopy_in_equations(eqs) + for eq in eqs + if has_homotopy(eq.lhs) || has_homotopy(eq.rhs) + return true + end + end + return false +end + +""" + rewrite_trivial_in_equations(eqs) + +Return a new vector of `Equation`s with every `homotopy(a, s)` replaced by `a` +on both lhs and rhs. Original `eqs` not mutated; the system caller is +responsible for swapping the equation vector into the system. +""" +function rewrite_trivial_in_equations(eqs) + return [Equation(rewrite_trivial(eq.lhs), rewrite_trivial(eq.rhs)) + for eq in eqs] +end diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 92be761255..5f79052e5d 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -65,4 +65,28 @@ using Symbolics @test !has_homotopy(rewritten[i]) end end + + @testset "Integration: homotopy in NonlinearSystem init" begin + using ModelingToolkitBase: System, mtkcompile + using ModelingToolkitBase: InitializationProblem + using NonlinearSolve: NewtonRaphson, solve + using SciMLBase + + @variables x + @parameters p + # Equation: homotopy(actual = x^2 - p, simplified = x - 1) = 0 + # L0 trivial must solve actual: x^2 = p ⇒ x = ±√p. + # If rewrite is broken and simplified got selected: x = 1. + # At p = 9: actual root x = 3 (or -3); simplified root x = 1. + # The x^2 ≈ p assertion distinguishes the two outcomes. + eqs = [0 ~ homotopy(x^2 - p, x - 1)] + @named sys = System(eqs, [x], [p]) + sys = mtkcompile(sys) + + prob = InitializationProblem{false}(sys, nothing, Dict(p => 9.0, x => 2.5)) + sol = solve(prob, NewtonRaphson(); abstol = 1e-10, reltol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation solved (NOT simplified, which would give x=1) + @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root + end end From a47c6fefc3f941084255b59cfa2854ae1a2cd800 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:06:15 +0800 Subject: [PATCH 21/24] test: refine RED integration test to ODESystem with algebraic homotopy The previous NonlinearSystem variant never exercised the rewrite hook because generate_initializesystem_timeindependent does not migrate the parent system's algebraic equations into the initialization system, so the failure mode was an empty u0 unrelated to homotopy. The new design uses an ODESystem with an algebraic constraint `0 ~ homotopy(y^2 - p, y - 1)` and a free init guess for `y`. This routes through generate_initializesystem_timevarying, which preserves the homotopy-bearing equation in the init system. The structural assertion `!has_homotopy_in_equations(equations(prob.f.sys))` distinguishes hooked from unhooked: without the rewrite, the init system carries the opaque homotopy() symbolic call even though runtime dispatch happens to give a numerically correct answer via the homotopy(::Real, ::Real) fallback. --- lib/ModelingToolkitBase/test/homotopy_init.jl | 33 ++++++++++++------- 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index 5f79052e5d..f2f28dbd5a 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -66,27 +66,36 @@ using Symbolics end end - @testset "Integration: homotopy in NonlinearSystem init" begin + @testset "Integration: homotopy in ODESystem init" begin using ModelingToolkitBase: System, mtkcompile using ModelingToolkitBase: InitializationProblem + using ModelingToolkitBase: t_nounits as t, D_nounits as D + using ModelingToolkitBase: has_homotopy_in_equations, equations using NonlinearSolve: NewtonRaphson, solve using SciMLBase - @variables x + @variables x(t) y(t) @parameters p - # Equation: homotopy(actual = x^2 - p, simplified = x - 1) = 0 - # L0 trivial must solve actual: x^2 = p ⇒ x = ±√p. - # If rewrite is broken and simplified got selected: x = 1. - # At p = 9: actual root x = 3 (or -3); simplified root x = 1. - # The x^2 ≈ p assertion distinguishes the two outcomes. - eqs = [0 ~ homotopy(x^2 - p, x - 1)] - @named sys = System(eqs, [x], [p]) + # `y` is algebraic, constrained by homotopy(actual = y^2 - p, simplified = y - 1). + # L0 trivial rewrite must replace the constraint with `y^2 - p = 0`, so init solves + # to y = ±√p. At p = 9 with guess y = 2.5, init converges to y ≈ 3. If the rewrite + # never fired, the init system would carry the opaque `homotopy(...)` symbolic call + # — which is the structural assertion below. + eqs = [D(x) ~ -x, + 0 ~ homotopy(y^2 - p, y - 1)] + @named sys = System(eqs, t; guesses = [y => 2.5]) sys = mtkcompile(sys) - prob = InitializationProblem{false}(sys, nothing, Dict(p => 9.0, x => 2.5)) + prob = InitializationProblem{false}(sys, 0.0, Dict(x => 1.0, p => 9.0)) + + # Structural: after the init pipeline applies rewrite_trivial, no equation in + # the wrapped initialization system should still contain a `homotopy(...)` node. + @test !has_homotopy_in_equations(equations(prob.f.sys)) + + # Numerical: init must converge to the actual root. sol = solve(prob, NewtonRaphson(); abstol = 1e-10, reltol = 1e-10) @test SciMLBase.successful_retcode(sol) - @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation solved (NOT simplified, which would give x=1) - @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root + @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation y^2 = p satisfied + @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root (y = 1) end end From 3cca63680ee93126119d3b3831e91508482152b7 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:06:31 +0800 Subject: [PATCH 22/24] feat: rewrite homotopy to actual inside InitializationProblem Detects homotopy() in initialization-system equations and applies rewrite_trivial. Gated by has_homotopy_in_equations so non-homotopy systems pay zero traversal cost. Solves the L0 trivial requirement of Modelica spec 3.7.4. The hook lands before mtkcompile(isys) so that downstream structural transforms operate on the actual equations rather than opaque homotopy() calls; future PRs that lower into a parameter sweep will replace this branch when a HomotopySweep algorithm is in use. --- .../src/problems/initializationproblem.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl index e24346cfea..7056591350 100644 --- a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl @@ -82,6 +82,17 @@ All other keyword arguments are forwarded to the wrapped problem constructor. binds[get_iv(sys)::SymbolicT] = Symbolics.COMMON_ZERO @set! isys.bindings = ROSymmapT(binds) end + # L0 trivial homotopy rewrite (Modelica spec 3.7.4 trivial form). + # Only runs when `homotopy(...)` is actually present, so non-homotopy + # systems pay zero overhead. PR2 (L1 parameter sweep) will replace this + # branch with λ-lowering when an explicit HomotopySweep algorithm is in use. + # Applied before `mtkcompile` so that downstream structural transforms + # operate on the actual equations, not on opaque `homotopy(...)` calls. + if has_homotopy_in_equations(equations(isys)) + new_eqs = rewrite_trivial_in_equations(equations(isys)) + @set! isys.eqs = new_eqs + end + if simplify_system isys = mtkcompile(isys; fully_determined, split = is_split(sys), initsys_mtkcompile_kwargs...) end From 743fd9e254ad7cc613e04b159c5668da47eabc5a Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:06:47 +0800 Subject: [PATCH 23/24] test: register homotopy_init.jl in Initialization group The new homotopy_init.jl covers rewrite_trivial (unit) plus an ODESystem InitializationProblem integration check. It belongs alongside the other Initialization-group entries (guess_propagation, initializationsystem, initial_values) rather than the broader InterfaceI bucket. --- lib/ModelingToolkitBase/test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/ModelingToolkitBase/test/runtests.jl b/lib/ModelingToolkitBase/test/runtests.jl index e024cc0200..732d61dfe5 100644 --- a/lib/ModelingToolkitBase/test/runtests.jl +++ b/lib/ModelingToolkitBase/test/runtests.jl @@ -63,6 +63,7 @@ end if GROUP == "All" || GROUP == "Initialization" @safetestset "Guess Propagation" include("guess_propagation.jl") @safetestset "InitializationSystem Test" include("initializationsystem.jl") + @safetestset "Homotopy Operator L0 Trivial" include("homotopy_init.jl") @safetestset "Initial Values Test" include("initial_values.jl") end From e2b82a34730d28972508404ce876d2e191924015 Mon Sep 17 00:00:00 2001 From: lf Date: Thu, 21 May 2026 15:45:27 +0800 Subject: [PATCH 24/24] test: Buildings PressureDrop-style fixture exercises homotopy init Mirrors the structure of Buildings/Fluid/FixedResistances/PressureDrop.mo (from_dp=true branch): nonlinear actual (sqrt-based turbulent law) + linear simplified (nominal-point scaling). PR1 asserts (a) hook fired (no residual homotopy in init equations), (b) init converged, (c) m_flow numerically matches the actual root sqrt(5), not the simplified root m_flow_nominal*dp/dp_nominal = 1.0. Numerical OMC parity deferred to PR2. --- lib/ModelingToolkitBase/test/homotopy_init.jl | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/lib/ModelingToolkitBase/test/homotopy_init.jl b/lib/ModelingToolkitBase/test/homotopy_init.jl index f2f28dbd5a..3f1bb048cb 100644 --- a/lib/ModelingToolkitBase/test/homotopy_init.jl +++ b/lib/ModelingToolkitBase/test/homotopy_init.jl @@ -98,4 +98,62 @@ using Symbolics @test abs(sol.u[1]^2 - 9.0) < 1e-6 # actual equation y^2 = p satisfied @test abs(sol.u[1] - 1.0) > 0.5 # not the simplified root (y = 1) end + + @testset "Buildings PressureDrop fixture (PR1: init doesn't break)" begin + # Modelica Buildings Fluid/FixedResistances/PressureDrop.mo from_dp=true branch: + # m_flow = homotopy( + # actual = basicFlowFunction_dp(dp, k, m_flow_turbulent), + # simplified = m_flow_nominal_pos * dp / dp_nominal_pos) + # PR1 scope = "init runs and converges to the *actual* equation, not + # the simplified one". Numerical OMC parity is deferred to PR2. + # Uses the ODESystem-with-algebraic-constraint pattern (same as the + # Integration testset) so the hook in InitializationProblem fires + # symbolically; NonlinearSystem direct construction does not migrate + # parent algebraic eqs into init system. + using ModelingToolkitBase: System, mtkcompile + using ModelingToolkitBase: t_nounits as t, D_nounits as D + using ModelingToolkitBase: has_homotopy_in_equations, equations + using OrdinaryDiffEqRosenbrock: Rodas5P + using OrdinaryDiffEqNonlinearSolve # needed for DAE init nlsolve + using SciMLBase + + @variables x(t) m_flow(t) + @parameters dp k m_flow_nominal dp_nominal + + # Turbulent sqrt-law actual, linear-nominal simplified. + basic_flow(dp_val, k_val) = sign(dp_val) * sqrt(abs(dp_val)) * k_val + + # `dp` is a parameter (boundary fixed externally); `m_flow` is + # algebraic, pinned by the homotopy constraint. `x` is a dummy + # differential state so the system is a well-posed DAE (mirrors + # the Integration testset's structure: 1 ODE + 1 algebraic eq). + eqs = [ + D(x) ~ -x, + 0 ~ m_flow - homotopy( + basic_flow(dp, k), + m_flow_nominal * dp / dp_nominal, + ), + ] + @named sys = System(eqs, t; guesses = [m_flow => 2.0]) + sys = mtkcompile(sys) + + prob = ODEProblem(sys, + Dict(x => 1.0, dp => 5.0, k => 1.0, + m_flow_nominal => 1.0, dp_nominal => 5.0), + (0.0, 1.0)) + sol = solve(prob, Rodas5P(); abstol = 1e-10, reltol = 1e-10) + + # Hook fired: no residual homotopy in init equations + @test !has_homotopy_in_equations(equations(prob.f.initialization_data.initializeprob.f.sys)) + # Solver succeeded + @test SciMLBase.successful_retcode(sol) + # m_flow at t=0 should solve the *actual* equation: + # m_flow = sign(dp) * sqrt(abs(dp)) * k at dp=5, k=1 + # m_flow = sqrt(5) ≈ 2.236 + # If the rewrite had failed and simplified were active: + # m_flow = m_flow_nominal * dp / dp_nominal = 1.0 * 5/5 = 1.0 + # The sqrt(5) result distinguishes the two outcomes. + m_flow_at_t0 = sol[m_flow][1] + @test abs(m_flow_at_t0 - sqrt(5.0)) < 1e-6 # actual root + end end