Skip to content

Commit f287cb3

Browse files
ChrisZYJChrisZYJ
committed
merge upstream
1 parent da282e0 commit f287cb3

File tree

1 file changed

+73
-88
lines changed

1 file changed

+73
-88
lines changed

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 73 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -178,46 +178,22 @@ contains
178178

179179
type(int_bounds_info), intent(IN) :: ix, iy, iz
180180

181-
if (riemann_solver == 1) then
182-
call s_hll_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
183-
dqL_prim_dy_vf, &
184-
dqL_prim_dz_vf, &
185-
qL_prim_vf, &
186-
qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
187-
dqR_prim_dy_vf, &
188-
dqR_prim_dz_vf, &
189-
qR_prim_vf, &
190-
q_prim_vf, &
191-
flux_vf, flux_src_vf, &
192-
flux_gsrc_vf, &
193-
norm_dir, ix, iy, iz)
194-
elseif (riemann_solver == 2) then
195-
call s_hllc_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
196-
dqL_prim_dy_vf, &
197-
dqL_prim_dz_vf, &
198-
qL_prim_vf, &
199-
qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
200-
dqR_prim_dy_vf, &
201-
dqR_prim_dz_vf, &
202-
qR_prim_vf, &
203-
q_prim_vf, &
204-
flux_vf, flux_src_vf, &
205-
flux_gsrc_vf, &
206-
norm_dir, ix, iy, iz)
207-
elseif (riemann_solver == 4) then
208-
call s_hlld_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
209-
dqL_prim_dy_vf, &
210-
dqL_prim_dz_vf, &
211-
qL_prim_vf, &
212-
qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
213-
dqR_prim_dy_vf, &
214-
dqR_prim_dz_vf, &
215-
qR_prim_vf, &
216-
q_prim_vf, &
217-
flux_vf, flux_src_vf, &
218-
flux_gsrc_vf, &
219-
norm_dir, ix, iy, iz)
220-
end if
181+
#:for NAME, NUM in [('hll', 1), ('hllc', 2), ('hlld', 4)]
182+
if (riemann_solver == ${NUM}$) then
183+
call s_${NAME}$_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
184+
dqL_prim_dy_vf, &
185+
dqL_prim_dz_vf, &
186+
qL_prim_vf, &
187+
qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
188+
dqR_prim_dy_vf, &
189+
dqR_prim_dz_vf, &
190+
qR_prim_vf, &
191+
q_prim_vf, &
192+
flux_vf, flux_src_vf, &
193+
flux_gsrc_vf, &
194+
norm_dir, ix, iy, iz)
195+
end if
196+
#:endfor
221197

222198
end subroutine s_riemann_solver
223199

