@@ -2995,19 +2995,15 @@ contains
29952995 ! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
29962996 if (mhd) then
29972997 if (n == 0 ) then ! 1D : constant Bx; By, Bz as variables; only in x so not permutated
2998- B%L(1 ) = Bx0
2999- B%R(1 ) = Bx0
3000- B%L(2 ) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg)
3001- B%R(2 ) = qR_prim_rs${XYZ}$_vf(j + 1 , k, l, B_idx%beg)
3002- B%L(3 ) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1 )
3003- B%R(3 ) = qR_prim_rs${XYZ}$_vf(j + 1 , k, l, B_idx%beg + 1 )
2998+ 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 )]
2999+ 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 )]
30043000 else ! 2D / 3D : Bx, By, Bz as variables
3005- B%L( 1 ) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1 ) - 1 )
3006- B%R( 1 ) = qR_prim_rs ${XYZ}$_vf(j + 1 , k, l, B_idx%beg + dir_idx(1 ) - 1 )
3007- B%L( 2 ) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2 ) - 1 )
3008- B%R( 2 ) = qR_prim_rs${XYZ}$_vf(j + 1 , k, l, B_idx%beg + dir_idx(2 ) - 1 )
3009- B%L( 3 ) = qL_prim_rs ${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3 ) - 1 )
3010- B%R( 3 ) = qR_prim_rs${XYZ}$_vf(j + 1 , k, l, B_idx%beg + dir_idx(3 ) - 1 )
3001+ B%L = [ qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1 ) - 1 ), &
3002+ qL_prim_rs ${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2 ) - 1 ), &
3003+ qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3 ) - 1 )]
3004+ B%R = [ qR_prim_rs${XYZ}$_vf(j + 1 , k, l, B_idx%beg + dir_idx(1 ) - 1 ), &
3005+ qR_prim_rs ${XYZ}$_vf(j + 1 , k, l, B_idx%beg + dir_idx(2 ) - 1 ), &
3006+ qR_prim_rs${XYZ}$_vf(j + 1 , k, l, B_idx%beg + dir_idx(3 ) - 1 )]
30113007 end if
30123008 end if
30133009
@@ -3058,74 +3054,38 @@ contains
30583054 E_starL = ((s_L - vel%L(1 ))* E%L - pTot_L* vel%L(1 ) + p_star* s_M)/ (s_L - s_M)
30593055 E_starR = ((s_R - vel%R(1 ))* E%R - pTot_R* vel%R(1 ) + p_star* s_M)/ (s_R - s_M)
30603056
3061- ! (5 ) Compute the left/ right conserved state vectors
3062- U_L(1 ) = rho%L
3063- U_L(2 ) = rho%L* vel%L(1 )
3064- U_L(3 ) = rho%L* vel%L(2 )
3065- U_L(4 ) = rho%L* vel%L(3 )
3066- U_L(5 ) = B%L(2 )
3067- U_L(6 ) = B%L(3 )
3068- U_L(7 ) = E%L
3069-
3070- U_R(1 ) = rho%R
3071- U_R(2 ) = rho%R* vel%R(1 )
3072- U_R(3 ) = rho%R* vel%R(2 )
3073- U_R(4 ) = rho%R* vel%R(3 )
3074- U_R(5 ) = B%R(2 )
3075- U_R(6 ) = B%R(3 )
3076- U_R(7 ) = E%R
3077-
3078- ! (6 ) Compute the left/ right star state vectors
3079- U_starL(1 ) = rhoL_star
3080- U_starL(2 ) = rhoL_star* s_M
3081- U_starL(3 ) = rhoL_star* vel%L(2 )
3082- U_starL(4 ) = rhoL_star* vel%L(3 )
3083- U_starL(5 ) = B%L(2 )
3084- U_starL(6 ) = B%L(3 )
3085- U_starL(7 ) = E_starL
3086-
3087- U_starR(1 ) = rhoR_star
3088- U_starR(2 ) = rhoR_star* s_M
3089- U_starR(3 ) = rhoR_star* vel%R(2 )
3090- U_starR(4 ) = rhoR_star* vel%R(3 )
3091- U_starR(5 ) = B%R(2 )
3092- U_starR(6 ) = B%R(3 )
3093- U_starR(7 ) = E_starR
3094-
3095- ! (7 ) Compute the left/ right fluxes
3096- F_L(1 ) = rho%L* vel%L(1 )
3097- F_L(2 ) = rho%L* vel%L(1 )* vel%L(1 ) - B%L(1 )* B%L(1 ) + pTot_L
3098- F_L(3 ) = rho%L* vel%L(1 )* vel%L(2 ) - B%L(1 )* B%L(2 )
3099- F_L(4 ) = rho%L* vel%L(1 )* vel%L(3 ) - B%L(1 )* B%L(3 )
3100- F_L(5 ) = vel%L(1 )* B%L(2 ) - vel%L(2 )* B%L(1 )
3101- F_L(6 ) = vel%L(1 )* B%L(3 ) - vel%L(3 )* B%L(1 )
3057+ ! (5 ) Compute left/ right state vectors and fluxes
3058+ U_L = [rho%L, rho%L* vel%L(1 :3 ), B%L(2 :3 ), E%L]
3059+ U_starL = [rhoL_star, rhoL_star* s_M, rhoL_star* vel%L(2 :3 ), B%L(2 :3 ), E_starL]
3060+ U_R = [rho%R, rho%R* vel%R(1 :3 ), B%R(2 :3 ), E%R]
3061+ U_starR = [rhoR_star, rhoR_star* s_M, rhoR_star* vel%R(2 :3 ), B%R(2 :3 ), E_starR]
3062+
3063+ ! Compute the left/ right fluxes
3064+ F_L(1 ) = U_L(2 )
3065+ F_L(2 ) = U_L(2 )* vel%L(1 ) - B%L(1 )* B%L(1 ) + pTot_L
3066+ F_L(3 :4 ) = U_L(2 )* vel%L(2 :3 ) - B%L(1 )* B%L(2 :3 )
3067+ F_L(5 :6 ) = vel%L(1 )* B%L(2 :3 ) - vel%L(2 :3 )* B%L(1 )
31023068 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 ))
31033069
3104- F_R(1 ) = rho%R* vel%R(1 )
3105- F_R(2 ) = rho%R* vel%R(1 )* vel%R(1 ) - B%R(1 )* B%R(1 ) + pTot_R
3106- F_R(3 ) = rho%R* vel%R(1 )* vel%R(2 ) - B%R(1 )* B%R(2 )
3107- F_R(4 ) = rho%R* vel%R(1 )* vel%R(3 ) - B%R(1 )* B%R(3 )
3108- F_R(5 ) = vel%R(1 )* B%R(2 ) - vel%R(2 )* B%R(1 )
3109- F_R(6 ) = vel%R(1 )* B%R(3 ) - vel%R(3 )* B%R(1 )
3070+ F_R(1 ) = U_R(2 )
3071+ F_R(2 ) = U_R(2 )* vel%R(1 ) - B%R(1 )* B%R(1 ) + pTot_R
3072+ F_R(3 :4 ) = U_R(2 )* vel%R(2 :3 ) - B%R(1 )* B%R(2 :3 )
3073+ F_R(5 :6 ) = vel%R(1 )* B%R(2 :3 ) - vel%R(2 :3 )* B%R(1 )
31103074 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 ))
3111-
3112- ! (8 ) Compute the left/ right star fluxes (note array operations)
3075+ ! Compute the star flux using HLL relation
31133076 F_starL = F_L + s_L* (U_starL - U_L)
31143077 F_starR = F_R + s_R* (U_starR - U_R)
3115-
3116- ! (9 ) Compute the rotational (Alfvén) speeds
3078+ ! Compute the rotational (Alfvén) speeds
31173079 s_starL = s_M - abs (B%L(1 ))/ sqrt (rhoL_star)
31183080 s_starR = s_M + abs (B%L(1 ))/ sqrt (rhoR_star)
3081+ ! Compute the double–star states [Miyoshi Eqns. (59 )- (62 )]
3082+ sqrt_rhoL_star = sqrt (rhoL_star); sqrt_rhoR_star = sqrt (rhoR_star)
3083+ vL_star = vel%L(2 ); wL_star = vel%L(3 )
3084+ vR_star = vel%R(2 ); wR_star = vel%R(3 )
31193085
3120- ! (10 ) Compute the double–star states [Miyoshi Eqns. (59 )- (62 )]
3121- sqrt_rhoL_star = sqrt (rhoL_star)
3122- sqrt_rhoR_star = sqrt (rhoR_star)
3086+ ! (6 ) Compute the double–star states [Miyoshi Eqns. (59 )- (62 )]
31233087 denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
31243088 sign_Bx = sign (1._wp , B%L(1 ))
3125- vL_star = vel%L(2 )
3126- wL_star = vel%L(3 )
3127- vR_star = vel%R(2 )
3128- wR_star = vel%R(3 )
31293089 v_double = (sqrt_rhoL_star* vL_star + sqrt_rhoR_star* vR_star + (B%R(2 ) - B%L(2 ))* sign_Bx)/ denom_ds
31303090 w_double = (sqrt_rhoL_star* wL_star + sqrt_rhoR_star* wR_star + (B%R(3 ) - B%L(3 ))* sign_Bx)/ denom_ds
31313091 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
@@ -3135,21 +3095,8 @@ contains
31353095 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
31363096 E_double = 0.5_wp * (E_doubleL + E_doubleR)
31373097
3138- U_doubleL(1 ) = rhoL_star
3139- U_doubleL(2 ) = rhoL_star* s_M
3140- U_doubleL(3 ) = rhoL_star* v_double
3141- U_doubleL(4 ) = rhoL_star* w_double
3142- U_doubleL(5 ) = By_double
3143- U_doubleL(6 ) = Bz_double
3144- U_doubleL(7 ) = E_double
3145-
3146- U_doubleR(1 ) = rhoR_star
3147- U_doubleR(2 ) = rhoR_star* s_M
3148- U_doubleR(3 ) = rhoR_star* v_double
3149- U_doubleR(4 ) = rhoR_star* w_double
3150- U_doubleR(5 ) = By_double
3151- U_doubleR(6 ) = Bz_double
3152- U_doubleR(7 ) = E_double
3098+ U_doubleL = [rhoL_star, rhoL_star* s_M, rhoL_star* v_double, rhoL_star* w_double, By_double, Bz_double, E_double]
3099+ U_doubleR = [rhoR_star, rhoR_star* s_M, rhoR_star* v_double, rhoR_star* w_double, By_double, Bz_double, E_double]
31533100
31543101 ! (11 ) Choose HLLD flux based on wave- speed regions
31553102 if (0.0_wp <= s_L) then
@@ -3170,16 +3117,12 @@ contains
31703117 ! Mass
31713118 flux_rs${XYZ}$_vf(j, k, l, 1 ) = F_hlld(1 ) ! TODO multi- component
31723119 ! Momentum
3173- flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1 )) = F_hlld(2 )
3174- flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(2 )) = F_hlld(3 )
3175- flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(3 )) = F_hlld(4 )
3120+ flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1 ), contxe + dir_idx(2 ), contxe + dir_idx(3 )]) = F_hlld([2 , 3 , 4 ])
31763121 ! Magnetic field
31773122 if (n == 0 ) then
3178- flux_rs${XYZ}$_vf(j, k, l, B_idx%beg) = F_hlld(5 )
3179- flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1 ) = F_hlld(6 )
3123+ flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg, B_idx%beg + 1 ]) = F_hlld([5 , 6 ])
31803124 else
3181- flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2 ) - 1 ) = F_hlld(5 )
3182- flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3 ) - 1 ) = F_hlld(6 )
3125+ 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 ])
31833126 end if
31843127 ! Energy
31853128 flux_rs${XYZ}$_vf(j, k, l, E_idx) = F_hlld(7 )
0 commit comments