diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 85f67069e6..f965f8ef3c 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -1601,12 +1601,12 @@ contains end subroutine s_finalize_variables_conversion_module #ifndef MFC_PRE_PROCESS - subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c) + subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv) $:GPU_ROUTINE(function_name='s_compute_speed_of_sound', & & parallelism='[seq]', cray_inline=True) real(wp), intent(in) :: pres - real(wp), intent(in) :: rho, gamma, pi_inf + real(wp), intent(in) :: rho, gamma, pi_inf, qv real(wp), intent(in) :: H real(wp), dimension(num_fluids), intent(in) :: adv real(wp), intent(in) :: vel_sum @@ -1642,7 +1642,7 @@ contains end do c = c/rho elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles_euler))) then - ! Sound speed for bubble mmixture to order O(\alpha) + ! Sound speed for bubble mixture to order O(\alpha) if (mpp_lim .and. (num_fluids > 1)) then c = (1._wp/gamma + 1._wp)* & @@ -1654,7 +1654,7 @@ contains (rho*(1._wp - adv(num_fluids))) end if else - c = ((H - 5.e-1*vel_sum)/gamma) + c = (H - 5.e-1*vel_sum - qv/rho)/gamma end if if (mixture_err .and. c < 0._wp) then diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index ea901fa438..41363136ca 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -1582,7 +1582,7 @@ contains impure subroutine s_write_energy_data_file(q_prim_vf, q_cons_vf) type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf, q_cons_vf real(wp) :: Elk, Egk, Elp, Egint, Vb, Vl, pres_av, Et - real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H + real(wp) :: rho, pres, dV, tmp, gamma, pi_inf, MaxMa, MaxMa_glb, maxvel, c, Ma, H, qv real(wp), dimension(num_vels) :: vel real(wp), dimension(num_fluids) :: adv integer :: i, j, k, l, s !looping indices @@ -1610,6 +1610,7 @@ contains rho = 0._wp gamma = 0._wp pi_inf = 0._wp + qv = 0._wp pres = q_prim_vf(E_idx)%sf(i, j, k) Egint = Egint + q_prim_vf(E_idx + 2)%sf(i, j, k)*(fluid_pp(2)%gamma*pres)*dV do s = 1, num_vels @@ -1625,13 +1626,14 @@ contains gamma = gamma + adv(l)*fluid_pp(l)%gamma pi_inf = pi_inf + adv(l)*fluid_pp(l)%pi_inf rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k) + qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*fluid_pp(l)%qv end do - H = ((gamma + 1._wp)*pres + pi_inf)/rho + H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho call s_compute_speed_of_sound(pres, rho, & gamma, pi_inf, & - H, adv, 0._wp, 0._wp, c) + H, adv, 0._wp, 0._wp, c, qv) Ma = maxvel/c if (Ma > MaxMa .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 151f25ccee..bfa0b1ae9e 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -700,7 +700,7 @@ contains call s_compute_speed_of_sound(pres, rho_sf(i, j, k), & gamma_sf(i, j, k), pi_inf_sf(i, j, k), & - H, adv, 0._wp, 0._wp, c) + H, adv, 0._wp, 0._wp, c, qv_sf(i, j, k)) q_sf(i, j, k) = c end do diff --git a/src/simulation/include/inline_riemann.fpp b/src/simulation/include/inline_riemann.fpp index 84e70aa7e7..3d1095324a 100644 --- a/src/simulation/include/inline_riemann.fpp +++ b/src/simulation/include/inline_riemann.fpp @@ -8,6 +8,7 @@ H_avg = 5.e-1_wp*(H_L + H_R) gamma_avg = 5.e-1_wp*(gamma_L + gamma_R) + qv_avg = 5.e-1_wp*(qv_L + qv_R) #:enddef arithmetic_avg @@ -32,6 +33,9 @@ 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 + qv_avg = (sqrt(rho_L)*qv_L + sqrt(rho_R)*qv_R)/ & + (sqrt(rho_L) + sqrt(rho_R)) + if (chemistry) then eps = 0.001_wp call get_species_enthalpies_rt(T_L, h_iL) diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 94c47f75b4..93bf4ae64a 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -854,7 +854,7 @@ contains H = (E + pres)/rho ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv) ! First-Order Spatial Derivatives of Primitive Variables diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index d3de561a35..1c9f2e2933 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -273,19 +273,20 @@ contains real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction real(wp) :: gamma !< Cell-avg. sp. heat ratio real(wp) :: pi_inf !< Cell-avg. liquid stiffness function + real(wp) :: qv !< Cell-avg. internal energy reference value real(wp) :: c !< Cell-avg. sound speed real(wp) :: H !< Cell-avg. enthalpy real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers integer :: j, k, l ! Computing Stability Criteria at Current Time-step - $:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l,vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l,vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv]') do l = 0, p do k = 0, n do j = 0, m - call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv) if (viscous) then call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) @@ -1292,7 +1293,7 @@ contains ! Compute mixture sound Speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, & - ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c) + ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv) accel = accel_mag(j - 2, k, l) end if @@ -1382,7 +1383,7 @@ contains end if ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, & - ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c) + ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv) end if end if @@ -1447,7 +1448,7 @@ contains ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, & - ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c) + ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, 0._wp, c, qv) accel = accel_mag(j - 2, k - 2, l - 2) end if diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 77c8d4e495..13fe9b35cc 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -317,6 +317,7 @@ contains real(wp) :: rho_avg real(wp) :: H_avg + real(wp) :: qv_avg real(wp) :: gamma_avg real(wp) :: c_avg @@ -358,7 +359,7 @@ contains #:for NORM_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')] if (norm_dir == ${NORM_DIR}$) then - $:GPU_PARALLEL_LOOP(collapse=3, private='[i,j,k,l,q,alpha_rho_L,alpha_rho_R,vel_L,vel_R,alpha_L,alpha_R,tau_e_L,tau_e_R,Re_L,Re_R,s_L,s_R,s_S,Ys_L,Ys_R,xi_field_L, xi_field_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, zcoef, vel_L_tmp, vel_R_tmp, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R]', copyin='[norm_dir]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[i,j,k,l,q,alpha_rho_L,alpha_rho_R,vel_L,vel_R,alpha_L,alpha_R,tau_e_L,tau_e_R,Re_L,Re_R,s_L,s_R,s_S,Ys_L,Ys_R,xi_field_L, xi_field_R, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, zcoef, vel_L_tmp, vel_R_tmp, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R]', copyin='[norm_dir]') do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end @@ -629,16 +630,16 @@ contains @:compute_average_state() call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, & - vel_L_rms, 0._wp, c_L) + vel_L_rms, 0._wp, c_L, qv_L) call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, & - vel_R_rms, 0._wp, c_R) + vel_R_rms, 0._wp, c_R, qv_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & - vel_avg_rms, c_sum_Yi_Phi, c_avg) + vel_avg_rms, c_sum_Yi_Phi, c_avg, qv_avg) if (mhd) then call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L) @@ -1356,10 +1357,10 @@ contains end if call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, & - vel_L_rms, 0._wp, c_L) + vel_L_rms, 0._wp, c_L, qv_L) call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, & - vel_R_rms, 0._wp, c_R) + vel_R_rms, 0._wp, c_R, qv_R) if (mhd) then call s_compute_fast_magnetosonic_speed(rho_L, c_L, B%L, norm_dir, c_fast%L, H_L) @@ -1937,6 +1938,7 @@ contains real(wp) :: rho_avg real(wp) :: H_avg real(wp) :: gamma_avg + real(wp) :: qv_avg real(wp) :: c_avg real(wp) :: s_L, s_R, s_M, s_P, s_S @@ -1998,7 +2000,7 @@ contains ! 6-EQUATION MODEL WITH HLLC if (model_eqns == 3) then !ME3 - $:GPU_PARALLEL_LOOP(collapse=3, private='[i,j,k,l, q, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, tau_e_L, tau_e_R, flux_ene_e, xi_field_L, xi_field_R, pcorr, zcoef, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, 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, xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP, idx1, idxi]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[i,j,k,l, q, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, tau_e_L, tau_e_R, flux_ene_e, xi_field_L, xi_field_R, pcorr, zcoef, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg,c_L, c_R, G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, 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, xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP, idx1, idxi]') do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end @@ -2156,15 +2158,15 @@ contains @:compute_average_state() call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, & - vel_L_rms, 0._wp, c_L) + vel_L_rms, 0._wp, c_L, qv_L) call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, & - vel_R_rms, 0._wp, c_R) + vel_R_rms, 0._wp, c_R, qv_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & - vel_avg_rms, 0._wp, c_avg) + vel_avg_rms, 0._wp, c_avg, qv_avg) if (viscous) then $:GPU_LOOP(parallelism='[seq]') @@ -2411,7 +2413,7 @@ contains elseif (model_eqns == 4) then !ME4 - $:GPU_PARALLEL_LOOP(collapse=3, private='[i, q, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, nbub_L, nbub_R, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, 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, xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[i, q, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, nbub_L, nbub_R, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg,c_L, c_R, G_L, G_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, 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, xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP]') do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end @@ -2472,16 +2474,16 @@ contains @:compute_average_state() call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, & - vel_L_rms, 0._wp, c_L) + vel_L_rms, 0._wp, c_L, qv_L) call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, & - vel_R_rms, 0._wp, c_R) + vel_R_rms, 0._wp, c_R, qv_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & - vel_avg_rms, 0._wp, c_avg) + vel_avg_rms, 0._wp, c_avg, qv_avg) if (wave_speeds == 1) then s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) @@ -2655,7 +2657,7 @@ contains $:END_GPU_PARALLEL_LOOP() elseif (model_eqns == 2 .and. bubbles_euler) then - $:GPU_PARALLEL_LOOP(collapse=3, private='[i, q, R0_L, R0_R, V0_L, V0_R, P0_L, P0_R, pbw_L, pbw_R, vel_L, vel_R, rho_avg, alpha_L, alpha_R, h_avg, gamma_avg, nbub_L, nbub_R, Re_L, Re_R, pcorr, zcoef, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, c_L, c_R, c_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, 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, xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP, nbub_L_denom, nbub_R_denom, Pbwr3Lbar, PbwR3Rbar, R3Lbar, R3Rbar, R3V2Lbar, R3V2Rbar]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[i, q, R0_L, R0_R, V0_L, V0_R, P0_L, P0_R, pbw_L, pbw_R, vel_L, vel_R, rho_avg, alpha_L, alpha_R, h_avg, gamma_avg, nbub_L, nbub_R, Re_L, Re_R, pcorr, zcoef, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg,c_L, c_R, c_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, 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, xi_M, xi_P, xi_L, xi_R, xi_MP, xi_PP, nbub_L_denom, nbub_R_denom, Pbwr3Lbar, PbwR3Rbar, R3Lbar, R3Rbar, R3V2Lbar, R3V2Rbar]') do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end @@ -2848,6 +2850,7 @@ contains rho_avg = 5.e-1_wp*(rho_L + rho_R) H_avg = 5.e-1_wp*(H_L + H_R) gamma_avg = 5.e-1_wp*(gamma_L + gamma_R) + qv_avg = 5.e-1_wp*(qv_L + qv_R) vel_avg_rms = 0._wp $:GPU_LOOP(parallelism='[seq]') @@ -2858,15 +2861,15 @@ contains end if call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, & - vel_L_rms, 0._wp, c_L) + vel_L_rms, 0._wp, c_L, qv_L) call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, & - vel_R_rms, 0._wp, c_R) + vel_R_rms, 0._wp, c_R, qv_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & - vel_avg_rms, 0._wp, c_avg) + vel_avg_rms, 0._wp, c_avg, qv_avg) if (viscous) then $:GPU_LOOP(parallelism='[seq]') @@ -3090,7 +3093,7 @@ contains $:END_GPU_PARALLEL_LOOP() else ! 5-EQUATION MODEL WITH HLLC - $:GPU_PARALLEL_LOOP(collapse=3, private='[Re_max, i, q, T_L, T_R, idx1, idxi, vel_L_rms, vel_R_rms, pres_L, pres_R, rho_L, gamma_L, pi_inf_L, qv_L, rho_R, gamma_R, pi_inf_R, qv_R, alpha_L_sum, alpha_R_sum, E_L, E_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, Y_L, Y_R, H_L, H_R, rho_avg, gamma_avg, H_avg, c_L, c_R, c_avg, s_P, s_M, xi_P, xi_M, xi_L, xi_R, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, s_L, s_R, s_S, vel_avg_rms, pcorr, zcoef, vel_L_tmp, vel_R_tmp, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, tau_e_L, tau_e_R, xi_field_L, xi_field_R, Yi_avg,Phi_avg, h_iL, h_iR, h_avg_2, G_L, G_R]', copyin='[is1, is2, is3]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[Re_max, i, q, T_L, T_R, idx1, idxi, vel_L_rms, vel_R_rms, pres_L, pres_R, rho_L, gamma_L, pi_inf_L, qv_L, rho_R, gamma_R, pi_inf_R, qv_R, alpha_L_sum, alpha_R_sum, E_L, E_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, Y_L, Y_R, H_L, H_R, qv_avg, rho_avg, gamma_avg, H_avg, c_L, c_R, c_avg, s_P, s_M, xi_P, xi_M, xi_L, xi_R, Ms_L, Ms_R, pres_SL, pres_SR, vel_L, vel_R, Re_L, Re_R, alpha_L, alpha_R, s_L, s_R, s_S, vel_avg_rms, pcorr, zcoef, vel_L_tmp, vel_R_tmp, Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, tau_e_L, tau_e_R, xi_field_L, xi_field_R, Yi_avg,Phi_avg, h_iL, h_iR, h_avg_2, G_L, G_R]', copyin='[is1, is2, is3]') do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end @@ -3295,15 +3298,15 @@ contains @:compute_average_state() call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, & - vel_L_rms, 0._wp, c_L) + vel_L_rms, 0._wp, c_L, qv_L) call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, & - vel_R_rms, 0._wp, c_R) + vel_R_rms, 0._wp, c_R, qv_R) !> The computation of c_avg does not require all the variables, and therefore the non '_avg' ! variables are placeholders to call the subroutine. call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & - vel_avg_rms, c_sum_Yi_Phi, c_avg) + vel_avg_rms, c_sum_Yi_Phi, c_avg, qv_avg) if (viscous) then if (chemistry) then @@ -3750,8 +3753,8 @@ contains H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) ! (2) Compute fast wave speeds - call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, 0._wp, c%L) - call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, 0._wp, c%R) + call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, 0._wp, c%L, qv%L) + call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, 0._wp, c%R, qv%R) call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L) call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R) diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 89f0b82dc9..3c81726b6b 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -92,7 +92,7 @@ contains !! @param j x index !! @param k y index !! @param l z index - subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) $:GPU_ROUTINE(function_name='s_compute_enthalpy',parallelism='[seq]', & & cray_inline=True) @@ -100,11 +100,12 @@ contains real(wp), intent(inout), dimension(num_fluids) :: alpha real(wp), intent(inout), dimension(num_vels) :: vel real(wp), intent(inout) :: rho, gamma, pi_inf, vel_sum, H, pres + real(wp), intent(out) :: qv integer, intent(in) :: j, k, l real(wp), dimension(2), intent(inout) :: Re real(wp), dimension(num_fluids) :: alpha_rho, Gs - real(wp) :: qv, E, G_local + real(wp) :: E, G_local integer :: i diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 8a852bb40b..f77666f81f 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -709,6 +709,7 @@ contains real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction real(wp) :: gamma !< Cell-avg. sp. heat ratio real(wp) :: pi_inf !< Cell-avg. liquid stiffness function + real(wp) :: qv !< Cell-avg. fluid reference energy real(wp) :: c !< Cell-avg. sound speed real(wp) :: H !< Cell-avg. enthalpy real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers @@ -725,18 +726,18 @@ contains idwint) end if - $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv]') do l = 0, p do k = 0, n do j = 0, m if (igr) then - call s_compute_enthalpy(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + call s_compute_enthalpy(q_cons_ts(1)%vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) else - call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) end if ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv) call s_compute_dt_from_cfl(vel, c, max_dt, rho, Re, j, k, l) end do