Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md with its PR number.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

Something in SciML broke again. Seems like a problem with a SplitODEFunction and implicit time stepping. My best guess would be that it is related to SciML/DiffEqBase.jl#1188 or SciML/SciMLBase.jl#1106 because they were also about split problems and differentiation and were released today. @ChrisRackauckas have you seen this kind of error lately or have an idea where it comes from (probably within the last few hours because CI was successful five hours ago)?

Sorry, thanks for pointing that out. Should be fixed now.

Though...

Seems like a problem with a SplitODEFunction and implicit time stepping.

How were you doing that before? It has been a known issue that doing this was not autodiff compatible since it was built. The PRs referenced make it so that it works with autodiff. Were you testing this just with finite differencing or IMEX? Since before the fallback cache wouldn't resize for autodiff automatically.

Thanks! We used finite differences, but good to know that AD should work now👍

Change to forwarddiff now!

The main thing caught here is we need to better test resizing in the SciML interface. Resizing has tests in OrdinaryDiffEq but isn't actually a listed part of the interface, and there's no trait stuff for "is resizable" and what that all means. We should flesh out that interface a bit and add some earlier testing that resize! either works or passes through.

Should be fixed now.

Sorry for bothering again @ChrisRackauckas. Can you point me to a PR that fixed this issue? We still have the same issue in CI after rerunning.

SciML/PreallocationTools.jl#135

resize! is defined. I guess the issue is that u0 is not something that's compatible with DiffCache... For now the right thing might be to restrict this to only when typeof(u0) <: AbstractArray{<:Number}

I did a hotfix to SciMLBase that will make u0 <: VectorOfArray not take this path. Not the greatest hotfix, but it'll get you something today. This also means you cannot ForwardDiff still 😢. Try the new SciMLBase patch.

I think the core of your problems might be that VectorOfArray in its current deprecated form might not be the best type to represent what you need. It might be good to get on a call and figure out exactly what you need there, and if the change to linear indexing will do it. I think with SciML/RecursiveArrayTools.jl#454 and SciML/RecursiveArrayTools.jl#453 it seems you need something that will make ragged "just work" in any solver, though ragged always has some odd issues if you also try to allow matrix indexing. In Julia's abstract arrays you cannot have size be like (4, (1,4,3)), but then certain operations act on size, and so then they assume all vectors in there are the same size. I think we need to then consider what would work best there for your use case:

1. Never assume this has a matrix structure and always assume it's an AbstractVector? If you go with 1, then DiffCache support is easy BTW
2. Assume it has some multi-dimensional stuff but is not an AbstractArray, but then that interface should get detailed and solidified.
3. Assume it's an abstract array, where size uses the maximum array size in there, and if you index beyond it gives... some sentinel like NaN?

I assume (1) might be what you want, which is where the issues are actually all stemming from. And this may be my issue for telling @ranocha this is the path we've got to go 😅 . I still think we need to make sure the u0 satisfies some interface, and in modern SciML we have now gotten much better at documenting and testing that interface:

https://docs.sciml.ai/SciMLBase/stable/interfaces/Array_and_Number/

(though it's not complete, we need to define things around resizability etc.)

But one of the things we have also learned from this exercise of getting better at interfaces is that AbstractVectorOfArray will need to change.

https://docs.sciml.ai/RecursiveArrayTools/stable/AbstractVectorOfArrayInterface/#RecursiveArrayTools.AbstractVectorOfArray

I'm not sure if the linear index fix is enough to make this the right type for you though, since it would be more like (3). Right now it's more like (2) but the ragged interface is not fully complete.

Thanks again. This fixed the tests. Just as a small clarification: SciML/RecursiveArrayTools.jl#454 and SciML/RecursiveArrayTools.jl#453 are not required for Trixi.jl right now, but these issues arise from another of my projects.

Yes, that should be possible.

Great! Could you please add such a test (that fails with the previous versions), add it here, and ping me for a review afterward?

Yes.

I added some tests also testing properties of the spectra, which failed with the previous two implementations. This is ready for a review, @ranocha.

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.71%. Comparing base (d9236b6) to head (c0fb7bb).
⚠️ Report is 1 commits behind head on main.

@@            Coverage Diff             @@
##             main    #2522      +/-   ##
==========================================
+ Coverage   95.20%   96.71%   +1.51%     
==========================================
  Files         512      512              
  Lines       42446    42448       +2     
==========================================
+ Hits        40410    41051     +641     
+ Misses       2036     1397     -639     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@ChrisRackauckas: Thanks a lot for your help! I cannot comment on the needs of @JoshuaLampert in the other project he mentioned that triggered SciML/RecursiveArrayTools.jl#454 and SciML/RecursiveArrayTools.jl#453, but here are some thoughts on Trixi.jl:

In Trixi.jl, we have a few solvers with different memory layouts. One class uses plain arrays of numbers (like Float64). To make these arrays resize!able in DiffEq, we use Vectors that we wrap in multi-dimensional arrays in Trixi.jl to access the data like we need. This part of the code should be fine since we do not need any AbstractVectorOfArray handling.

Another solver class used multi-dimensional arrays of SVectors (all with the same size). We wrapped them in a VectorOfArray to make it work with DiffEq. Inside of Trixi.jl, we extract the underlying multi-dimensional array of SVectors to work with it.

Based on your list in #2522 (comment), option 1 is fine with us as long as we can extract underlying data or view them in the multi-dimensional way we do right now.