Refac helpers

[sbryngelson]

User description

Refactors some helper functions.

---

PR Type

Enhancement

---

Description

• Extract common CFL calculation logic into reusable function

• Consolidate repeated code patterns in stability and dt calculations

• Simplify viscous CFL calculations with inline expressions

• Improve code maintainability and reduce duplication

---

Changes walkthrough 📝

	Relevant files
Enhancement	m_sim_helpers.f90 Refactor CFL calculations with helper function                      src/simulation/m_sim_helpers.f90 Added f_compute_cfl_terms function to consolidate CFL calculations Refactored s_compute_stability_from_dt to use new helper function Refactored s_compute_dt_from_cfl to use new helper function Simplified viscous CFL calculations with inline expressions +97/-117

---

Need help?

Type /help how to ... in the comments thread for any questions about Qodo Merge usage.

Check out the documentation for more information.

[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 Logic Bug The 1D case handling in the new helper function differs from the original code. Original code uses multiplication for stability calculation, but the helper function uses division with reciprocal logic that may not be equivalent. !1D - special handling for different calculation types if (for_dt_calc) then cfl_terms = dx(j)/(abs(vel(1)) + c) else cfl_terms = (1._wp/dx(j))*(abs(vel(1)) + c) end if end if Possible Issue Complex inline expressions for cylindrical geometry calculations are hard to verify and may introduce calculation errors. The fltr_dtheta logic is embedded in long expressions that differ from the original extracted variable approach. vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k), & 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)))**2._wp Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), dy(k)*(abs(vel(2)) + c), & 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp))*(abs(vel(3)) + c))/maxval(1._wp/Re_l) Code Duplication The cylindrical geometry filtering calculation is duplicated in both the helper function and the viscous calculation sections, defeating the purpose of the refactoring. if (grid_geometry == 3) then if (k == 0) then fltr_dtheta = 2._wp*pi*y_cb(0)/3._wp elseif (k <= fourier_rings) then Nfq = min(floor(2._wp*real(k, wp)*pi), (p + 1)/2 + 1) fltr_dtheta = 2._wp*pi*y_cb(k - 1)/real(Nfq, wp) else fltr_dtheta = y_cb(k - 1)*dz(l) end if end if

[qodo-merge-pro]

PR Code Suggestions ✨

Latest suggestions up to b007e50

Category	Suggestion	Impact
Possible issue	Prevent division by zero errors The function can cause division by zero when abs(vel(i)) + c equals zero, leading to undefined behavior. Add a small epsilon value to the denominator to prevent division by zero and ensure numerical stability. src/simulation/m_sim_helpers.f90 [48-74] pure function f_compute_multidim_cfl_terms(vel, c, j, k, l) result(cfl_terms) !$acc routine seq real(wp), dimension(num_vels), intent(in) :: vel real(wp), intent(in) :: c integer, intent(in) :: j, k, l real(wp) :: cfl_terms real(wp) :: fltr_dtheta + real(wp), parameter :: eps = 1.0e-12_wp fltr_dtheta = f_compute_filtered_dtheta(k, l) if (p > 0) then !3D if (grid_geometry == 3) then - cfl_terms = min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - fltr_dtheta/(abs(vel(3)) + c)) + cfl_terms = min(dx(j)/(abs(vel(1)) + c + eps), & + dy(k)/(abs(vel(2)) + c + eps), & + fltr_dtheta/(abs(vel(3)) + c + eps)) else - cfl_terms = min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - dz(l)/(abs(vel(3)) + c)) + cfl_terms = min(dx(j)/(abs(vel(1)) + c + eps), & + dy(k)/(abs(vel(2)) + c + eps), & + dz(l)/(abs(vel(3)) + c + eps)) end if else !2D - cfl_terms = min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c)) + cfl_terms = min(dx(j)/(abs(vel(1)) + c + eps), & + dy(k)/(abs(vel(2)) + c + eps)) end if end function f_compute_multidim_cfl_terms Suggestion importance[1-10]: 8 __ Why: The suggestion correctly identifies a potential division-by-zero issue if abs(vel(i)) + c is zero. Since c (speed of sound) and vel (velocity) are physical quantities that can be zero, adding a small epsilon to the denominator is a robust way to prevent runtime errors and improve numerical stability.	Medium
Update

