diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 1379081ba..77c2b1071 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -942,7 +942,7 @@ contains if (model_eqns /= 4) then qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & /rho_K - dyn_pres_K = dyn_pres_K + 5e-1_wp*qK_cons_vf(i)%sf(j, k, l) & + dyn_pres_K = dyn_pres_K + 0.5_wp*qK_cons_vf(i)%sf(j, k, l) & *qK_prim_vf(i)%sf(j, k, l) else qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & @@ -1341,7 +1341,7 @@ contains ! Computing the energy from the pressure E_K = gamma_K*pres_K + pi_inf_K & - + 5e-1_wp*rho_K*vel_K_sum + qv_K + + 0.5_wp*rho_K*vel_K_sum + qv_K ! mass flux, this should be \alpha_i \rho_i u_i !$acc loop seq @@ -1466,7 +1466,7 @@ contains (rho*(1._wp - adv(num_fluids))) end if else - c = ((H - 5e-1*vel_sum)/gamma) + c = ((H - 0.5_wp*vel_sum)/gamma) end if if (mixture_err .and. c < 0._wp) then diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index c8e37a7c4..0d611ff08 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -522,13 +522,13 @@ contains call get_mixture_energy_mass(T_L, Ys_L, E_L) call get_mixture_energy_mass(T_R, Ys_R, E_R) - E_L = rho_L*E_L + 5e-1*rho_L*vel_L_rms - E_R = rho_R*E_R + 5e-1*rho_R*vel_R_rms + E_L = rho_L*E_L + 0.5_wp*rho_L*vel_L_rms + E_R = rho_R*E_R + 0.5_wp*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R else - E_L = gamma_L*pres_L + pi_inf_L + 5e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5e-1*rho_R*vel_R_rms + qv_R + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R end if @@ -568,7 +568,7 @@ contains xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) end do - G_L = 0_wp; G_R = 0_wp; + G_L = 0._wp; G_R = 0._wp; !$acc loop seq do i = 1, num_fluids ! Mixture left and right shear modulus @@ -644,35 +644,35 @@ contains rho_R*(s_R - vel_R(dir_idx(1)))) elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - xi_M = (5e-1_wp + sign(5e-1_wp, s_L)) & - + (5e-1_wp - sign(5e-1_wp, s_L)) & - *(5e-1_wp + sign(5e-1_wp, s_R)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_R)) & - + (5e-1_wp - sign(5e-1_wp, s_L)) & - *(5e-1_wp + sign(5e-1_wp, s_R)) + xi_M = (0.5_wp + sign(0.5_wp, s_L)) & + + (0.5_wp - sign(0.5_wp, s_L)) & + *(0.5_wp + sign(0.5_wp, s_R)) + xi_P = (0.5_wp - sign(0.5_wp, s_R)) & + + (0.5_wp - sign(0.5_wp, s_L)) & + *(0.5_wp + sign(0.5_wp, s_R)) ! Mass !$acc loop seq @@ -1193,7 +1193,7 @@ contains tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) end do - G_L = 0_wp; G_R = 0_wp + G_L = 0._wp; G_R = 0._wp !$acc loop seq do i = 1, num_fluids G_L = G_L + alpha_L(i)*Gs(i) @@ -1221,7 +1221,7 @@ contains xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) end do - G_L = 0_wp; G_R = 0_wp; + G_L = 0._wp; G_R = 0._wp; !$acc loop seq do i = 1, num_fluids ! Mixture left and right shear modulus @@ -1285,25 +1285,25 @@ contains end if elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -1317,8 +1317,8 @@ contains ! goes with numerical star velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(0.5_wp, s_S)) - xi_P = (5e-1_wp - sign(0.5_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) ! goes with the numerical velocity in x/y/z directions ! xi_P/M (pressure) = min/max(0. sgn(1,sL/sR)) @@ -1335,7 +1335,7 @@ contains rho_Star = xi_M*(rho_L*(xi_MP*xi_L + 1._wp - xi_MP)) + & xi_P*(rho_R*(xi_PP*xi_R + 1._wp - xi_PP)) - vel_K_Star = vel_L(idx1)*(1_wp - xi_MP) + xi_MP*vel_R(idx1) + & + vel_K_Star = vel_L(idx1)*(1._wp - xi_MP) + xi_MP*vel_R(idx1) + & xi_MP*xi_PP*(s_S - vel_R(idx1)) ! computing fluxes @@ -1354,7 +1354,7 @@ contains do i = 1, num_dims idxi = dir_idx(i) flux_rs${XYZ}$_vf(j, k, l, contxe + idxi) = rho_Star*vel_K_Star* & - (dir_flg(idxi)*vel_K_Star + (1_wp - dir_flg(idxi))*(xi_M*vel_L(idxi) + xi_P*vel_R(idxi))) + dir_flg(idxi)*p_Star + (dir_flg(idxi)*vel_K_Star + (1._wp - dir_flg(idxi))*(xi_M*vel_L(idxi) + xi_P*vel_R(idxi))) + dir_flg(idxi)*p_Star end do ! energy flux. @@ -1363,7 +1363,7 @@ contains ! elasticity. elastic shear stress additions for the momentum and energy flux if (elasticity) then - flux_ene_e = 0_wp; + flux_ene_e = 0._wp; !$acc loop seq do i = 1, num_dims idxi = dir_idx(i) @@ -1402,10 +1402,10 @@ contains ! K-th pressure and velocity in preparation for the internal energy flux !$acc loop seq do i = 1, num_fluids - p_K_Star = xi_M*(xi_MP*((pres_L + pi_infs(i)/(1_wp + gammas(i)))* & - xi_L**(1_wp/gammas(i) + 1_wp) - pi_infs(i)/(1_wp + gammas(i)) - pres_L) + pres_L) + & - xi_P*(xi_PP*((pres_R + pi_infs(i)/(1_wp + gammas(i)))* & - xi_R**(1_wp/gammas(i) + 1_wp) - pi_infs(i)/(1_wp + gammas(i)) - pres_R) + pres_R) + p_K_Star = xi_M*(xi_MP*((pres_L + pi_infs(i)/(1._wp + gammas(i)))* & + xi_L**(1._wp/gammas(i) + 1._wp) - pi_infs(i)/(1._wp + gammas(i)) - pres_L) + pres_L) + & + xi_P*(xi_PP*((pres_R + pi_infs(i)/(1._wp + gammas(i)))* & + xi_R**(1._wp/gammas(i) + 1._wp) - pi_infs(i)/(1._wp + gammas(i)) - pres_R) + pres_R) flux_rs${XYZ}$_vf(j, k, l, i + intxb - 1) = & ((xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1))* & @@ -1463,7 +1463,7 @@ contains ! Geometrical source of the void fraction(s) is zero !$acc loop seq do i = advxb, advxe - flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0_wp + flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp end do end if #:endif @@ -1471,7 +1471,7 @@ contains if (grid_geometry == 3) then !$acc loop seq do i = 1, sys_size - flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0_wp + flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp end do flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) = & flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) - p_Star @@ -1542,9 +1542,9 @@ contains qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do - E_L = gamma_L*pres_L + pi_inf_L + 5e-1_wp*rho_L*vel_L_rms + qv_L + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5e-1_wp*rho_R*vel_R_rms + qv_R + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -1574,25 +1574,25 @@ contains /(rho_L*(s_L - vel_L(dir_idx(1))) - & rho_R*(s_R - vel_R(dir_idx(1)))) elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -1606,8 +1606,8 @@ contains ! goes with numerical velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(5e-1_wp, s_S)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) !$acc loop seq do i = 1, contxe @@ -1850,8 +1850,8 @@ contains end if end if - E_L = gamma_L*pres_L + pi_inf_L + 5e-1_wp*rho_L*vel_L_rms - E_R = gamma_R*pres_R + pi_inf_R + 5e-1_wp*rho_R*vel_R_rms + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -1954,14 +1954,14 @@ contains if ((ptilde_L /= ptilde_L) .or. (ptilde_R /= ptilde_R)) then end if - rho_avg = 5e-1_wp*(rho_L + rho_R) - H_avg = 5e-1_wp*(H_L + H_R) - gamma_avg = 5e-1_wp*(gamma_L + gamma_R) + rho_avg = 0.5_wp*(rho_L + rho_R) + H_avg = 0.5_wp*(H_L + H_R) + gamma_avg = 0.5_wp*(gamma_L + gamma_R) vel_avg_rms = 0._wp !$acc loop seq do i = 1, num_dims - vel_avg_rms = vel_avg_rms + (5e-1_wp*(vel_L(i) + vel_R(i)))**2._wp + vel_avg_rms = vel_avg_rms + (0.5_wp*(vel_L(i) + vel_R(i)))**2._wp end do end if @@ -1999,25 +1999,25 @@ contains /(rho_L*(s_L - vel_L(dir_idx(1))) - & rho_R*(s_R - vel_R(dir_idx(1)))) elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -2031,8 +2031,8 @@ contains ! goes with numerical velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(5e-1_wp, s_S)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) if (low_Mach == 1) then @:compute_low_Mach_correction() @@ -2359,14 +2359,13 @@ contains call get_mixture_energy_mass(T_L, Ys_L, E_L) call get_mixture_energy_mass(T_R, Ys_R, E_R) - E_L = rho_L*E_L + 5e-1*rho_L*vel_L_rms - E_R = rho_R*E_R + 5e-1*rho_R*vel_R_rms + E_L = rho_L*E_L + 0.5_wp*rho_L*vel_L_rms + E_R = rho_R*E_R + 0.5_wp*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R else - E_L = gamma_L*pres_L + pi_inf_L + 5e-1*rho_L*vel_L_rms + qv_L - - E_R = gamma_R*pres_R + pi_inf_R + 5e-1*rho_R*vel_R_rms + qv_R + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -2476,25 +2475,25 @@ contains end if elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(idx1) - & - vel_R(idx1))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(idx1) - & + vel_R(idx1))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(idx1) - c_L*Ms_L s_R = vel_R(idx1) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(idx1) + vel_R(idx1)) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(idx1) + vel_R(idx1)) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -2508,8 +2507,8 @@ contains ! goes with numerical velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(5e-1_wp, s_S)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) ! Computing the hllc fluxes @@ -3442,8 +3441,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = (4._wp/3._wp)*dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -3468,8 +3467,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -3496,18 +3495,18 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j + 1, k, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j + 1, k, l)) !$acc loop seq do i = 1, 2 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(2) = 5e-1_wp*(dvelL_dx_vf(2)%sf(j, k, l) & - + dvelR_dx_vf(2)%sf(j + 1, k, l)) + dvel_avg_dx(2) = 0.5_wp*(dvelL_dx_vf(2)%sf(j, k, l) & + + dvelR_dx_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*(dvel_avg_dy(2) + & avg_vel(2)/y_cc(k))/ & @@ -3538,11 +3537,11 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j + 1, k, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j + 1, k, l)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j + 1, k, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = (dvel_avg_dy(2) + & avg_vel(2)/y_cc(k))/ & @@ -3573,12 +3572,12 @@ contains !$acc loop seq do i = 1, 3, 2 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(3) = 5e-1_wp*(dvelL_dx_vf(3)%sf(j, k, l) & - + dvelR_dx_vf(3)%sf(j + 1, k, l)) + dvel_avg_dx(3) = 0.5_wp*(dvelL_dx_vf(3)%sf(j, k, l) & + + dvelR_dx_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dz(3)/y_cc(k)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -3611,8 +3610,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j + 1, k, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dz(3)/y_cc(k)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -3642,19 +3641,19 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k + 1, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j, k + 1, l)) !$acc loop seq do i = 1, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k + 1, l)) dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k + 1, l)) end do @@ -3691,14 +3690,14 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k + 1, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j, k + 1, l)) - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k + 1, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k + 1, l)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k + 1, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k + 1, l)) tau_Re(2, 2) = (dvel_avg_dx(1) + dvel_avg_dy(2) + & avg_vel(2)/y_cb(k))/ & @@ -3726,18 +3725,18 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(3) = 5e-1_wp*(velL_vf(3)%sf(j, k, l) & - + velR_vf(3)%sf(j, k + 1, l)) + avg_vel(3) = 0.5_wp*(velL_vf(3)%sf(j, k, l) & + + velR_vf(3)%sf(j, k + 1, l)) !$acc loop seq do i = 2, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k + 1, l)) end do - dvel_avg_dy(3) = 5e-1_wp*(dvelL_dy_vf(3)%sf(j, k, l) & - + dvelR_dy_vf(3)%sf(j, k + 1, l)) + dvel_avg_dy(3) = 0.5_wp*(dvelL_dy_vf(3)%sf(j, k, l) & + + dvelR_dy_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = -(2._wp/3._wp)*dvel_avg_dz(3)/y_cb(k)/ & Re_avg_rsy_vf(k, j, l, 1) @@ -3771,8 +3770,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k + 1, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = dvel_avg_dz(3)/y_cb(k)/ & Re_avg_rsy_vf(k, j, l, 2) @@ -3803,28 +3802,28 @@ contains !$acc loop seq do i = 2, 3 - avg_vel(i) = 5e-1_wp*(velL_vf(i)%sf(j, k, l) & - + velR_vf(i)%sf(j, k, l + 1)) + avg_vel(i) = 0.5_wp*(velL_vf(i)%sf(j, k, l) & + + velR_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 1, 3, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k, l + 1)) end do do i = 2, 3 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 1, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k, l + 1)) end do tau_Re(3, 1) = (dvel_avg_dz(1)/y_cc(k) + dvel_avg_dx(3))/ & @@ -3866,17 +3865,17 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k, l + 1)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j, k, l + 1)) - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k, l + 1)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k, l + 1)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k, l + 1)) - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k, l + 1)) tau_Re(3, 3) = (dvel_avg_dx(1) & + dvel_avg_dy(2) & @@ -3966,8 +3965,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = (4._wp/3._wp)*dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -3992,8 +3991,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4023,12 +4022,12 @@ contains !$acc loop seq do i = 1, 2 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(2) = 5e-1_wp*(dvelL_dx_vf(2)%sf(j, k, l) & - + dvelR_dx_vf(2)%sf(j + 1, k, l)) + dvel_avg_dx(2) = 0.5_wp*(dvelL_dx_vf(2)%sf(j, k, l) & + + dvelR_dx_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dy(2)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -4061,8 +4060,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j + 1, k, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dy(2)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4092,12 +4091,12 @@ contains !$acc loop seq do i = 1, 3, 2 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(3) = 5e-1_wp*(dvelL_dx_vf(3)%sf(j, k, l) & - + dvelR_dx_vf(3)%sf(j + 1, k, l)) + dvel_avg_dx(3) = 0.5_wp*(dvelL_dx_vf(3)%sf(j, k, l) & + + dvelR_dx_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dz(3)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -4129,8 +4128,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j + 1, k, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dz(3)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4163,12 +4162,12 @@ contains do i = 1, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k + 1, l)) dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k + 1, l)) end do @@ -4204,11 +4203,11 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k + 1, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k + 1, l)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k + 1, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k + 1, l)) tau_Re(2, 2) = (dvel_avg_dx(1) + dvel_avg_dy(2))/ & Re_avg_rsy_vf(k, j, l, 2) @@ -4238,12 +4237,12 @@ contains !$acc loop seq do i = 2, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k + 1, l)) end do - dvel_avg_dy(3) = 5e-1_wp*(dvelL_dy_vf(3)%sf(j, k, l) & - + dvelR_dy_vf(3)%sf(j, k + 1, l)) + dvel_avg_dy(3) = 0.5_wp*(dvelL_dy_vf(3)%sf(j, k, l) & + + dvelR_dy_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = -(2._wp/3._wp)*dvel_avg_dz(3)/ & Re_avg_rsy_vf(k, j, l, 1) @@ -4276,8 +4275,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k + 1, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = dvel_avg_dz(3)/ & Re_avg_rsy_vf(k, j, l, 2) @@ -4309,22 +4308,22 @@ contains !$acc loop seq do i = 1, 3, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 2, 3 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 1, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k, l + 1)) end do tau_Re(3, 1) = (dvel_avg_dz(1) + dvel_avg_dx(3))/ & @@ -4363,14 +4362,14 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k, l + 1)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k, l + 1)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k, l + 1)) - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k, l + 1)) tau_Re(3, 3) = (dvel_avg_dx(1) & + dvel_avg_dy(2) & diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index 2fc787cf5..2a4371111 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -75,7 +75,7 @@ subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, pres = q_prim_vf(E_idx)%sf(j, k, l) - E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_sum + qv + E = gamma*pres + pi_inf + 0.5_wp*rho*vel_sum + qv ! energy adjustments for hyperelastic energy if (hyperelasticity) then @@ -246,7 +246,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) dy(k)/(abs(vel(2)) + c)) if (viscous) then - vcfl_dt = cfl_target*(min(dx(j), dy(k))**2._wp)/maxval((1/Re_l)/rho) + vcfl_dt = cfl_target*(min(dx(j), dy(k))**2._wp)/maxval((1._wp/Re_l)/rho) end if else @@ -254,7 +254,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) if (viscous) then - vcfl_dt = cfl_target*(dx(j)**2._wp)/minval(1/(rho*Re_l)) + vcfl_dt = cfl_target*(dx(j)**2._wp)/minval(1._wp/(rho*Re_l)) end if end if