Skip to content

Commit aca89c0

Browse files
sbryngelsonsbryngelson
committed
refactor: extend riemann_states to HLL, LF, HLLC solvers
Replace flat rho_L/rho_R, pres_L/pres_R, E_L/E_R, H_L/H_R, gamma_L/R, pi_inf_L/R, qv_L/R, c_L/R, T_L/R, MW_L/R, R_gas_L/R, Cp_L/R, Cv_L/R, Gamm_L/R, G_L/R, Y_L/R, vel_L_rms/R_rms, vel_L_tmp/R_tmp, nbub_L/R, ptilde_L/R with type(riemann_states) :: rho, pres, E, H, ... accessed as rho%L / rho%R. Replace vel_L/vel_R with type(riemann_states_vec3) :: vel accessed as vel%L(i) / vel%R(i). HLLD already used this pattern; HLL, LF, and HLLC now match. GPU private lists updated to name the struct once instead of _L and _R.
1 parent 31f8421 commit aca89c0

2 files changed

Lines changed: 913 additions & 970 deletions

File tree

‎src/simulation/include/inline_riemann.fpp‎

Lines changed: 41 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -1,77 +1,77 @@
11
#:def arithmetic_avg()
2-
rho_avg = 5.e-1_wp*(rho_L + rho_R)
2+
rho_avg = 5.e-1_wp*(rho%L + rho%R)
33
vel_avg_rms = 0._wp
44
$:GPU_LOOP(parallelism='[seq]')
55
do i = 1, num_vels
6-
vel_avg_rms = vel_avg_rms + (5.e-1_wp*(vel_L(i) + vel_R(i)))**2._wp
6+
vel_avg_rms = vel_avg_rms + (5.e-1_wp*(vel%L(i) + vel%R(i)))**2._wp
77
end do
88

9-
H_avg = 5.e-1_wp*(H_L + H_R)
10-
gamma_avg = 5.e-1_wp*(gamma_L + gamma_R)
11-
qv_avg = 5.e-1_wp*(qv_L + qv_R)
9+
H_avg = 5.e-1_wp*(H%L + H%R)
10+
gamma_avg = 5.e-1_wp*(gamma%L + gamma%R)
11+
qv_avg = 5.e-1_wp*(qv%L + qv%R)
1212
#:enddef arithmetic_avg
1313

1414
#:def roe_avg()
15-
rho_avg = sqrt(rho_L*rho_R)
15+
rho_avg = sqrt(rho%L*rho%R)
1616

1717
vel_avg_rms = 0._wp
1818

1919
$:GPU_LOOP(parallelism='[seq]')
2020
do i = 1, num_vels
21-
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp
21+
vel_avg_rms = vel_avg_rms + (sqrt(rho%L)*vel%L(i) + sqrt(rho%R)*vel%R(i))**2._wp/(sqrt(rho%L) + sqrt(rho%R))**2._wp
2222
end do
2323

24-
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/(sqrt(rho_L) + sqrt(rho_R))
24+
H_avg = (sqrt(rho%L)*H%L + sqrt(rho%R)*H%R)/(sqrt(rho%L) + sqrt(rho%R))
2525

26-
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/(sqrt(rho_L) + sqrt(rho_R))
26+
gamma_avg = (sqrt(rho%L)*gamma%L + sqrt(rho%R)*gamma%R)/(sqrt(rho%L) + sqrt(rho%R))
2727

28-
vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2._wp/(sqrt(rho_L) + sqrt(rho_R))**2._wp
28+
vel_avg_rms = (sqrt(rho%L)*vel%L(1) + sqrt(rho%R)*vel%R(1))**2._wp/(sqrt(rho%L) + sqrt(rho%R))**2._wp
2929

30-
qv_avg = (sqrt(rho_L)*qv_L + sqrt(rho_R)*qv_R)/(sqrt(rho_L) + sqrt(rho_R))
30+
qv_avg = (sqrt(rho%L)*qv%L + sqrt(rho%R)*qv%R)/(sqrt(rho%L) + sqrt(rho%R))
3131

