Refactor m_riemann_solvers Module (HLLD Solver)

[Malmahrouqi3]

User description

Description

Reduced PR from (#855),

Essentially, I refactored math-critical items in s_hlld_riemann_solver subroutine.

---

PR Type

Enhancement

---

Description

• Refactored HLLD Riemann solver for improved code clarity

• Replaced individual array assignments with array constructors

• Simplified magnetic field and state vector computations

• Consolidated flux calculations using vectorized operations

---

Changes diagram

flowchart LR
  A["Individual Array Assignments"] --> B["Array Constructors"]
  C["Separate Flux Calculations"] --> D["Vectorized Operations"]
  E["Complex State Vector Logic"] --> F["Simplified Computations"]
  B --> G["Cleaner HLLD Solver"]
  D --> G
  F --> G

Loading

---

Changes walkthrough 📝

	Relevant files
Enhancement	m_riemann_solvers.fpp HLLD solver refactoring with array operations                        src/simulation/m_riemann_solvers.fpp Replaced individual magnetic field assignments with array constructors Simplified state vector computations using vectorized operations Consolidated flux calculations with array slicing syntax Streamlined double-star state vector assignments +35/-92

---

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-code-review]

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

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

🧪 No relevant tests

🔒 No security concerns identified

⚡ Recommended focus areas for review Array Slicing The new array slicing syntax for flux assignments may not be compatible with all Fortran compilers or may have different behavior than individual assignments. This should be tested to ensure numerical equivalence. flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1), contxe + dir_idx(2), contxe + dir_idx(3)]) = F_hlld([2, 3, 4]) ! Magnetic field if (n == 0) then flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg, B_idx%beg + 1]) = F_hlld([5, 6]) else flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg + dir_idx(2) - 1, B_idx%beg + dir_idx(3) - 1]) = F_hlld([5, 6]) Missing Tests This refactoring of a critical HLLD Riemann solver lacks accompanying tests to verify that the mathematical operations produce identical results to the original implementation. B%L = [Bx0, qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg), qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1)] B%R = [Bx0, qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg), qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + 1)] else ! 2D/3D: Bx, By, Bz as variables B%L = [qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1) - 1), & qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1), & qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1)] B%R = [qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(1) - 1), & qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(2) - 1), & qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(3) - 1)] end if end if ! Sum properties of all fluid components rho%L = 0._wp; gamma%L = 0._wp; pi_inf%L = 0._wp; qv%L = 0._wp rho%R = 0._wp; gamma%R = 0._wp; pi_inf%R = 0._wp; qv%R = 0._wp !$acc loop seq do i = 1, num_fluids rho%L = rho%L + alpha_rho_L(i) gamma%L = gamma%L + alpha_L(i)*gammas(i) pi_inf%L = pi_inf%L + alpha_L(i)*pi_infs(i) qv%L = qv%L + alpha_rho_L(i)*qvs(i) rho%R = rho%R + alpha_rho_R(i) gamma%R = gamma%R + alpha_R(i)*gammas(i) pi_inf%R = pi_inf%R + alpha_R(i)*pi_infs(i) qv%R = qv%R + alpha_rho_R(i)*qvs(i) end do pres_mag%L = 0.5_wp*sum(B%L**2._wp) pres_mag%R = 0.5_wp*sum(B%R**2._wp) E%L = gamma%L*pres%L + pi_inf%L + 0.5_wp*rho%L*vel_rms%L + qv%L + pres_mag%L E%R = gamma%R*pres%R + pi_inf%R + 0.5_wp*rho%R*vel_rms%R + qv%R + pres_mag%R ! includes magnetic energy H_no_mag%L = (E%L + pres%L - pres_mag%L)/rho%L H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) ! (2) Compute fast wave speeds call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, 0._wp, c%L) call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, 0._wp, c%R) call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L) call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R) ! (3) Compute contact speed s_M [Miyoshi Equ. (38)] s_L = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R) s_R = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L) pTot_L = pres%L + pres_mag%L pTot_R = pres%R + pres_mag%R s_M = (((s_R - vel%R(1))*rho%R*vel%R(1) - & (s_L - vel%L(1))*rho%L*vel%L(1) - pTot_R + pTot_L)/ & ((s_R - vel%R(1))*rho%R - (s_L - vel%L(1))*rho%L)) ! (4) Compute star state variables rhoL_star = rho%L*(s_L - vel%L(1))/(s_L - s_M) rhoR_star = rho%R*(s_R - vel%R(1))/(s_R - s_M) p_star = pTot_L + rho%L*(s_L - vel%L(1))*(s_M - vel%L(1))/(s_L - s_M) E_starL = ((s_L - vel%L(1))*E%L - pTot_L*vel%L(1) + p_star*s_M)/(s_L - s_M) E_starR = ((s_R - vel%R(1))*E%R - pTot_R*vel%R(1) + p_star*s_M)/(s_R - s_M) ! (5) Compute left/right state vectors and fluxes U_L = [rho%L, rho%L*vel%L(1:3), B%L(2:3), E%L] U_starL = [rhoL_star, rhoL_star*s_M, rhoL_star*vel%L(2:3), B%L(2:3), E_starL] U_R = [rho%R, rho%R*vel%R(1:3), B%R(2:3), E%R] U_starR = [rhoR_star, rhoR_star*s_M, rhoR_star*vel%R(2:3), B%R(2:3), E_starR] ! Compute the left/right fluxes F_L(1) = U_L(2) F_L(2) = U_L(2)*vel%L(1) - B%L(1)*B%L(1) + pTot_L F_L(3:4) = U_L(2)*vel%L(2:3) - B%L(1)*B%L(2:3) F_L(5:6) = vel%L(1)*B%L(2:3) - vel%L(2:3)*B%L(1) F_L(7) = (E%L + pTot_L)*vel%L(1) - B%L(1)*(vel%L(1)*B%L(1) + vel%L(2)*B%L(2) + vel%L(3)*B%L(3)) F_R(1) = U_R(2) F_R(2) = U_R(2)*vel%R(1) - B%R(1)*B%R(1) + pTot_R F_R(3:4) = U_R(2)*vel%R(2:3) - B%R(1)*B%R(2:3) F_R(5:6) = vel%R(1)*B%R(2:3) - vel%R(2:3)*B%R(1) F_R(7) = (E%R + pTot_R)*vel%R(1) - B%R(1)*(vel%R(1)*B%R(1) + vel%R(2)*B%R(2) + vel%R(3)*B%R(3)) ! Compute the star flux using HLL relation F_starL = F_L + s_L*(U_starL - U_L) F_starR = F_R + s_R*(U_starR - U_R) ! Compute the rotational (Alfvén) speeds s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star) s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star) ! Compute the double–star states [Miyoshi Eqns. (59)-(62)] sqrt_rhoL_star = sqrt(rhoL_star); sqrt_rhoR_star = sqrt(rhoR_star) vL_star = vel%L(2); wL_star = vel%L(3) vR_star = vel%R(2); wR_star = vel%R(3) ! (6) Compute the double–star states [Miyoshi Eqns. (59)-(62)] denom_ds = sqrt_rhoL_star + sqrt_rhoR_star sign_Bx = sign(1._wp, B%L(1)) v_double = (sqrt_rhoL_star*vL_star + sqrt_rhoR_star*vR_star + (B%R(2) - B%L(2))*sign_Bx)/denom_ds w_double = (sqrt_rhoL_star*wL_star + sqrt_rhoR_star*wR_star + (B%R(3) - B%L(3))*sign_Bx)/denom_ds By_double = (sqrt_rhoL_star*B%R(2) + sqrt_rhoR_star*B%L(2) + sqrt_rhoL_star*sqrt_rhoR_star*(vR_star - vL_star)*sign_Bx)/denom_ds Bz_double = (sqrt_rhoL_star*B%R(3) + sqrt_rhoR_star*B%L(3) + sqrt_rhoL_star*sqrt_rhoR_star*(wR_star - wL_star)*sign_Bx)/denom_ds E_doubleL = E_starL - sqrt_rhoL_star*((vL_star*B%L(2) + wL_star*B%L(3)) - (v_double*By_double + w_double*Bz_double))*sign_Bx E_doubleR = E_starR + sqrt_rhoR_star*((vR_star*B%R(2) + wR_star*B%R(3)) - (v_double*By_double + w_double*Bz_double))*sign_Bx E_double = 0.5_wp*(E_doubleL + E_doubleR) U_doubleL = [rhoL_star, rhoL_star*s_M, rhoL_star*v_double, rhoL_star*w_double, By_double, Bz_double, E_double] U_doubleR = [rhoR_star, rhoR_star*s_M, rhoR_star*v_double, rhoR_star*w_double, By_double, Bz_double, E_double] ! (11) Choose HLLD flux based on wave-speed regions if (0.0_wp <= s_L) then F_hlld = F_L else if (0.0_wp <= s_starL) then F_hlld = F_L + s_L*(U_starL - U_L) else if (0.0_wp <= s_M) then F_hlld = F_starL + s_starL*(U_doubleL - U_starL) else if (0.0_wp <= s_starR) then F_hlld = F_starR + s_starR*(U_doubleR - U_starR) else if (0.0_wp <= s_R) then F_hlld = F_R + s_R*(U_starR - U_R) else F_hlld = F_R end if ! (12) Reorder and write temporary variables to the flux array ! Mass flux_rs${XYZ}$_vf(j, k, l, 1) = F_hlld(1) ! TODO multi-component ! Momentum flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1), contxe + dir_idx(2), contxe + dir_idx(3)]) = F_hlld([2, 3, 4]) ! Magnetic field if (n == 0) then flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg, B_idx%beg + 1]) = F_hlld([5, 6]) else flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg + dir_idx(2) - 1, B_idx%beg + dir_idx(3) - 1]) = F_hlld([5, 6]) end if ! Energy flux_rs${XYZ}$_vf(j, k, l, E_idx) = F_hlld(7)

