@jlchan mentioned some progress in SciML/RecursiveArrayTools.jl#359 and SciML/RecursiveArrayTools.jl#357

from trixi.jl.

Thanks a lot @jlchan for working on this! No worries, take your time. Get well soon!

from trixi.jl.

Is there any update regarding this? This is getting quite annoying because the restriction to DiffEqBase.jl prevents us from getting newer versions of OrdinaryDiffEq.jl.

from trixi.jl.

Thanks for the references. Good to hear there is some progress!

from trixi.jl.

Yes - we should be able to wrap a array of arrays u with VectorOfArray(u) to get around OrdinaryDiffEq.jl dropping support for arrays of arrays. We can unwrap u via parent(VectorOfArray(u)), which should be possible inside wrap_array.

I need to get a PR together for this, but I've got a bit of a backlog, esp with COVID this week and then travel next week. Apologies for the delay - I hope things will clear up around the second week of March...

from trixi.jl.

There's one more issue that needs to be resolved upstream related to broadcasting over VectorOfArray with multidimensional parent arrays SciML/RecursiveArrayTools.jl#373

from trixi.jl.

I've started working on having DGMulti solvers use VectorOfArray. Some more upstream issues: SciML/RecursiveArrayTools.jl#378

from trixi.jl.

Thanks a lot!

from trixi.jl.

I'm consider another alternative since we're running into upstream issues with RecursiveArrayTools.jl.

Proposal: switch DGMulti to Matrix{<:SVector} storage as discussed previously in #1240.

What DGMulti currently uses

DGMulti assumes that u[index, element] is an SVector, so it uses two array of SVector:

1. Matrix{<:SVector}, since the data layout provides efficient memory accesses
2. StructArray{<:SVector{nvariables}, 2}, since the memory layout (lazy zipping of nvariables Matrix arrays) enabled fast matrix-vector products using Octavian.jl (since simplices tend to have larger dense matrix operations).

What is changing in Julia v1.11

Octavian.jl is being deprecated in Julia v1.11 (see #1906). On an unrelated note, OrdinaryDiffEq.jl no longer supports Array{<:SVector}, so ideally DGMulti would use a flat (or easily flattened) array format.

What I'm proposing

I don't think there is much reason to support StructArray{<:SVector} for DGMulti solutions anymore since Octavian is being deprecated and mul! with Matrix{<:SVector} is about as fast (or faster!) than mul! using Matrix{Float64} now (with the exception of an odd allocation for the 5-arg mul!).

Instead, I think we can do wrap_array with the SVector reinterpret trick, or if it is inefficient, maybe we could use wrap_array to do an unsafe conversion.

I'll start a PR for #1906 first to a staging branch.

from trixi.jl.

Thanks! This sounds good to me - the wrap_array infrastructure we use in Trixi.jl should help with that. Please let me know if I can help you or should review something

from trixi.jl.