Add velocity functions for different equations #2086

Arpit-Babbar · 2024-09-23T13:30:36Z

Trixi.jl has functions to compute density and pressure for all equations, but most equations do not have those for velocity. This PR adds those functions to CompressibleEuler, MultiComponentEuler, PolytropicEuler, IdealMHD, ShallowWater equations.

github-actions · 2024-09-23T13:30:52Z

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-09-23T14:58:46Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 96.35%. Comparing base (1672c2a) to head (82976ac).
Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2086      +/-   ##
==========================================
+ Coverage   96.34%   96.35%   +0.01%     
==========================================
  Files         470      470              
  Lines       37501    37588      +87     
==========================================
+ Hits        36129    36216      +87     
  Misses       1372     1372

Flag	Coverage Δ
unittests	`96.35% <100.00%> (+0.01%)`	⬆️

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.

ranocha

Thanks a lot! Could you please

export velocity in the main file
write a docstring describing velocity in general - maybe in equations.jl since it is a rather general interface?
add some unit tests checking whether the correct values are returned in https://github.com/trixi-framework/Trixi.jl/blob/main/test/test_unit.jl?

…e 3D for all dims

ranocha

Thanks a lot! I just have some minor comments.

src/Trixi.jl

test/test_unit.jl

ranocha · 2024-09-26T03:42:05Z

Two of the GLM MHD 2D tests are failing: https://github.com/trixi-framework/Trixi.jl/actions/runs/11035734678/job/30652521365?pr=2086#step:7:1405

src/equations/ideal_glm_mhd_2d.jl

src/equations/ideal_glm_mhd_1d.jl

andrewwinters5000 · 2024-09-26T09:09:22Z

The ShallowWaterEquations in 1D and 2D already have a velocity function, e.g., here that returns each Cartesian directions velocity. The normal direction velocity is then computed within the function if necessary.

andrewwinters5000 · 2024-09-26T09:12:56Z

src/equations/compressible_euler_2d.jl

