fix: SDE has no tearing state

Author: vyudu

src/systems/diffeqs/sdesystem.jl

@@ -642,20 +642,16 @@ function DiffEqBase.SDEFunction{iip, specialize}(sys::SDESystem, dvs = unknowns(
         _Wfact, _Wfact_t = nothing, nothing
     end
 
+    M = calculate_massmatrix(sys)
     if sparse
         uElType = u0 === nothing ? Float64 : eltype(u0)
-        if jac
-            jac_prototype = similar(calculate_jacobian(sys; sparse), uElType)
-        else
-            jac_prototype = similar(jacobian_sparsity(sys), uElType)
-        end
-        W_prototype = similar(W_sparsity(sys), uElType)
+        jac_prototype = similar(calculate_jacobian(sys; sparse), uElType)
+        W_prototype = similar(jac_prototype .+ M, uElType)
     else
         jac_prototype = nothing
         W_prototype = nothing
     end
 
-    M = calculate_massmatrix(sys)
     _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0', M)
 
     observedfun = ObservedFunctionCache(
@@ -742,15 +738,14 @@ function SDEFunctionExpr{iip}(sys::SDESystem, dvs = unknowns(sys),
         _jac = :nothing
     end
 
-    jac_prototype = if sparse
+    M = calculate_massmatrix(sys)
+    if sparse
         uElType = u0 === nothing ? Float64 : eltype(u0)
-        if jac
-            similar(calculate_jacobian(sys, sparse = sparse), uElType)
-        else
-            similar(jacobian_sparsity(sys), uElType)
-        end
+        jac_prototype = similar(calculate_jacobian(sys; sparse), uElType)
+        W_prototype = similar(jac_prototype .+ M, uElType)
     else
-        nothing
+        jac_prototype = nothing
+        W_prototype = nothing
     end
 
     if Wfact
@@ -763,8 +758,6 @@ function SDEFunctionExpr{iip}(sys::SDESystem, dvs = unknowns(sys),
         _Wfact, _Wfact_t = :nothing, :nothing
     end
 
-    M = calculate_massmatrix(sys)
-
     _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0', M)
 
     ex = quote

test/jacobiansparsity.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, SparseArrays#, OrdinaryDiffEq
+using ModelingToolkit, SparseArrays, OrdinaryDiffEq
 
 N = 3
 xyd_brusselator = range(0, stop = 1, length = N)
@@ -100,6 +100,7 @@ prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = true)
     W_prototype = ModelingToolkit.W_sparsity(pend)
     @test nnz(W_prototype) == nnz(jac_prototype) + 2
 
+    # jac_prototype should be the same as W_prototype
     @test findnz(prob.f.jac_prototype)[1:2] == findnz(W_prototype)[1:2]
 
     u = zeros(5)

test/nonlinearsystem.jl

@@ -30,7 +30,7 @@ eqs = [0 ~ σ * (y - x) * h,
 @test eval(toexpr(ns)) == ns
 test_nlsys_inference("standard", ns, (x, y, z), (σ, ρ, β))
 @test begin
-    f = eval(generate_function(ns, [x, y, z], [σ, ρ, β])[2])
+    f = generate_function(ns, [x, y, z], [σ, ρ, β], expression = Val{false})[2]
     du = [0.0, 0.0, 0.0]
     f(du, [1, 2, 3], [1, 2, 3])
     du ≈ [1, -3, -7]