You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
This PR links LAPACK to post_process and adds Liutex to the outputs of post_process using LAPACK subroutines.
Liutex (formerly named Rortex) is relatively recently invented vortex visualization method which provides more physical interpretation of identified vortical structures. See Xu et al. (2019) and Liu et al. (2018) for details.
New feature (non-breaking change which adds functionality)
Scope
This PR comprises a set of related changes with a common goal
How Has This Been Tested?
Below is the Liutex visualization of examples/3D_turb_mixing. liutex_mag = 1 and colored by pressure.
Test Configuration:
What computers and compilers did you use to test this: Carpenter GNU
Checklist
I have added comments for the new code
I added Doxygen docstrings to the new code
I have made corresponding changes to the documentation (docs/)
I have added regression tests to the test suite so that people can verify in the future that the feature is behaving as expected
I have added example cases in examples/ that demonstrate my new feature performing as expected.
They run to completion and demonstrate "interesting physics"
I ran ./mfc.sh format before committing my code
New and existing tests pass locally with my changes, including with GPU capability enabled (both NVIDIA hardware with NVHPC compilers and AMD hardware with CRAY compilers) and disabled
This PR does not introduce any repeated code (it follows the DRY principle)
I cannot think of a way to condense this code and reduce any introduced additional line count
PR Type
Enhancement
Description
Add LAPACK dependency integration to build system
Implement Liutex vortex visualization method in post_process
Add input validation and parameter handling for Liutex
Update documentation with Liutex configuration options
Diagram Walkthrough
flowchart LR
A["LAPACK dependency"] --> B["Build system"]
B --> C["post_process module"]
C --> D["Liutex computation"]
D --> E["Vortex visualization"]
F["Input validation"] --> D
G["Documentation"] --> D
The eigenvalue selection logic in the Liutex computation may be incorrect. The code finds the eigenvector with the smallest imaginary part but then uses a different eigenvalue index for computing the imaginary part of the complex eigenvalue, which could lead to incorrect results.
idx =1do r =2, 3if (abs(li(r)) < abs(li(idx))) then
idx = r
end ifend do
eigvec = vr(:, idx)
! Normalize real eigenvector if it is effectively non-zero
eigvec_mag =sqrt(eigvec(1)**2._wp &
+ eigvec(2)**2._wp &
+ eigvec(3)**2._wp)
if (eigvec_mag /=0._wp) then
eigvec = eigvec/eigvec_mag
end if
! Compute vorticity projected on the eigenvector
omega_proj = omega(1)*eigvec(1) &
+ omega(2)*eigvec(2) &
+ omega(3)*eigvec(3)
! As eigenvector can have +/- signs, we can choose the sign
! so that omega_proj is positive
if (omega_proj < 0._wp) then
eigvec =-eigvec
omega_proj =-omega_proj
end if
! Find imaginary part of complex eigenvalue
lci = li(mod(idx, 3) +1)
The Liutex computation involves expensive LAPACK eigenvalue decomposition at every grid point in a triple nested loop, which could significantly impact performance for large grids. Consider optimization strategies or parallelization.
do l =-offset_z%beg, p + offset_z%end
do k =-offset_y%beg, n + offset_y%end
do j =-offset_x%beg, m + offset_x%end
! Get velocity gradient tensor (VGT)
vgt(:, :) =0._wpdo r =-fd_number, fd_number
do i =1, 3
! d()/dx
vgt(i, 1) = &
vgt(i, 1) + &
fd_coeff_x(r, j)* &
q_prim_vf(mom_idx%beg + i -1)%sf(r + j, k, l)
! d()/dy
vgt(i, 2) = &
vgt(i, 2) + &
fd_coeff_y(r, k)* &
q_prim_vf(mom_idx%beg + i -1)%sf(j, r + k, l)
! d()/dz
vgt(i, 3) = &
vgt(i, 3) + &
fd_coeff_z(r, l)* &
q_prim_vf(mom_idx%beg + i -1)%sf(j, k, r + l)
end doend do
! Compute vorticity
omega(1) = vgt(3, 2) - vgt(2, 3)
omega(2) = vgt(1, 3) - vgt(3, 1)
omega(3) = vgt(2, 1) - vgt(1, 2)
! Call appropriate LAPACK routine based on precision
#ifdef MFC_SINGLE_PRECISION
call cgeev(ivl, ivr, nm, vgt, nm, lr, li, vl, nm, vr, nm, work, lwork, info)
#elsecall dgeev(ivl, ivr, nm, vgt, nm, lr, li, vl, nm, vr, nm, work, lwork, info)
#endif
! Find real eigenvector
idx =1do r =2, 3if (abs(li(r)) < abs(li(idx))) then
idx = r
end ifend do
eigvec = vr(:, idx)
! Normalize real eigenvector if it is effectively non-zero
eigvec_mag =sqrt(eigvec(1)**2._wp &
+ eigvec(2)**2._wp &
+ eigvec(3)**2._wp)
if (eigvec_mag /=0._wp) then
eigvec = eigvec/eigvec_mag
end if
! Compute vorticity projected on the eigenvector
omega_proj = omega(1)*eigvec(1) &
+ omega(2)*eigvec(2) &
+ omega(3)*eigvec(3)
! As eigenvector can have +/- signs, we can choose the sign
! so that omega_proj is positive
if (omega_proj < 0._wp) then
eigvec =-eigvec
omega_proj =-omega_proj
end if
! Find imaginary part of complex eigenvalue
lci = li(mod(idx, 3) +1)
! Compute Liutex magnitude
alpha = omega_proj**2._wp-4._wp*lci**2._wp ! (2*alpha)^2if (alpha > 0._wp) then
liutex_mag(j, k, l) = omega_proj -sqrt(alpha)
else
liutex_mag(j, k, l) = omega_proj
end if
! Compute Liutex axis
liutex_axis(j, k, l, 1) = eigvec(1)
liutex_axis(j, k, l, 2) = eigvec(2)
liutex_axis(j, k, l, 3) = eigvec(3)
end doend doend do
The velocity gradient tensor computation code is duplicated from vorticity calculation. This violates DRY principle and should be refactored into a shared utility function.
! Get velocity gradient tensor (VGT)
vgt(:, :) =0._wpdo r =-fd_number, fd_number
do i =1, 3
! d()/dx
vgt(i, 1) = &
vgt(i, 1) + &
fd_coeff_x(r, j)* &
q_prim_vf(mom_idx%beg + i -1)%sf(r + j, k, l)
! d()/dy
vgt(i, 2) = &
vgt(i, 2) + &
fd_coeff_y(r, k)* &
q_prim_vf(mom_idx%beg + i -1)%sf(j, r + k, l)
! d()/dz
vgt(i, 3) = &
vgt(i, 3) + &
fd_coeff_z(r, l)* &
q_prim_vf(mom_idx%beg + i -1)%sf(j, k, r + l)
end doend do
! Compute vorticity
omega(1) = vgt(3, 2) - vgt(2, 3)
omega(2) = vgt(1, 3) - vgt(3, 1)
omega(3) = vgt(2, 1) - vgt(1, 2)
The calculation of lci using mod(idx, 3) + 1 is incorrect and will not give the proper complex eigenvalue pair. For a 3x3 matrix, complex eigenvalues come in conjugate pairs, so you need to identify which eigenvalues form the complex pair rather than using a modular arithmetic approach.
! Find imaginary part of complex eigenvalue
-lci = li(mod(idx, 3) + 1)+! Complex eigenvalues come in conjugate pairs+if (idx == 1) then+ lci = max(abs(li(2)), abs(li(3)))+elseif (idx == 2) then+ lci = max(abs(li(1)), abs(li(3)))+else+ lci = max(abs(li(1)), abs(li(2)))+end if
Suggestion importance[1-10]: 8
__
Why: The suggestion correctly identifies a bug in the logic for finding the complex eigenvalue, as the mod function can fail, and provides a robust alternative that correctly handles all cases.
Medium
Validate real eigenvalue selection logic
The logic for finding the real eigenvector is incorrect. It should find the eigenvalue with the smallest imaginary part, but the current implementation finds the index with the smallest absolute imaginary part, which could still be complex. Add a check to ensure the selected eigenvalue is actually real (imaginary part near zero).
! Find real eigenvector
idx = 1
do r = 2, 3
if (abs(li(r)) < abs(li(idx))) then
idx = r
end if
end do
+! Ensure we have a real eigenvalue (imaginary part should be near zero)+if (abs(li(idx)) > 1.0e-12_wp) then+ ! If no real eigenvalue found, use the one with smallest imaginary part+ ! but this indicates a potential issue with the computation+end if
eigvec = vr(:, idx)
Suggestion importance[1-10]: 6
__
Why: The suggestion correctly points out a potential robustness issue by proposing a check to ensure the identified eigenvalue is truly real, which is a valid defensive programming practice.
Low
General
✅ Use tolerance for eigenvector normalizationSuggestion Impact:The suggestion was implemented by replacing the exact equality check (eigvec_mag /= 0._wp) with a tolerance-based comparison (eigvec_mag .gt. sgm_eps) and adding an else clause to set eigvec to zero when the magnitude is below the tolerance
code diff:
- if (eigvec_mag /= 0._wp) then+ if (eigvec_mag .gt. sgm_eps) then
eigvec = eigvec/eigvec_mag
+ else+ eigvec = 0._wp
end if
Using exact equality comparison with floating-point zero is unreliable and can lead to division by very small numbers. Use a tolerance-based comparison to avoid numerical instability when the eigenvector magnitude is very small but not exactly zero.
! Normalize real eigenvector if it is effectively non-zero
eigvec_mag = sqrt(eigvec(1)**2._wp &
+ eigvec(2)**2._wp &
+ eigvec(3)**2._wp)
-if (eigvec_mag /= 0._wp) then+if (eigvec_mag > 1.0e-12_wp) then
eigvec = eigvec/eigvec_mag
+else+ eigvec = 0._wp
end if
Suggestion importance[1-10]: 7
__
Why: The suggestion correctly advises using a tolerance for floating-point comparison to prevent potential numerical instability, which is a standard and important practice in numerical code.
❌ Patch coverage is 20.89552% with 53 lines in your changes missing coverage. Please review.
✅ Project coverage is 43.57%. Comparing base (a1d1576) to head (51b37d4). ⚠️ Report is 3 commits behind head on master.
@sbryngelson I have no idea why Frontier test fails with the Could NOT find LAPACK error after merging master branch. Before merging it, all tests passed (at the commit 1686f50). I reviewed the change, but unclear what's going on. Do you have any thoughts?
If you look at the printed logs from the earlier CI runs that passed, it looks like it never actually built LAPACK on the Frontier cases... which is strange but is certainly why they passed.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
User description
Description
This PR links LAPACK to post_process and adds Liutex to the outputs of post_process using LAPACK subroutines.
Liutex (formerly named Rortex) is relatively recently invented vortex visualization method which provides more physical interpretation of identified vortical structures. See Xu et al. (2019) and Liu et al. (2018) for details.
LAPACK implementation is based on #939
Type of change
Scope
How Has This Been Tested?
Below is the Liutex visualization of
examples/3D_turb_mixing.liutex_mag = 1and colored by pressure.Test Configuration:
Checklist
docs/)examples/that demonstrate my new feature performing as expected.They run to completion and demonstrate "interesting physics"
./mfc.sh formatbefore committing my codePR Type
Enhancement
Description
Add LAPACK dependency integration to build system
Implement Liutex vortex visualization method in post_process
Add input validation and parameter handling for Liutex
Update documentation with Liutex configuration options
Diagram Walkthrough
File Walkthrough
11 files
Add Liutex input validation checksImplement Liutex computation using LAPACK eigenvalue solverAdd `liutex_wrt` global parameterAdd `liutex_wrt` to MPI broadcast variablesIntegrate Liutex output in data processing pipelineEnable Liutex output in turbulent mixing exampleAdd LAPACK as dependency for post_process targetAdd `liutex_wrt` parameter to case dictionaryCreate LAPACK CMake finder with Cray system supportIntegrate LAPACK linking in MFC build systemAdd LAPACK external dependency build option2 files
Document `liutex_wrt` parameter and usageAdd Liutex method reference citation