Use `VectorOfArray` in `wrap_array` for `DGMulti` solvers #2150

jlchan · 2024-11-07T15:29:42Z

@huiyuxie FYI this PR goes towards addressing #1789.

d

github-actions · 2024-11-07T15:30:00Z

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.

codecov · 2024-11-07T16:12:17Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 95.44%. Comparing base (91eaaf6) to head (9b21e12).
Report is 2 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2150      +/-   ##
==========================================
- Coverage   96.36%   95.44%   -0.92%     
==========================================
  Files         480      480              
  Lines       38028    38025       -3     
==========================================
- Hits        36645    36292     -353     
- Misses       1383     1733     +350

Flag	Coverage Δ
unittests	`95.44% <ø> (-0.92%)`	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

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

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

jlchan · 2024-11-07T23:09:17Z

@DanielDoehring if I remember correctly, you implemented @allocated tests, right?

I am running into some @allocated failures (for example https://github.com/trixi-framework/Trixi.jl/actions/runs/11732466598/job/32684751742?pr=2150)

elixir_euler_weakform.jl (SBP, EC): Test Failed at /home/runner/work/Trixi.jl/Trixi.jl/test/test_threaded.jl:388
  Expression: #= /home/runner/work/Trixi.jl/Trixi.jl/test/test_threaded.jl:388 =# @allocated(Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
   Evaluated: 118704 < 5000

However, at least locally, this appears to only happen the first time Trixi.rhs! is run. I tried running Trixi.rhs! to avoid these allocations in CI, but it doesn't seem to work.

Did you run into this issue before?

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Update test/test_dgmulti_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

…ework/Trixi.jl into jc/wrap_VectorOfArray

DanielDoehring · 2024-11-08T07:36:13Z

However, at least locally, this appears to only happen the first time Trixi.rhs! is run. I tried running Trixi.rhs! to avoid these allocations in CI, but it doesn't seem to work.

Did you run into this issue before?

Hm, I have not encountered this. Can you maybe track down the allocations using the SummaryCallback ? I guess there must be some type instability. But I have no idea why the behaviour after executing rhs! once is different depending on the machine.

jlchan · 2024-11-08T21:07:48Z

Hm, I have not encountered this. Can you maybe track down the allocations using the SummaryCallback ? I guess there must be some type instability. But I have no idea why the behaviour after executing rhs! once is different depending on the machine.

There doesn't seem to be a type instability. There's a weird GC call stack that @profview_allocs shows, I have no idea what's going on. I'll bring it up at the next Trixi meeting.

DanielDoehring · 2024-11-09T12:15:00Z

There doesn't seem to be a type instability. There's a weird GC call stack that @profview_allocs shows, I have no idea what's going on. I'll bring it up at the next Trixi meeting.

Oh wow, that sounds strange. So the allocations are only present in the test "environment"?

jlchan · 2024-11-10T16:05:22Z

There doesn't seem to be a type instability. There's a weird GC call stack that @profview_allocs shows, I have no idea what's going on. I'll bring it up at the next Trixi meeting.

Oh wow, that sounds strange. So the allocations are only present in the test "environment"?

Not exactly - on my machine, the allocations disappear after running rhs! once. In the test environment (or on CI, I can't tell) this isn't the case.

Here's the output from Profile.Allocs.@profile sample_rate=0.001 Trixi.rhs!(du, ode.u0, semi, 0.0) printed to a txt file: profile_allocs.txt.

Project.toml

src/auxiliary/precompile.jl

test/test_dgmulti_1d.jl

huiyuxie · 2024-11-11T03:12:23Z

https://github.com/JuliaLang/julia/blob/3318941e585db632423366b8b703ea55a6ba8421/base/timing.jl#L479-L489
If the allocation is less in the second call compared to the first during a single session, the GC might have been triggered between these two gc_bytes calls. But I also suspect it could be caused by type instability introduced by RecursiveArrayTools.jl. Why do you think it doesn’t seem to be a type instability issue?

huiyuxie · 2024-11-11T03:12:48Z

Here's the output from Profile.Allocs.@Profile sample_rate=0.001 Trixi.rhs!(du, ode.u0, semi, 0.0) printed to a txt file: profile_allocs.txt.

The text file only shows half of the information on each line. Could you provide it in a different way?

jlchan · 2024-11-11T03:37:27Z

Why do you think it doesn’t seem to be a type instability issue?

I checked @code_warntype on the first call; nothing is red and I don't see an Any instance anywhere. If you have any suggestions for where else to check, please let me know.

The text file only shows half of the information on each line. Could you provide it in a different way?

Thanks - how is this?
profile_allocs.txt

huiyuxie · 2024-11-11T03:41:38Z

Same issue - how about telling me how to reproduce it? Use your latest commit and then run Profile.Allocs.@profile sample_rate=0.001 Trixi.rhs!(du, ode.u0, semi, 0.0)?

…ework/Trixi.jl into jc/wrap_VectorOfArray

test/test_visualization.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

test/test_dgmulti_3d.jl

huiyuxie · 2024-11-15T02:55:19Z

src/auxiliary/precompile.jl

+    # @assert Base.precompile(Tuple{DiscreteCallback{typeof(Trixi.summary_callback),
+    #                                                typeof(Trixi.summary_callback),
+    #                                                typeof(Trixi.initialize_summary_callback),
+    #                                                typeof(SciMLBase.FINALIZE_DEFAULT)}})


Will the CI fail if you uncomment these lines?

Tried in fdf1ff5; lets see

huiyuxie · 2024-11-15T03:25:58Z

One thing that might be wired for me is that I tried the allocation profiling (use your MWE) for both main branch and this PR - but found that the memory allocation bytes are similar

Total bytes: 199336 - from main
Total bytes: 191504 - from this PR

Can you confirm this on your side @jlchan (no rush)? Or is there anything we're doing wrong...

jlchan · 2024-11-15T03:29:08Z

One thing that might be wired for me is that I tried the allocation profiling (use your MWE) for both main branch and this PR - but found that the memory allocation bytes are similar

Total bytes: 199336 - from main Total bytes: 191504 - from this PR

Can you confirm this on your side @jlchan (no rush)? Or is there anything we're doing wrong...

Thanks for checking; I can confirm tomorrow.

If you run Trixi.rhs! a second time, do the total number of allocations decrease?

jlchan · 2024-11-15T03:49:30Z

@huiyuxie I see something very weird. If I print out typeof(du) in this function, then for the failing elixir I see typeof(du) = Matrix{SVector{...}}.

However, for the passing elixir I see typeof(du) = StructArray{SVector{. However, u is the same type for both...

EDIT: aha, it is the time-stepper. If I use CarpenterKennedy with a fixed time-step size, it works. If I use SSPRK43 with adaptivity turned on, it fails. In general, turning adaptivity off fixes the failing CI.

huiyuxie · 2024-11-15T03:59:36Z

But if you said that u is the same type for both, I highly suspect there is a type instability problem with this package...

jlchan · 2024-11-15T04:15:42Z

But if you said that u is the same type for both, I highly suspect there is a type instability problem with this package...

No, I think it's a weird subtlety with how OrdinaryDiffEq.jl initializes the du container for specifically the SSPRK43 time-stepper. Here's a MWE:

using Trixi, OrdinaryDiffEq

dg = DGMulti(polydeg = 3, element_type = Line(),
             surface_integral = SurfaceIntegralWeakForm(FluxLaxFriedrichs()),
             volume_integral = VolumeIntegralWeakForm())

equations = CompressibleEulerEquations1D(1.4)
initial_condition = initial_condition_constant

mesh = DGMultiMesh(dg, (2, ), periodicity = true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
ode = semidiscretize(semi, (0.0, .01))

integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 0.01) # works
integrator = init(ode, Tsit5()) # works
integrator = init(ode, RDPK3SpFSAL49()) # works
integrator = init(ode, SSPRK43(), dt = 0.01) # works

integrator = init(ode, SSPRK43()) # fails?!?

The error for the last line isn't very helpful, but if I run @show typeof(du) in Trixi.rhs!, I get typeof(du) = Matrix{SVector{3, Float64}}. It should be StructArray{SVector{3, Float64}, ...} instead, so I think OrdinaryDiffEq.jl is initializing the integrator.k field (that gets passed into Trixi.rhs! as du) incorrectly.

@ranocha - you know the SSPRK time-steppers in OrdinaryDiffEq.jl better than I do. Any idea what might be causing this?

ranocha · 2024-11-20T19:19:30Z

Which versions of the packages are you using?

julia> using Pkg; Pkg.activate(temp = true); Pkg.add(["OrdinaryDiffEq", "StaticArrays", "StructArrays"])
[...]
⌃ [1dea7af3] + OrdinaryDiffEq v6.89.0
  [90137ffa] + StaticArrays v1.9.8
⌃ [09ab397b] + StructArrays v0.6.18
[...]

julia> using OrdinaryDiffEq, StructArrays, StaticArrays

julia> u = StructArray([SVector(1.0, 1.0), SVector(2.0, 2.0)])
2-element StructArray(::Vector{Float64}, ::Vector{Float64}) with eltype SVector{2, Float64}:
 [1.0, 1.0]
 [2.0, 2.0]

julia> ode = ODEProblem(u, (0.0, 1.0)) do du, u, p, t
           @show typeof(du)
           du .= u
           return nothing
       end

julia> solve(ode, Tsit5())
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
[...]

julia> solve(ode, SSPRK43())
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
[...]

I can't see the difference here that you describe 🤷

jlchan · 2024-11-20T20:47:20Z

Which versions of the packages are you using?

julia> using Pkg; Pkg.activate(temp = true); Pkg.add(["OrdinaryDiffEq", "StaticArrays", "StructArrays"])
[...]
⌃ [1dea7af3] + OrdinaryDiffEq v6.89.0
  [90137ffa] + StaticArrays v1.9.8
⌃ [09ab397b] + StructArrays v0.6.18
[...]

julia> using OrdinaryDiffEq, StructArrays, StaticArrays

julia> u = StructArray([SVector(1.0, 1.0), SVector(2.0, 2.0)])
2-element StructArray(::Vector{Float64}, ::Vector{Float64}) with eltype SVector{2, Float64}:
 [1.0, 1.0]
 [2.0, 2.0]

julia> ode = ODEProblem(u, (0.0, 1.0)) do du, u, p, t
           @show typeof(du)
           du .= u
           return nothing
       end

julia> solve(ode, Tsit5())
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
[...]

julia> solve(ode, SSPRK43())
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
typeof(du) = StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}
[...]

I can't see the difference here that you describe 🤷

Odd! I'll check tomorrow (Thursdays are my Trixi PR day).

jlchan · 2024-11-20T20:48:19Z

@ranocha I misread your MWE; I couldn't reproduce this with a simplified MWE either. It only occurs in Trixi.jl.

I'll try to isolate it into a MWE tomorrow.

huiyuxie · 2024-11-22T00:46:42Z

Project.toml

@@ -94,7 +94,7 @@ PrecompileTools = "1.1"
 Preferences = "1.3"
 Printf = "1"
 RecipesBase = "1.1"
-RecursiveArrayTools = "2.38.10"
+RecursiveArrayTools = "3.27.1"


You are setting this compatibility version so high that SciMLBase will also adopt a higher version, forcing you to use the new version of the DiscreteCallback struct.

That's why even if you're using Julia 1.10 but it still failed here https://github.com/trixi-framework/Trixi.jl/actions/runs/11871869411/job/33084910314?pr=2150#step:7:24

Either change this to a lower version or adopt the new DiscreteCallback struct in precompile.

JoshuaLampert · 2024-12-03T16:08:14Z

Is there any update on this? I also need Trixi.jl being compatible with newer versions of RecursiveArrayTools.jl and OrdinaryDiffEq.jl.

jlchan · 2024-12-03T19:13:22Z

Is there any update on this? I also need Trixi.jl being compatible with newer versions of RecursiveArrayTools.jl and OrdinaryDiffEq.jl.

I've got a conference this week so I won't be able to return to this until next Monday. However, there aren't any major numerical issues to address IMO. The main failing tests are related to allocation (e.g., #2150 (comment)) and precompilation (#2150 (comment)). I believe @huiyuxie has identified some of these issues already.

If you have any time to look at those before I get back to it next week I'd welcome the help.

huiyuxie · 2024-12-04T00:05:36Z

I will open a PR to address this but I'm not familiar with new DiscreteCallback struct I may need help from @ranocha

ranocha · 2024-12-04T05:44:46Z

Please let me know if you have questions

jlchan added 2 commits November 7, 2024 09:28

update bound on RecursiveArrayTools to have access to 3.27.3

3b02770

d

comment out broken precompile statements

9b21e12

jlchan and others added 4 commits November 7, 2024 16:31

adding VoA for DGMulti

e66bcd9

Apply suggestions from code review

8da6bfa

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Merge branch 'main' into jc/wrap_VectorOfArray

d997c8f

run Trixi.rhs! twice to try to avoid excessive allocations

a20dc52

jlchan and others added 9 commits November 7, 2024 17:17

bump lower compat of RecursiveArrayTools

781a7d8

run Trixi.rhs! twice to reduce allocation count

1d04471

update l2, linf errors

c10b421

Apply suggestions from code review

88d08a1

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Update test/test_dgmulti_2d.jl

595a75e

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Apply suggestions from code review

b5c218d

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

Merge branch 'jc/wrap_VectorOfArray' of https://github.com/trixi-fram…

fea32d2

…ework/Trixi.jl into jc/wrap_VectorOfArray

bump LinearAlgebra lower compat for Downgrade CI

d69dbc7

jlchan added 2 commits November 8, 2024 09:28

Merge branch 'main' into jc/wrap_VectorOfArray

823caf1

Merge branch 'main' into jc/wrap_VectorOfArray

bd16969

Merge branch 'main' into jc/wrap_VectorOfArray

0ab4c2e

huiyuxie reviewed Nov 11, 2024

View reviewed changes

Project.toml Show resolved Hide resolved

src/auxiliary/precompile.jl Outdated Show resolved Hide resolved

test/test_dgmulti_1d.jl Show resolved Hide resolved

jlchan added 2 commits November 14, 2024 14:06

bump Julia compat and ci.yml to v1.10

181ac56

Merge branch 'jc/wrap_VectorOfArray' of https://github.com/trixi-fram…

1a89762

…ework/Trixi.jl into jc/wrap_VectorOfArray

github-actions bot reviewed Nov 14, 2024

View reviewed changes

test/test_visualization.jl Outdated Show resolved Hide resolved

jlchan and others added 2 commits November 14, 2024 14:09

Update test/test_visualization.jl

27729cc

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

fix VoA in test_visualization

025d968

jlchan mentioned this pull request Nov 14, 2024

Broadcasting over VectorOfArray does not preserve type of parent array SciML/RecursiveArrayTools.jl#410

Closed

jlchan commented Nov 14, 2024

View reviewed changes

test/test_dgmulti_3d.jl Outdated Show resolved Hide resolved

jlchan commented Nov 14, 2024

View reviewed changes

test/test_dgmulti_3d.jl Outdated Show resolved Hide resolved

jlchan added 2 commits November 14, 2024 19:54

Update test/test_dgmulti_3d.jl

d7b70a6

Update test/test_dgmulti_3d.jl

8481581

huiyuxie reviewed Nov 15, 2024

View reviewed changes

uncommenting precompile statement

fdf1ff5

jlchan mentioned this pull request Nov 15, 2024

Fix sign error in 3D DGMultiFDSBP elixir #2164

Merged

Merge branch 'main' into jc/wrap_VectorOfArray

31a85a5

huiyuxie mentioned this pull request Nov 18, 2024

Incorrect memory allocation measurements due to global variables #2166

Closed

huiyuxie suggested changes Nov 22, 2024

View reviewed changes

JoshuaLampert mentioned this pull request Dec 2, 2024

Use OrdinaryDiffEqTsit5.jl and OrdinaryDiffEqLowStorageRK.jl JoshuaLampert/DispersiveShallowWater.jl#163

Merged

huiyuxie mentioned this pull request Dec 8, 2024

Update RecursiveArrayTools.jl compatibility to version 3 #2194

Draft

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Use `VectorOfArray` in `wrap_array` for `DGMulti` solvers #2150

Use `VectorOfArray` in `wrap_array` for `DGMulti` solvers #2150

jlchan commented Nov 7, 2024 •

edited

Loading

github-actions bot commented Nov 7, 2024

codecov bot commented Nov 7, 2024 •

edited

Loading

jlchan commented Nov 7, 2024 •

edited

Loading

DanielDoehring commented Nov 8, 2024

jlchan commented Nov 8, 2024

DanielDoehring commented Nov 9, 2024

jlchan commented Nov 10, 2024 •

edited

Loading

huiyuxie commented Nov 11, 2024

huiyuxie commented Nov 11, 2024

jlchan commented Nov 11, 2024

huiyuxie commented Nov 11, 2024

huiyuxie Nov 15, 2024

jlchan Nov 15, 2024

huiyuxie Nov 15, 2024

huiyuxie commented Nov 15, 2024 •

edited

Loading

jlchan commented Nov 15, 2024

jlchan commented Nov 15, 2024 •

edited

Loading

huiyuxie commented Nov 15, 2024 •

edited

Loading

jlchan commented Nov 15, 2024

ranocha commented Nov 20, 2024

jlchan commented Nov 20, 2024

jlchan commented Nov 20, 2024

huiyuxie Nov 22, 2024

JoshuaLampert commented Dec 3, 2024

jlchan commented Dec 3, 2024

huiyuxie commented Dec 4, 2024

ranocha commented Dec 4, 2024

Use VectorOfArray in wrap_array for DGMulti solvers #2150

Are you sure you want to change the base?

Use VectorOfArray in wrap_array for DGMulti solvers #2150

Conversation

jlchan commented Nov 7, 2024 • edited Loading

github-actions bot commented Nov 7, 2024

Review checklist

Purpose and scope

Code quality

Documentation

Testing

Performance

Verification

codecov bot commented Nov 7, 2024 • edited Loading

Codecov Report

jlchan commented Nov 7, 2024 • edited Loading

DanielDoehring commented Nov 8, 2024

jlchan commented Nov 8, 2024

DanielDoehring commented Nov 9, 2024

jlchan commented Nov 10, 2024 • edited Loading

huiyuxie commented Nov 11, 2024

huiyuxie commented Nov 11, 2024

jlchan commented Nov 11, 2024

huiyuxie commented Nov 11, 2024

huiyuxie Nov 15, 2024

Choose a reason for hiding this comment

jlchan Nov 15, 2024

Choose a reason for hiding this comment

huiyuxie Nov 15, 2024

Choose a reason for hiding this comment

huiyuxie commented Nov 15, 2024 • edited Loading

jlchan commented Nov 15, 2024

jlchan commented Nov 15, 2024 • edited Loading

huiyuxie commented Nov 15, 2024 • edited Loading

jlchan commented Nov 15, 2024

ranocha commented Nov 20, 2024

jlchan commented Nov 20, 2024

jlchan commented Nov 20, 2024

huiyuxie Nov 22, 2024

Choose a reason for hiding this comment

JoshuaLampert commented Dec 3, 2024

jlchan commented Dec 3, 2024

huiyuxie commented Dec 4, 2024

ranocha commented Dec 4, 2024

Use `VectorOfArray` in `wrap_array` for `DGMulti` solvers #2150

Use `VectorOfArray` in `wrap_array` for `DGMulti` solvers #2150

jlchan commented Nov 7, 2024 •

edited

Loading

codecov bot commented Nov 7, 2024 •

edited

Loading

jlchan commented Nov 7, 2024 •

edited

Loading

jlchan commented Nov 10, 2024 •

edited

Loading

huiyuxie commented Nov 15, 2024 •

edited

Loading

jlchan commented Nov 15, 2024 •

edited

Loading

huiyuxie commented Nov 15, 2024 •

edited

Loading