76 rho_L = std::max(small,rho_L);
77 rho_R = std::max(small,rho_R);
83 Set::Scalar T_L=0, T_R=0, p_L=0, p_R=0, gamma_L=0, gamma_R=0;
87 p_L = gas.
ComputeP(rho_L, T_L,
X, i-1, j, k);
88 gamma_L = gas.
gamma(T_L,
X, i-1, j, k);
90 p_R = gas.
ComputeP(rho_R, T_R,
X, i, j, k);
91 gamma_R = gas.
gamma(T_R,
X, i, j, k);
95 p_L = gas.
ComputeP(rho_L, T_L,
X, i, j, k);
96 gamma_L = gas.
gamma(T_L,
X, i, j, k);
98 p_R = gas.
ComputeP(rho_R, T_R,
X, i+1, j, k);
99 gamma_R = gas.
gamma(T_R,
X, i+1, j, k);
103 p_L = gas.
ComputeP(rho_L, T_L,
X, i, j-1, k);
104 gamma_L = gas.
gamma(T_L,
X, i, j-1, k);
106 p_R = gas.
ComputeP(rho_R, T_R,
X, i, j, k);
107 gamma_R = gas.
gamma(T_R,
X, i, j, k);
111 p_L = gas.
ComputeP(rho_L, T_L,
X, i, j, k);
112 gamma_L = gas.
gamma(T_L,
X, i, j, k);
114 p_R = gas.
ComputeP(rho_R, T_R,
X, i, j+1, k);
115 gamma_R = gas.
gamma(T_R,
X, i, j+1, k);
119 p_L = gas.
ComputeP(rho_L, T_L,
X, i, j, k-1);
120 gamma_L = gas.
gamma(T_L,
X, i, j, k-1);
122 p_R = gas.
ComputeP(rho_R, T_R,
X, i, j, k);
123 gamma_R = gas.
gamma(T_R,
X, i, j, k);
127 p_L = gas.
ComputeP(rho_L, T_L,
X, i, j, k);
128 gamma_L = gas.
gamma(T_L,
X, i, j, k);
130 p_R = gas.
ComputeP(rho_R, T_R,
X, i, j, k+1);
131 gamma_R = gas.
gamma(T_R,
X, i, j, k+1);
138 Set::Scalar ke_L = 0.5 * (Mn_L * Mn_L ) / (rho_L+small);
140 Set::Scalar h_TL = (ke_L + ue_L + p_L) / (rho_L+small);
142 Set::Scalar ke_R = 0.5 * (Mn_R * Mn_R ) / (rho_R+small);
144 Set::Scalar h_TR = (ke_R + ue_R + p_R) / (rho_R+small);
146 Set::Scalar u_L = Mn_L/(rho_L+small), u_R = Mn_R/(rho_R+small);
147 Set::Scalar v_L = Mt_L/(rho_L+small), v_R = Mt_R/(rho_R+small);
153 Set::Scalar u_RL = (std::sqrt(rho_L) * u_L + std::sqrt(rho_R) * u_R ) / (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
154 Set::Scalar h_RL = (std::sqrt(rho_L) * h_TL + std::sqrt(rho_R) * h_TR) / (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
157 Set::Scalar gamma_RL = (std::sqrt(rho_L) * gamma_L + std::sqrt(rho_R) * gamma_R) / (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
158 Set::Scalar a_RL_sq = std::max(0.0,(gamma_RL - 1.0) * (h_RL - 0.5 * u_RL * u_RL));
161 if (
verbose && ((a_RL_sq<0) || (a_RL_sq!=a_RL_sq)))
202 Set::Scalar v_RL = (std::sqrt(rho_L) * v_L + std::sqrt(rho_R) * v_R ) /
203 (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
204 Set::Scalar Ma = (std::abs(u_RL) + std::abs(v_RL)) / (a_RL + small);
205 Ma = std::min(Ma, 1.0);
209 Set::Scalar deltav1 = deltarho - deltap / (a_RL_sq + small);
210 Set::Scalar deltav2 = deltau + deltap / (rho_RL * a_RL + small);
211 Set::Scalar deltav3 = deltau - deltap / (rho_RL * a_RL + small);
220 Set::Scalar r22 = 0.5*rho_RL/a_RL * ( u_RL + a_RL );
221 Set::Scalar r23 = 0.5*rho_RL/a_RL * ( h_RL + a_RL*u_RL );
223 Set::Scalar r32 = -0.5*rho_RL/a_RL * ( u_RL - a_RL );
224 Set::Scalar r33 = -0.5*rho_RL/a_RL * ( h_RL - a_RL*u_RL );
234 lambda1 = fabs(lambda1);
235 lambda2 = fabs(lambda2);
236 lambda3 = fabs(lambda3);
237 if ( lambda1 < deltau ) lambda1 = 0.5*(lambda1*lambda1 + deltau*deltau)/deltau;
238 if ( lambda2 < deltau ) lambda2 = 0.5*(lambda2*lambda2 + deltau*deltau)/deltau;
239 if ( lambda3 < deltau ) lambda3 = 0.5*(lambda3*lambda3 + deltau*deltau)/deltau;
242 Set::Scalar a_L = std::sqrt(gamma_L * p_L / (rho_L + small));
243 Set::Scalar a_R = std::sqrt(gamma_R * p_R / (rho_R + small));
247 Set::Scalar fix1 = std::max(0.0, std::max(lambda1 - lambda1_L, lambda1_R - lambda1));
248 Set::Scalar fix2 = std::max(0.0, std::max(lambda2 - lambda2_L, lambda2_R - lambda2));
249 Set::Scalar fix3 = std::max(0.0, std::max(lambda3 - lambda3_L, lambda3_R - lambda3));
250 if ( lambda1 < fix1 ) lambda1 = fix1;
251 if ( lambda2 < fix2 ) lambda2 = fix2;
252 if ( lambda3 < fix3 ) lambda3 = fix3;
260 fl.
mass = (0.5*(rho_L*u_L + rho_R*u_R) - 0.5*(
261 r11*fabs(lambda1)*deltav1 +
262 r21*fabs(lambda2)*deltav2 +
263 r31*fabs(lambda3)*deltav3)
290 fl.
momentum_normal = ( 0.5*(rho_L*u_L*u_L + p_L + rho_R*u_R*u_R + p_R) - 0.5*(
291 r12*fabs(lambda1)*deltav1 +
292 r22*fabs(lambda2)*deltav2 +
293 r32*fabs(lambda3)*deltav3)
297 fl.
energy = ( 0.5*(u_L*(ke_L + p_L + ue_L) + u_R*(ke_R + p_R + ue_R)) - 0.5*
299 r13*fabs(lambda1)*deltav1 +
300 r23*fabs(lambda2)*deltav2 +
301 r33*fabs(lambda3)*deltav3)