`maybe_zeros` check does not look inside a division, so tearing divides by a parameter expression that is zero (NaN in initialization)

Start in lib/ModelingToolkitTearing/src/stateselection_interface.jl, at MaybeZeroQueryFn and maybe_zeros_descend, and at how query_maybe_zeros is called from _check_allow_symbolic_parameter. Run the reproducer with and without maybe_zeros to confirm that a BSImpl.Div root returns false. Done means a quotient whose numerator may be zero is flagged, the trim returns uc = [0, 0, 0] with Success when maybe_zeros is declared, and a regression test covers the Div case.

_check_allow_symbolic_parameter relies on query_maybe_zeros to refuse dividing by a coefficient that may be zero. Neither maybe_zeros_descend nor MaybeZeroQueryFn handles BSImpl.Div, so a coefficient of the form num / den is always treated as safe, even when every term of num is a product of maybe_zeros parameters. Tearing then solves for a variable through that coefficient, and the generated code divides by zero at run time. The result is NaN with no warning.

I hit this when trimming a rigid body by initialization. Angular acceleration is written D(w) ~ inv(J) * M, with inv(J) spelled out as adjugate / det so that the ODE stays explicit (the mass-matrix form runs into SciML/ModelingToolkit.jl#5279). The moments M are parameters set to missing, and the trim condition is D(w) ~ 0. When the products of inertia Ixy, Ixz and Iyz are zero, the initialization system solves D(w[1]) = 0 for uc[3] and divides by (Ixy*Iyz + Ixz*Iyy), which is 0. Declaring maybe_zeros = [Ixy, Ixz, Iyz] does not help, because the coefficient is (Ixy*Iyz + Ixz*Iyy) / det(J).

Reproducer

using ModelingToolkit, OrdinaryDiffEqTsit5, SciMLBase, LinearAlgebra
using SimpleNonlinearSolve: SimpleNewtonRaphson
using ModelingToolkit: t_nounits as t, D_nounits as D

# D(w) ~ inv(J) * uc, with inv(J) written out as adjugate / det (so the ODE stays explicit).
# Trim: find the parameters uc such that D(w) = 0 (answer: uc = 0).
function model(; mz = false)
    @parameters Ixx Iyy Izz Ixy Ixz Iyz uc[1:3] = zeros(3)
    @variables w(t)[1:3] = zeros(3)
    J = [Ixx -Ixy -Ixz; -Ixy Iyy -Iyz; -Ixz -Iyz Izz]
    Jinv = [J[2,2]*J[3,3]-J[2,3]*J[3,2]  J[1,3]*J[3,2]-J[1,2]*J[3,3]  J[1,2]*J[2,3]-J[1,3]*J[2,2];
            J[2,3]*J[3,1]-J[2,1]*J[3,3]  J[1,1]*J[3,3]-J[1,3]*J[3,1]  J[1,3]*J[2,1]-J[1,1]*J[2,3];
            J[2,1]*J[3,2]-J[2,2]*J[3,1]  J[1,2]*J[3,1]-J[1,1]*J[3,2]  J[1,1]*J[2,2]-J[1,2]*J[2,1]] ./ det(J)
    kw = mz ? (; maybe_zeros = [Ixy, Ixz, Iyz]) : (;)
    System([D(w[i]) ~ sum(Jinv[i, j] * uc[j] for j in 1:3) for i in 1:3], t; name = :m, kw...)
end

for mz in (false, true)
    s = mtkcompile(model(; mz))
    vals = [s.Ixx => 3.4, s.Iyy => 8.2, s.Izz => 11.4, s.Ixy => 0.0, s.Ixz => 0.0, s.Iyz => 0.0]
    prob = ODEProblem(s, [vals; s.uc => missing], (0.0, 1.0); guesses = [s.uc => ones(3)],
                      initialization_eqs = collect(D(s.w) .~ 0))
    isys = prob.f.initialization_data.initializeprob.f.sys
    torn = filter(eq -> any(isequal(eq.lhs), collect(s.uc)), observed(isys))
    sol = solve(prob, Tsit5(); initializealg = SciMLBase.OverrideInit(nlsolve = SimpleNewtonRaphson()))
    println("maybe_zeros = ", mz ? "[Ixy, Ixz, Iyz]" : "none")
    println("  torn in the initialization system: ", join(string.(torn), "; "))
    println("  retcode = ", sol.retcode, ", uc = ", sol.ps[s.uc])
end

Output

maybe_zeros = none
  torn in the initialization system: uc[3] ~ (<det>*(-(((Ixy*Izz + Ixz*Iyz)*uc[2]) / <det>) - (((Iyy*Izz - (Iyz^2))*uc[1]) / <det>))) / (Ixy*Iyz + Ixz*Iyy)
  retcode = DtNaN, uc = [NaN, NaN, NaN]
maybe_zeros = [Ixy, Ixz, Iyz]
  torn in the initialization system: uc[3] ~ (<det>*(...)) / (Ixy*Iyz + Ixz*Iyy)
  retcode = DtNaN, uc = [NaN, NaN, NaN]

(<det> abbreviates Ixx*(Iyy*Izz - (Iyz^2)) + Ixy*(-Ixy*Izz - Ixz*Iyz) - (Ixy*Iyz + Ixz*Iyy)*Ixz.)

With the same model, all three products non-zero (Ixy = 0.15, Ixz = 0.2, Iyz = -0.1) works, and so does Ixz alone. Which zero pattern fails depends on the pivot that tearing picks: in a larger model (the full 6-DOF body with forces) Ixz alone was the failing case.

Cause

find_eq_solvables! gets the coefficient a = (Ixy*Iyz + Ixz*Iyy) / <det> for uc[3] in D(w[1]) ~ 0. It is not constant, so _check_allow_symbolic_parameter asks query_maybe_zeros. That runs SU.query(MaybeZeroQueryFn(...), a; recurse = maybe_zeros_descend). The root node is a BSImpl.Div: the functor falls through to _ => return false, and maybe_zeros_descend also returns false for Div, so the numerator is never inspected. The ADD rule for the numerator would have returned true, because both of its terms are products containing a maybe_zeros parameter.

Suggested fix

Treat a quotient as possibly zero wherever its numerator may be zero, in MaybeZeroQueryFn (lib/ModelingToolkitTearing/src/stateselection_interface.jl), in the same way #79 added max, min and ifelse:

        BSImpl.Div(; num) => return query_maybe_zeros(qf, num)

I checked this locally with an equivalent override of query_maybe_zeros (if operation(ex) === /, query the numerator) on ModelingToolkitTearing 1.20.7. With maybe_zeros = [Ixy, Ixz, Iyz] the initialization system then pivots on (Ixx*Iyy - Ixy^2) and the trim gives uc = [0, 0, 0], retcode Success. Without maybe_zeros it still fails, as expected.

Related

• A model compiled with non-zero defaults for Ixy, Ixz and Iyz, then run with them set to 0, fails in the same way even with the fix, unless maybe_zeros is declared. Values given only to ODEProblem are not used by discover_maybe_zeros; the initialization system is built at ODEProblem time and inherits maybe_zeros(sys) from the compiled system. I'm not sure whether that is intended. The general case (a pivot that is zero at run time with nothing to flag it) is SciML/ModelingToolkit.jl#5212.
• mtkcompile(sys; allow_parameter = false) does not change the initialization system's choice either. That looks expected, since the initialization system is compiled separately, but there seems to be no way to pass it through.
• SciML/ModelingToolkit.jl#5279: same physical case (an inertia tensor with zero products), different code path (inline_linear_sccs).

Versions: Julia 1.12.5, ModelingToolkit 11.45.4, ModelingToolkitBase 1.77.4, ModelingToolkitTearing 1.20.7, StateSelection 1.11.1, SymbolicUtils 4.51.0. maybe_zeros_descend and MaybeZeroQueryFn on StateSelection.jl main (checked 2026-10-16) are the same as in 1.20.7.