Refactor m_riemann_solvers Module (Non-Solver Subroutines)

[Malmahrouqi3]

User description

Description

Reduced PR from (#855),

Essentially, I reduced non-math-critical subroutines slightly. In addition, I added s_compute_wave_speed as a comprehensive subroutine to tackle its calcs for all Riemann Solvers (HLL, HLLC, HLLD).

---

PR Type

Enhancement

---

Description

• Extracted wave speed computation into centralized s_compute_wave_speed subroutine

• Consolidated duplicate Riemann buffer population logic across directions

• Refactored output data reshaping to reduce code duplication

• Added debug validation for wave speed calculations

---

Diagram Walkthrough

flowchart LR
  A["Duplicate wave speed code"] --> B["s_compute_wave_speed subroutine"]
  C["Direction-specific buffer logic"] --> D["Unified buffer population"]
  E["Repeated output reshaping"] --> F["Consolidated reshaping loops"]
  B --> G["All Riemann solvers"]
  D --> G
  F --> G

Loading

File Walkthrough

Relevant files

Enhancement

m_variables_conversion.fpp Add centralized wave speed computation subroutine                src/common/m_variables_conversion.fpp Added new s_compute_wave_speed subroutine with comprehensive wave speed calculations Supports all Riemann solver types (HLL, HLLC, HLLD) and physics models Includes debug validation for wave speed consistency Added function export to module interface +115/-0  m_riemann_solvers.fpp Refactor Riemann solvers with unified computation patterns src/simulation/m_riemann_solvers.fpp Replaced duplicate wave speed calculations with calls to s_compute_wave_speed Consolidated direction-specific buffer population into unified logic Refactored output data reshaping to eliminate code duplication Updated variable declarations and intent specifications Added pointer-based approach for direction-agnostic operations +177/-697

---

[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 new s_compute_wave_speed subroutine has complex conditional logic with multiple physics models (elasticity, mhd, hypoelasticity, hyperelasticity) but lacks comprehensive validation. The wave speed calculations involve intricate mathematical formulations that could introduce numerical errors or incorrect physics handling. subroutine s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & c_L, c_R, c_avg, c_fast_L, c_fast_R, G_L, G_R, & tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, idx, idx_tau) ! Computes the wave speeds for the Riemann solver #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_wave_speed #else !$acc routine seq #endif ! Input parameters integer, intent(in) :: wave_speeds integer, intent(in) :: idx, idx_tau real(wp), intent(in) :: rho_L, rho_R real(wp), dimension(:), intent(in) :: vel_L, vel_R, tau_e_L, tau_e_R real(wp), intent(in) :: pres_L, pres_R, c_L, c_R real(wp), intent(in) :: gamma_L, gamma_R, pi_inf_L, pi_inf_R real(wp), intent(in) :: rho_avg, c_avg real(wp), intent(in) :: c_fast_L, c_fast_R real(wp), intent(in) :: G_L, G_R ! Local variables real(wp) :: pres_SL, pres_SR, Ms_L, Ms_R ! Output parameters real(wp), intent(out) :: s_L, s_R, s_S, s_M, s_P if (wave_speeds == 1) then if (elasticity) then s_L = min(vel_L(idx) - sqrt(c_L*c_L + & (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L), vel_R(idx) - sqrt(c_R*c_R + & (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R)) s_R = max(vel_R(idx) + sqrt(c_R*c_R + & (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R), vel_L(idx) + sqrt(c_L*c_L + & (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L)) s_S = (pres_R - tau_e_R(idx_tau) - pres_L + & tau_e_L(idx_tau) + rho_L*vel_L(idx)*(s_L - vel_L(idx)) - & rho_R*vel_R(idx)*(s_R - vel_R(idx)))/(rho_L*(s_L - vel_L(idx)) - & rho_R*(s_R - vel_R(idx))) else if (mhd) then s_L = min(vel_L(idx) - c_fast_L, vel_R(idx) - c_fast_R) s_R = max(vel_R(idx) + c_fast_R, vel_L(idx) + c_fast_L) s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) else if (hypoelasticity) then s_L = min(vel_L(idx) - sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + & tau_e_L(idx_tau))/rho_L) & , vel_R(idx) - sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + & tau_e_R(idx_tau))/rho_R)) s_R = max(vel_R(idx) + sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + & tau_e_R(idx_tau))/rho_R) & , vel_L(idx) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + & tau_e_L(idx_tau))/rho_L)) s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) else if (hyperelasticity) then s_L = min(vel_L(idx) - sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L) & , vel_R(idx) - sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R)) s_R = max(vel_R(idx) + sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R) & , vel_L(idx) + sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L)) s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) else s_L = min(vel_L(idx) - c_L, vel_R(idx) - c_R) s_R = max(vel_R(idx) + c_R, vel_L(idx) + c_L) s_S = (pres_R - pres_L + rho_L*vel_L(idx)* & (s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) & /(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx))) end if else if (wave_speeds == 2) then pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg*(vel_L(idx) - vel_R(idx))) pres_SR = pres_SL Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(idx) - c_L*Ms_L s_R = vel_R(idx) + c_R*Ms_R s_S = 5e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg)) end if ! ! follows Einfeldt et al. ! s_M/P = min/max(0.,s_L/R) s_M = min(0._wp, s_L) s_P = max(0._wp, s_R) #ifdef DEBUG ! Check for potential issues in wave speed calculation if (s_R <= s_L) then print *, 'WARNING: Wave speed issue detected in s_compute_wave_speed' print *, 'Left wave speed >= Right wave speed:', s_L, s_R print *, 'Input velocities :', vel_L(idx), vel_R(idx) print *, 'Sound speeds:', c_L, c_R print *, 'Densities:', rho_L, rho_R print *, 'Pressures:', pres_L, pres_R print *, 'Wave speeds method:', wave_speeds if (elasticity .or. hypoelasticity .or. hyperelasticity) then print *, 'Shear moduli:', G_L, G_R end if call s_mpi_abort('Error: Invalid wave speeds in s_compute_wave_speed') end if #endif end subroutine s_compute_wave_speed #endif Code Duplication After the refactoring, there are still redundant wave speed calculations present. Line 690 shows s_M and s_P being recalculated immediately after calling s_compute_wave_speed, which already computes these values. This suggests incomplete refactoring and potential inconsistencies. call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) Logic Error In the buffer population logic, there are undefined variables like 'j' being used in conditional blocks without proper declaration or loop context. This could lead to compilation errors or undefined behavior during runtime. dqL_prim_dx_vf(i)%sf(j, -1, l) = dqR_prim_dx_vf(i)%sf(j, 0, l) dqL_prim_dy_vf(i)%sf(j, -1, l) = dqR_prim_dy_vf(i)%sf(j, 0, l) if (p > 0) then dqL_prim_dz_vf(i)%sf(j, -1, l) = dqR_prim_dz_vf(i)%sf(j, 0, l) end if else dqL_prim_dx_vf(i)%sf(j, k, -1) = dqR_prim_dx_vf(i)%sf(j, k, 0) dqL_prim_dy_vf(i)%sf(j, k, -1) = dqR_prim_dy_vf(i)%sf(j, k, 0) dqL_prim_dz_vf(i)%sf(j, k, -1) = dqR_prim_dz_vf(i)%sf(j, k, 0) end if

[Copilot]

Pull Request Overview

This PR refactors the m_riemann_solvers module by extracting all inline wave‐speed logic into a new s_compute_wave_speed subroutine, and cleans up buffer initializations with pointers and target attributes.

• Adds s_compute_wave_speed (in m_variables_conversion) to centralize all HLL/HLLC/HLLD wave‐speed formulas.
• Updates m_riemann_solvers.fpp to call s_compute_wave_speed instead of duplicating logic, and refactors boundary‐buffer routines using pointer/target.
• Adjusts variable declarations (intent, allocatable, pointer, new bc_beg, bc_end, end_val) for directional BC handling.

Reviewed Changes

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

File	Description
src/simulation/m_riemann_solvers.fpp	Replaced duplicated wave‐speed blocks with calls to s_compute_wave_speed; refactored Riemann buffer loops with pointers/targets and unified BC handling
src/common/m_variables_conversion.fpp	Exported and implemented new s_compute_wave_speed subroutine under #ifndef MFC_PRE_PROCESS

Comments suppressed due to low confidence (2)

src/common/m_variables_conversion.fpp:1714

• There’s no !> docstring for s_compute_wave_speed. Add a header explaining the purpose, inputs, and outputs, so future readers can quickly understand the routine.

    subroutine s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, &

src/common/m_variables_conversion.fpp:1714

• The newly introduced s_compute_wave_speed contains complex branching for multiple wave‐speed methods but no dedicated unit tests. Consider adding focused tests covering both wave_speeds==1 and wave_speeds==2 paths.

    subroutine s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, &

[qodo-merge-pro]

PR Code Suggestions ✨

Latest suggestions up to dd158d0

Category	Suggestion	Impact
Possible issue	Fix wave speed validation condition The condition s_R <= s_L should be s_R < s_L to properly detect invalid wave speeds. Wave speeds can be equal in degenerate cases without being invalid. The current condition may trigger false positives when s_R = s_L. src/common/m_variables_conversion.fpp [1809-1821] -if (s_R <= s_L) then +if (s_R < s_L) then print *, 'WARNING: Wave speed issue detected in s_compute_wave_speed' - print *, 'Left wave speed >= Right wave speed:', s_L, s_R + print *, 'Left wave speed > Right wave speed:', s_L, s_R print *, 'Input velocities :', vel_L(idx), vel_R(idx) print *, 'Sound speeds:', c_L, c_R print *, 'Densities:', rho_L, rho_R print *, 'Pressures:', pres_L, pres_R print *, 'Wave speeds method:', wave_speeds if (elasticity .or. hypoelasticity .or. hyperelasticity) then print *, 'Shear moduli:', G_L, G_R end if call s_mpi_abort('Error: Invalid wave speeds in s_compute_wave_speed') end if Suggestion importance[1-10]: 6 __ Why: The suggestion correctly identifies that the check s_R <= s_L could cause a false positive abort when s_R == s_L, which can be a valid physical state, thus improving the robustness of the debug validation.	Low
Possible issue	Prevent division by zero error Add validation to prevent division by zero in the s_S calculation. The denominator rho_L(s_L - vel_L(idx)) - rho_R(s_R - vel_R(idx)) could become zero in degenerate cases, leading to numerical instability.** src/common/m_variables_conversion.fpp [1743-1754] if (wave_speeds == 1) then if (elasticity) then s_L = min(vel_L(idx) - sqrt(c_L*c_L + & (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L), vel_R(idx) - sqrt(c_R*c_R + & (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R)) s_R = max(vel_R(idx) + sqrt(c_R*c_R + & (((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R), vel_L(idx) + sqrt(c_L*c_L + & (((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L)) - s_S = (pres_R - tau_e_R(idx_tau) - pres_L + & - tau_e_L(idx_tau) + rho_L*vel_L(idx)*(s_L - vel_L(idx)) - & - rho_R*vel_R(idx)*(s_R - vel_R(idx)))/(rho_L*(s_L - vel_L(idx)) - & - rho_R*(s_R - vel_R(idx))) + + real(wp) :: denom + denom = rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx)) + if (abs(denom) < epsilon(1.0_wp)) then + s_S = 0.5_wp * (vel_L(idx) + vel_R(idx)) + else + s_S = (pres_R - tau_e_R(idx_tau) - pres_L + & + tau_e_L(idx_tau) + rho_L*vel_L(idx)*(s_L - vel_L(idx)) - & + rho_R*vel_R(idx)*(s_R - vel_R(idx)))/denom + end if Suggestion importance[1-10]: 8 __ Why: This is a valid and important suggestion to prevent a potential division-by-zero runtime error, which could cause a crash, by adding a check on the denominator.	Medium
	Add input validation for norm_dir Add validation to ensure norm_dir is within valid range (1-3) before using it in the conditional logic. Invalid values could cause undefined behavior or incorrect pointer assignments. src/simulation/m_riemann_solvers.fpp [3177-3204] +if (norm_dir < 1 .or. norm_dir > 3) then + call s_mpi_abort('Error: Invalid norm_dir value in s_populate_riemann_states_variables_buffers') +end if + if (norm_dir == 1) then is1 = ix; is2 = iy; is3 = iz dir_idx = (/1, 2, 3/); dir_flg = (/1._wp, 0._wp, 0._wp/) bc_beg = bc_x%beg; bc_end = bc_x%end end_val = m qL_prim_rs_vf => qL_prim_rsx_vf qR_prim_rs_vf => qR_prim_rsx_vf dqL_prim_d_vf => dqL_prim_dx_vf dqR_prim_d_vf => dqR_prim_dx_vf else if (norm_dir == 2) then is1 = iy; is2 = ix; is3 = iz dir_idx = (/2, 1, 3/); dir_flg = (/0._wp, 1._wp, 0._wp/) bc_beg = bc_y%beg; bc_end = bc_y%end end_val = n qL_prim_rs_vf => qL_prim_rsy_vf qR_prim_rs_vf => qR_prim_rsy_vf dqL_prim_d_vf => dqL_prim_dy_vf dqR_prim_d_vf => dqR_prim_dy_vf else is1 = iz; is2 = iy; is3 = ix dir_idx = (/3, 1, 2/); dir_flg = (/0._wp, 0._wp, 1._wp/) bc_beg = bc_z%beg; bc_end = bc_z%end end_val = p qL_prim_rs_vf => qL_prim_rsz_vf qR_prim_rs_vf => qR_prim_rsz_vf dqL_prim_d_vf => dqL_prim_dz_vf dqR_prim_d_vf => dqR_prim_dz_vf end if Suggestion importance[1-10]: 7 __ Why: The suggestion correctly identifies that an invalid norm_dir would fall into the else block, causing incorrect pointer assignments; adding a check improves code robustness.	Medium
Update

---

Previous suggestions

✅ Suggestions up to commit b684a4e

Category	Suggestion	Impact
Possible issue	Fix wave speed validation condition The debug check uses s_R <= s_L which will trigger for valid cases where wave speeds are equal. This condition should be s_R < s_L to only catch truly invalid configurations where the right wave speed is strictly less than the left wave speed. src/common/m_variables_conversion.fpp [1809-1822] -if (s_R <= s_L) then +if (s_R < s_L) then print *, 'WARNING: Wave speed issue detected in s_compute_wave_speed' - print *, 'Left wave speed >= Right wave speed:', s_L, s_R + print *, 'Left wave speed > Right wave speed:', s_L, s_R print *, 'Input velocities :', vel_L(idx), vel_R(idx) print *, 'Sound speeds:', c_L, c_R print *, 'Densities:', rho_L, rho_R print *, 'Pressures:', pres_L, pres_R print *, 'Wave speeds method:', wave_speeds if (elasticity .or. hypoelasticity .or. hyperelasticity) then print *, 'Shear moduli:', G_L, G_R end if call s_mpi_abort('Error: Invalid wave speeds in s_compute_wave_speed') end if Suggestion importance[1-10]: 7 __ Why: The suggestion correctly identifies that the debug check s_R <= s_L is too strict and would cause an abort in valid cases where s_R == s_L, improving the robustness of the validation logic.	Medium
General	✅ Remove redundant wave speed calculations Suggestion Impact: The commit removed multiple instances of redundant s_M and s_P calculations throughout the code, exactly as suggested. The lines "s_M = min(0._wp, s_L); s_P = max(0._wp, s_R)" were removed from several locations. code diff: - - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_L)) & + (5.e-1_wp - sign(5.e-1_wp, s_L)) & @@ -1445,10 +1443,6 @@ tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - ! goes with q_star_L/R = xi_L/R * (variable) ! xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) xi_L = (s_L - vel_L(idx1))/(s_L - s_S) @@ -1717,10 +1711,6 @@ c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) ! xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) @@ -2120,10 +2110,6 @@ tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - ! goes with q_star_L/R = xi_L/R * (variable) ! xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) @@ -2562,10 +2548,6 @@ tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1)) - ! follows Einfeldt et al. - ! s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) This line duplicates the wave speed calculation that is already performed inside s_compute_wave_speed. The s_M and s_P values are already computed and returned by the subroutine, making this redundant calculation unnecessary. src/simulation/m_riemann_solvers.fpp [690] -s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) +! s_M and s_P are already computed in s_compute_wave_speed Suggestion importance[1-10]: 6 __ Why: The suggestion correctly points out that the calculation of s_M and s_P is redundant, as these values are now computed and returned by the new s_compute_wave_speed subroutine.	Low