@@ -1112,6 +1088,19 @@ contains
11121088
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
11131089
end do
11141090
end if
1091+
1092+
if (cyl_coord .and. hypoelasticity) then
1093+
! += tau_sigmasigma using HLL
1094+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) = &
1095+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) + &
1096+
(s_M*tau_e_R(4) - s_P*tau_e_L(4)) &
1097+
/(s_M - s_P)
1098+
1099+
!$acc loop seq
1100+
do i = strxb, strxe
1101+
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
1102+
end do
1103+
end if
11151104
#:endif
11161105
11171106
end do
@@ -3012,7 +3001,7 @@ contains
30123001
30133002
! Local variables:
30143003
real(wp), dimension(num_fluids) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R
3015-
type(riemann_states), dimension(num_vels) :: vel
3004+
type(riemann_states_vec3) :: vel
30163005
type(riemann_states) :: rho, pres, E, H_no_mag
30173006
type(riemann_states) :: gamma, pi_inf, qv
30183007
type(riemann_states) :: vel_rms
@@ -3065,16 +3054,12 @@ contains
30653054
30663055
! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
30673056
do i = 1, num_vels
3068-
vel(i)%L = qL_prim_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i))
3069-
vel(i)%R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, contxe + dir_idx(i))
3057+
vel%L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i))
3058+
vel%R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, contxe + dir_idx(i))
30703059
end do
30713060
3072-
vel_rms%L = 0._wp; vel_rms%R = 0._wp
3073-
3074-
do i = 1, num_vels
3075-
vel_rms%L = vel_rms%L + vel(i)%L**2._wp
3076-
vel_rms%R = vel_rms%R + vel(i)%R**2._wp
3077-
end do
3061+
vel_rms%L = sum(vel%L**2._wp)
3062+
vel_rms%R = sum(vel%R**2._wp)
30783063
30793064
do i = 1, num_fluids
30803065
alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)
@@ -3133,73 +3118,73 @@ contains
31333118
call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R)
31343119
31353120
! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
3136-
s_L = min(vel(1)%L - c_fast%L, vel(1)%R - c_fast%R)
3137-
s_R = max(vel(1)%R + c_fast%R, vel(1)%L + c_fast%L)
3121+
s_L = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
3122+
s_R = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
31383123
31393124
pTot_L = pres%L + pres_mag%L
31403125
pTot_R = pres%R + pres_mag%R
31413126
3142-
s_M = (((s_R - vel(1)%R)*rho%R*vel(1)%R - &
3143-
(s_L - vel(1)%L)*rho%L*vel(1)%L - pTot_R + pTot_L)/ &
3144-
((s_R - vel(1)%R)*rho%R - (s_L - vel(1)%L)*rho%L))
3127+
s_M = (((s_R - vel%R(1))*rho%R*vel%R(1) - &
3128+
(s_L - vel%L(1))*rho%L*vel%L(1) - pTot_R + pTot_L)/ &
3129+
((s_R - vel%R(1))*rho%R - (s_L - vel%L(1))*rho%L))
31453130
31463131
! (4) Compute star state variables
3147-
rhoL_star = rho%L*(s_L - vel(1)%L)/(s_L - s_M)
3148-
rhoR_star = rho%R*(s_R - vel(1)%R)/(s_R - s_M)
3149-
p_star = pTot_L + rho%L*(s_L - vel(1)%L)*(s_M - vel(1)%L)/(s_L - s_M)
3150-
E_starL = ((s_L - vel(1)%L)*E%L - pTot_L*vel(1)%L + p_star*s_M)/(s_L - s_M)
3151-
E_starR = ((s_R - vel(1)%R)*E%R - pTot_R*vel(1)%R + p_star*s_M)/(s_R - s_M)
3132+
rhoL_star = rho%L*(s_L - vel%L(1))/(s_L - s_M)
3133+
rhoR_star = rho%R*(s_R - vel%R(1))/(s_R - s_M)
3134+
p_star = pTot_L + rho%L*(s_L - vel%L(1))*(s_M - vel%L(1))/(s_L - s_M)
3135+
E_starL = ((s_L - vel%L(1))*E%L - pTot_L*vel%L(1) + p_star*s_M)/(s_L - s_M)
3136+
E_starR = ((s_R - vel%R(1))*E%R - pTot_R*vel%R(1) + p_star*s_M)/(s_R - s_M)
31523137
31533138
! (5) Compute the left/right conserved state vectors
31543139
U_L(1) = rho%L
3155-
U_L(2) = rho%L*vel(1)%L
3156-
U_L(3) = rho%L*vel(2)%L
3157-
U_L(4) = rho%L*vel(3)%L
3140+
U_L(2) = rho%L*vel%L(1)
3141+
U_L(3) = rho%L*vel%L(2)
3142+
U_L(4) = rho%L*vel%L(3)
31583143
U_L(5) = B%L(2)
31593144
U_L(6) = B%L(3)
31603145
U_L(7) = E%L
31613146
31623147
U_R(1) = rho%R
3163-
U_R(2) = rho%R*vel(1)%R
3164-
U_R(3) = rho%R*vel(2)%R
3165-
U_R(4) = rho%R*vel(3)%R
3148+
U_R(2) = rho%R*vel%R(1)
3149+
U_R(3) = rho%R*vel%R(2)
3150+
U_R(4) = rho%R*vel%R(3)
31663151
U_R(5) = B%R(2)
31673152
U_R(6) = B%R(3)
31683153
U_R(7) = E%R
31693154
31703155
! (6) Compute the left/right star state vectors
31713156
U_starL(1) = rhoL_star
31723157
U_starL(2) = rhoL_star*s_M
3173-
U_starL(3) = rhoL_star*vel(2)%L
3174-
U_starL(4) = rhoL_star*vel(3)%L
3158+
U_starL(3) = rhoL_star*vel%L(2)
3159+
U_starL(4) = rhoL_star*vel%L(3)
31753160
U_starL(5) = B%L(2)
31763161
U_starL(6) = B%L(3)
31773162
U_starL(7) = E_starL
31783163
31793164
U_starR(1) = rhoR_star
31803165
U_starR(2) = rhoR_star*s_M
3181-
U_starR(3) = rhoR_star*vel(2)%R
3182-
U_starR(4) = rhoR_star*vel(3)%R
3166+
U_starR(3) = rhoR_star*vel%R(2)
3167+
U_starR(4) = rhoR_star*vel%R(3)
31833168
U_starR(5) = B%R(2)
31843169
U_starR(6) = B%R(3)
31853170
U_starR(7) = E_starR
31863171
31873172
! (7) Compute the left/right fluxes
3188-
F_L(1) = rho%L*vel(1)%L
3189-
F_L(2) = rho%L*vel(1)%L*vel(1)%L - B%L(1)*B%L(1) + pTot_L
3190-
F_L(3) = rho%L*vel(1)%L*vel(2)%L - B%L(1)*B%L(2)
3191-
F_L(4) = rho%L*vel(1)%L*vel(3)%L - B%L(1)*B%L(3)
3192-
F_L(5) = vel(1)%L*B%L(2) - vel(2)%L*B%L(1)
3193-
F_L(6) = vel(1)%L*B%L(3) - vel(3)%L*B%L(1)
3194-
F_L(7) = (E%L + pTot_L)*vel(1)%L - B%L(1)*(vel(1)%L*B%L(1) + vel(2)%L*B%L(2) + vel(3)%L*B%L(3))
3195-
3196-
F_R(1) = rho%R*vel(1)%R
3197-
F_R(2) = rho%R*vel(1)%R*vel(1)%R - B%R(1)*B%R(1) + pTot_R
3198-
F_R(3) = rho%R*vel(1)%R*vel(2)%R - B%R(1)*B%R(2)
3199-
F_R(4) = rho%R*vel(1)%R*vel(3)%R - B%R(1)*B%R(3)
3200-
F_R(5) = vel(1)%R*B%R(2) - vel(2)%R*B%R(1)
3201-
F_R(6) = vel(1)%R*B%R(3) - vel(3)%R*B%R(1)
3202-
F_R(7) = (E%R + pTot_R)*vel(1)%R - B%R(1)*(vel(1)%R*B%R(1) + vel(2)%R*B%R(2) + vel(3)%R*B%R(3))
3173+
F_L(1) = rho%L*vel%L(1)
3174+
F_L(2) = rho%L*vel%L(1)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
3175+
F_L(3) = rho%L*vel%L(1)*vel%L(2) - B%L(1)*B%L(2)
3176+
F_L(4) = rho%L*vel%L(1)*vel%L(3) - B%L(1)*B%L(3)
3177+
F_L(5) = vel%L(1)*B%L(2) - vel%L(2)*B%L(1)
3178+
F_L(6) = vel%L(1)*B%L(3) - vel%L(3)*B%L(1)
3179+
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))
3180+
3181+
F_R(1) = rho%R*vel%R(1)
3182+
F_R(2) = rho%R*vel%R(1)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
3183+
F_R(3) = rho%R*vel%R(1)*vel%R(2) - B%R(1)*B%R(2)
3184+
F_R(4) = rho%R*vel%R(1)*vel%R(3) - B%R(1)*B%R(3)
3185+
F_R(5) = vel%R(1)*B%R(2) - vel%R(2)*B%R(1)
3186+
F_R(6) = vel%R(1)*B%R(3) - vel%R(3)*B%R(1)
3187+
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))
32033188
32043189
! (8) Compute the left/right star fluxes (note array operations)
32053190
F_starL = F_L + s_L*(U_starL - U_L)
@@ -3214,10 +3199,10 @@ contains
32143199
sqrt_rhoR_star = sqrt(rhoR_star)
32153200
denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
32163201
sign_Bx = sign(1._wp, B%L(1))
3217-
vL_star = vel(2)%L
3218-
wL_star = vel(3)%L
3219-
vR_star = vel(2)%R
3220-
wR_star = vel(3)%R
3202+
vL_star = vel%L(2)
3203+
wL_star = vel%L(3)
3204+
vR_star = vel%R(2)
3205+
wR_star = vel%R(3)
32213206
v_double = (sqrt_rhoL_star*vL_star + sqrt_rhoR_star*vR_star + (B%R(2) - B%L(2))*sign_Bx)/denom_ds
32223207
w_double = (sqrt_rhoL_star*wL_star + sqrt_rhoR_star*wR_star + (B%R(3) - B%L(3))*sign_Bx)/denom_ds
32233208
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

0 commit comments

Comments
 (0)