WIP: refactor: use ImplicitDiscreteSystem to implement affects in callback systems

Project.toml

@@ -28,6 +28,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+ImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
@@ -40,6 +41,7 @@ NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -156,8 +158,8 @@ julia = "1.9"
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
 BoundaryValueDiffEqAscher = "7227322d-7511-4e07-9247-ad6ff830280e"
+BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"

src/ModelingToolkit.jl

@@ -54,6 +54,7 @@ import Moshi
 using Moshi.Data: @data
 using NonlinearSolve
 import SCCNonlinearSolve
+using ImplicitDiscreteSolve
 using Reexport
 using RecursiveArrayTools
 import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
@@ -154,7 +155,6 @@ include("systems/model_parsing.jl")
 include("systems/connectors.jl")
 include("systems/analysis_points.jl")
 include("systems/imperative_affect.jl")
-include("systems/callbacks.jl")
 include("systems/codegen_utils.jl")
 include("systems/problem_utils.jl")
 include("linearization.jl")
@@ -164,19 +164,20 @@ include("systems/optimization/optimizationsystem.jl")
 include("systems/optimization/modelingtoolkitize.jl")
 
 include("systems/nonlinear/nonlinearsystem.jl")
-include("systems/nonlinear/homotopy_continuation.jl")
+include("systems/discrete_system/discrete_system.jl")
+include("systems/discrete_system/implicit_discrete_system.jl")
+include("systems/callbacks.jl")
+
 include("systems/diffeqs/odesystem.jl")
 include("systems/diffeqs/sdesystem.jl")
 include("systems/diffeqs/abstractodesystem.jl")
+include("systems/nonlinear/homotopy_continuation.jl")
 include("systems/nonlinear/modelingtoolkitize.jl")
 include("systems/nonlinear/initializesystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
 include("systems/diffeqs/basic_transformations.jl")
 
-include("systems/discrete_system/discrete_system.jl")
-include("systems/discrete_system/implicit_discrete_system.jl")
-
 include("systems/jumps/jumpsystem.jl")
 
 include("systems/pde/pdesystem.jl")
@@ -299,6 +300,7 @@ export initialization_equations, guesses, defaults, parameter_dependencies, hier
 export structural_simplify, expand_connections, linearize, linearization_function,
        LinearizationProblem
 export solve
+export Pre
 
 export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function
 export calculate_control_jacobian, generate_control_jacobian

src/linearization.jl

@@ -438,7 +438,6 @@ function linearize_symbolic(sys::AbstractSystem, inputs,
         if !iszero(Bs)
             if !allow_input_derivatives
                 der_inds = findall(vec(any(!iszero, Bs, dims = 1)))
-                @show typeof(der_inds)
                 error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(ModelingToolkit.inputs(sys)[der_inds]). Call `linearize_symbolic` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
             end
             B = [B [zeros(nx, nu); Bs]]

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -96,7 +96,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
     ieqs = Int[]
     filtered_vars = BitSet()
     free_eqs = free_equations(graph, var_sccs, var_eq_matching, varfilter)
-    is_overdetemined = !isempty(free_eqs)
+    is_overdetermined = !isempty(free_eqs)
     for vars in var_sccs
         for var in vars
             if varfilter(var)
@@ -112,7 +112,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
             filtered_vars, isder)
         # If the systems is overdetemined, we cannot assume the free equations
         # will not form algebraic loops with equations in the sccs.
-        if !is_overdetemined
+        if !is_overdetermined
             vargraph.ne = 0
             for var in vars
                 vargraph.matching[var] = unassigned
@@ -121,7 +121,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
         empty!(ieqs)
         empty!(filtered_vars)
     end
-    if is_overdetemined
+    if is_overdetermined
         free_vars = findall(x -> !(x isa Int), var_eq_matching)
         tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, free_eqs,
             BitSet(free_vars), isder)

src/structural_transformation/utils.jl

@@ -218,6 +218,7 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
     all_int_vars = true
     coeffs === nothing || empty!(coeffs)
     empty!(to_rm)
+
     for j in 𝑠neighbors(graph, ieq)
         var = fullvars[j]
         isirreducible(var) && (all_int_vars = false; continue)