ArrayPartition + Implicit + OOP Fails

[ChrisRackauckas]

using DiffEqOperators, OrdinaryDiffEq, Plots

ω = 1; k = 1; c = (ω/k)^2
number_of_spatial_cells = 100
dx = 1/(number_of_spatial_cells+1)
ord_spacial_deriv = 2
u_0 = 0.2.*sin.(range(0, 2*pi, length=number_of_spatial_cells))
u_dot_0 = zeros(number_of_spatial_cells)
A_x = CenteredDifference{1}(2, 2, dx, number_of_spatial_cells)
Q = Dirichlet0BC(Float64)

# Inplace works
u_dot_dot_iip(out, du, u, p, t) = out .= c*A_x*Q*u
prob = SecondOrderODEProblem(u_dot_dot_iip, u_dot_0, u_0, (0.0, 1.0))
sol = solve(prob, Rosenbrock23())

# Out of place fails
u_dot_dot(du, u, p, t) = c*A_x*Q*u
prob = SecondOrderODEProblem(u_dot_dot, u_dot_0, u_0, (0.0, 1.0))
sol = solve(prob, Rosenbrock23())

[tmistele]

I think I hit the same error but ended up with a more minimal reproducer:

using OrdinaryDiffEq

calc_ddu(du,u,p,t) = [-u[1]]
calc_ddu_inplace(ddu,du,u,p,t) = ddu[1] = -u[1]

# These three work
solve(SecondOrderODEProblem(calc_ddu_inplace, [0.0], [1.0], (0.0, 10.0)), Tsit5())
solve(SecondOrderODEProblem(calc_ddu_inplace, [0.0], [1.0], (0.0, 10.0)), Rosenbrock23())
solve(SecondOrderODEProblem(calc_ddu, [0.0], [1.0], (0.0, 10.0)), Tsit5())

# This one dies with `type Array has no field x`
solve(SecondOrderODEProblem(calc_ddu, [0.0], [1.0], (0.0, 10.0)), Rosenbrock23())

@ChrisRackauckas can you confirm this has the same root cause? If not I can open a separate issue.

fwiw, I used OrdinaryDiffEq.jl v6.30.1

[gerlero]

Same issue as far as I can tell

[ChrisRackauckas]

Yes same issue. The workaround is to just use in-place for now.

[gerlero]

Does rewriting the problem as a first-order one also count as a workaround, or does DifferentialEquations.jl do anything different when it sees a SecondOrderODE?

I'm asking because, for small problems, I've found that rewriting as a generic ODEProblem while keeping the functions out-of-place is faster.

[ChrisRackauckas]

Does rewriting the problem as a first-order one also count as a workaround, or does DifferentialEquations.jl do anything different when it sees a SecondOrderODE?

It does a little bit. I mean, there are some solvers which are only defined for second order ODEs (Nystrom methods) and they are pretty efficient. But that's a really special case, and honestly they need some optimizations (not the implementation but the methods themselves, this needs some research: I don't like the current error estimators). So writing it as a first order one yourself is quite fine.

That said, this could get fixed with some linear algebra overloads. It's not too hard but I just haven't found the time for it yet.

[gerlero]

Thank you!