diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 0e7b88091..d7339ceaa 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -314,6 +314,9 @@ module m_derived_types logical :: diffusion logical :: reactions + !> Method of determining gamma. + !> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. + !> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats. integer :: gamma_method end type chemistry_parameters diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 32753aadb..7dfc4c07b 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -942,7 +942,7 @@ contains qK_prim_vf(i)%sf(j, k, l) = max(0d0, qK_cons_vf(i)%sf(j, k, l)/rho_K) end do - qK_prim_vf(tempxb)%sf(j, k, l) = qK_cons_vf(tempxb)%sf(j, k, l) + qK_prim_vf(T_idx)%sf(j, k, l) = qK_cons_vf(T_idx)%sf(j, k, l) else !$acc loop seq do i = 1, contxe @@ -980,7 +980,7 @@ contains qK_prim_vf(E_idx)%sf(j, k, l) = pres if (chemistry) then - qK_prim_vf(tempxb)%sf(j, k, l) = T + qK_prim_vf(T_idx)%sf(j, k, l) = T end if if (bubbles) then @@ -1133,13 +1133,13 @@ contains end do call get_mixture_molecular_weight(Ys, mix_mol_weight) - T = q_prim_vf(tempxb)%sf(j, k, l) + T = q_prim_vf(T_idx)%sf(j, k, l) call get_mixture_energy_mass(T, Ys, e_mix) q_cons_vf(E_idx)%sf(j, k, l) = & dyn_pres + rho*e_mix - q_cons_vf(tempxb)%sf(j, k, l) = q_prim_vf(tempxb)%sf(j, k, l) + q_cons_vf(T_idx)%sf(j, k, l) = q_prim_vf(T_idx)%sf(j, k, l) #:else ! Computing the energy from the pressure if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then @@ -1398,7 +1398,7 @@ 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_avggg, c) + subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c) #ifdef CRAY_ACC_WAR !DIR$ INLINEALWAYS s_compute_speed_of_sound #else @@ -1409,16 +1409,14 @@ contains real(kind(0d0)), intent(in) :: H real(kind(0d0)), dimension(num_fluids), intent(in) :: adv real(kind(0d0)), intent(in) :: vel_sum - real(kind(0d0)), intent(in) :: c_avggg + real(kind(0d0)), intent(in) :: c_c real(kind(0d0)), intent(out) :: c real(kind(0d0)) :: blkmod1, blkmod2 - real(kind(0d0)) :: Tolerance, c_c + real(kind(0d0)) :: Tolerance integer :: q if (chemistry) then - c_c = c_avggg - if (avg_state == 1 .and. abs(c_c) > Tolerance) then c = sqrt(c_c - (gamma - 1.0d0)*(vel_sum - H)) else diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index d822761dc..e74f179ce 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -129,14 +129,16 @@ module m_global_parameters type(int_bounds_info) :: stress_idx !< Indices of elastic stresses integer :: c_idx !< Index of color function type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns. - type(int_bounds_info) :: temperature_idx !< Indexes of first & last temperature eqns. + integer :: T_idx !< Index of temperature eqn. !> @} ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). + ! Stands for "InDices With BUFFer". type(int_bounds_info) :: idwint(1:3) ! Cell Indices for the entire (local) domain. In simulation, this includes ! the buffer region. idwbuff and idwint are the same otherwise. + ! Stands for "InDices With BUFFer". type(int_bounds_info) :: idwbuff(1:3) !> @name Boundary conditions in the x-, y- and z-coordinate directions @@ -284,7 +286,6 @@ module m_global_parameters integer :: bubxb, bubxe integer :: strxb, strxe integer :: chemxb, chemxe - integer :: tempxb, tempxe !> @} contains @@ -624,14 +625,12 @@ contains species_idx%end = sys_size + num_species sys_size = species_idx%end - temperature_idx%beg = sys_size + 1 - temperature_idx%end = sys_size + 1 - sys_size = temperature_idx%end + T_idx = sys_size + 1 + sys_size = T_idx else species_idx%beg = 1 species_idx%end = 1 - temperature_idx%beg = 1 - temperature_idx%end = 1 + T_idx = 1 end if momxb = mom_idx%beg @@ -648,8 +647,6 @@ contains intxe = internalEnergies_idx%end chemxb = species_idx%beg chemxe = species_idx%end - tempxb = temperature_idx%beg - tempxe = temperature_idx%end #ifdef MFC_MPI allocate (MPI_IO_DATA%view(1:sys_size)) diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index d67dcf9d8..946d6cfc5 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -312,9 +312,9 @@ subroutine s_save_data(t_step, varname, pres, c, H) end do if (chem_wrt_T) then - q_sf = q_prim_vf(tempxb)%sf(-offset_x%beg:m + offset_x%end, & - -offset_y%beg:n + offset_y%end, & - -offset_z%beg:p + offset_z%end) + q_sf = q_prim_vf(T_idx)%sf(-offset_x%beg:m + offset_x%end, & + -offset_y%beg:n + offset_y%end, & + -offset_z%beg:p + offset_z%end) write (varname, '(A)') 'T' call s_write_variable_to_formatted_database_file(varname, t_step) diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index 5aa19e6a3..f9280d86d 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -191,7 +191,7 @@ contains end block call get_mixture_molecular_weight(Ys, mean_molecular_weight) - q_prim_vf(tempxb)%sf(j, k, l) = & + q_prim_vf(T_idx)%sf(j, k, l) = & q_prim_vf(E_idx)%sf(j, k, l)*mean_molecular_weight & /(gas_constant*q_prim_vf(1)%sf(j, k, l)) #:endif @@ -570,7 +570,7 @@ contains end block call get_mixture_molecular_weight(Ys, mean_molecular_weight) - q_prim_vf(tempxb)%sf(j, k, l) = & + q_prim_vf(T_idx)%sf(j, k, l) = & q_prim_vf(E_idx)%sf(j, k, l)*mean_molecular_weight/(gas_constant*q_prim_vf(1)%sf(j, k, l)) #:endif diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 50b09604f..0ddc3deaa 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -261,7 +261,7 @@ contains .or. & ((i >= chemxb) .and. (i <= chemxe)) & .or. & - ((i == tempxb)) & + ((i == T_idx)) & ) then write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0) else if (i == mom_idx%beg) then !u @@ -825,7 +825,7 @@ contains write (1, '(I3,A20,A20)') chemxb + i - 1, "Y_{"//trim(species_names(i))//"} \rho", "Y_{"//trim(species_names(i))//"}" end do - write (1, '(I3,A20,A20)') tempxb, "T", "T" + write (1, '(I3,A20,A20)') T_idx, "T", "T" end if write (1, '(A)') "" @@ -837,7 +837,7 @@ contains if (strxb /= 0) write (1, '("[",I2,",",I2,"]",A)') strxb, strxe, " Stress" if (intxb /= 0) write (1, '("[",I2,",",I2,"]",A)') intxb, intxe, " Internal Energies" if (chemxb /= 0) write (1, '("[",I2,",",I2,"]",A)') chemxb, chemxe, " Chemistry" - if (tempxb /= 0) write (1, '("[",I2,",",I2,"]",A)') tempxb, tempxe, " Temperature" + if (T_idx /= 0) write (1, '("[",I2,",",I2,"]",A)') T_idx, T_idx, " Temperature" close (1) diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 1c032bba7..4ce08a0c4 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -106,13 +106,15 @@ module m_global_parameters type(int_bounds_info) :: stress_idx !< Indexes of elastic shear stress eqns. integer :: c_idx !< Index of the color function type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns. - type(int_bounds_info) :: temperature_idx !< Indexes of first & last temperature eqns. + integer :: T_idx !< Index of temperature eqn. ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). + ! Stands for "InDices With BUFFer". type(int_bounds_info) :: idwint(1:3) ! Cell Indices for the entire (local) domain. In simulation and post_process, ! this includes the buffer region. idwbuff and idwint are the same otherwise. + ! Stands for "InDices With BUFFer". type(int_bounds_info) :: idwbuff(1:3) type(int_bounds_info) :: bc_x, bc_y, bc_z !< @@ -237,7 +239,6 @@ module m_global_parameters integer :: bubxb, bubxe integer :: strxb, strxe integer :: chemxb, chemxe - integer :: tempxb, tempxe !> @} integer, allocatable, dimension(:, :, :) :: logic_grid @@ -707,9 +708,8 @@ contains species_idx%end = sys_size + num_species sys_size = species_idx%end - temperature_idx%beg = sys_size + 1 - temperature_idx%end = sys_size + 1 - sys_size = temperature_idx%end + T_idx = sys_size + 1 + sys_size = T_idx end if momxb = mom_idx%beg @@ -726,8 +726,6 @@ contains intxe = internalEnergies_idx%end chemxb = species_idx%beg chemxe = species_idx%end - tempxb = temperature_idx%beg - tempxe = temperature_idx%end ! Configuring Coordinate Direction Indexes ========================= idwint(1)%beg = 0; idwint(2)%beg = 0; idwint(3)%beg = 0 diff --git a/src/simulation/include/inline_riemann.fpp b/src/simulation/include/inline_riemann.fpp index 79759dd76..c0805d082 100644 --- a/src/simulation/include/inline_riemann.fpp +++ b/src/simulation/include/inline_riemann.fpp @@ -58,7 +58,7 @@ gamma_avg = Cp_avg/Cv_avg Phi_avg(:) = (gamma_avg - 1.d0)*(vel_avg_rms/2.0d0 - h_avg_2(:)) + gamma_avg*gas_constant/mol_weights(:)*T_avg - c_avggg = sum(Yi_avg(:)*Phi_avg(:)) + c_sum_Yi_Phi = sum(Yi_avg(:)*Phi_avg(:)) #:endif #:enddef roe_avg diff --git a/src/simulation/m_chemistry.fpp b/src/simulation/m_chemistry.fpp index f8fddde9a..063f38a2e 100644 --- a/src/simulation/m_chemistry.fpp +++ b/src/simulation/m_chemistry.fpp @@ -88,7 +88,7 @@ contains rhs_vf(eqn)%sf(x, y, z) = flux_x + flux_y + flux_z end do - rhs_vf(tempxb)%sf(x, y, z) = 0d0 + rhs_vf(T_idx)%sf(x, y, z) = 0d0 end do end do end do @@ -130,7 +130,7 @@ contains end do rho = q_cons_qp(contxe)%sf(x, y, z) - T = q_prim_qp(tempxb)%sf(x, y, z) + T = q_prim_qp(T_idx)%sf(x, y, z) call get_net_production_rates(rho, T, Ys, omega) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 755ce1775..dd5808ebb 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -239,17 +239,19 @@ module m_global_parameters type(int_bounds_info) :: stress_idx !< Indexes of first and last shear stress eqns. integer :: c_idx ! Index of the color function type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns. - type(int_bounds_info) :: temperature_idx !< Indexes of first & last temperature eqns. + integer :: T_idx !< Index of the temperature equation !> @} !$acc declare create(bub_idx) ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). + ! Stands for "InDices With INTerior". type(int_bounds_info) :: idwint(1:3) !$acc declare create(idwint) ! Cell Indices for the entire (local) domain. In simulation and post_process, ! this includes the buffer region. idwbuff and idwint are the same otherwise. + ! Stands for "InDices With BUFFer". type(int_bounds_info) :: idwbuff(1:3) !$acc declare create(idwbuff) @@ -301,7 +303,7 @@ module m_global_parameters integer :: startx, starty, startz - !$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, species_idx) + !$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, T_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, species_idx) ! END: Simulation Algorithm Parameters ===================================== @@ -467,9 +469,8 @@ module m_global_parameters integer :: bubxb, bubxe integer :: strxb, strxe integer :: chemxb, chemxe - integer :: tempxb, tempxe -!$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe, chemxb, chemxe, tempxb, tempxe) +!$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe, chemxb, chemxe) #ifdef CRAY_ACC_WAR @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps) @@ -1004,9 +1005,8 @@ contains species_idx%end = sys_size + num_species sys_size = species_idx%end - temperature_idx%beg = sys_size + 1 - temperature_idx%end = sys_size + 1 - sys_size = temperature_idx%end + T_idx = sys_size + 1 + sys_size = T_idx end if if (qbmm .and. .not. polytropic) then @@ -1118,10 +1118,8 @@ contains intxe = internalEnergies_idx%end chemxb = species_idx%beg chemxe = species_idx%end - tempxb = temperature_idx%beg - tempxe = temperature_idx%end - !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe, tempxb, tempxe) + !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, T_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe) !$acc update device(species_idx) !$acc update device(cfl_target, m, n, p) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index a0fa59729..095cdc0a1 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -106,6 +106,8 @@ contains @:ALLOCATE_GLOBAL(ghost_points(num_gps)) @:ALLOCATE_GLOBAL(inner_points(num_inner_gps)) + print *, "num_gps was ", num_gps + !$acc enter data copyin(ghost_points, inner_points) call s_find_ghost_points(ghost_points, inner_points) @@ -451,6 +453,8 @@ contains :: subsection_3D integer :: i, j, k, l, q !< Iterator variables + num_gps = 0 + do i = 0, m do j = 0, n if (p == 0) then diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 02a81f448..6e5123746 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -308,7 +308,7 @@ contains real(kind(0d0)), dimension(num_species) :: Ys_L, Ys_R real(kind(0d0)), dimension(num_species) :: Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR real(kind(0d0)), dimension(num_species) :: Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2 - real(kind(0d0)) :: Cp_avg, Cv_avg, gamma_avggg, T_avg, eps, c_avggg + real(kind(0d0)) :: Cp_avg, Cv_avg, gamma_avggg, T_avg, eps, c_sum_Yi_Phi real(kind(0d0)) :: T_L, T_R real(kind(0d0)) :: Y_L, Y_R real(kind(0d0)) :: MW_L, MW_R @@ -506,12 +506,14 @@ contains call get_species_specific_heats_r(T_R, Cp_iR) if (chem_params%gamma_method == 1) then + ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. Gamma_iL = Cp_iL/(Cp_iL - 1.0d0) Gamma_iR = Cp_iR/(Cp_iR - 1.0d0) gamma_L = sum(Xs_L(:)/(Gamma_iL(:) - 1.0d0)) gamma_R = sum(Xs_R(:)/(Gamma_iR(:) - 1.0d0)) 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_cp_mass(T_L, Ys_L, Cp_L) call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R) call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L) @@ -580,7 +582,7 @@ contains ! 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_avggg, c_avg) + vel_avg_rms, c_sum_Yi_Phi, c_avg) if (any(Re_size > 0)) then !$acc loop seq @@ -924,7 +926,7 @@ contains real(kind(0d0)), dimension(num_fluids) :: alpha_L, alpha_R real(kind(0d0)), dimension(num_species) :: Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR real(kind(0d0)), dimension(num_species) :: Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2 - real(kind(0d0)) :: Cp_avg, Cv_avg, T_avg, c_avggg, eps + real(kind(0d0)) :: Cp_avg, Cv_avg, T_avg, c_sum_Yi_Phi, eps real(kind(0d0)) :: T_L, T_R real(kind(0d0)) :: MW_L, MW_R real(kind(0d0)) :: R_gas_L, R_gas_R @@ -2181,7 +2183,7 @@ contains end if #:if chemistry - c_avggg = 0.0d0 + c_sum_Yi_Phi = 0.0d0 !$acc loop seq do i = chemxb, chemxe Ys_L(i - chemxb + 1) = qL_prim_rs${XYZ}$_vf(j, k, l, i) @@ -2204,12 +2206,14 @@ contains call get_species_specific_heats_r(T_R, Cp_iR) if (chem_params%gamma_method == 1) then + !> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97. Gamma_iL = Cp_iL/(Cp_iL - 1.0d0) Gamma_iR = Cp_iR/(Cp_iR - 1.0d0) gamma_L = sum(Xs_L(:)/(Gamma_iL(:) - 1.0d0)) gamma_R = sum(Xs_R(:)/(Gamma_iR(:) - 1.0d0)) 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_cp_mass(T_L, Ys_L, Cp_L) call get_mixture_specific_heat_cp_mass(T_R, Ys_R, Cp_R) call get_mixture_specific_heat_cv_mass(T_L, Ys_L, Cv_L) @@ -2238,7 +2242,6 @@ contains #:endif @:compute_average_state() - !print *, c_avggg call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, & vel_L_rms, 0d0, c_L) @@ -2250,7 +2253,7 @@ contains ! 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_avggg, c_avg) + vel_avg_rms, c_sum_Yi_Phi, c_avg) if (any(Re_size > 0)) then !$acc loop seq diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 079d65863..97cc9fdac 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -204,10 +204,10 @@ contains @:ACC_SETUP_SFs(q_prim_vf(i)) end do - @:ALLOCATE(q_prim_vf(tempxb)%sf(idwbuff(1)%beg:idwbuff(1)%end, & + @:ALLOCATE(q_prim_vf(T_idx)%sf(idwbuff(1)%beg:idwbuff(1)%end, & idwbuff(2)%beg:idwbuff(2)%end, & idwbuff(3)%beg:idwbuff(3)%end)) - @:ACC_SETUP_SFs(q_prim_vf(tempxb)) + @:ACC_SETUP_SFs(q_prim_vf(T_idx)) end if @:ALLOCATE_GLOBAL(pb_ts(1:2))