Switch `DGMulti` to `Matrix{<:SVector}` solution storage #1952

jlchan · 2024-05-23T16:22:19Z

Also bumps StartUpDG.jl compat to 1.0.

* remove Octavian and matmul! * removing some glue code introduced for matmul! * actually remove Octavian this time * Revert "removing some glue code introduced for matmul!" This reverts commit 1fd1897.

github-actions · 2024-05-23T16:22:35Z

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-05-23T17:04:52Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 82.03%. Comparing base (c2513e2) to head (c7b38bb).

Additional details and impacted files

@@             Coverage Diff             @@
##             main    #1952       +/-   ##
===========================================
- Coverage   96.09%   82.03%   -14.06%     
===========================================
  Files         457      457               
  Lines       36734    36661       -73     
===========================================
- Hits        35298    30074     -5224     
- Misses       1436     6587     +5151

Flag	Coverage Δ
unittests	`82.03% <100.00%> (-14.06%)`	⬇️

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.

jlchan · 2024-05-23T19:50:15Z

I am running into issues due to Polyester.@batch turning arrays into StrideArraysCore.PtrArrays". For a Navier-Stokes solve, this dispatches to Octavian.matmul! via https://github.com/JuliaSIMD/StrideArrays.jl/blob/d31fb6cf28b0374162202da8c81694f2b704e0f3/src/blas.jl#L6.

For now, I'll just focus on the single-threaded case, but dealing with multi-threading will require figuring out if Polyester and StrideArrays will continue after stopping support of LoopVec.jl and Octavian.jl in Julia v1.11 (cc @chriselrod)

chriselrod · 2024-05-23T23:42:19Z

For a Navier-Stokes solve, this dispatches to Octavian.matmul! via https://github.com/JuliaSIMD/StrideArrays.jl/blob/d31fb6cf28b0374162202da8c81694f2b704e0f3/src/blas.jl#L6.

Feel free to PR it to drop Octavian as a dependency.

For now, I'll just focus on the single-threaded case, but dealing with multi-threading will require figuring out if Polyester and StrideArrays will continue after stopping support of LoopVec.jl and Octavian.jl in Julia v1.11

Regardless of the official status, you don't have to rush to make changes. There are no test failures on 1.11.

chriselrod · 2024-05-23T23:48:48Z

src/solvers/dgmulti/dg.jl

-# Allocate nested array type for DGMulti solution storage.
 function allocate_nested_array(uEltype, nvars, array_dimensions, dg)
-    # store components as separate arrays, combine via StructArrays
-    return StructArray{SVector{nvars, uEltype}}(ntuple(_ -> zeros(uEltype,
-                                                                  array_dimensions...),
-                                                       nvars))


Out of curiosity, what's the motivation for the change and how does it perform?
The StructArray version should be SIMD-friendlier, while the Matrix{<:SVector} would be better if SIMD isn't possible.

The motivation is due mainly to #1789. This impacts my specific Trixi solver DGMulti and locks Trixi.jl to a specific version of OrdinaryDiffEq. I looked at the array wrappers in RecursiveArrayTools.jl to try to address this, but it takes a lot more energy for me to figure out customizing broadcasting over custom arrays than just using Base.Array. I'm also concerned about changes like SciML/RecursiveArrayTools.jl#343.

My plan is to swap between Matrix{<:SVector{N, T}} and Array{T, 3} using unsafe_wrap, which is currently what Trixi does for other solvers too. My guess is that it won't be nearly as performant once we make these changes, since we have a few pieces of code which are accelerated by Octavian.jl and LoopVec.jl.

However, I'm less worried about efficiency at the moment and more about maintainability, which I think this would help with.

jlchan · 2024-05-24T04:41:51Z

For a Navier-Stokes solve, this dispatches to Octavian.matmul! via https://github.com/JuliaSIMD/StrideArrays.jl/blob/d31fb6cf28b0374162202da8c81694f2b704e0f3/src/blas.jl#L6.

Feel free to PR it to drop Octavian as a dependency.

For now, I'll just focus on the single-threaded case, but dealing with multi-threading will require figuring out if Polyester and StrideArrays will continue after stopping support of LoopVec.jl and Octavian.jl in Julia v1.11

Regardless of the official status, you don't have to rush to make changes. There are no test failures on 1.11.

Thanks - I removed the LinearAlgebra.mul! overloads in https://github.com/jlchan/StrideArrays.jl/blob/d31fb6cf28b0374162202da8c81694f2b704e0f3/src/blas.jl#L2-L36. It runs, but the solution gets NaNs. When I switch @batch to Threads.@threads, I get an error that seems related to @turbo:

ERROR: TaskFailedException

    nested task error: UndefRefError: access to undefined reference

so I think I'm going to have to dig a little deeper beyond removing Octavian dependencies to figure out what's going wrong...

huiyuxie · 2024-10-17T02:11:08Z

Can you @jlchan try replace ::StructArray with ::Matrix for some functions like below. If you change allocate_nested_array to make arrays from cache not store in StructArray (now in Matrix) then function like

Trixi.jl/src/solvers/dgmulti/flux_differencing.jl

Lines 401 to 407 in 00ef308

    
           @inline function cons2entropy!(entropy_var_values::StructArray, 
        
                                          u_values::StructArray, 
        
                                          equations) 
        
               @threaded for i in eachindex(u_values) 
        
                   entropy_var_values[i] = cons2entropy(u_values[i], equations) 
        
               end 
        
           end

does not work since Matrix is not a sub type of StructArray.

I am still trying to fix but I am not so familiar with dgmulti so things look a little bit messy for me now.

huiyuxie · 2024-10-17T02:13:30Z

