Skip to content

Commit 967c128

Browse files
Malmahrouqi3Malmahrouqi3
committed
s_compute_wave_speed
1 parent 4864d36 commit 967c128

File tree

2 files changed

+135
-198
lines changed

2 files changed

+135
-198
lines changed

‎src/common/m_variables_conversion.fpp‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@ module m_variables_conversion
4444
#ifndef MFC_PRE_PROCESS
4545
s_compute_speed_of_sound, &
4646
s_compute_fast_magnetosonic_speed, &
47+
s_compute_wave_speed, &
4748
#endif
4849
s_finalize_variables_conversion_module
4950

@@ -1709,4 +1710,118 @@ contains
17091710
end subroutine s_compute_fast_magnetosonic_speed
17101711
#endif
17111712

1713+
#ifndef MFC_PRE_PROCESS
1714+
subroutine s_compute_wave_speed(wave_speeds, vel_L, vel_R, pres_L, pres_R, rho_L, rho_R, rho_avg, &
1715+
c_L, c_R, c_avg, c_fast_L, c_fast_R, G_L, G_R, &
1716+
tau_e_L, tau_e_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, &
1717+
s_L, s_R, s_S, s_M, s_P, idx, idx_tau)
1718+
1719+
! Computes the wave speeds for the Riemann solver
1720+
#ifdef _CRAYFTN
1721+
!DIR$ INLINEALWAYS s_compute_wave_speed
1722+
#else
1723+
!$acc routine seq
1724+
#endif
1725+
1726+
! Input parameters
1727+
integer, intent(in) :: wave_speeds
1728+
integer, intent(in) :: idx, idx_tau
1729+
real(wp), intent(in) :: rho_L, rho_R
1730+
real(wp), dimension(:), intent(in) :: vel_L, vel_R, tau_e_L, tau_e_R
1731+
real(wp), intent(in) :: pres_L, pres_R, c_L, c_R
1732+
real(wp), intent(in) :: gamma_L, gamma_R, pi_inf_L, pi_inf_R
1733+
real(wp), intent(in) :: rho_avg, c_avg
1734+
real(wp), intent(in) :: c_fast_L, c_fast_R
1735+
real(wp), intent(in) :: G_L, G_R
1736+
1737+
! Local variables
1738+
real(wp) :: pres_SL, pres_SR, Ms_L, Ms_R
1739+
1740+
! Output parameters
1741+
real(wp), intent(out) :: s_L, s_R, s_S, s_M, s_P
1742+
1743+
if (wave_speeds == 1) then
1744+
if (elasticity) then
1745+
s_L = min(vel_L(idx) - sqrt(c_L*c_L + &
1746+
(((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L), vel_R(idx) - sqrt(c_R*c_R + &
1747+
(((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R))
1748+
s_R = max(vel_R(idx) + sqrt(c_R*c_R + &
1749+
(((4_wp*G_R)/3_wp) + tau_e_R(idx_tau))/rho_R), vel_L(idx) + sqrt(c_L*c_L + &
1750+
(((4_wp*G_L)/3_wp) + tau_e_L(idx_tau))/rho_L))
1751+
s_S = (pres_R - tau_e_R(idx_tau) - pres_L + &
1752+
tau_e_L(idx_tau) + rho_L*vel_L(idx)*(s_L - vel_L(idx)) - &
1753+
rho_R*vel_R(idx)*(s_R - vel_R(idx)))/(rho_L*(s_L - vel_L(idx)) - &
1754+
rho_R*(s_R - vel_R(idx)))
1755+
else if (mhd) then
1756+
s_L = min(vel_L(idx) - c_fast_L, vel_R(idx) - c_fast_R)
1757+
s_R = max(vel_R(idx) + c_fast_R, vel_L(idx) + c_fast_L)
1758+
s_S = (pres_R - pres_L + rho_L*vel_L(idx)* &
1759+
(s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) &
1760+
/(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx)))
1761+
else if (hypoelasticity) then
1762+
s_L = min(vel_L(idx) - sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + &
1763+
tau_e_L(idx_tau))/rho_L) &
1764+
, vel_R(idx) - sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + &
1765+
tau_e_R(idx_tau))/rho_R))
1766+
s_R = max(vel_R(idx) + sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + &
1767+
tau_e_R(idx_tau))/rho_R) &
1768+
, vel_L(idx) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + &
1769+
tau_e_L(idx_tau))/rho_L))
1770+
s_S = (pres_R - pres_L + rho_L*vel_L(idx)* &
1771+
(s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) &
1772+
/(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx)))
1773+
else if (hyperelasticity) then
1774+
s_L = min(vel_L(idx) - sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L) &
1775+
, vel_R(idx) - sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R))
1776+
s_R = max(vel_R(idx) + sqrt(c_R*c_R + (4._wp*G_R/3._wp)/rho_R) &
1777+
, vel_L(idx) + sqrt(c_L*c_L + (4._wp*G_L/3._wp)/rho_L))
1778+
s_S = (pres_R - pres_L + rho_L*vel_L(idx)* &
1779+
(s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) &
1780+
/(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx)))
1781+
else
1782+
s_L = min(vel_L(idx) - c_L, vel_R(idx) - c_R)
1783+
s_R = max(vel_R(idx) + c_R, vel_L(idx) + c_L)
1784+
s_S = (pres_R - pres_L + rho_L*vel_L(idx)* &
1785+
(s_L - vel_L(idx)) - rho_R*vel_R(idx)*(s_R - vel_R(idx))) &
1786+
/(rho_L*(s_L - vel_L(idx)) - rho_R*(s_R - vel_R(idx)))
1787+
end if
1788+
else if (wave_speeds == 2) then
1789+
pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg*(vel_L(idx) - vel_R(idx)))
1790+
pres_SR = pres_SL
1791+
Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* &
1792+
(pres_SL/pres_L - 1._wp)*pres_L/ &
1793+
((pres_L + pi_inf_L/(1._wp + gamma_L)))))
1794+
Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* &
1795+
(pres_SR/pres_R - 1._wp)*pres_R/ &
1796+
((pres_R + pi_inf_R/(1._wp + gamma_R)))))
1797+
s_L = vel_L(idx) - c_L*Ms_L
1798+
s_R = vel_R(idx) + c_R*Ms_R
1799+
s_S = 5e-1_wp*((vel_L(idx) + vel_R(idx)) + (pres_L - pres_R)/(rho_avg*c_avg))
1800+
end if
1801+
1802+
! ! follows Einfeldt et al.
1803+
! s_M/P = min/max(0.,s_L/R)
1804+
s_M = min(0._wp, s_L)
1805+
s_P = max(0._wp, s_R)
1806+
1807+
#ifdef DEBUG
1808+
! Check for potential issues in wave speed calculation
1809+
if (s_R <= s_L) then
1810+
print *, 'WARNING: Wave speed issue detected in s_compute_wave_speed'
1811+
print *, 'Left wave speed >= Right wave speed:', s_L, s_R
1812+
print *, 'Input velocities :', vel_L(idx), vel_R(idx)
1813+
print *, 'Sound speeds:', c_L, c_R
1814+
print *, 'Densities:', rho_L, rho_R
1815+
print *, 'Pressures:', pres_L, pres_R
1816+
print *, 'Wave speeds method:', wave_speeds
1817+
if (elasticity .or. hypoelasticity .or. hyperelasticity) then
1818+
print *, 'Shear moduli:', G_L, G_R
1819+
end if
1820+
call s_mpi_abort('Error: Invalid wave speeds in s_compute_wave_speed')
1821+
end if
1822+
#endif
1823+
1824+
end subroutine s_compute_wave_speed
1825+
#endif
1826+
17121827
end module m_variables_conversion

0 commit comments

Comments
 (0)