@@ -2995,19 +2995,15 @@ contains
2995
2995
! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
2996
2996
if (mhd) then
2997
2997
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 )]
3004
3000
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 )]
3011
3007
end if
3012
3008
end if
3013
3009
@@ -3058,74 +3054,38 @@ contains
3058
3054
E_starL = ((s_L - vel%L(1 ))* E%L - pTot_L* vel%L(1 ) + p_star* s_M)/ (s_L - s_M)
3059
3055
E_starR = ((s_R - vel%R(1 ))* E%R - pTot_R* vel%R(1 ) + p_star* s_M)/ (s_R - s_M)
3060
3056
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 )
3102
3068
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 ))
3103
3069
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 )
3110
3074
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
3113
3076
F_starL = F_L + s_L* (U_starL - U_L)
3114
3077
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
3117
3079
s_starL = s_M - abs (B%L(1 ))/ sqrt (rhoL_star)
3118
3080
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 )
3119
3085
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 )]
3123
3087
denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
3124
3088
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 )
3129
3089
v_double = (sqrt_rhoL_star* vL_star + sqrt_rhoR_star* vR_star + (B%R(2 ) - B%L(2 ))* sign_Bx)/ denom_ds
3130
3090
w_double = (sqrt_rhoL_star* wL_star + sqrt_rhoR_star* wR_star + (B%R(3 ) - B%L(3 ))* sign_Bx)/ denom_ds
3131
3091
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
3135
3095
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
3136
3096
E_double = 0.5_wp * (E_doubleL + E_doubleR)
3137
3097
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]
3153
3100
3154
3101
! (11 ) Choose HLLD flux based on wave- speed regions
3155
3102
if (0.0_wp <= s_L) then
@@ -3170,16 +3117,12 @@ contains
3170
3117
! Mass
3171
3118
flux_rs${XYZ}$_vf(j, k, l, 1 ) = F_hlld(1 ) ! TODO multi- component
3172
3119
! 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 ])
3176
3121
! Magnetic field
3177
3122
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 ])
3180
3124
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 ])
3183
3126
end if
3184
3127
! Energy
3185
3128
flux_rs${XYZ}$_vf(j, k, l, E_idx) = F_hlld(7 )
0 commit comments