diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index dbbf4457e1..c3a62deeab 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -1467,6 +1467,8 @@ contains real(wp) :: qv_K real(wp), dimension(2) :: Re_K real(wp) :: G_K + real(wp), dimension(num_species) :: Y_K + real(wp) :: T_K, mix_mol_weight, R_gas integer :: i, j, k, l !< Generic loop iterators @@ -1479,7 +1481,7 @@ contains ! Computing the flux variables from the primitive variables, without ! accounting for the contribution of either viscosity or capillarity #ifdef MFC_SIMULATION - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_K, vel_K, alpha_K, Re_K) + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_K, vel_K, alpha_K, Re_K, Y_K) do l = is3b, is3e do k = is2b, is2e do j = is1b, is1e @@ -1518,8 +1520,23 @@ contains end if ! 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 + + if (chemistry) then + !$acc loop seq + do i = chemxb, chemxe + Y_K(i - chemxb + 1) = qK_prim_vf(j, k, l, i) + end do + !Computing the energy from the internal energy of the mixture + call get_mixture_molecular_weight(Y_k, mix_mol_weight) + R_gas = gas_constant/mix_mol_weight + T_K = pres_K/rho_K/R_gas + call get_mixture_energy_mass(T_K, Y_K, E_K) + E_K = rho_K*E_K + 5e-1_wp*rho_K*vel_K_sum + else + ! 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 + end if ! mass flux, this should be \alpha_i \rho_i u_i !$acc loop seq @@ -1538,6 +1555,14 @@ contains ! energy flux, u(E+p) FK_vf(j, k, l, E_idx) = vel_K(dir_idx(1))*(E_K + pres_K) + ! Species advection Flux, \rho*u*Y + if (chemistry) then + !$acc loop seq + do i = 1, num_species + FK_vf(j, k, l, i - 1 + chemxb) = vel_K(dir_idx(1))*(rho_K*Y_K(i)) + end do + end if + if (riemann_solver == 1 .or. riemann_solver == 4) then !$acc loop seq do i = advxb, advxe diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 869e333bb6..f90d25ebcd 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -30,6 +30,13 @@ module m_cbc use m_compute_cbc + use m_thermochem, only: & + get_mixture_energy_mass, get_mixture_specific_heat_cv_mass, & + get_mixture_specific_heat_cp_mass, gas_constant, & + get_mixture_molecular_weight, get_species_enthalpies_rt, & + molecular_weights, get_species_specific_heats_r, & + get_mole_fractions, get_species_specific_heats_r + implicit none private; public :: s_initialize_cbc_module, s_cbc, s_finalize_cbc_module @@ -96,7 +103,8 @@ module m_cbc integer :: dj integer :: bcxb, bcxe, bcyb, bcye, bczb, bcze integer :: cbc_dir, cbc_loc - !$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc) + integer :: flux_cbc_index + !$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc, flux_cbc_index) !! GRCBC inputs for subsonic inflow and outflow conditions consisting of !! inflow velocities, pressure, density and void fraction as well as @@ -124,6 +132,13 @@ contains integer :: i logical :: is_cbc + if (chemistry) then + flux_cbc_index = sys_size + else + flux_cbc_index = adv_idx%end + end if + !$acc update device(flux_cbc_index) + call s_any_cbc_boundaries(is_cbc) if (is_cbc .eqv. .false.) return @@ -153,7 +168,7 @@ contains @:ALLOCATE(F_rsx_vf(0:buff_size, & is2%beg:is2%end, & - is3%beg:is3%end, 1:adv_idx%end)) + is3%beg:is3%end, 1:flux_cbc_index)) @:ALLOCATE(F_src_rsx_vf(0:buff_size, & is2%beg:is2%end, & @@ -163,7 +178,7 @@ contains @:ALLOCATE(flux_rsx_vf_l(-1:buff_size, & is2%beg:is2%end, & - is3%beg:is3%end, 1:adv_idx%end)) + is3%beg:is3%end, 1:flux_cbc_index)) @:ALLOCATE(flux_src_rsx_vf_l(-1:buff_size, & is2%beg:is2%end, & @@ -196,7 +211,7 @@ contains @:ALLOCATE(F_rsy_vf(0:buff_size, & is2%beg:is2%end, & - is3%beg:is3%end, 1:adv_idx%end)) + is3%beg:is3%end, 1:flux_cbc_index)) @:ALLOCATE(F_src_rsy_vf(0:buff_size, & is2%beg:is2%end, & @@ -206,7 +221,7 @@ contains @:ALLOCATE(flux_rsy_vf_l(-1:buff_size, & is2%beg:is2%end, & - is3%beg:is3%end, 1:adv_idx%end)) + is3%beg:is3%end, 1:flux_cbc_index)) @:ALLOCATE(flux_src_rsy_vf_l(-1:buff_size, & is2%beg:is2%end, & @@ -241,7 +256,7 @@ contains @:ALLOCATE(F_rsz_vf(0:buff_size, & is2%beg:is2%end, & - is3%beg:is3%end, 1:adv_idx%end)) + is3%beg:is3%end, 1:flux_cbc_index)) @:ALLOCATE(F_src_rsz_vf(0:buff_size, & is2%beg:is2%end, & @@ -251,7 +266,7 @@ contains @:ALLOCATE(flux_rsz_vf_l(-1:buff_size, & is2%beg:is2%end, & - is3%beg:is3%end, 1:adv_idx%end)) + is3%beg:is3%end, 1:flux_cbc_index)) @:ALLOCATE(flux_src_rsz_vf_l(-1:buff_size, & is2%beg:is2%end, & @@ -651,6 +666,9 @@ contains real(wp) :: qv !< Cell averaged fluid reference energy real(wp) :: c real(wp) :: Ma + real(wp) :: T, sum_Enthalpies + real(wp) :: Cv, Cp, e_mix, Mw, R_gas + real(wp), dimension(num_species) :: Ys, h_k, dYs_dt, dYs_ds, Xs, Gamma_i, Cp_i real(wp) :: vel_K_sum, vel_dv_dt_sum @@ -684,7 +702,7 @@ contains is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg) !$acc parallel loop collapse(3) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do r = is3%beg, is3%end do k = is2%beg, is2%end flux_rs${XYZ}$_vf_l(0, k, r, i) = F_rs${XYZ}$_vf(0, k, r, i) & @@ -715,7 +733,7 @@ contains is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg) !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do j = 0, 1 do r = is3%beg, is3%end do k = is2%beg, is2%end @@ -757,7 +775,7 @@ contains end if ! FD2 or FD4 of RHS at j = 0 - !$acc parallel loop collapse(2) gang vector default(present) private(alpha_rho, vel, adv, mf, dvel_ds, dadv_ds, Re_cbc, dalpha_rho_ds,dvel_dt, dadv_dt, dalpha_rho_dt,L, lambda) + !$acc parallel loop collapse(2) gang vector default(present) private(alpha_rho, vel, adv, mf, dvel_ds, dadv_ds, Re_cbc, dalpha_rho_ds,dvel_dt, dadv_dt, dalpha_rho_dt,L, lambda,Ys,dYs_dt,dYs_ds,h_k,Cp_i,Gamma_i,Xs) do r = is3%beg, is3%end do k = is2%beg, is2%end @@ -796,7 +814,33 @@ contains mf(i) = alpha_rho(i)/rho end do - E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_K_sum + if (chemistry) then + !$acc loop seq + do i = chemxb, chemxe + Ys(i - chemxb + 1) = q_prim_rs${XYZ}$_vf(0, k, r, i) + end do + + call get_mixture_molecular_weight(Ys, Mw) + R_gas = gas_constant/Mw + T = pres/rho/R_gas + call get_mixture_specific_heat_cp_mass(T, Ys, Cp) + call get_mixture_energy_mass(T, Ys, e_mix) + E = rho*e_mix + 5e-1_wp*rho*vel_K_sum + if (chem_params%gamma_method == 1) then + !> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. + call get_mole_fractions(Mw, Ys, Xs) + call get_species_specific_heats_r(T, Cp_i) + Gamma_i = Cp_i/(Cp_i - 1.0_wp) + gamma = sum(Xs(:)/(Gamma_i(:) - 1.0_wp)) + else if (chem_params%gamma_method == 2) then + !> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats. + call get_mixture_specific_heat_cv_mass(T, Ys, Cv) + gamma = 1.0_wp/(Cp/Cv - 1.0_wp) + end if + else + E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_K_sum + end if + H = (E + pres)/rho ! Compute mixture sound speed @@ -820,6 +864,13 @@ contains dadv_ds(i) = 0._wp end do + if (chemistry) then + !$acc loop seq + do i = 1, num_species + dYs_ds(i) = 0._wp + end do + end if + !$acc loop seq do j = 0, buff_size @@ -845,6 +896,15 @@ contains fd_coef_${XYZ}$ (j, cbc_loc) + & dadv_ds(i) end do + + if (chemistry) then + !$acc loop seq + do i = 1, num_species + dYs_ds(i) = q_prim_rs${XYZ}$_vf(j, k, r, chemxb - 1 + i)* & + fd_coef_${XYZ}$ (j, cbc_loc) + & + dYs_ds(i) + end do + end if end do ! First-Order Temporal Derivatives of Primitive Variables @@ -859,7 +919,7 @@ contains call s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_BUFFER) .or. & (cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_BUFFER)) then - call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) + call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_INFLOW) .or. & (cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_INFLOW)) then call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) @@ -883,7 +943,7 @@ contains end if else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_OUTFLOW) .or. & (cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_OUTFLOW)) then - call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) + call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) ! Add GRCBC for Subsonic Outflow (Pressure) if (bc_${XYZ}$%grcbc_out) then L(advxe) = c*(1._wp - Ma)*(pres - pres_out(${CBC_DIR}$))/Del_out(${CBC_DIR}$) @@ -904,7 +964,7 @@ contains call s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_SUP_OUTFLOW) .or. & (cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_SUP_OUTFLOW)) then - call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) + call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) end if ! Be careful about the cylindrical coordinate! @@ -935,6 +995,13 @@ contains vel_dv_dt_sum = vel_dv_dt_sum + vel(i)*dvel_dt(i) end do + if (chemistry) then + !$acc loop seq + do i = 1, num_species + dYs_dt(i) = -1._wp*L(chemxb + i - 1) + end do + end if + ! The treatment of void fraction source is unclear if (cyl_coord .and. cbc_dir == 2 .and. cbc_loc == 1) then !$acc loop seq @@ -978,13 +1045,31 @@ contains + rho*dvel_dt(i - contxe)) end do - flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) & - + ds(0)*(pres*dgamma_dt & - + gamma*dpres_dt & - + dpi_inf_dt & - + dqv_dt & - + rho*vel_dv_dt_sum & - + 5e-1_wp*drho_dt*vel_K_sum) + if (chemistry) then + ! Evolution of LODI equation of energy for real gases adjusted to perfect gas, doi:10.1006/jcph.2002.6990 + call get_species_enthalpies_rt(T, h_k) + sum_Enthalpies = 0._wp + !$acc loop seq + do i = 1, num_species + h_k(i) = h_k(i)*gas_constant/molecular_weights(i)*T + sum_Enthalpies = sum_Enthalpies + (rho*h_k(i) - pres*Mw/molecular_weights(i)*Cp/R_gas)*dYs_dt(i) + end do + flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) & + + ds(0)*((E/rho + pres/rho)*drho_dt + rho*vel_dv_dt_sum + Cp*T*L(2)/(c*c) + sum_Enthalpies) + !$acc loop seq + do i = 1, num_species + flux_rs${XYZ}$_vf_l(-1, k, r, i - 1 + chemxb) = flux_rs${XYZ}$_vf_l(0, k, r, chemxb + i - 1) & + + ds(0)*(drho_dt*Ys(i) + rho*dYs_dt(i)) + end do + else + flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) & + + ds(0)*(pres*dgamma_dt & + + gamma*dpres_dt & + + dpi_inf_dt & + + dqv_dt & + + rho*vel_dv_dt_sum & + + 5e-1_wp*drho_dt*vel_K_sum) + end if if (riemann_solver == 1) then !$acc loop seq @@ -1106,7 +1191,7 @@ contains end do !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do r = is3%beg, is3%end do k = is2%beg, is2%end do j = -1, buff_size @@ -1182,7 +1267,7 @@ contains end do !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do r = is3%beg, is3%end do k = is2%beg, is2%end do j = -1, buff_size @@ -1258,7 +1343,7 @@ contains end do !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do r = is3%beg, is3%end do k = is2%beg, is2%end do j = -1, buff_size @@ -1339,7 +1424,7 @@ contains if (cbc_dir == 1) then !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do r = is3%beg, is3%end do k = is2%beg, is2%end do j = -1, buff_size @@ -1390,7 +1475,7 @@ contains elseif (cbc_dir == 2) then !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do r = is3%beg, is3%end do k = is2%beg, is2%end do j = -1, buff_size @@ -1443,7 +1528,7 @@ contains else !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, advxe + do i = 1, flux_cbc_index do r = is3%beg, is3%end do k = is2%beg, is2%end do j = -1, buff_size diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp index 8f4e90d852..dc0339ec17 100644 --- a/src/simulation/m_compute_cbc.fpp +++ b/src/simulation/m_compute_cbc.fpp @@ -53,7 +53,7 @@ contains !! see pg. 13 of Thompson (1987). The nonreflecting subsonic !! buffer reduces the amplitude of any reflections caused by !! outgoing waves. - subroutine s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) + subroutine s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_buffer_L #else @@ -66,6 +66,7 @@ contains real(wp), intent(in) :: dpres_ds real(wp), dimension(num_dims), intent(in) :: dvel_ds real(wp), dimension(num_fluids), intent(in) :: dadv_ds + real(wp), dimension(num_species), intent(in) :: dYs_ds integer :: i !< Generic loop iterator @@ -90,6 +91,13 @@ contains L(advxe) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(3)))*lambda(3) & *(dpres_ds + rho*c*dvel_ds(dir_idx(1))) + if (chemistry) then + do i = chemxb, chemxe + L(i) = (5e-1_wp - 5e-1_wp*sign(1._wp, lambda(2)))*lambda(2) & + *(dYs_ds(i - chemxb + 1)) + end do + end if + end subroutine s_compute_nonreflecting_subsonic_buffer_L !> The L variables for the nonreflecting subsonic inflow CBC !! see pg. 455, Thompson (1990). This nonreflecting subsonic @@ -117,13 +125,19 @@ contains L(i) = 0._wp end do + if (chemistry) then + do i = chemxb, chemxe + L(i) = 0._wp + end do + end if + end subroutine s_compute_nonreflecting_subsonic_inflow_L !> The L variables for the nonreflecting subsonic outflow !! CBC see pg. 454 of Thompson (1990). This nonreflecting !! subsonic CBC presumes an outgoing flow and reduces the !! amplitude of any reflections caused by outgoing waves. - subroutine s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) + subroutine s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_outflow_L #else @@ -136,6 +150,7 @@ contains real(wp), intent(in) :: dpres_ds real(wp), dimension(num_dims), intent(in) :: dvel_ds real(wp), dimension(num_fluids), intent(in) :: dadv_ds + real(wp), dimension(num_species), intent(in) :: dYs_ds integer :: i !> Generic loop iterator @@ -156,6 +171,12 @@ contains ! bubble index L(advxe) = 0._wp + if (chemistry) then + do i = chemxb, chemxe + L(i) = lambda(2)*dYs_ds(i - chemxb + 1) + end do + end if + end subroutine s_compute_nonreflecting_subsonic_outflow_L !> The L variables for the force-free subsonic outflow CBC, @@ -261,13 +282,19 @@ contains L(i) = 0._wp end do + if (chemistry) then + do i = chemxb, chemxe + L(i) = 0._wp + end do + end if + end subroutine s_compute_supersonic_inflow_L !> L variables for the supersonic outflow CBC, see pg. 453 !! of Thompson (1990). For the supersonic outflow CBC, the !! flow evolution at the boundary is determined completely !! by the interior data. - subroutine s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) + subroutine s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds) #ifdef _CRAYFTN !DIR$ INLINEALWAYS s_compute_supersonic_outflow_L #else @@ -280,7 +307,7 @@ contains real(wp), intent(in) :: dpres_ds real(wp), dimension(num_dims), intent(in) :: dvel_ds real(wp), dimension(num_fluids), intent(in) :: dadv_ds - + real(wp), dimension(num_species), intent(in) :: dYs_ds integer :: i !< Generic loop iterator L(1) = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) @@ -299,6 +326,12 @@ contains L(advxe) = lambda(3)*(dpres_ds + rho*c*dvel_ds(dir_idx(1))) + if (chemistry) then + do i = chemxb, chemxe + L(i) = lambda(2)*dYs_ds(i - chemxb + 1) + end do + end if + end subroutine s_compute_supersonic_outflow_L end module m_compute_cbc