3232
if (chemistry) then
3333
eps = 0.001_wp
34-
call get_species_enthalpies_rt(T_L, h_iL)
35-
call get_species_enthalpies_rt(T_R, h_iR)
34+
call get_species_enthalpies_rt(T%L, h_iL)
35+
call get_species_enthalpies_rt(T%R, h_iR)
3636
#:if USING_AMD
37-
h_iL = h_iL*gas_constant/molecular_weights_nonparameter*T_L
38-
h_iR = h_iR*gas_constant/molecular_weights_nonparameter*T_R
37+
h_iL = h_iL*gas_constant/molecular_weights_nonparameter*T%L
38+
h_iR = h_iR*gas_constant/molecular_weights_nonparameter*T%R
3939
#:else
40-
h_iL = h_iL*gas_constant/molecular_weights*T_L
41-
h_iR = h_iR*gas_constant/molecular_weights*T_R
40+
h_iL = h_iL*gas_constant/molecular_weights*T%L
41+
h_iR = h_iR*gas_constant/molecular_weights*T%R
4242
#:endif
43-
call get_species_specific_heats_r(T_L, Cp_iL)
44-
call get_species_specific_heats_r(T_R, Cp_iR)
43+
call get_species_specific_heats_r(T%L, Cp_iL)
44+
call get_species_specific_heats_r(T%R, Cp_iR)
4545

46-
h_avg_2 = (sqrt(rho_L)*h_iL + sqrt(rho_R)*h_iR)/(sqrt(rho_L) + sqrt(rho_R))
47-
Yi_avg = (sqrt(rho_L)*Ys_L + sqrt(rho_R)*Ys_R)/(sqrt(rho_L) + sqrt(rho_R))
48-
T_avg = (sqrt(rho_L)*T_L + sqrt(rho_R)*T_R)/(sqrt(rho_L) + sqrt(rho_R))
46+
h_avg_2 = (sqrt(rho%L)*h_iL + sqrt(rho%R)*h_iR)/(sqrt(rho%L) + sqrt(rho%R))
47+
Yi_avg = (sqrt(rho%L)*Ys_L + sqrt(rho%R)*Ys_R)/(sqrt(rho%L) + sqrt(rho%R))
48+
T_avg = (sqrt(rho%L)*T%L + sqrt(rho%R)*T%R)/(sqrt(rho%L) + sqrt(rho%R))
4949
#:if USING_AMD
50-
if (abs(T_L - T_R) < eps) then
51-
! Case when T_L and T_R are very close
50+
if (abs(T%L - T%R) < eps) then
51+
! Case when T%L and T%R are very close
5252
Cp_avg = sum(Yi_avg(:)*(0.5_wp*Cp_iL(:) + 0.5_wp*Cp_iR(:))*gas_constant/molecular_weights_nonparameter(:))
5353
Cv_avg = sum(Yi_avg(:)*((0.5_wp*Cp_iL(:) + 0.5_wp*Cp_iR(:))*gas_constant/molecular_weights_nonparameter(:) &
5454
& - gas_constant/molecular_weights_nonparameter(:)))
5555
else
56-
! Normal calculation when T_L and T_R are sufficiently different
57-
Cp_avg = sum(Yi_avg(:)*(h_iR(:) - h_iL(:))/(T_R - T_L))
58-
Cv_avg = sum(Yi_avg(:)*((h_iR(:) - h_iL(:))/(T_R - T_L) - gas_constant/molecular_weights_nonparameter(:)))
56+
! Normal calculation when T%L and T%R are sufficiently different
57+
Cp_avg = sum(Yi_avg(:)*(h_iR(:) - h_iL(:))/(T%R - T%L))
58+
Cv_avg = sum(Yi_avg(:)*((h_iR(:) - h_iL(:))/(T%R - T%L) - gas_constant/molecular_weights_nonparameter(:)))
5959
end if
6060
gamma_avg = Cp_avg/Cv_avg
6161