---

Previous suggestions

✅ Suggestions up to commit 47fa3d7

Category	Suggestion	Impact
General	✅ Extract complex inline expression Suggestion Impact: The suggestion was implemented by creating a dedicated function f_compute_filtered_dtheta that extracts the complex dtheta calculation. The commit goes beyond the suggestion by creating a separate function rather than just a local variable, but achieves the same goal of improving code readability and eliminating duplication. code diff: + !> Computes the modified dtheta for Fourier filtering in azimuthal direction + !! @param k y coordinate index !! @param l z coordinate index - !! @param for_dt_calc logical flag: true for dt calculation, false for stability calculation - !! @return cfl_terms computed CFL terms (with appropriate scaling applied) - pure function f_compute_cfl_terms(vel, c, j, k, l, for_dt_calc) result(cfl_terms) - !$acc routine seq - real(wp), dimension(num_vels), intent(in) :: vel - real(wp), intent(in) :: c - integer, intent(in) :: j, k, l - logical, intent(in) :: for_dt_calc - real(wp) :: cfl_terms + !! @return fltr_dtheta Modified dtheta value for cylindrical coordinates + pure function f_compute_filtered_dtheta(k, l) result(fltr_dtheta) + !$acc routine seq + integer, intent(in) :: k, l real(wp) :: fltr_dtheta integer :: Nfq - ! Compute filtered dtheta for cylindrical coordinates if (grid_geometry == 3) then if (k == 0) then fltr_dtheta = 2._wp*pi*y_cb(0)/3._wp @@ -46,38 +33,45 @@ else fltr_dtheta = y_cb(k - 1)*dz(l) end if - end if - - ! Compute CFL terms based on dimensionality + else + fltr_dtheta = 0._wp + end if + end function f_compute_filtered_dtheta The inline expression for filtered dtheta is complex and duplicated. Extract this calculation to a local variable for better readability and maintainability. src/simulation/m_sim_helpers.f90 [181-182] -vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k), & - 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)))**2._wp +fltr_dtheta = 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)) +vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k), fltr_dtheta)**2._wp Suggestion importance[1-10]: 7 __ Why: The suggestion correctly identifies a complex, hard-to-read inline expression for fltr_dtheta. Extracting this logic into a local variable, as proposed, significantly improves code readability and maintainability, which aligns with the PR's refactoring goal.	Medium
	✅ Extract duplicated complex calculation Suggestion Impact: The suggestion was implemented by creating a dedicated function f_compute_filtered_dtheta that extracts the complex filtered dtheta calculation. This function is then called in multiple places to replace the duplicated code, exactly addressing the maintainability concern raised in the suggestion. code diff: + !> Computes the modified dtheta for Fourier filtering in azimuthal direction + !! @param k y coordinate index !! @param l z coordinate index - !! @param for_dt_calc logical flag: true for dt calculation, false for stability calculation - !! @return cfl_terms computed CFL terms (with appropriate scaling applied) - pure function f_compute_cfl_terms(vel, c, j, k, l, for_dt_calc) result(cfl_terms) - !$acc routine seq - real(wp), dimension(num_vels), intent(in) :: vel - real(wp), intent(in) :: c - integer, intent(in) :: j, k, l - logical, intent(in) :: for_dt_calc - real(wp) :: cfl_terms + !! @return fltr_dtheta Modified dtheta value for cylindrical coordinates + pure function f_compute_filtered_dtheta(k, l) result(fltr_dtheta) + !$acc routine seq + integer, intent(in) :: k, l real(wp) :: fltr_dtheta integer :: Nfq - ! Compute filtered dtheta for cylindrical coordinates if (grid_geometry == 3) then if (k == 0) then fltr_dtheta = 2._wp*pi*y_cb(0)/3._wp @@ -46,38 +33,45 @@ else fltr_dtheta = y_cb(k - 1)*dz(l) end if - end if - - ! Compute CFL terms based on dimensionality + else + fltr_dtheta = 0._wp + end if + end function f_compute_filtered_dtheta + + !> Computes inviscid CFL terms for multi-dimensional cases (2D/3D only) + !! @param vel directional velocities + !! @param c mixture speed of sound + !! @param j x coordinate index + !! @param k y coordinate index + !! @param l z coordinate index + !! @return cfl_terms computed CFL terms for 2D/3D cases + pure function f_compute_multidim_cfl_terms(vel, c, j, k, l) result(cfl_terms) + !$acc routine seq + real(wp), dimension(num_vels), intent(in) :: vel + real(wp), intent(in) :: c + integer, intent(in) :: j, k, l + real(wp) :: cfl_terms + real(wp) :: fltr_dtheta + + fltr_dtheta = f_compute_filtered_dtheta(k, l) + if (p > 0) then !3D if (grid_geometry == 3) then cfl_terms = min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - fltr_dtheta/(abs(vel(3)) + c)) + dy(k)/(abs(vel(2)) + c), & + fltr_dtheta/(abs(vel(3)) + c)) else cfl_terms = min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - dz(l)/(abs(vel(3)) + c)) - end if - elseif (n > 0) then + dy(k)/(abs(vel(2)) + c), & + dz(l)/(abs(vel(3)) + c)) + end if + else !2D cfl_terms = min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c)) - else - !1D - special handling for different calculation types - if (for_dt_calc) then - cfl_terms = dx(j)/(abs(vel(1)) + c) - else - cfl_terms = (1._wp/dx(j))*(abs(vel(1)) + c) - end if - end if - - ! Apply appropriate factor for stability vs dt calculation - if (.not. for_dt_calc .and. (p > 0 .or. n > 0)) then - cfl_terms = 1._wp/cfl_terms - end if - end function f_compute_cfl_terms + dy(k)/(abs(vel(2)) + c)) + end if + end function f_compute_multidim_cfl_terms !> Computes enthalpy !! @param q_prim_vf cell centered primitive variables @@ -170,26 +164,44 @@ real(wp), dimension(2), intent(in) :: Re_l integer, intent(in) :: j, k, l + real(wp) :: fltr_dtheta + ! Inviscid CFL calculation - icfl_sf(j, k, l) = dt*f_compute_cfl_terms(vel, c, j, k, l, .false.) - - ! Viscous calculations (simplified with common patterns) + if (p > 0 .or. n > 0) then + ! 2D/3D cases + icfl_sf(j, k, l) = dt/f_compute_multidim_cfl_terms(vel, c, j, k, l) + else + ! 1D case - exact original logic + icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c) + end if + + ! Viscous calculations if (viscous) then + fltr_dtheta = f_compute_filtered_dtheta(k, l) + if (p > 0) then !3D if (grid_geometry == 3) then - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k), & - 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)))**2._wp - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), dy(k)*(abs(vel(2)) + c), & - 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp))*(abs(vel(3)) + c))/maxval(1._wp/Re_l) + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & + /min(dx(j), dy(k), fltr_dtheta)**2._wp + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c), & + fltr_dtheta*(abs(vel(3)) + c)) & + /maxval(1._wp/Re_l) else - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k), dz(l))**2._wp - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), dy(k)*(abs(vel(2)) + c), dz(l)*(abs(vel(3)) + c))/maxval(1._wp/Re_l) + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & + /min(dx(j), dy(k), dz(l))**2._wp + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c), & + dz(l)*(abs(vel(3)) + c)) & + /maxval(1._wp/Re_l) end if elseif (n > 0) then !2D vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2._wp - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), dy(k)*(abs(vel(2)) + c))/maxval(1._wp/Re_l) + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c)) & + /maxval(1._wp/Re_l) else !1D vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2._wp @@ -216,26 +228,36 @@ integer, intent(in) :: j, k, l real(wp) :: icfl_dt, vcfl_dt + real(wp) :: fltr_dtheta ! Inviscid CFL calculation - icfl_dt = cfl_target*f_compute_cfl_terms(vel, c, j, k, l, .true.) - - ! Viscous calculations (simplified with common patterns) + if (p > 0 .or. n > 0) then + ! 2D/3D cases + icfl_dt = cfl_target*f_compute_multidim_cfl_terms(vel, c, j, k, l) + else + ! 1D case - exact original logic + icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) + end if + + ! Viscous calculations if (viscous) then + fltr_dtheta = f_compute_filtered_dtheta(k, l) + if (p > 0) then !3D if (grid_geometry == 3) then - vcfl_dt = cfl_target*min(dx(j), dy(k), & - 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)))**2._wp/minval(1/(rho*Re_l)) + vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2._wp) & + /minval(1/(rho*Re_l)) The same complex filtered dtheta calculation is duplicated here. Extract this to a local variable to improve code maintainability and reduce duplication. src/simulation/m_sim_helpers.f90 [228-229] -vcfl_dt = cfl_target*min(dx(j), dy(k), & - 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)))**2._wp/minval(1/(rho*Re_l)) +fltr_dtheta = 2._wp*pi*y_cb(max(0,k-1))/max(3._wp, real(min(floor(2._wp*real(max(1,k), wp)*pi), (p + 1)/2 + 1), wp)) +vcfl_dt = cfl_target*min(dx(j), dy(k), fltr_dtheta)**2._wp/minval(1/(rho*Re_l)) Suggestion importance[1-10]: 7 __ Why: This suggestion correctly points out that the same complex expression from the previous suggestion is duplicated in another subroutine. Extracting it to a local variable improves readability and maintainability by reducing code duplication.	Medium