Trixi.jl/src/solvers/dgmulti/flux_differencing.jl

Lines 384 to 399 in 00ef308

    
           function entropy_projection!(cache, u, mesh::DGMultiMesh, equations, dg::DGMulti) 
        
               rd = dg.basis 
        
               @unpack Vq = rd 
        
               @unpack VhP, entropy_var_values, u_values = cache 
        
               @unpack projected_entropy_var_values, entropy_projected_u_values = cache 
        
               apply_to_each_field(mul_by!(Vq), u_values, u) 
        
               cons2entropy!(entropy_var_values, u_values, equations) 
        
               # "VhP" fuses the projection "P" with interpolation to volume and face quadrature "Vh" 
        
               apply_to_each_field(mul_by!(VhP), projected_entropy_var_values, entropy_var_values) 
        
               entropy2cons!(entropy_projected_u_values, projected_entropy_var_values, equations) 
        
               return nothing 
        
           end

For example the cons2entropy! above will not work

jlchan · 2024-10-17T02:33:56Z

@huiyuxie thanks for taking a look. I don't think I end up using those functions in this PR. Since I switched the underlying solution type to Matrix{<:SVector}, I believe I am calling some other versions of those functions (I don't remember their precise location at the moment). This change should allow me to delete those functions; I just haven't gotten around to removing those yet since the CI is still failing.

FYI @huiyuxie, I am discussing how best to clean up DGMulti with the other core Trixi PIs. It is hard to maintain right now because there is a lot of specialized code, but I do not think I will use much of it. Therefore, we are considering removing significant parts of DGMulti, but we have to have some discussions with some of the users of DGMulti first.

Since there will likely be significant changes to DGMulti in the near future, I don't want to ask you to try to digest/modify DGMulti code if it will ultimately be wasted. What do you think @huiyuxie?

huiyuxie · 2024-10-17T03:30:24Z

I see thanks for letting me know @jlchan!

Does this mean that I don't have to keep fixing this issue and can wait for you to resolve it? I am also curious why DGMulti is the only part that uses StructArray, i.e., why not use Matrix directly during the development stage?

Since there will likely be significant changes to DGMulti in the near future, I don't want to ask you to try to digest/modify DGMulti code if it will ultimately be wasted.

Will it take long? If you need help implementing the new version, please let me know; I am willing to help.

jlchan · 2024-10-17T03:37:52Z

I see thanks for letting me know @jlchan! Does this mean that I don't have to keep fixing this issue and can wait for you to resolve it?

Yes, I think that would be best. I will try to keep you in the loop, however - and will ask for your help (if you are still interested at the time) once we have planned out how best to fix things.

I am also curious why DGMulti is the only part that uses StructArray, i.e., why not use Matrix directly during the development stage?

Good question; it is because in DGMulti, certain solvers need to do matrix-matrix multiplication for each variable. Suppose I have 2 variables;StructArrays combines 2 separate matrices A, B together lazily into one matrix whose entries are SVector{A[i,j], B[i,j]}, which is convenient for certain parts of our solvers. However, I can still access the original A and B arrays and multiply them by other matrices, which is more efficient than using a Matrix{<:SVector} layout.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Switch `DGMulti` to `Matrix{<:SVector}` solution storage #1952

Switch `DGMulti` to `Matrix{<:SVector}` solution storage #1952

jlchan commented May 23, 2024

github-actions bot commented May 23, 2024

codecov bot commented May 23, 2024 •

edited

Loading

jlchan commented May 23, 2024

chriselrod commented May 23, 2024

chriselrod May 23, 2024

jlchan May 24, 2024 •

edited

Loading

jlchan commented May 24, 2024 •

edited

Loading

huiyuxie commented Oct 17, 2024

huiyuxie commented Oct 17, 2024 •

edited

Loading

jlchan commented Oct 17, 2024

huiyuxie commented Oct 17, 2024 •

edited

Loading

jlchan commented Oct 17, 2024 •

edited

Loading

jlchan commented Oct 17, 2024

huiyuxie commented Oct 17, 2024

Switch DGMulti to Matrix{<:SVector} solution storage #1952

Are you sure you want to change the base?

Switch DGMulti to Matrix{<:SVector} solution storage #1952

Conversation

jlchan commented May 23, 2024

github-actions bot commented May 23, 2024

Review checklist

Purpose and scope

Code quality

Documentation

Testing

Performance

Verification

codecov bot commented May 23, 2024 • edited Loading

Codecov Report

jlchan commented May 23, 2024

chriselrod commented May 23, 2024

chriselrod May 23, 2024

Choose a reason for hiding this comment

jlchan May 24, 2024 • edited Loading

Choose a reason for hiding this comment

jlchan commented May 24, 2024 • edited Loading

huiyuxie commented Oct 17, 2024

huiyuxie commented Oct 17, 2024 • edited Loading

jlchan commented Oct 17, 2024

huiyuxie commented Oct 17, 2024 • edited Loading

jlchan commented Oct 17, 2024 • edited Loading

jlchan commented Oct 17, 2024

huiyuxie commented Oct 17, 2024

Switch `DGMulti` to `Matrix{<:SVector}` solution storage #1952

Switch `DGMulti` to `Matrix{<:SVector}` solution storage #1952

codecov bot commented May 23, 2024 •

edited

Loading

jlchan May 24, 2024 •

edited

Loading

jlchan commented May 24, 2024 •

edited

Loading

huiyuxie commented Oct 17, 2024 •

edited

Loading

huiyuxie commented Oct 17, 2024 •

edited

Loading

jlchan commented Oct 17, 2024 •

edited

Loading