Skip to content

Commit b652700

Browse files
sbryngelsonsbryngelson
committed
get missed avg state
1 parent a077bd6 commit b652700

File tree

2 files changed

+76
-36
lines changed

2 files changed

+76
-36
lines changed

‎src/simulation/inline_riemann.fpp‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,12 @@
88

99
H_avg = 5d-1*(H_L + H_R)
1010
gamma_avg = 5d-1*(gamma_L + gamma_R)
11+
12+
13+
1114
#:enddef arithmetic_avg
1215

16+
1317
#:def roe_avg()
1418
rho_avg = sqrt(rho_L*rho_R)
1519
vel_avg_rms = 0d0
@@ -24,6 +28,11 @@
2428

2529
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
2630
(sqrt(rho_L) + sqrt(rho_R))
31+
32+
rho_avg = sqrt(rho_L*rho_R)
33+
vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2d0/ &
34+
(sqrt(rho_L) + sqrt(rho_R))**2d0
35+
2736
#:enddef roe_avg
2837

2938
#:def compute_average_state()
@@ -40,6 +49,71 @@ end if
4049

4150
! #:def compute_wave_speeds()
4251

52+
! if (wave_speeds == 1) then
53+
! s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R)
54+
! s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L)
55+
56+
! s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* &
57+
! (s_L - vel_L(dir_idx(1))) - &
58+
! rho_R*vel_R(dir_idx(1))* &
59+
! (s_R - vel_R(dir_idx(1)))) &
60+
! /(rho_L*(s_L - vel_L(dir_idx(1))) - &
61+
! rho_R*(s_R - vel_R(dir_idx(1))))
62+
! elseif (wave_speeds == 2) then
63+
! pres_SL = 5d-1*(pres_L + pres_R + rho_avg*c_avg* &
64+
! (vel_L(dir_idx(1)) - &
65+
! vel_R(dir_idx(1))))
66+
67+
! pres_SR = pres_SL
68+
69+
! Ms_L = max(1d0, sqrt(1d0 + ((5d-1 + gamma_L)/(1d0 + gamma_L))* &
70+
! (pres_SL/pres_L - 1d0)*pres_L/ &
71+
! ((pres_L + pi_inf_L/(1d0 + gamma_L)))))
72+
! Ms_R = max(1d0, sqrt(1d0 + ((5d-1 + gamma_R)/(1d0 + gamma_R))* &
73+
! (pres_SR/pres_R - 1d0)*pres_R/ &
74+
! ((pres_R + pi_inf_R/(1d0 + gamma_R)))))
75+
76+
! s_L = vel_L(dir_idx(1)) - c_L*Ms_L
77+
! s_R = vel_R(dir_idx(1)) + c_R*Ms_R
78+
79+
! s_S = 5d-1*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + &
80+
! (pres_L - pres_R)/ &
81+
! (rho_avg*c_avg))
82+
! end if
83+
84+
85+
! if (wave_speeds == 1) then
86+
! s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R)
87+
! s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L)
88+
89+
! s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))* &
90+
! (s_L - vel_L(dir_idx(1))) - &
91+
! rho_R*vel_R(dir_idx(1))* &
92+
! (s_R - vel_R(dir_idx(1)))) &
93+
! /(rho_L*(s_L - vel_L(dir_idx(1))) - &
94+
! rho_R*(s_R - vel_R(dir_idx(1))))
95+
! elseif (wave_speeds == 2) then
96+
! pres_SL = 5d-1*(pres_L + pres_R + rho_avg*c_avg* &
97+
! (vel_L(dir_idx(1)) - &
98+
! vel_R(dir_idx(1))))
99+
100+
! pres_SR = pres_SL
101+
102+
! Ms_L = max(1d0, sqrt(1d0 + ((5d-1 + gamma_L)/(1d0 + gamma_L))* &
103+
! (pres_SL/pres_L - 1d0)*pres_L/ &
104+
! ((pres_L + pi_inf_L/(1d0 + gamma_L)))))
105+
! Ms_R = max(1d0, sqrt(1d0 + ((5d-1 + gamma_R)/(1d0 + gamma_R))* &
106+
! (pres_SR/pres_R - 1d0)*pres_R/ &
107+
! ((pres_R + pi_inf_R/(1d0 + gamma_R)))))
108+
109+
! s_L = vel_L(dir_idx(1)) - c_L*Ms_L
110+
! s_R = vel_R(dir_idx(1)) + c_R*Ms_R
111+
112+
! s_S = 5d-1*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + &
113+
! (pres_L - pres_R)/ &
114+
! (rho_avg*c_avg))
115+
! end if
116+
43117
! #:enddef compute_wave_speeds
44118

45119

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 2 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -1560,6 +1560,7 @@ contains
15601560
(pres_L - pres_R)/ &
15611561
(rho_avg*c_avg))
15621562
end if
1563+
15631564
! follows Einfeldt et al.
15641565
! s_M/P = min/max(0.,s_L/R)
15651566
s_M = min(0d0, s_L); s_P = max(0d0, s_R)
@@ -2245,43 +2246,8 @@ contains
22452246
22462247
H_L = (E_L + pres_L)/rho_L
22472248
H_R = (E_R + pres_R)/rho_R
2248-
if (avg_state == 2) then
2249-
2250-
rho_avg = 5d-1*(rho_L + rho_R)
2251-
vel_avg_rms = (5d-1*(vel_L(1) + vel_R(1)))**2d0
2252-
if (num_dims >= 2) then
2253-
vel_avg_rms = vel_avg_rms + (5d-1*(vel_L(2) + vel_R(2)))**2d0
2254-
end if
2255-
if (num_dims == 3) then
2256-
vel_avg_rms = vel_avg_rms + (5d-1*(vel_L(3) + vel_R(3)))**2d0
2257-
end if
2258-
2259-
H_avg = 5d-1*(H_L + H_R)
2260-
2261-
gamma_avg = 5d-1*(gamma_L + gamma_R)
22622249
2263-
elseif (avg_state == 1) then
2264-
2265-
rho_avg = sqrt(rho_L*rho_R)
2266-
vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2d0/ &
2267-
(sqrt(rho_L) + sqrt(rho_R))**2d0
2268-
2269-
if (num_dims >= 2) then
2270-
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(2) + sqrt(rho_R)*vel_R(2))**2d0/ &
2271-
(sqrt(rho_L) + sqrt(rho_R))**2d0
2272-
end if
2273-
2274-
if (num_dims == 3) then
2275-
vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(3) + sqrt(rho_R)*vel_R(3))**2d0/ &
2276-
(sqrt(rho_L) + sqrt(rho_R))**2d0
2277-
end if
2278-
2279-
H_avg = (sqrt(rho_L)*H_L + sqrt(rho_R)*H_R)/ &
2280-
(sqrt(rho_L) + sqrt(rho_R))
2281-
2282-
gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
2283-
(sqrt(rho_L) + sqrt(rho_R))
2284-
end if
2250+
@:compute_average_state()
22852251
22862252
if (mixture_err) then
22872253
if ((H_avg - 5d-1*vel_avg_rms) < 0d0) then

0 commit comments

Comments
 (0)