|
| 1 | +function tearing_sub(expr, dict, s) |
| 2 | + expr = ModelingToolkit.fixpoint_sub(expr, dict) |
| 3 | + s ? simplify(expr) : expr |
| 4 | +end |
| 5 | + |
| 6 | +function tearing_reassemble(sys; simplify=false) |
| 7 | + s = structure(sys) |
| 8 | + @unpack fullvars, partitions, var_eq_matching, graph, scc = s |
| 9 | + eqs = equations(sys) |
| 10 | + |
| 11 | + ### extract partition information |
| 12 | + rhss = [] |
| 13 | + solvars = [] |
| 14 | + ns, nd = nsrcs(graph), ndsts(graph) |
| 15 | + active_eqs = trues(ns) |
| 16 | + active_vars = trues(nd) |
| 17 | + rvar2reqs = Vector{Vector{Int}}(undef, nd) |
| 18 | + for (ith_scc, partition) in enumerate(partitions) |
| 19 | + @unpack e_solved, v_solved, e_residual, v_residual = partition |
| 20 | + for ii in eachindex(e_solved) |
| 21 | + ieq = e_solved[ii]; ns -= 1 |
| 22 | + iv = v_solved[ii]; nd -= 1 |
| 23 | + rvar2reqs[iv] = e_solved |
| 24 | + |
| 25 | + active_eqs[ieq] = false |
| 26 | + active_vars[iv] = false |
| 27 | + |
| 28 | + eq = eqs[ieq] |
| 29 | + var = fullvars[iv] |
| 30 | + rhs = value(solve_for(eq, var; simplify=simplify, check=false)) |
| 31 | + # if we don't simplify the rhs and the `eq` is not solved properly |
| 32 | + (!simplify && occursin(rhs, var)) && (rhs = SymbolicUtils.polynormalize(rhs)) |
| 33 | + # Since we know `eq` is linear wrt `var`, so the round off must be a |
| 34 | + # linear term. We can correct the round off error by a linear |
| 35 | + # correction. |
| 36 | + rhs -= expand_derivatives(Differential(var)(rhs))*var |
| 37 | + @assert !(var in vars(rhs)) """ |
| 38 | + When solving |
| 39 | + $eq |
| 40 | + $var remainded in |
| 41 | + $rhs. |
| 42 | + """ |
| 43 | + push!(rhss, rhs) |
| 44 | + push!(solvars, var) |
| 45 | + end |
| 46 | + # DEBUG: |
| 47 | + #@show ith_scc solvars .~ rhss |
| 48 | + #Main._nlsys[] = eqs[e_solved], fullvars[v_solved] |
| 49 | + #ModelingToolkit.topsort_equations(solvars .~ rhss, fullvars) |
| 50 | + #empty!(solvars); empty!(rhss) |
| 51 | + end |
| 52 | + |
| 53 | + ### update SCC |
| 54 | + eq_reidx = Vector{Int}(undef, nsrcs(graph)) |
| 55 | + idx = 0 |
| 56 | + for (i, active) in enumerate(active_eqs) |
| 57 | + eq_reidx[i] = active ? (idx += 1) : -1 |
| 58 | + end |
| 59 | + |
| 60 | + rmidxs = Int[] |
| 61 | + newscc = Vector{Int}[]; sizehint!(newscc, length(scc)) |
| 62 | + for component′ in newscc |
| 63 | + component = copy(component′) |
| 64 | + for (idx, eq) in enumerate(component) |
| 65 | + if active_eqs[eq] |
| 66 | + component[idx] = eq_reidx[eq] |
| 67 | + else |
| 68 | + push!(rmidxs, idx) |
| 69 | + end |
| 70 | + end |
| 71 | + push!(newscc, component) |
| 72 | + deleteat!(component, rmidxs) |
| 73 | + empty!(rmidxs) |
| 74 | + end |
| 75 | + |
| 76 | + ### update graph |
| 77 | + var_reidx = Vector{Int}(undef, ndsts(graph)) |
| 78 | + idx = 0 |
| 79 | + for (i, active) in enumerate(active_vars) |
| 80 | + var_reidx[i] = active ? (idx += 1) : -1 |
| 81 | + end |
| 82 | + |
| 83 | + newgraph = BipartiteGraph(ns, nd, Val(false)) |
| 84 | + |
| 85 | + |
| 86 | + ### update equations |
| 87 | + odestats = [] |
| 88 | + for idx in eachindex(fullvars); isdervar(s, idx) && continue |
| 89 | + push!(odestats, fullvars[idx]) |
| 90 | + end |
| 91 | + newstates = setdiff(odestats, solvars) |
| 92 | + varidxmap = Dict(newstates .=> 1:length(newstates)) |
| 93 | + neweqs = Vector{Equation}(undef, ns) |
| 94 | + newalgeqs = falses(ns) |
| 95 | + |
| 96 | + dict = Dict(value.(solvars) .=> value.(rhss)) |
| 97 | + |
| 98 | + visited = falses(ndsts(graph)) |
| 99 | + for ieq in Iterators.flatten(scc); active_eqs[ieq] || continue |
| 100 | + eq = eqs[ieq] |
| 101 | + ridx = eq_reidx[ieq] |
| 102 | + |
| 103 | + fill!(visited, false) |
| 104 | + compact_graph!(newgraph, graph, visited, ieq, ridx, rvar2reqs, var_reidx, active_vars) |
| 105 | + |
| 106 | + if isdiffeq(eq) |
| 107 | + neweqs[ridx] = eq.lhs ~ tearing_sub(eq.rhs, dict, simplify) |
| 108 | + else |
| 109 | + newalgeqs[ridx] = true |
| 110 | + if !(eq.lhs isa Number && eq.lhs != 0) |
| 111 | + eq = 0 ~ eq.rhs - eq.lhs |
| 112 | + end |
| 113 | + rhs = tearing_sub(eq.rhs, dict, simplify) |
| 114 | + if rhs isa Symbolic |
| 115 | + neweqs[ridx] = 0 ~ rhs |
| 116 | + else # a number |
| 117 | + if abs(rhs) > 100eps(float(rhs)) |
| 118 | + @warn "The equation $eq is not consistent. It simplifed to 0 == $rhs." |
| 119 | + end |
| 120 | + neweqs[ridx] = 0 ~ fullvars[invview(var_eq_matching)[ieq]] |
| 121 | + end |
| 122 | + end |
| 123 | + end |
| 124 | + |
| 125 | + ### update partitions |
| 126 | + newpartitions = similar(partitions, 0) |
| 127 | + emptyintvec = Int[] |
| 128 | + for (ii, partition) in enumerate(partitions) |
| 129 | + @unpack e_residual, v_residual = partition |
| 130 | + isempty(v_residual) && continue |
| 131 | + new_e_residual = similar(e_residual) |
| 132 | + new_v_residual = similar(v_residual) |
| 133 | + for ii in eachindex(e_residual) |
| 134 | + new_e_residual[ii] = eq_reidx[ e_residual[ii]] |
| 135 | + new_v_residual[ii] = var_reidx[v_residual[ii]] |
| 136 | + end |
| 137 | + # `emptyintvec` is aliased to save memory |
| 138 | + # We need them for type stability |
| 139 | + newpart = SystemPartition(emptyintvec, emptyintvec, new_e_residual, new_v_residual) |
| 140 | + push!(newpartitions, newpart) |
| 141 | + end |
| 142 | + |
| 143 | + obseqs = solvars .~ rhss |
| 144 | + |
| 145 | + @set! s.graph = newgraph |
| 146 | + @set! s.scc = newscc |
| 147 | + @set! s.fullvars = fullvars[active_vars] |
| 148 | + @set! s.vartype = s.vartype[active_vars] |
| 149 | + @set! s.partitions = newpartitions |
| 150 | + @set! s.algeqs = newalgeqs |
| 151 | + |
| 152 | + @set! sys.structure = s |
| 153 | + @set! sys.eqs = neweqs |
| 154 | + @set! sys.states = newstates |
| 155 | + @set! sys.observed = [observed(sys); obseqs] |
| 156 | + return sys |
| 157 | +end |
| 158 | + |
| 159 | +# removes the solved equations and variables |
| 160 | +function compact_graph!(newgraph, graph, visited, eq, req, rvar2reqs, var_reidx, active_vars) |
| 161 | + for ivar in 𝑠neighbors(graph, eq) |
| 162 | + # Note that we need to check `ii` against the rhs states to make |
| 163 | + # sure we don't run in circles. |
| 164 | + visited[ivar] && continue |
| 165 | + visited[ivar] = true |
| 166 | + |
| 167 | + if active_vars[ivar] |
| 168 | + add_edge!(newgraph, req, var_reidx[ivar]) |
| 169 | + else |
| 170 | + # If a state is reduced, then we go to the rhs and collect |
| 171 | + # its states. |
| 172 | + for ieq in rvar2reqs[ivar] |
| 173 | + compact_graph!(newgraph, graph, visited, ieq, req, rvar2reqs, var_reidx, active_vars) |
| 174 | + end |
| 175 | + end |
| 176 | + end |
| 177 | + return nothing |
| 178 | +end |
| 179 | + |
| 180 | +""" |
| 181 | + tearing(sys; simplify=false) |
| 182 | +
|
| 183 | +Tear the nonlinear equations in system. When `simplify=true`, we simplify the |
| 184 | +new residual residual equations after tearing. End users are encouraged to call [`structural_simplify`](@ref) |
| 185 | +instead, which calls this function internally. |
| 186 | +""" |
| 187 | +tearing(sys; simplify=false) = tearing_reassemble(tear_graph(algebraic_equations_scc(sys)); simplify=simplify) |
0 commit comments