[Malmahrouqi3]

Code Duplication
After the refactoring, there are still redundant wave speed calculations present. Line 690 shows s_M and s_P being recalculated immediately after calling s_compute_wave_speed, which already computes these values. This suggests incomplete refactoring and potential inconsistencies.

Good catch, I forgot about it honestly when replacing wave speed calcs with the subroutine call. Thnx.

[codecov]

Codecov Report

❌ Patch coverage is 40.50633% with 94 lines in your changes missing coverage. Please review.
✅ Project coverage is 40.99%. Comparing base (5e8f116) to head (6e7281b).
⚠️ Report is 2 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_riemann_solvers.fpp	43.13%	47 Missing and 11 partials ⚠️
src/common/m_variables_conversion.fpp	35.71%	34 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #908      +/-   ##
==========================================
+ Coverage   40.91%   40.99%   +0.08%     
==========================================
  Files          70       70              
  Lines       20288    20147     -141     
  Branches     2517     2500      -17     
==========================================
- Hits         8301     8260      -41     
+ Misses      10450    10366      -84     
+ Partials     1537     1521      -16     

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

[sbryngelson]

be sure to run tests on frontier and phoenix (amd and nvidia gpu) before pushing. helps you debug faster and saves the runners to do other things.

[Malmahrouqi3]

