feat: improve error message when system contains array equations

[AayushSabharwal]

Close #2597

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[ChrisRackauckas]

Why would this be required? The codegen doesn't require it if it's a true ODE.

[AayushSabharwal]

It doesn't really work very reliably.

@variables x(t)[1:2]
@named sys = ODESystem(D(x) ~ x, t)
sys = complete(sys)
prob = ODEProblem(sys, ...)

This errors in calculate_massmatrix because that can't handle array equations and needs them scalarized. If I hack around that and comment out check_eqs_u0 this errors in generate_function because check_lhs finds that the variable inside the derivative D(x) isn't in the list of unknowns.

I could make these cases work, but all of the fixes boil down to scalarizing the equations and/or unknowns.

[AayushSabharwal]

And if there are multiple array unknowns or a mix of array/scalars then codegen also won't work without scalarizing the equations.

[ChrisRackauckas]

calculate_massmatrix needs to semantically scalarize and show it's constant.

And if there are multiple array unknowns or a mix of array/scalars then codegen also won't work without scalarizing the equations.

Is that really the case? Symbolics.jl handled this awhile back.

[ChrisRackauckas]

I see no issue with a temporary error message that told people as a temporary workaround this is required, though I think we should look at some of the codegen here since Symbolics should already handle these cases.

[AayushSabharwal]

The problem with the codegen is that this system:

julia> @named sys = ODESystem([D(x) ~ x, D(y) ~ t], t, [x..., y])

Generates this code:

RuntimeGeneratedFunction(#=in ModelingToolkit=#, #=using ModelingToolkit=#, :((ˍ₋out, ˍ₋arg1, t)->begin
          #= /Users/aayush/.julia/packages/SymbolicUtils/ij6YM/src/code.jl:385 =#
          #= /Users/aayush/.julia/packages/SymbolicUtils/ij6YM/src/code.jl:386 =#
          #= /Users/aayush/.julia/packages/SymbolicUtils/ij6YM/src/code.jl:387 =#
          begin
              begin
                  var"x(t)" = view(ˍ₋arg1, 1:2)
                  begin
                      begin
                          #= /Users/aayush/.julia/packages/Symbolics/XnDVB/src/build_function.jl:546 =#
                          #= /Users/aayush/.julia/packages/SymbolicUtils/ij6YM/src/code.jl:434 =# @inbounds begin
                                  #= /Users/aayush/.julia/packages/SymbolicUtils/ij6YM/src/code.jl:430 =#
                                  ˍ₋out[1] = var"x(t)"
                                  ˍ₋out[2] = t
                                  #= /Users/aayush/.julia/packages/SymbolicUtils/ij6YM/src/code.jl:432 =#
                                  nothing
                              end
                      end
                  end
              end
          end
      end))

Which assumes the first element of u is a Vector and the second a float.

Which again points to the fact that the system needs to be scalarized

[AayushSabharwal]

Also without structural_simplify it's up to the user to ensure the equations and unknowns are in the correct order.

[ChrisRackauckas]

Okay we should fix that, it's fixable. For now, mention that a temporary workaround is to scalarize or structural simplify.

[ChrisRackauckas]

@AayushSabharwal some test greening https://github.com/SciML/ModelingToolkit.jl/actions/runs/11185499006/job/31098502328#step:6:2267 ?