-
Notifications
You must be signed in to change notification settings - Fork 108
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.
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
base: master
Are you sure you want to change the base?
Changes from all commits
967c128
89280ad
3bfa9f4
ddaa1b6
12eabf4
4e2a439
b684a4e
6fb8c65
dd158d0
46de417
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -44,6 +44,7 @@ | |
#ifndef MFC_PRE_PROCESS | ||
s_compute_speed_of_sound, & | ||
s_compute_fast_magnetosonic_speed, & | ||
s_compute_wave_speed, & | ||
#endif | ||
s_finalize_variables_conversion_module | ||
|
||
|
@@ -1709,4 +1710,118 @@ | |
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, & | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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? There was a problem hiding this comment. Choose a reason for hiding this commentThe 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
|
||
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 + & | ||
sbryngelson marked this conversation as resolved.
Show resolved
Hide resolved
|
||
(((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 |
Uh oh!
There was an error while loading. Please reload this page.