Add W indices to sparse J to avoid iszero check

Author: hersle

src/systems/codegen.jl

@@ -186,16 +186,14 @@ function calculate_jacobian(sys::System;
     if sparse
         jac = sparsejacobian(rhs, dvs; simplify)
         if get_iv(sys) !== nothing
-            W_s = W_sparsity(sys)
-            (Is, Js, Vs) = findnz(W_s)
-            # Add nonzeros of W as non-structural zeros of the Jacobian (to ensure equal
-            # results for oop and iip Jacobian)
-            for (i, j) in zip(Is, Js)
-                iszero(jac[i, j]) && begin
-                    jac[i, j] = 1
-                    jac[i, j] = 0
-                end
-            end
+            # Add nonzeros of W as non-structural zeros of the Jacobian
+            # (to ensure equal results for oop and iip Jacobian)
+            JIs, JJs, JVs = findnz(jac)
+            WIs, WJs, _ = findnz(W_sparsity(sys))
+            append!(JIs, WIs) # explicitly put all W's indices also in J,
+            append!(JJs, WJs) # even if it duplicates some indices
+            append!(JVs, zeros(eltype(JVs), length(WIs))) # add zero
+            jac = SparseArrays.sparse(JIs, JJs, JVs) # values at duplicate indices are summed; not overwritten
         end
     else
         jac = jacobian(rhs, dvs; simplify)