[sbryngelson]

/improve

[codecov]

Codecov Report

Attention: Patch coverage is 40.00000% with 18 lines in your changes missing coverage. Please review.

Project coverage is 45.82%. Comparing base (a34f0ca) to head (6132923).
Report is 3 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_sim_helpers.f90	40.00%	9 Missing and 9 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #901      +/-   ##
==========================================
+ Coverage   45.77%   45.82%   +0.05%     
==========================================
  Files          68       68              
  Lines       18664    18656       -8     
  Branches     2251     2247       -4     
==========================================
+ Hits         8543     8550       +7     
+ Misses       8763     8753      -10     
+ Partials     1358     1353       -5     

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

[Copilot]

Pull Request Overview

This PR refactors helper routines for stability, dt, and CFL calculations to reduce code duplication and improve maintainability. Key changes include the extraction of common CFL logic into the new function f_compute_multidim_cfl_terms, the introduction of f_compute_filtered_dtheta for Fourier filtering adjustments, and updates to inline expressions in viscous CFL calculations.

Reviewed Changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 1 comment.

File	Description
src/simulation/m_sim_helpers.f90	Refactors CFL calculation code by introducing helper functions and consolidating repeated logic in stability and dt routines.
.github/workflows/bench.yml	Updates the conditional for running benchmarks to include additional criteria for pull request events.

Comments suppressed due to low confidence (1)

src/simulation/m_sim_helpers.f90:170

• [nitpick] The check 'if (p > 0 .or. n > 0)' differentiates the 2D/3D and 1D cases but may be unclear to the reader. Consider adding a comment to explain the meaning and origin of the variables 'p' and 'n' in this context.

        if (p > 0 .or. n > 0) then