Workaround for SciML dropping support of arrays of SVectors

[ranocha]

See for example SciML/OrdinaryDiffEq.jl#2087
This has just happened and lets CI fail for us. Some related PRs with hotfixes are #1788, #1785, #1784. However, we should get a proper fix for it. It only affects some DGMulti solvers.

#1788 is merged. It includes a restriction of DiffEqBase.jl that has to be removed when a proper fix is implemented.

[JoshuaLampert]

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.

[ranocha]

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

[JoshuaLampert]

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

[jlchan]

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...

[JoshuaLampert]

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

[jlchan]

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

[jlchan]

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

[ranocha]

Thanks a lot!

[jlchan]

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.

[ranocha]

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

[ranocha]

Notes from another discussion:

• LoopVectorization.jl supports Julia v1.11
• DGMulti could switch to VectorOfArray or Matrix{SVector}. Both options have some unresolved issues that need to be fixed. It is not clear which one is better.