Skip to content

Refactor m_riemann_solvers Module (Non-Solver Subroutines) #908

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 25 commits into
base: master
Choose a base branch
from
Open

Refactor m_riemann_solvers Module (Non-Solver Subroutines) #908

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
967c128
s_compute_wave_speed
Malmahrouqi3 Jun 29, 2025
89280ad
low-stakes refactoring & lint/format
Malmahrouqi3 Jun 29, 2025
3bfa9f4
fixed few errors
Malmahrouqi3 Jun 29, 2025
ddaa1b6
quick fix
Malmahrouqi3 Jun 29, 2025
12eabf4
simple fix
Malmahrouqi3 Jun 29, 2025
4e2a439
removed &
Malmahrouqi3 Jun 29, 2025
b684a4e
remove allocatable from pointers
Malmahrouqi3 Jun 29, 2025
6fb8c65
m_riemann_solvers cleanup & format/lint
Malmahrouqi3 Jun 29, 2025
dd158d0
Merge branch 'master' into refactor
sbryngelson Jul 1, 2025
46de417
Merge branch 'master' into refactor
sbryngelson Jul 3, 2025
dcfc3e0
lint-fix
Malmahrouqi3 Jul 11, 2025
7f0c19d
Merge remote-tracking branch 'upstream/master' into refactor
Malmahrouqi3 Aug 14, 2025
6e7281b
format
Aug 14, 2025
3ad2ddb
lint/format + few fixes
Aug 18, 2025
dc8019f
removed raw directive in s_compute_wave_speed
Aug 18, 2025
f3cde3e
corrected something
Aug 18, 2025
156ec55
Merge branch 'MFlowCode:master' into refactor
Malmahrouqi3 Aug 18, 2025
b6a352e
Merge branch 'master' into refactor
Malmahrouqi3 Sep 1, 2025
3c4b4e1
Merge branch 'master' into refactor
Malmahrouqi3 Sep 6, 2025
797c40f
Merge branch 'master' into refactor
Malmahrouqi3 Sep 11, 2025
676082e
undid all changes, retained s_compute_wave_speed
Sep 15, 2025
f5e8c52
minor changes
Sep 15, 2025
e9bc89c
Merge branch 'master' into refactor
Malmahrouqi3 Sep 15, 2025
89064af
Merge branch 'master' into refactor
Malmahrouqi3 Sep 17, 2025
5ba2e99
Merge branch 'master' into refactor
sbryngelson Oct 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
112 changes: 112 additions & 0 deletions src/common/m_variables_conversion.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ module m_variables_conversion
#ifndef MFC_PRE_PROCESS
s_compute_speed_of_sound, &
s_compute_fast_magnetosonic_speed, &
s_compute_wave_speed, &
#endif
s_finalize_variables_conversion_module

Expand Down Expand Up @@ -1726,4 +1727,115 @@ contains
end subroutine s_compute_fast_magnetosonic_speed
#endif

#ifndef MFC_PRE_PROCESS
subroutine s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, &
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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))

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
$:GPU_ROUTINE(function_name='s_compute_wave_speed', &
& parallelism='[seq]', cray_inline=True)

! 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))
Comment on lines +1803 to +1813
Copy link

Copilot AI Sep 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Inconsistent use of floating point literal notation. Line 1804, 1806, 1809, and 1814 use 5e-1_wp while other parts of the codebase use 5.e-1_wp (with decimal point). Should use 5.e-1_wp for consistency with established patterns in the codebase.

Suggested change
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))
pres_SL = 5.e-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 + ((5.e-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 + ((5.e-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 = 5.e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg))

Copilot uses AI. Check for mistakes.
Comment on lines +1803 to +1813
Copy link

Copilot AI Sep 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Inconsistent use of floating point literal notation. Line 1804, 1806, 1809, and 1814 use 5e-1_wp while other parts of the codebase use 5.e-1_wp (with decimal point). Should use 5.e-1_wp for consistency with established patterns in the codebase.

Suggested change
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))
pres_SL = 5.e-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 + ((5.e-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 + ((5.e-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 = 5.e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg))

Copilot uses AI. Check for mistakes.
Comment on lines +1803 to +1813
Copy link

Copilot AI Sep 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Inconsistent use of floating point literal notation. Line 1804, 1806, 1809, and 1814 use 5e-1_wp while other parts of the codebase use 5.e-1_wp (with decimal point). Should use 5.e-1_wp for consistency with established patterns in the codebase.

Suggested change
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))
pres_SL = 5.e-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 + ((5.e-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 + ((5.e-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 = 5.e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg))

Copilot uses AI. Check for mistakes.
Comment on lines +1803 to +1813
Copy link

Copilot AI Sep 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Inconsistent use of floating point literal notation. Line 1804, 1806, 1809, and 1814 use 5e-1_wp while other parts of the codebase use 5.e-1_wp (with decimal point). Should use 5.e-1_wp for consistency with established patterns in the codebase.

Suggested change
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))
pres_SL = 5.e-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 + ((5.e-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 + ((5.e-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 = 5.e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg))

Copilot uses AI. Check for mistakes.
end if

! ! follows Einfeldt et al.
Copy link

Copilot AI Sep 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Double comment marker ! ! is inconsistent with standard single comment style used elsewhere in the codebase. Should use single ! for consistency.

Suggested change
! ! follows Einfeldt et al.
! follows Einfeldt et al.

Copilot uses AI. Check for mistakes.
! 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

end module m_variables_conversion
Loading
Loading