Roger that

[sbryngelson]

surely some of these arguments should be optional since not every riemann solver uses them all?

[Malmahrouqi3]

yup yup, not all of them will be used in every scenario. I will look into a way to reduce that as the subroutine call looks the same in every solver yet not neat-looking

                            
call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, &
                     c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, &
                     tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, &
                     s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1))

[sbryngelson]

Code Duplication
After the refactoring, there are still redundant wave speed calculations present. Line 690 shows s_M and s_P being recalculated immediately after calling s_compute_wave_speed, which already computes these values. This suggests incomplete refactoring and potential inconsistencies.

Good catch, I forgot about it honestly when replacing wave speed calcs with the subroutine call. Thnx.

@Malmahrouqi3
You don't have to reply to every Copilot or Qodo PR comment. They're just supposed to be handy to look at. I look too, but ignore most of them. But sometimes a useful comment or suggestion pops up.

[Malmahrouqi3]

Not frankly sure what is happening. I tried interactive where all tests passed vs. slurm where all tests failed for phoenix gpu. No specific error messages e.g. tolerance issue or such, yet prolly due to the lint issue. I will fix that and re-run the slurm job.

[sbryngelson]

The error is only on GPU cases and is because you did something with OpenACC incorrectly:

FATAL ERROR: data in PRESENT clause was not found on device 1: name=c_fast host:0x3d77820
 file:/storage/scratch1/6/sbryngelson3/mfc-runners/actions-runner-03/_work/MFC/MFC/src/simulation/m_riemann_solvers.fpp s_hllc_riemann_solver line:2287

