Skip to content

Commit e60e9d7

Browse files
Malmahrouqi3Malmahrouqi3
committed
s_hlld_riemann_solver
1 parent 4426684 commit e60e9d7

File tree

1 file changed

+44
-99
lines changed

1 file changed

+44
-99
lines changed

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 44 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -3090,19 +3090,15 @@ contains
30903090
! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
30913091
if (mhd) then
30923092
if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
3093-
B%L(1) = Bx0
3094-
B%R(1) = Bx0
3095-
B%L(2) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg)
3096-
B%R(2) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg)
3097-
B%L(3) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1)
3098-
B%R(3) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + 1)
3093+
B%L = [Bx0, qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg), qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1)]
3094+
B%R = [Bx0, qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg), qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + 1)]
30993095
else ! 2D/3D: Bx, By, Bz as variables
3100-
B%L(1) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1) - 1)
3101-
B%R(1) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(1) - 1)
3102-
B%L(2) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1)
3103-
B%R(2) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(2) - 1)
3104-
B%L(3) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1)
3105-
B%R(3) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(3) - 1)
3096+
B%L = [qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1) - 1), &
3097+
qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1), &
3098+
qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1)]
3099+
B%R = [qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(1) - 1), &
3100+
qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(2) - 1), &
3101+
qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(3) - 1)]
31063102
end if
31073103
end if
31083104
@@ -3153,74 +3149,37 @@ contains
31533149
E_starL = ((s_L - vel%L(1))*E%L - pTot_L*vel%L(1) + p_star*s_M)/(s_L - s_M)
31543150
E_starR = ((s_R - vel%R(1))*E%R - pTot_R*vel%R(1) + p_star*s_M)/(s_R - s_M)
31553151
3156-
! (5) Compute the left/right conserved state vectors
3157-
U_L(1) = rho%L
3158-
U_L(2) = rho%L*vel%L(1)
3159-
U_L(3) = rho%L*vel%L(2)
3160-
U_L(4) = rho%L*vel%L(3)
3161-
U_L(5) = B%L(2)
3162-
U_L(6) = B%L(3)
3163-
U_L(7) = E%L
3164-
3165-
U_R(1) = rho%R
3166-
U_R(2) = rho%R*vel%R(1)
3167-
U_R(3) = rho%R*vel%R(2)
3168-
U_R(4) = rho%R*vel%R(3)
3169-
U_R(5) = B%R(2)
3170-
U_R(6) = B%R(3)
3171-
U_R(7) = E%R
3172-
3173-
! (6) Compute the left/right star state vectors
3174-
U_starL(1) = rhoL_star
3175-
U_starL(2) = rhoL_star*s_M
3176-
U_starL(3) = rhoL_star*vel%L(2)
3177-
U_starL(4) = rhoL_star*vel%L(3)
3178-
U_starL(5) = B%L(2)
3179-
U_starL(6) = B%L(3)
3180-
U_starL(7) = E_starL
3181-
3182-
U_starR(1) = rhoR_star
3183-
U_starR(2) = rhoR_star*s_M
3184-
U_starR(3) = rhoR_star*vel%R(2)
3185-
U_starR(4) = rhoR_star*vel%R(3)
3186-
U_starR(5) = B%R(2)
3187-
U_starR(6) = B%R(3)
3188-
U_starR(7) = E_starR
3189-
3190-
! (7) Compute the left/right fluxes
3191-
F_L(1) = rho%L*vel%L(1)
3192-
F_L(2) = rho%L*vel%L(1)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
3193-
F_L(3) = rho%L*vel%L(1)*vel%L(2) - B%L(1)*B%L(2)
3194-
F_L(4) = rho%L*vel%L(1)*vel%L(3) - B%L(1)*B%L(3)
3195-
F_L(5) = vel%L(1)*B%L(2) - vel%L(2)*B%L(1)
3196-
F_L(6) = vel%L(1)*B%L(3) - vel%L(3)*B%L(1)
3197-
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))
3152+
! (5) Compute left/right state vectors and fluxes
3153+
U_L = [rho%L, rho%L*vel%L(1:3), B%L(2:3), E%L]
3154+
U_starL = [rhoL_star, rhoL_star*s_M, rhoL_star*vel%L(2:3), B%L(2:3), E_starL]
3155+
U_R = [rho%R, rho%R*vel%R(1:3), B%R(2:3), E%R]
3156+
U_starR = [rhoR_star, rhoR_star*s_M, rhoR_star*vel%R(2:3), B%R(2:3), E_starR]
31983157
3199-
F_R(1) = rho%R*vel%R(1)
3200-
F_R(2) = rho%R*vel%R(1)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
3201-
F_R(3) = rho%R*vel%R(1)*vel%R(2) - B%R(1)*B%R(2)
3202-
F_R(4) = rho%R*vel%R(1)*vel%R(3) - B%R(1)*B%R(3)
3203-
F_R(5) = vel%R(1)*B%R(2) - vel%R(2)*B%R(1)
3204-
F_R(6) = vel%R(1)*B%R(3) - vel%R(3)*B%R(1)
3158+
F_L(1) = U_L(2)
3159+
F_L(2) = U_L(2)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
3160+
F_L(3:4) = U_L(2)*vel%L(2:3) - B%L(1)*B%L(2:3)
3161+
F_L(5:6) = vel%L(1)*B%L(2:3) - vel%L(2:3)*B%L(1)
3162+
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))
3163+
3164+
F_R(1) = U_R(2)
3165+
F_R(2) = U_R(2)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
3166+
F_R(3:4) = U_R(2)*vel%R(2:3) - B%R(1)*B%R(2:3)
3167+
F_R(5:6) = vel%R(1)*B%R(2:3) - vel%R(2:3)*B%R(1)
32053168
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))
3206-
3207-
! (8) Compute the left/right star fluxes (note array operations)
3208-
F_starL = F_L + s_L*(U_starL - U_L)
3209-
F_starR = F_R + s_R*(U_starR - U_R)
3210-
3211-
! (9) Compute the rotational (Alfvén) speeds
3169+
! Compute the star flux using HLL relation
3170+
F_starL = F_L + s_M*(U_starL - U_L)
3171+
F_starR = F_R + s_M*(U_starR - U_R)
3172+
! Compute the rotational (Alfvén) speeds
32123173
s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star)
32133174
s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star)
3175+
! Compute the double–star states [Miyoshi Eqns. (59)-(62)]
3176+
sqrt_rhoL_star = sqrt(rhoL_star); sqrt_rhoR_star = sqrt(rhoR_star)
3177+
vL_star = vel%L(2); wL_star = vel%L(3)
3178+
vR_star = vel%R(2); wR_star = vel%R(3)
32143179
3215-
! (10) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
3216-
sqrt_rhoL_star = sqrt(rhoL_star)
3217-
sqrt_rhoR_star = sqrt(rhoR_star)
3180+
! (6) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
32183181
denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
32193182
sign_Bx = sign(1._wp, B%L(1))
3220-
vL_star = vel%L(2)
3221-
wL_star = vel%L(3)
3222-
vR_star = vel%R(2)
3223-
wR_star = vel%R(3)
32243183
v_double = (sqrt_rhoL_star*vL_star + sqrt_rhoR_star*vR_star + (B%R(2) - B%L(2))*sign_Bx)/denom_ds
32253184
w_double = (sqrt_rhoL_star*wL_star + sqrt_rhoR_star*wR_star + (B%R(3) - B%L(3))*sign_Bx)/denom_ds
32263185
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
@@ -3230,23 +3189,14 @@ contains
32303189
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
32313190
E_double = 0.5_wp*(E_doubleL + E_doubleR)
32323191
3233-
U_doubleL(1) = rhoL_star
3234-
U_doubleL(2) = rhoL_star*s_M
3235-
U_doubleL(3) = rhoL_star*v_double
3236-
U_doubleL(4) = rhoL_star*w_double
3237-
U_doubleL(5) = By_double
3238-
U_doubleL(6) = Bz_double
3239-
U_doubleL(7) = E_double
3240-
3241-
U_doubleR(1) = rhoR_star
3242-
U_doubleR(2) = rhoR_star*s_M
3243-
U_doubleR(3) = rhoR_star*v_double
3244-
U_doubleR(4) = rhoR_star*w_double
3245-
U_doubleR(5) = By_double
3246-
U_doubleR(6) = Bz_double
3247-
U_doubleR(7) = E_double
3248-
3249-
! (11) Choose HLLD flux based on wave-speed regions
3192+
U_doubleL = [rhoL_star, rhoL_star*s_M, rhoL_star*v_double, rhoL_star*w_double, By_double, Bz_double, E_double]
3193+
U_doubleR = [rhoR_star, rhoR_star*s_M, rhoR_star*w_double, rhoR_star*w_double, By_double, Bz_double, E_double]
3194+
3195+
! (7) Compute the rotational (Alfvén) speeds
3196+
s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star)
3197+
s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star)
3198+
3199+
! (8) Choose HLLD flux based on wave-speed regions
32503200
if (0.0_wp <= s_L) then
32513201
F_hlld = F_L
32523202
else if (0.0_wp <= s_starL) then
@@ -3261,20 +3211,16 @@ contains
32613211
F_hlld = F_R
32623212
end if
32633213
3264-
! (12) Reorder and write temporary variables to the flux array
3214+
! (9) Reorder and write temporary variables to the flux array
32653215
! Mass
32663216
flux_rs${XYZ}$_vf(j, k, l, 1) = F_hlld(1) ! TODO multi-component
32673217
! Momentum
3268-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1)) = F_hlld(2)
3269-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(2)) = F_hlld(3)
3270-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(3)) = F_hlld(4)
3218+
flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1), contxe + dir_idx(2), contxe + dir_idx(3)]) = F_hlld([2, 3, 4])
32713219
! Magnetic field
32723220
if (n == 0) then
3273-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg) = F_hlld(5)
3274-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1) = F_hlld(6)
3221+
flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg, B_idx%beg + 1]) = F_hlld([5, 6])
32753222
else
3276-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1) = F_hlld(5)
3277-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1) = F_hlld(6)
3223+
flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg + dir_idx(2) - 1, B_idx%beg + dir_idx(3) - 1]) = F_hlld([5, 6])
32783224
end if
32793225
! Energy
32803226
flux_rs${XYZ}$_vf(j, k, l, E_idx) = F_hlld(7)
@@ -3283,7 +3229,6 @@ contains
32833229
do i = advxb, advxe
32843230
flux_rs${XYZ}$_vf(j, k, l, i) = 0._wp ! TODO multi-component (zero for now)
32853231
end do
3286-
32873232
flux_src_rs${XYZ}$_vf(j, k, l, advxb) = 0._wp
32883233
end do
32893234
end do

0 commit comments

Comments
 (0)