Skip to content

Commit 4b85649

Browse files
sbryngelsonsbryngelson
committed
fix: extend xi_L_m1/xi_R_m1 cancellation fix to all four HLLC blocks
Blocks 1 (model_eqns==3), 2 (model_eqns==4), and 3 (model_eqns==2 with bubbles_euler) still computed xi_L/R - 1 by subtraction, losing bits near a sonic contact where xi_L/R approx 1. Add xi_L_m1/xi_R_m1 = (s_S - u_L/R)/(s_L/R - s_S) to all three blocks and replace xi_L - 1 in mass, advection, bubble, and vel_src flux loops.
1 parent 47ac8cf commit 4b85649

1 file changed

Lines changed: 35 additions & 30 deletions

File tree

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 35 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -1815,7 +1815,7 @@ contains
18151815
& rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, &
18161816
& vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, &
18171817
& alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, vel_K_star, s_L, s_R, s_M, s_P, s_S, &
1818-
& xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP]')
1818+
& xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, xi_PP]')
18191819
do l = is3%beg, is3%end
18201820
do k = is2%beg, is2%end
18211821
do j = is1%beg, is1%end
@@ -2043,6 +2043,8 @@ contains
20432043
! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) )
20442044
xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
20452045
xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
2046+
xi_L_m1 = (s_S - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
2047+
xi_R_m1 = (s_S - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
20462048

20472049
! goes with numerical star velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star)
20482050
xi_M = (5.e-1_wp + sign(0.5_wp, s_S))
@@ -2074,8 +2076,8 @@ contains
20742076
$:GPU_LOOP(parallelism='[seq]')
20752077
do i = 1, eqn_idx%cont%end
20762078
flux_rs${XYZ}$_vf(j, k, l, i) = xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, &
2077-
& i)*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) + xi_P*qR_prim_rs${XYZ}$_vf(j &
2078-
& + 1, k, l, i)*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2079+
& i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, &
2080+
& k, l, i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
20792081
end do
20802082

20812083
! MOMENTUM FLUX. f = \rho u u - \sigma, q = \rho u, q_star = \xi * \rho*(s_star, v, w)
@@ -2122,8 +2124,8 @@ contains
21222124
do i = 1, num_dims
21232125
vel_src_rs${XYZ}$_vf(j, k, l, &
21242126
& dir_idx(i)) = xi_M*(vel_L(dir_idx(i)) + dir_flg(dir_idx(i)) &
2125-
& *(s_S*(xi_MP*(xi_L - 1) + 1) - vel_L(dir_idx(i)))) &
2126-
& + xi_P*(vel_R(dir_idx(i)) + dir_flg(dir_idx(i))*(s_S*(xi_PP*(xi_R - 1) &
2127+
& *(s_S*(xi_MP*xi_L_m1 + 1) - vel_L(dir_idx(i)))) &
2128+
& + xi_P*(vel_R(dir_idx(i)) + dir_flg(dir_idx(i))*(s_S*(xi_PP*xi_R_m1 &
21272129
& + 1) - vel_R(dir_idx(i))))
21282130
end do
21292131

@@ -2230,8 +2232,8 @@ contains
22302232
& G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, &
22312233
& vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, &
22322234
& alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, vel_K_star, s_L, s_R, s_M, s_P, s_S, &
2233-
& xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, &
2234-
& Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2]')
2235+
& xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, xi_PP, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, &
2236+
& Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2]')
22352237
do l = is3%beg, is3%end
22362238
do k = is2%beg, is2%end
22372239
do j = is1%beg, is1%end
@@ -2334,6 +2336,8 @@ contains
23342336
! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) )
23352337
xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
23362338
xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
2339+
xi_L_m1 = (s_S - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
2340+
xi_R_m1 = (s_S - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
23372341

23382342
! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star)
23392343
xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S))
@@ -2342,8 +2346,8 @@ contains
23422346
$:GPU_LOOP(parallelism='[seq]')
23432347
do i = 1, eqn_idx%cont%end
23442348
flux_rs${XYZ}$_vf(j, k, l, &
2345-
& i) = xi_M*alpha_rho_L(i)*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) &
2346-
& + xi_P*alpha_rho_R(i)*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2349+
& i) = xi_M*alpha_rho_L(i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) &
2350+
& + xi_P*alpha_rho_R(i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
23472351
end do
23482352

23492353
! Momentum flux. f = \rho u u + p I, q = \rho u, q_star = \xi * \rho*(s_star, v, w)
@@ -2374,8 +2378,8 @@ contains
23742378
$:GPU_LOOP(parallelism='[seq]')
23752379
do i = eqn_idx%alf, eqn_idx%alf ! only advect the void fraction
23762380
flux_rs${XYZ}$_vf(j, k, l, i) = xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, &
2377-
& i)*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) + xi_P*qR_prim_rs${XYZ}$_vf(j &
2378-
& + 1, k, l, i)*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2381+
& i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, &
2382+
& k, l, i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
23792383
end do
23802384

23812385
! Advection velocity source: interface velocity for volume fraction transport
@@ -2392,9 +2396,9 @@ contains
23922396
$:GPU_LOOP(parallelism='[seq]')
23932397
do i = eqn_idx%bub%beg, eqn_idx%bub%end
23942398
flux_rs${XYZ}$_vf(j, k, l, i) = xi_M*nbub_L*qL_prim_rs${XYZ}$_vf(j, k, l, &
2395-
& i)*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) &
2399+
& i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) &
23962400
& + xi_P*nbub_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, &
2397-
& i)*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2401+
& i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
23982402
end do
23992403
end if
24002404