[sbryngelson]

i don't think this is right

[Malmahrouqi3]

I added this initially like a while back as it gave me errors and refused to complete the build. More likely the reason tho, I will look into this.

[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 new debug check in s_compute_wave_speed unconditionally calls s_mpi_abort when s_R <= s_L, which may be too strict in borderline states or due to round-off, potentially causing unintended simulation aborts. Consider tolerances or logging-only behavior in DEBUG. #ifdef DEBUG ! Check for potential issues in wave speed calculation if (s_R <= s_L) then print *, 'WARNING: Wave speed issue detected in s_compute_wave_speed' print *, 'Left wave speed >= Right wave speed:', s_L, s_R print *, 'Input velocities :', vel_L(idx), vel_R(idx) print *, 'Sound speeds:', c_L, c_R print *, 'Densities:', rho_L, rho_R print *, 'Pressures:', pres_L, pres_R print *, 'Wave speeds method:', wave_speeds if (elasticity .or. hypoelasticity .or. hyperelasticity) then print *, 'Shear moduli:', G_L, G_R end if call s_mpi_abort('Error: Invalid wave speeds in s_compute_wave_speed') end if Boundary Indexing In unified buffer population, derivative copies use j in conditions for norm_dir 2 and 3 without ensuring j loops are active in the surrounding parallel region; verify that j is defined and iterates over correct bounds in those branches to avoid out-of-bounds or uninitialized index usage. $:GPU_PARALLEL_LOOP(collapse=3) do i = momxb, momxe do l = isz%beg, isz%end do k = isy%beg, isy%end if (norm_dir == 1) then dqL_prim_dx_vf(i)%sf(-1, k, l) = dqR_prim_dx_vf(i)%sf(0, k, l) if (n > 0) then dqL_prim_dy_vf(i)%sf(-1, k, l) = dqR_prim_dy_vf(i)%sf(0, k, l) if (p > 0) then dqL_prim_dz_vf(i)%sf(-1, k, l) = dqR_prim_dz_vf(i)%sf(0, k, l) end if end if else if (norm_dir == 2) then dqL_prim_dx_vf(i)%sf(j, -1, l) = dqR_prim_dx_vf(i)%sf(j, 0, l) dqL_prim_dy_vf(i)%sf(j, -1, l) = dqR_prim_dy_vf(i)%sf(j, 0, l) if (p > 0) then dqL_prim_dz_vf(i)%sf(j, -1, l) = dqR_prim_dz_vf(i)%sf(j, 0, l) end if else dqL_prim_dx_vf(i)%sf(j, k, -1) = dqR_prim_dx_vf(i)%sf(j, k, 0) dqL_prim_dy_vf(i)%sf(j, k, -1) = dqR_prim_dy_vf(i)%sf(j, k, 0) dqL_prim_dz_vf(i)%sf(j, k, -1) = dqR_prim_dz_vf(i)%sf(j, k, 0) end if end do end do end do end if end if if (bc_end == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at end $:GPU_PARALLEL_LOOP(collapse=3) do i = 1, sys_size do l = is3%beg, is3%end do k = is2%beg, is2%end qR_prim_rs_vf(end_val + 1, k, l, i) = qL_prim_rs_vf(end_val, k, l, i) end do end do end do if (viscous) then $:GPU_PARALLEL_LOOP(collapse=3) do i = momxb, momxe do l = isz%beg, isz%end do k = isy%beg, isy%end if (norm_dir == 1) then dqR_prim_dx_vf(i)%sf(end_val + 1, k, l) = dqL_prim_dx_vf(i)%sf(end_val, k, l) if (n > 0) then dqR_prim_dy_vf(i)%sf(end_val + 1, k, l) = dqL_prim_dy_vf(i)%sf(end_val, k, l) if (p > 0) then dqR_prim_dz_vf(i)%sf(end_val + 1, k, l) = dqL_prim_dz_vf(i)%sf(end_val, k, l) end if end if else if (norm_dir == 2) then dqR_prim_dx_vf(i)%sf(j, end_val + 1, l) = dqL_prim_dx_vf(i)%sf(j, end_val, l) dqR_prim_dy_vf(i)%sf(j, end_val + 1, l) = dqL_prim_dy_vf(i)%sf(j, end_val, l) if (p > 0) then dqR_prim_dz_vf(i)%sf(j, end_val + 1, l) = dqL_prim_dz_vf(i)%sf(j, end_val, l) end if else dqR_prim_dx_vf(i)%sf(j, k, end_val + 1) = dqL_prim_dx_vf(i)%sf(j, k, end_val) dqR_prim_dy_vf(i)%sf(j, k, end_val + 1) = dqL_prim_dy_vf(i)%sf(j, k, end_val) dqR_prim_dz_vf(i)%sf(j, k, end_val + 1) = dqL_prim_dz_vf(i)%sf(j, k, end_val) end if API Consistency s_compute_wave_speed signature requires many inputs (elasticity, mhd, hypoelasticity flags used inside). Ensure these globals are consistently available across call sites and that idx/idx_tau map correctly for all directions; otherwise speeds could be computed with wrong components. end if ! Low Mach correction if (low_Mach == 2) then @:compute_low_Mach_correction() end if ! COMPUTING THE DIRECT WAVE SPEEDS call s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, & c_L, c_R, c_avg, c_fast%L, c_fast%R, G_L, G_R, & tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, & s_L, s_R, s_S, s_M, s_P, dir_idx(1), dir_idx_tau(1))

[qodo-merge-pro]

PR Code Suggestions ✨

Explore these optional code suggestions:

Category	Suggestion	Impact
Possible issue	Fix undefined loop indices Variables j and k are used in the norm_dir 2/3 branches without being defined in the surrounding loops, which will lead to undefined indexing. Add appropriate loop indices for j/k in those branches to mirror the original per-direction population logic. src/simulation/m_riemann_solvers.fpp [3175-3212] -if (bc_beg == BC_RIEMANN_EXTRAP) then ! Riemann state extrap. BC at beginning +if (bc_beg == BC_RIEMANN_EXTRAP) then $:GPU_PARALLEL_LOOP(collapse=3) do i = 1, sys_size do l = is3%beg, is3%end do k = is2%beg, is2%end qL_prim_rs_vf(-1, k, l, i) = qR_prim_rs_vf(0, k, l, i) end do end do end do if (viscous) then - $:GPU_PARALLEL_LOOP(collapse=3) - do i = momxb, momxe - do l = isz%beg, isz%end - do k = isy%beg, isy%end - if (norm_dir == 1) then + if (norm_dir == 1) then + $:GPU_PARALLEL_LOOP(collapse=3) + do i = momxb, momxe + do l = isz%beg, isz%end + do k = isy%beg, isy%end dqL_prim_dx_vf(i)%sf(-1, k, l) = dqR_prim_dx_vf(i)%sf(0, k, l) - if (n > 0) then - dqL_prim_dy_vf(i)%sf(-1, k, l) = dqR_prim_dy_vf(i)%sf(0, k, l) - if (p > 0) then - dqL_prim_dz_vf(i)%sf(-1, k, l) = dqR_prim_dz_vf(i)%sf(0, k, l) - end if - end if - else if (norm_dir == 2) then + if (n > 0) dqL_prim_dy_vf(i)%sf(-1, k, l) = dqR_prim_dy_vf(i)%sf(0, k, l) + if (p > 0) dqL_prim_dz_vf(i)%sf(-1, k, l) = dqR_prim_dz_vf(i)%sf(0, k, l) + end do + end do + end do + else if (norm_dir == 2) then + $:GPU_PARALLEL_LOOP(collapse=3) + do i = momxb, momxe + do l = isz%beg, isz%end + do j = isx%beg, isx%end dqL_prim_dx_vf(i)%sf(j, -1, l) = dqR_prim_dx_vf(i)%sf(j, 0, l) dqL_prim_dy_vf(i)%sf(j, -1, l) = dqR_prim_dy_vf(i)%sf(j, 0, l) - if (p > 0) then - dqL_prim_dz_vf(i)%sf(j, -1, l) = dqR_prim_dz_vf(i)%sf(j, 0, l) - end if - else + if (p > 0) dqL_prim_dz_vf(i)%sf(j, -1, l) = dqR_prim_dz_vf(i)%sf(j, 0, l) + end do + end do + end do + else + $:GPU_PARALLEL_LOOP(collapse=3) + do i = momxb, momxe + do k = isy%beg, isy%end + do j = isx%beg, isx%end dqL_prim_dx_vf(i)%sf(j, k, -1) = dqR_prim_dx_vf(i)%sf(j, k, 0) dqL_prim_dy_vf(i)%sf(j, k, -1) = dqR_prim_dy_vf(i)%sf(j, k, 0) dqL_prim_dz_vf(i)%sf(j, k, -1) = dqR_prim_dz_vf(i)%sf(j, k, 0) - end if + end do end do end do - end do + end if end if end if Suggestion importance[1-10]: 10 __ Why: This suggestion identifies a critical bug where the loop indices j and k are used without being defined, which was introduced during refactoring and would lead to undefined behavior or a crash.	High
	Add safe-division floors Guard against divisions by zero or near-zero denominators in the pressure-based and contact wave calculations. Add small epsilon floors to pres_L, pres_R, and rho_avg*c_avg to prevent NaNs/Inf when pressures are tiny or averages degenerate. src/common/m_variables_conversion.fpp [1806-1818] else if (wave_speeds == 2) then + real(wp), parameter :: eps = 1.e-14_wp + real(wp) :: pres_L_eff, pres_R_eff, rc_eff + + pres_L_eff = max(pres_L, eps) + pres_R_eff = max(pres_R, eps) + rc_eff = max(rho_avg*c_avg, eps) + pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg*(vel_L(idx) - vel_R(idx))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & - (pres_SL/pres_L - 1._wp)*pres_L/ & - ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & - (pres_SR/pres_R - 1._wp)*pres_R/ & - ((pres_R + pi_inf_R/(1._wp + gamma_R))))) + + Ms_L = max(1._wp, sqrt(max(0._wp, 1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L)) * & + (pres_SL/pres_L_eff - 1._wp)*pres_L_eff / (pres_L_eff + pi_inf_L/(1._wp + gamma_L)) ))) + Ms_R = max(1._wp, sqrt(max(0._wp, 1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R)) * & + (pres_SR/pres_R_eff - 1._wp)*pres_R_eff / (pres_R_eff + pi_inf_R/(1._wp + gamma_R)) ))) + s_L = vel_L(idx) - c_L*Ms_L s_R = vel_R(idx) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg)) + s_S = 5e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/rc_eff) end if Suggestion importance[1-10]: 8 __ Why: The suggestion correctly identifies potential numerical instability from division by zero in the new s_compute_wave_speed subroutine and provides a robust fix, enhancing the reliability of the new code.	Medium
More