Skip to content

Commit b74ce95

Browse files
mohdsaid497566Malmahrouqi3
mohdsaid497566
authored and
Malmahrouqi3
committed
cleaned up s_hlld_riemann_solver redundancy
1 parent cce093a commit b74ce95

File tree

1 file changed

+60
-84
lines changed

1 file changed

+60
-84
lines changed

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 60 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -2849,7 +2849,7 @@ contains
28492849
28502850
real(wp), dimension(7) :: U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR
28512851
real(wp), dimension(7) :: F_L, F_R, F_starL, F_starR, F_hlld
2852-
2852+
real
28532853
! Indices for U and F: (rho, rho*vel(1), rho*vel(2), rho*vel(3), By, Bz, E)
28542854
! Note: vel and B are permutated, so vel(1) is the normal velocity, and x is the normal direction
28552855
! Note: Bx is omitted as the magnetic flux is always zero in the normal direction
@@ -2859,6 +2859,50 @@ contains
28592859
real(wp) :: v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double
28602860
28612861
integer :: i, j, k, l
2862+
contains
2863+
function s_compute_U_double(rho_star, s_M, v_double, w_double, By_double, Bz_double, E_double) result(U_double)
2864+
implicit none
2865+
real(wp), intent(in) :: rho_star, s_M, v_double, w_double, By_double, Bz_double, E_double
2866+
real(wp) :: U_double(7)
2867+
2868+
U_double(1) = rho_star
2869+
U_double(2) = rho_star*s_M
2870+
U_double(3) = rho_star*v_double
2871+
U_double(4) = rho_star*w_double
2872+
U_double(5) = By_double
2873+
U_double(6) = Bz_double
2874+
U_double(7) = E_double
2875+
end function s_compute_U_double
2876+
2877+
subroutine s_compute_hlld_state_variables (side, rho, vel, B, E, pTot, rho_star, s_M, E_star, s_wave, &
2878+
U, F, U_star, F_star, sqrt_rho_star, v_star, w_star)
2879+
implicit none
2880+
! Input parameters
2881+
character(len=1), intent(in) :: side ! dummy 'L' for left or 'R' for right
2882+
real(wp), intent(in) :: rho, pTot, rho_star, s_M, E_star, s_wave, E
2883+
real(wp), dimension(:), intent(in) :: vel, B
2884+
! Output parameters
2885+
real(wp), dimension(7), intent(out) :: U, F, U_star
2886+
real(wp), intent(out) :: sqrt_rho_star, v_star, w_star
2887+
real(wp), dimension(7), intent(out) :: F_star
2888+
! Compute the base state vector
2889+
U(1) = rho, U(2) = rho*vel(1), U(3) = rho*vel(2), U(4) = rho*vel(3)
2890+
U(5) = B(2), U(6) = B(3), U(7) = E
2891+
! Compute the flux vector
2892+
F(1) = U(2), F(2) = U(2)*vel(1) - B(1)*B(1) + pTot, F(3) = U(2)*vel(2) - B(1)*B(2)
2893+
F(4) = U(2)*vel(3) - B(1)*B(3), F(5) = vel(1)*B(2) - vel(2)*B(1)
2894+
F(6) = vel(1)*B(3) - vel(3)*B(1), F(7) = (E + pTot)*vel(1) - B(1)*(vel(1)*B(1) + vel(2)*B(2) + vel(3)*B(3))
2895+
! Compute the star state
2896+
U_star(1) = rho_star, U_star(2) = rho_star*s_M, U_star(3) = rho_star*vel(2)
2897+
U_star(4) = rho_star*vel(3), U_star(5) = B(2), U_star(6) = B(3)
2898+
U_star(7) = E_star
2899+
! Compute the star flux using HLL relation
2900+
F_star = F + s_wave*(U_star - U)
2901+
! Compute additional parameters needed for double-star states
2902+
sqrt_rho_star = sqrt(rho_star)
2903+
v_star = vel(2)
2904+
w_star = vel(3)
2905+
end subroutine s_compute_hlld_state_variables
28622906
28632907
call s_populate_riemann_states_variables_buffers( &
28642908
qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
@@ -2969,74 +3013,15 @@ contains
29693013
E_starL = ((s_L - vel%L(1))*E%L - pTot_L*vel%L(1) + p_star*s_M)/(s_L - s_M)
29703014
E_starR = ((s_R - vel%R(1))*E%R - pTot_R*vel%R(1) + p_star*s_M)/(s_R - s_M)
29713015
2972-
! (5) Compute the left/right conserved state vectors
2973-
U_L(1) = rho%L
2974-
U_L(2) = rho%L*vel%L(1)
2975-
U_L(3) = rho%L*vel%L(2)
2976-
U_L(4) = rho%L*vel%L(3)
2977-
U_L(5) = B%L(2)
2978-
U_L(6) = B%L(3)
2979-
U_L(7) = E%L
2980-
2981-
U_R(1) = rho%R
2982-
U_R(2) = rho%R*vel%R(1)
2983-
U_R(3) = rho%R*vel%R(2)
2984-
U_R(4) = rho%R*vel%R(3)
2985-
U_R(5) = B%R(2)
2986-
U_R(6) = B%R(3)
2987-
U_R(7) = E%R
2988-
2989-
! (6) Compute the left/right star state vectors
2990-
U_starL(1) = rhoL_star
2991-
U_starL(2) = rhoL_star*s_M
2992-
U_starL(3) = rhoL_star*vel%L(2)
2993-
U_starL(4) = rhoL_star*vel%L(3)
2994-
U_starL(5) = B%L(2)
2995-
U_starL(6) = B%L(3)
2996-
U_starL(7) = E_starL
2997-
2998-
U_starR(1) = rhoR_star
2999-
U_starR(2) = rhoR_star*s_M
3000-
U_starR(3) = rhoR_star*vel%R(2)
3001-
U_starR(4) = rhoR_star*vel%R(3)
3002-
U_starR(5) = B%R(2)
3003-
U_starR(6) = B%R(3)
3004-
U_starR(7) = E_starR
3005-
3006-
! (7) Compute the left/right fluxes
3007-
F_L(1) = rho%L*vel%L(1)
3008-
F_L(2) = rho%L*vel%L(1)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
3009-
F_L(3) = rho%L*vel%L(1)*vel%L(2) - B%L(1)*B%L(2)
3010-
F_L(4) = rho%L*vel%L(1)*vel%L(3) - B%L(1)*B%L(3)
3011-
F_L(5) = vel%L(1)*B%L(2) - vel%L(2)*B%L(1)
3012-
F_L(6) = vel%L(1)*B%L(3) - vel%L(3)*B%L(1)
3013-
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))
3014-
3015-
F_R(1) = rho%R*vel%R(1)
3016-
F_R(2) = rho%R*vel%R(1)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
3017-
F_R(3) = rho%R*vel%R(1)*vel%R(2) - B%R(1)*B%R(2)
3018-
F_R(4) = rho%R*vel%R(1)*vel%R(3) - B%R(1)*B%R(3)
3019-
F_R(5) = vel%R(1)*B%R(2) - vel%R(2)*B%R(1)
3020-
F_R(6) = vel%R(1)*B%R(3) - vel%R(3)*B%R(1)
3021-
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))
3022-
3023-
! (8) Compute the left/right star fluxes (note array operations)
3024-
F_starL = F_L + s_L*(U_starL - U_L)
3025-
F_starR = F_R + s_R*(U_starR - U_R)
3026-
3027-
! (9) Compute the rotational (Alfvén) speeds
3028-
s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star)
3029-
s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star)
3016+
! (5) Compute left/right state vectors and fluxes
3017+
call s_compute_hlld_state_variables('L', rho%L, vel%L, B%L, E%L, pTot_L, rhoL_star, s_M, E_starL, s_L, &
3018+
U_L, F_L, U_starL, F_starL, sqrt_rhoL_star, vL_star, wL_star)
3019+
call s_compute_hlld_state_variables('R', rho%R, vel%R, B%R, E%R, pTot_R, rhoR_star, s_M, E_starR, s_R, &
3020+
U_R, F_R, U_starR, F_starR, sqrt_rhoR_star, vR_star, wR_star)
30303021
3031-
! (10) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
3032-
sqrt_rhoL_star = sqrt(rhoL_star)
3033-
sqrt_rhoR_star = sqrt(rhoR_star)
3022+
! (6) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
30343023
denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
30353024
sign_Bx = sign(1._wp, B%L(1))
3036-
vL_star = vel%L(2)
3037-
wL_star = vel%L(3)
3038-
vR_star = vel%R(2)
3039-
wR_star = vel%R(3)
30403025
v_double = (sqrt_rhoL_star*vL_star + sqrt_rhoR_star*vR_star + (B%R(2) - B%L(2))*sign_Bx)/denom_ds
30413026
w_double = (sqrt_rhoL_star*wL_star + sqrt_rhoR_star*wR_star + (B%R(3) - B%L(3))*sign_Bx)/denom_ds
30423027
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
@@ -3046,23 +3031,14 @@ contains
30463031
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
30473032
E_double = 0.5_wp*(E_doubleL + E_doubleR)
30483033
3049-
U_doubleL(1) = rhoL_star
3050-
U_doubleL(2) = rhoL_star*s_M
3051-
U_doubleL(3) = rhoL_star*v_double
3052-
U_doubleL(4) = rhoL_star*w_double
3053-
U_doubleL(5) = By_double
3054-
U_doubleL(6) = Bz_double
3055-
U_doubleL(7) = E_double
3056-
3057-
U_doubleR(1) = rhoR_star
3058-
U_doubleR(2) = rhoR_star*s_M
3059-
U_doubleR(3) = rhoR_star*v_double
3060-
U_doubleR(4) = rhoR_star*w_double
3061-
U_doubleR(5) = By_double
3062-
U_doubleR(6) = Bz_double
3063-
U_doubleR(7) = E_double
3064-
3065-
! (11) Choose HLLD flux based on wave-speed regions
3034+
U_doubleL = s_compute_U_double(rhoL_star, s_M, v_double, w_double, By_double, Bz_double, E_double)
3035+
U_doubleR = s_compute_U_double(rhoR_star, s_M, v_double, w_double, By_double, Bz_double, E_double)
3036+
3037+
! (7) Compute the rotational (Alfvén) speeds
3038+
s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star)
3039+
s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star)
3040+
3041+
! (8) Choose HLLD flux based on wave-speed regions
30663042
if (0.0_wp <= s_L) then
30673043
F_hlld = F_L
30683044
else if (0.0_wp <= s_starL) then
@@ -3077,7 +3053,7 @@ contains
30773053
F_hlld = F_R
30783054
end if
30793055
3080-
! (12) Reorder and write temporary variables to the flux array
3056+
! (9) Reorder and write temporary variables to the flux array
30813057
! Mass
30823058
flux_rs${XYZ}$_vf(j, k, l, 1) = F_hlld(1) ! TODO multi-component
30833059
! Momentum

0 commit comments

Comments
 (0)