Add Liutex to post_process using LAPACK

[hyeoksu-lee]

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

• 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

Loading

File Walkthrough

	Relevant files
Enhancement	11 files m_checker.fpp Add Liutex input validation checks                                              +13/-2    m_derived_variables.fpp Implement Liutex computation using LAPACK eigenvalue solver +138/-3  m_global_parameters.fpp Add `liutex_wrt` global parameter                                                +3/-1      m_mpi_proxy.fpp Add `liutex_wrt` to MPI broadcast variables                            +1/-1      m_start_up.f90 Integrate Liutex output in data processing pipeline            +37/-5    case.py Enable Liutex output in turbulent mixing example                  +1/-0      build.py Add LAPACK as dependency for post_process target                  +3/-2      case_dicts.py Add `liutex_wrt` parameter to case dictionary                        +1/-0      FindLAPACK.cmake Create LAPACK CMake finder with Cray system support            +60/-0    CMakeLists.txt Integrate LAPACK linking in MFC build system                          +8/-2      CMakeLists.txt Add LAPACK external dependency build option                            +27/-0
Documentation	2 files case.md Document `liutex_wrt` parameter and usage                                +2/-1      references.md Add Liutex method reference citation                                          +2/-0

---

[qodo-merge-pro]

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪

🧪 No relevant tests

🔒 No security concerns identified

⚡ Recommended focus areas for review Possible Issue 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 = 1 do r = 2, 3 if (abs(li(r)) < abs(li(idx))) then idx = r end if end 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) Performance Concern 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._wp do 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 do end 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) #else call dgeev(ivl, ivr, nm, vgt, nm, lr, li, vl, nm, vr, nm, work, lwork, info) #endif ! Find real eigenvector idx = 1 do r = 2, 3 if (abs(li(r)) < abs(li(idx))) then idx = r end if end 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)^2 if (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 do end do end do Duplicate Code 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._wp do 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 do end 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)

[qodo-merge-pro]

PR Code Suggestions ✨

Explore these optional code suggestions:

Category	Suggestion	Impact
Possible issue	Fix complex eigenvalue pair identification 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. src/post_process/m_derived_variables.fpp [672-673] ! 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). src/post_process/m_derived_variables.fpp [643-650] ! 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 normalization Suggestion 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. src/post_process/m_derived_variables.fpp [652-658] ! 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.	Medium
Update

[codecov]

Codecov Report

❌ 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.

Files with missing lines	Patch %	Lines
src/post_process/m_derived_variables.fpp	2.32%	42 Missing ⚠️
src/post_process/m_start_up.f90	31.25%	10 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #970      +/-   ##
==========================================
- Coverage   43.67%   43.57%   -0.10%     
==========================================
  Files          70       70              
  Lines       19818    19881      +63     
  Branches     2473     2479       +6     
==========================================
+ Hits         8655     8663       +8     
- Misses       9644     9697      +53     
- Partials     1519     1521       +2     

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

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[hyeoksu-lee]

@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?

[sbryngelson]

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.

[sbryngelson]

i think i fixed it. can you include an example in the documentation? this is a cool visualization!