@@ -2450,9 +2454,9 @@ contains
24502454
& rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, &
24512455
& qv_R, qv_avg, c_L, c_R, c_avg, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, &
24522456
& Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, s_L, s_R, s_M, s_P, s_S, xi_M, &
2453-
& xi_P, xi_L, xi_R, xi_MP, xi_PP, nbub_L, nbub_R, PbwR3Lbar, PbwR3Rbar, R3Lbar, R3Rbar, &
2454-
& R3V2Lbar, R3V2Rbar, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Yi_avg, &
2455-
& Phi_avg, h_iL, h_iR, h_avg_2]')
2457+
& xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, xi_PP, nbub_L, nbub_R, PbwR3Lbar, PbwR3Rbar, &
2458+
& R3Lbar, R3Rbar, R3V2Lbar, R3V2Rbar, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, &
2459+
& Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2]')
24562460
do l = is3%beg, is3%end
24572461
do k = is2%beg, is2%end
24582462
do j = is1%beg, is1%end
@@ -2691,6 +2695,8 @@ contains
26912695
! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) )
26922696
xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
26932697
xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
2698+
xi_L_m1 = (s_S - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps)
2699+
xi_R_m1 = (s_S - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps)
26942700

26952701
! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star)
26962702
xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S))
@@ -2706,8 +2712,8 @@ contains
27062712
$:GPU_LOOP(parallelism='[seq]')
27072713
do i = 1, eqn_idx%cont%end
27082714
flux_rs${XYZ}$_vf(j, k, l, i) = xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, &
2709-
& i)*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) + xi_P*qR_prim_rs${XYZ}$_vf(j &
2710-
& + 1, k, l, i)*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2715+
& i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, &
2716+
& k, l, i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
27112717
end do
27122718

27132719
if (bubbles_euler .and. (num_fluids > 1)) then
@@ -2758,17 +2764,17 @@ contains
27582764
$:GPU_LOOP(parallelism='[seq]')
27592765
do i = eqn_idx%adv%beg, eqn_idx%adv%end
27602766
flux_rs${XYZ}$_vf(j, k, l, i) = xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, &
2761-
& i)*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) + xi_P*qR_prim_rs${XYZ}$_vf(j &
2762-
& + 1, k, l, i)*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2767+
& i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, &
2768+
& k, l, i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
27632769
end do
27642770

27652771
! Advection velocity source: interface velocity for volume fraction transport
27662772
$:GPU_LOOP(parallelism='[seq]')
27672773
do i = 1, num_dims
27682774
vel_src_rs${XYZ}$_vf(j, k, l, &
2769-
& dir_idx(i)) = xi_M*(vel_L(dir_idx(i)) + dir_flg(dir_idx(i))*s_M*(xi_L &
2770-
& - 1._wp)) + xi_P*(vel_R(dir_idx(i)) + dir_flg(dir_idx(i))*s_P*(xi_R &
2771-
& - 1._wp))
2775+
& dir_idx(i)) = xi_M*(vel_L(dir_idx(i)) + dir_flg(dir_idx(i)) &
2776+
& *s_M*xi_L_m1) + xi_P*(vel_R(dir_idx(i)) + dir_flg(dir_idx(i)) &
2777+
& *s_P*xi_R_m1)
27722778

27732779
! IF ( (model_eqns == 4) .or. (num_fluids==1) ) vel_src_rs_vf(dir_idx(i))%sf(j,k,l) = 0._wp
27742780
end do
@@ -2779,21 +2785,20 @@ contains
27792785
$:GPU_LOOP(parallelism='[seq]')
27802786
do i = eqn_idx%bub%beg, eqn_idx%bub%end
27812787
flux_rs${XYZ}$_vf(j, k, l, i) = xi_M*nbub_L*qL_prim_rs${XYZ}$_vf(j, k, l, &
2782-
& i)*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) &
2783-
& + xi_P*nbub_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, &
2784-
& i)*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2788+
& i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) + xi_P*nbub_R*qR_prim_rs${XYZ}$_vf(j &
2789+
& + 1, k, l, i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
27852790
end do
27862791

27872792
if (qbmm) then
27882793
flux_rs${XYZ}$_vf(j, k, l, &
2789-
& eqn_idx%bub%beg) = xi_M*nbub_L*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) &
2790-
& + xi_P*nbub_R*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2794+
& eqn_idx%bub%beg) = xi_M*nbub_L*(vel_L(dir_idx(1)) + s_M*xi_L_m1) &
2795+
& + xi_P*nbub_R*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
27912796
end if
27922797

27932798
if (adv_n) then
27942799
flux_rs${XYZ}$_vf(j, k, l, &
2795-
& eqn_idx%n) = xi_M*nbub_L*(vel_L(dir_idx(1)) + s_M*(xi_L - 1._wp)) &
2796-
& + xi_P*nbub_R*(vel_R(dir_idx(1)) + s_P*(xi_R - 1._wp))
2800+
& eqn_idx%n) = xi_M*nbub_L*(vel_L(dir_idx(1)) + s_M*xi_L_m1) &
2801+
& + xi_P*nbub_R*(vel_R(dir_idx(1)) + s_P*xi_R_m1)
27972802
end if
27982803

27992804
! Geometrical source flux for cylindrical coordinates

0 commit comments

Comments
 (0)