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 11 commits into
base: master
Choose a base branch
from
+292 −697
Open

Refactor m_riemann_solvers Module (Non-Solver Subroutines) #908

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 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
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
115 changes: 115 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 @@ -1709,4 +1710,118 @@ 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
#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

end module m_variables_conversion
Loading