+@inline function velocity(u, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v = v1 * normal_direction[1] + v2 * normal_direction[2]
+    return v
+end


Why do we need three versions of this functionality? Wouldn't it be sufficient to use the one that only dispatches on the equation set and returns the SVector of Cartesian fluxes. Then computing a particular normal direction, using orientation or a vector, could be done within the function that calls velocity.

Thank you for the observation!

I updated the PR to have a general function for the velocity in normal_direction.

However, having different versions for the orientation is slightly cheaper, as the general velocity function will compute velocity in all directions even though we need only one. Thus, there is a minor readability-performance trade-off and I would like your thoughts on it.

I can see a reason for this. Checking a simple example, the compiler seems to be smart enough to optimize away the case of constant orientations but not with variable ones, e.g.,

julia> using StaticArrays julia> @inline function velocity(u) rho = u[1] v1 = u[2] / rho v2 = u[3] / rho return SVector(v1, v2) end velocity (generic function with 1 method) julia> @inline function velocity(u, orientation) rho = u[1] v = u[1 + orientation] / rho return v end velocity (generic function with 2 methods) julia> foo1(u) = velocity(u)[1] foo1 (generic function with 1 method) julia> foo2(u) = velocity(u, 1) foo2 (generic function with 1 method) julia> u = SVector(1.0, 0.1, -0.2, 1.5) 4-element SVector{4, Float64} with indices SOneTo(4): 1.0 0.1 -0.2 1.5 julia> @code_llvm debuginfo=:none foo1(u) define double @julia_foo1_188([1 x [4 x double]]* nocapture noundef nonnull readonly align 8 dereferenceable(32) %0) #0 { top: %1 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 0 %2 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 1 %unbox = load double, double* %2, align 8 %unbox1 = load double, double* %1, align 8 %3 = fdiv double %unbox, %unbox1 ret double %3 } julia> @code_llvm debuginfo=:none foo2(u) define double @julia_foo2_201([1 x [4 x double]]* nocapture noundef nonnull readonly align 8 dereferenceable(32) %0) #0 { top: %1 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 0 %2 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 1 %unbox = load double, double* %2, align 8 %unbox1 = load double, double* %1, align 8 %3 = fdiv double %unbox, %unbox1 ret double %3 } julia> bar1(u, i) = velocity(u)[i] bar1 (generic function with 1 method) julia> bar2(u, i) = velocity(u, i) bar2 (generic function with 1 method) julia> @code_llvm debuginfo=:none bar1(u, 1) define double @julia_bar1_215([1 x [4 x double]]* nocapture noundef nonnull readonly align 8 dereferenceable(32) %0, i64 signext %1) #0 { top: %newstruct = alloca <2 x double>, align 16 %2 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 0 %3 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 1 %unbox1 = load double, double* %2, align 8 %4 = bitcast double* %3 to <2 x double>* %5 = load <2 x double>, <2 x double>* %4, align 8 %6 = insertelement <2 x double> poison, double %unbox1, i64 0 %7 = shufflevector <2 x double> %6, <2 x double> poison, <2 x i32> zeroinitializer %8 = fdiv <2 x double> %5, %7 store <2 x double> %8, <2 x double>* %newstruct, align 16 %9 = add i64 %1, -1 %boundscheck = icmp ult i64 %9, 2 br i1 %boundscheck, label %pass, label %fail fail: ; preds = %top %10 = bitcast <2 x double>* %newstruct to i8* call void @ijl_bounds_error_unboxed_int(i8* %10, {}* nonnull inttoptr (i64 4803813760 to {}*), i64 %1) unreachable pass: ; preds = %top %11 = getelementptr inbounds <2 x double>, <2 x double>* %newstruct, i64 0, i64 %9 %unbox4 = load double, double* %11, align 8 ret double %unbox4 } julia> @code_llvm debuginfo=:none bar2(u, 1) define double @julia_bar2_217([1 x [4 x double]]* nocapture noundef nonnull readonly align 8 dereferenceable(32) %0, i64 signext %1) #0 { top: %boundscheck = icmp ult i64 %1, 4 br i1 %boundscheck, label %pass, label %fail fail: ; preds = %top %2 = add i64 %1, 1 %3 = bitcast [1 x [4 x double]]* %0 to i8* call void @ijl_bounds_error_unboxed_int(i8* nonnull %3, {}* nonnull inttoptr (i64 4766129648 to {}*), i64 %2) unreachable pass: ; preds = %top %4 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 0 %5 = getelementptr inbounds [1 x [4 x double]], [1 x [4 x double]]* %0, i64 0, i64 0, i64 %1 %unbox = load double, double* %5, align 8 %unbox1 = load double, double* %4, align 8 %6 = fdiv double %unbox, %unbox1 ret double %6 }

What do you think, @andrewwinters5000?

ranocha

Are you fine with this, @andrewwinters5000?

Add velocity functions for different equations

217788a

Arpit-Babbar added 3 commits September 23, 2024 15:36

Formatter

9d2a0eb

Fix CI

7386197

Format

35b776d

Arpit-Babbar marked this pull request as ready for review September 24, 2024 08:56

ranocha requested changes Sep 24, 2024

View reviewed changes

Arpit-Babbar added 3 commits September 24, 2024 16:56

Export velocity, add docs

893accd

Export velocity, add docstring, add unit tests, make MHD velocity lik…

6f2ae51

…e 3D for all dims

Refractor bug fixgit push

7428094

ranocha requested changes Sep 26, 2024

View reviewed changes

src/Trixi.jl Outdated Show resolved Hide resolved

test/test_unit.jl Outdated Show resolved Hide resolved

ranocha requested changes Sep 26, 2024

View reviewed changes

src/equations/ideal_glm_mhd_2d.jl Outdated Show resolved Hide resolved

src/equations/ideal_glm_mhd_1d.jl Outdated Show resolved Hide resolved

Changes as per review

cdf2653

andrewwinters5000 reviewed Sep 26, 2024

View reviewed changes

Arpit-Babbar and others added 2 commits September 29, 2024 20:56

General velocity normal function

d29df2b

Merge branch 'main' into velocity

82976ac

ranocha requested a review from andrewwinters5000 September 30, 2024 18:24

ranocha approved these changes Oct 2, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add velocity functions for different equations #2086

Add velocity functions for different equations #2086

Arpit-Babbar commented Sep 23, 2024 •

edited

Loading

github-actions bot commented Sep 23, 2024

codecov bot commented Sep 23, 2024 •

edited

Loading

ranocha left a comment

ranocha left a comment •

edited

Loading

ranocha commented Sep 26, 2024

andrewwinters5000 commented Sep 26, 2024

andrewwinters5000 Sep 26, 2024

Arpit-Babbar Sep 29, 2024

ranocha Sep 30, 2024

ranocha Sep 30, 2024

ranocha left a comment

Add velocity functions for different equations #2086

Are you sure you want to change the base?

Add velocity functions for different equations #2086

Conversation

Arpit-Babbar commented Sep 23, 2024 • edited Loading

github-actions bot commented Sep 23, 2024

Review checklist

Purpose and scope

Code quality

Documentation

Testing

Performance

Verification

codecov bot commented Sep 23, 2024 • edited Loading

Codecov Report

ranocha left a comment

Choose a reason for hiding this comment

ranocha left a comment • edited Loading

Choose a reason for hiding this comment

ranocha commented Sep 26, 2024

andrewwinters5000 commented Sep 26, 2024

andrewwinters5000 Sep 26, 2024

Choose a reason for hiding this comment

Arpit-Babbar Sep 29, 2024

Choose a reason for hiding this comment

ranocha Sep 30, 2024

Choose a reason for hiding this comment

ranocha Sep 30, 2024

Choose a reason for hiding this comment

ranocha left a comment

Choose a reason for hiding this comment

Arpit-Babbar commented Sep 23, 2024 •

edited

Loading

codecov bot commented Sep 23, 2024 •

edited

Loading

ranocha left a comment •

edited

Loading