[Copilot]

Pull Request Overview

This PR refactors the HLLD Riemann solver implementation in the m_riemann_solvers module to make the code more concise and maintainable by replacing many element‐wise assignments with vectorized (array slice) assignments.

• Replaced individual assignments of magnetic fields and state vectors with array constructors and slicing.
• Consolidated duplicate flux computations into vectorized updates.
• Updated indexing for momentum and magnetic fluxes using concise slice notation.

Comments suppressed due to low confidence (1)

src/simulation/m_riemann_solvers.fpp:3217

• Verify that the vectorized indexing for the momentum flux assignment correctly corresponds to the intended physical indices and ordering.

                            flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1), contxe + dir_idx(2), contxe + dir_idx(3)]) = F_hlld([2, 3, 4])

[qodo-code-review]

PR Code Suggestions ✨

No code suggestions found for the PR.

[Malmahrouqi3]

Pushed after making sure it passed HLLD specific tests which initially failed, yet I was able to resolve the issues.

[Malmahrouqi3]

Pull Request Overview

This PR refactors the HLLD Riemann solver implementation in the m_riemann_solvers module to make the code more concise and maintainable by replacing many element‐wise assignments with vectorized (array slice) assignments.

• Replaced individual assignments of magnetic fields and state vectors with array constructors and slicing.
• Consolidated duplicate flux computations into vectorized updates.
• Updated indexing for momentum and magnetic fluxes using concise slice notation.

Comments suppressed due to low confidence (1)
src/simulation/m_riemann_solvers.fpp:3217

• Verify that the vectorized indexing for the momentum flux assignment correctly corresponds to the intended physical indices and ordering.

                            flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1), contxe + dir_idx(2), contxe + dir_idx(3)]) = F_hlld([2, 3, 4])

I will reverse the vectorization here if it turns to be error-prone as alleged.

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 25 lines in your changes missing coverage. Please review.

Project coverage is 44.15%. Comparing base (7114681) to head (93ad7f0).
Report is 2 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_riemann_solvers.fpp	0.00%	25 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #909      +/-   ##
==========================================
+ Coverage   44.03%   44.15%   +0.11%     
==========================================
  Files          68       68              
  Lines       18395    18347      -48     
  Branches     2227     2227              
==========================================
  Hits         8101     8101              
+ Misses       8991     8943      -48     
  Partials     1303     1303              

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