6262
Phi_avg(:) = (gamma_avg - 1._wp)*(vel_avg_rms/2.0_wp - h_avg_2(:)) &
6363
& + gamma_avg*gas_constant/molecular_weights_nonparameter(:)*T_avg
6464
c_sum_Yi_Phi = sum(Yi_avg(:)*Phi_avg(:))
6565
#:else
66-
if (abs(T_L - T_R) < eps) then
67-
! Case when T_L and T_R are very close
66+
if (abs(T%L - T%R) < eps) then
67+
! Case when T%L and T%R are very close
6868
Cp_avg = sum(Yi_avg(:)*(0.5_wp*Cp_iL(:) + 0.5_wp*Cp_iR(:))*gas_constant/molecular_weights(:))
6969
Cv_avg = sum(Yi_avg(:)*((0.5_wp*Cp_iL(:) + 0.5_wp*Cp_iR(:))*gas_constant/molecular_weights(:) &
7070
& - gas_constant/molecular_weights(:)))
7171
else
72-
! Normal calculation when T_L and T_R are sufficiently different
73-
Cp_avg = sum(Yi_avg(:)*(h_iR(:) - h_iL(:))/(T_R - T_L))
74-
Cv_avg = sum(Yi_avg(:)*((h_iR(:) - h_iL(:))/(T_R - T_L) - gas_constant/molecular_weights(:)))
72+
! Normal calculation when T%L and T%R are sufficiently different
73+
Cp_avg = sum(Yi_avg(:)*(h_iR(:) - h_iL(:))/(T%R - T%L))
74+
Cv_avg = sum(Yi_avg(:)*((h_iR(:) - h_iL(:))/(T%R - T%L) - gas_constant/molecular_weights(:)))
7575
end if
7676
gamma_avg = Cp_avg/Cv_avg
7777

@@ -93,24 +93,24 @@
9393

9494
#:def compute_low_Mach_correction()
9595
if (riemann_solver == 1 .or. riemann_solver == 5) then
96-
zcoef = min(1._wp, max(vel_L_rms**5.e-1_wp/c_L, vel_R_rms**5.e-1_wp/c_R))
96+
zcoef = min(1._wp, max(vel_rms%L**5.e-1_wp/c%L, vel_rms%R**5.e-1_wp/c%R))
9797
pcorr = 0._wp
9898

9999
if (low_Mach == 1) then
100-
pcorr = -(s_P - s_M)*(rho_L + rho_R)/8._wp*(zcoef - 1._wp)
100+
pcorr = -(s_P - s_M)*(rho%L + rho%R)/8._wp*(zcoef - 1._wp)
101101
end if
102102
else if (riemann_solver == 2) then
103-
zcoef = min(1._wp, max(vel_L_rms**5.e-1_wp/c_L, vel_R_rms**5.e-1_wp/c_R))
103+
zcoef = min(1._wp, max(vel_rms%L**5.e-1_wp/c%L, vel_rms%R**5.e-1_wp/c%R))
104104
pcorr = 0._wp
105105

106106
if (low_Mach == 1) then
107-
pcorr = rho_L*rho_R*(s_L - vel_L(dir_idx(1)))*(s_R - vel_R(dir_idx(1)))*(vel_R(dir_idx(1)) - vel_L(dir_idx(1))) &
108-
& /(rho_R*(s_R - vel_R(dir_idx(1))) - rho_L*(s_L - vel_L(dir_idx(1))))*(zcoef - 1._wp)
107+
pcorr = rho%L*rho%R*(s_L - vel%L(dir_idx(1)))*(s_R - vel%R(dir_idx(1)))*(vel%R(dir_idx(1)) - vel%L(dir_idx(1))) &
108+
& /(rho%R*(s_R - vel%R(dir_idx(1))) - rho%L*(s_L - vel%L(dir_idx(1))))*(zcoef - 1._wp)
109109
else if (low_Mach == 2) then
110-
vel_L_tmp = 5.e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_L(dir_idx(1)) - vel_R(dir_idx(1))))
111-
vel_R_tmp = 5.e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + zcoef*(vel_R(dir_idx(1)) - vel_L(dir_idx(1))))
112-
vel_L(dir_idx(1)) = vel_L_tmp
113-
vel_R(dir_idx(1)) = vel_R_tmp
110+
vel_tmp%L = 5.e-1_wp*((vel%L(dir_idx(1)) + vel%R(dir_idx(1))) + zcoef*(vel%L(dir_idx(1)) - vel%R(dir_idx(1))))
111+
vel_tmp%R = 5.e-1_wp*((vel%L(dir_idx(1)) + vel%R(dir_idx(1))) + zcoef*(vel%R(dir_idx(1)) - vel%L(dir_idx(1))))
112+
vel%L(dir_idx(1)) = vel_tmp%L
113+
vel%R(dir_idx(1)) = vel_tmp%R
114114
end if
115115
end if
116116
#:enddef compute_low_Mach_correction

0 commit comments

Comments
 (0)