DAE mass matrix with ArrayPartition

[gbene]

Question❓

Hi again!

I am trying to solve a DAE having as u₀ three NxN Arrays. The code that I am using to solve my DAE is a bit long for a MWE (has many parameters and different things going on) so I will write here a pseudo-problem that emulates the gist of it and reproduces my problem. The DAE makes 0 sense but again it is just to show the problem that I've encountered:

"""
Solve the following DAE (dummy example)

ẋ = z
ẏ = z*y/L
0 = x*L-z


with x,y and z as NxN arrays
"""

using OrdinaryDiffEq, RecursiveArrayTools, LinearAlgebra

function MatrixDAE!(du, u, L, t)

      x, y, z = u.x
  
      du.x[1] .= z
      du.x[2] .= (z*y)./L
      du.x[3] .= x*L-z
  
      nothing
end

N = 5

x = rand(N,N)
y = rand(N,N)
z = rand(N,N)

L=0.5
 
u₀ = ArrayPartition(x, y, z)

end_time = 1.0
tₛ = (0.0, end_time)


func = ODEFunction(MatrixDAE!, mass_matrix = Diagonal([1, 1, 0])) 

problem = ODEProblem(func, u₀, tₛ, L)

solver = Rodas4()

#4. Solve

solution = solve(problem, solver,save_everystep=false)

As per documentation I am using ArrayPartition to group the three arrays and give that as u₀. However by doing this, solve throws the following error:

Mass matrix size is incompatible with initial condition
sizing. The mass matrix must represent the `vec`
form of the initial condition `u0`, i.e.
`size(mm,1) == size(mm,2) == length(u)`

size(prob.f.mass_matrix,1): 3
length(u0): 75

This error makes sense because length(u₀) is 75 (3x5x5). However I thought that in this case the checked length should have been length(u₀.x) (that is 3). Is this an oversight or should I set the mass matrix in a different way (how)?

I am (again) not sure if this is an oversight or it is just a limitation of the solvers for DAEs. If it is the former I am sorry if I miss labelled the issue!

Best!

[ChrisRackauckas]

The mass matrix isn't able to be treated in block form, so the call that works is:

func = ODEFunction(MatrixDAE!, mass_matrix = Diagonal([ones(length(x)+length(y)); zeros(length(z))])) 

which is treating the whole vector. In the future, making it work as block diagonal would be cool, but we don't have that right now.