HLLC Roe speed of sound

Co-Authored-By: Henry Le Berre <[email protected]>

Author: DimAdam-01

src/common/m_derived_types.fpp

@@ -313,6 +313,8 @@ module m_derived_types
 
         logical :: diffusion
         logical :: reactions
+
+        integer :: gamma_method
     end type chemistry_parameters
 
 end module m_derived_types

src/common/m_thermochem.fpp

@@ -8,7 +8,7 @@ module m_thermochem
             get_temperature, get_net_production_rates, get_pressure, &
             get_mixture_molecular_weight, get_mixture_energy_mass, &
             get_mixture_specific_heat_cv_mass, get_mixture_specific_heat_cp_mass, &
-            get_species_enthalpies_rt
+            get_species_enthalpies_rt, get_species_specific_heats_r
     #:endif
 
     implicit none

src/common/m_variables_conversion.fpp

@@ -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)
+    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c, c_avggg)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_speed_of_sound
 #else
@@ -1412,44 +1412,57 @@ contains
         real(kind(0d0)), intent(out) :: c
 
         real(kind(0d0)) :: blkmod1, blkmod2
-
+        real(kind(0d0)) :: Tolerance, c_c
+        real(kind(0d0)), optional, intent(in) :: c_avggg
         integer :: q
+        if (chemistry) then
+            if (present(c_avggg)) then
+                c_c = c_avggg
+            else
+                c_c = 0.0d0  ! Default value if 'b' is absent
+            end if
 
-        if (alt_soundspeed) then
-            blkmod1 = ((gammas(1) + 1d0)*pres + &
-                       pi_infs(1))/gammas(1)
-            blkmod2 = ((gammas(2) + 1d0)*pres + &
-                       pi_infs(2))/gammas(2)
-            c = (1d0/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
-        elseif (model_eqns == 3) then
-            c = 0d0
-            !$acc loop seq
-            do q = 1, num_fluids
-                c = c + adv(q)*(1d0/gammas(q) + 1d0)* &
-                    (pres + pi_infs(q)/(gammas(q) + 1d0))
-            end do
-            c = c/rho
-
-        elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles))) then
-            ! Sound speed for bubble mmixture to order O(\alpha)
-
-            if (mpp_lim .and. (num_fluids > 1)) then
-                c = (1d0/gamma + 1d0)* &
-                    (pres + pi_inf/(gamma + 1d0))/rho
+            if (avg_state == 1 .and. abs(c_c) > Tolerance) then
+                c = sqrt(c_c - (gamma - 1.0d0)*(vel_sum - H))
             else
-                c = &
-                    (1d0/gamma + 1d0)* &
-                    (pres + pi_inf/(gamma + 1d0))/ &
-                    (rho*(1d0 - adv(num_fluids)))
+                c = sqrt((1.0d0 + 1.0d0/gamma)*pres/rho)
             end if
         else
-            c = ((H - 5d-1*vel_sum)/gamma)
-        end if
+            if (alt_soundspeed) then
+                blkmod1 = ((gammas(1) + 1d0)*pres + &
+                           pi_infs(1))/gammas(1)
+                blkmod2 = ((gammas(2) + 1d0)*pres + &
+                           pi_infs(2))/gammas(2)
+                c = (1d0/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2)))
+            elseif (model_eqns == 3) then
+                c = 0d0
+                !$acc loop seq
+                do q = 1, num_fluids
+                    c = c + adv(q)*(1d0/gammas(q) + 1d0)* &
+                        (pres + pi_infs(q)/(gammas(q) + 1d0))
+                end do
+                c = c/rho
+            elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles))) then
+                ! Sound speed for bubble mmixture to order O(\alpha)
+
+                if (mpp_lim .and. (num_fluids > 1)) then
+                    c = (1d0/gamma + 1d0)* &
+                        (pres + pi_inf/(gamma + 1d0))/rho
+                else
+                    c = &
+                        (1d0/gamma + 1d0)* &
+                        (pres + pi_inf/(gamma + 1d0))/ &
+                        (rho*(1d0 - adv(num_fluids)))
+                end if
+            else
+                c = ((H - 5d-1*vel_sum)/gamma)
+            end if
 
-        if (mixture_err .and. c < 0d0) then
-            c = 100.d0*sgm_eps
-        else
-            c = sqrt(c)
+            if (mixture_err .and. c < 0d0) then
+                c = 100.d0*sgm_eps
+            else
+                c = sqrt(c)
+            end if
         end if
     end subroutine s_compute_speed_of_sound
 #endif

src/post_process/m_global_parameters.fpp

@@ -112,6 +112,8 @@ module m_global_parameters
     logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
     !> @}
 
+    integer :: avg_state       !< Average state evaluation method
+
     !> @name Annotations of the structure, i.e. the organization, of the state vectors
     !> @{
     type(int_bounds_info) :: cont_idx              !< Indexes of first & last continuity eqns.
@@ -379,6 +381,7 @@ contains
         schlieren_alpha = dflt_real
 
         fd_order = dflt_int
+        avg_state = dflt_int
 
         ! Tait EOS
         rhoref = dflt_real

src/post_process/m_start_up.f90

@@ -70,7 +70,7 @@ subroutine s_read_input_file
             weno_order, bc_x, &
             bc_y, bc_z, fluid_pp, format, precision, &
             hypoelasticity, G, &
-            chem_wrt_Y, chem_wrt_T, &
+            chem_wrt_Y, chem_wrt_T, avg_state, &
             alpha_rho_wrt, rho_wrt, mom_wrt, vel_wrt, &
             E_wrt, pres_wrt, alpha_wrt, gamma_wrt, &
             heat_ratio_wrt, pi_inf_wrt, pres_inf_wrt, &

src/simulation/include/inline_riemann.fpp

@@ -12,8 +12,10 @@
 #:enddef arithmetic_avg
 
 #:def roe_avg()
+
     rho_avg = sqrt(rho_L*rho_R)
     vel_avg_rms = 0d0
+
     !$acc loop seq
     do i = 1, num_dims
         vel_avg_rms = vel_avg_rms + (sqrt(rho_L)*vel_L(i) + sqrt(rho_R)*vel_R(i))**2d0/ &
@@ -26,10 +28,45 @@
     gamma_avg = (sqrt(rho_L)*gamma_L + sqrt(rho_R)*gamma_R)/ &
                 (sqrt(rho_L) + sqrt(rho_R))
 
-    rho_avg = sqrt(rho_L*rho_R)
     vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2d0/ &
                   (sqrt(rho_L) + sqrt(rho_R))**2d0
 
+    #:if chemistry
+        eps = 0.001d0
+        call get_species_enthalpies_rt(T_L, h_iL)
+        call get_species_enthalpies_rt(T_R, h_iR)
+
+        h_iL = h_iL*gas_constant/mol_weights*T_L
+        h_iR = h_iR*gas_constant/mol_weights*T_R
+        call get_species_specific_heats_r(T_L, Cp_iL)
+        call get_species_specific_heats_r(T_R, Cp_iR)
+
+        h_avg_2 = (sqrt(rho_L)*h_iL + sqrt(rho_R)*h_iR)/(sqrt(rho_L) + sqrt(rho_R))
+        Yi_avg = (sqrt(rho_L)*Ys_L + sqrt(rho_R)*Ys_R)/(sqrt(rho_L) + sqrt(rho_R))
+        T_avg = (sqrt(rho_L)*T_L + sqrt(rho_R)*T_R)/(sqrt(rho_L) + sqrt(rho_R))
+        Cp_avg = 0.0d0
+        Cv_avg = 0.0d0
+
+        do i = 1, num_species
+            if (abs(T_L - T_R) < eps) then
+                ! Case when T_L and T_R are very close
+                Cp_avg = Cp_avg + Yi_avg(i)*(0.5d0*Cp_iL(i) + 0.5d0*Cp_iR(i))*gas_constant/mol_weights(i)
+                Cv_avg = Cv_avg + Yi_avg(i)*((0.5d0*Cp_iL(i) + 0.5d0*Cp_iR(i))*gas_constant/mol_weights(i) - gas_constant/mol_weights(i))
+            else
+                ! Normal calculation when T_L and T_R are sufficiently different
+                Cp_avg = Cp_avg + Yi_avg(i)*(h_iR(i) - h_iL(i))/(T_R - T_L)
+                Cv_avg = Cv_avg + Yi_avg(i)*((h_iR(i) - h_iL(i))/(T_R - T_L) - gas_constant/mol_weights(i))
+            end if
+        end do
+        gamma_avg = Cp_avg/Cv_avg
+        c_avggg = 0.0d0
+
+        do i = 1, num_species
+            Phi_avg(i) = (gamma_avg - 1.d0)*(vel_avg_rms/2.0d0 - h_avg_2(i)) + gamma_avg*gas_constant/mol_weights(i)*T_avg
+            c_avggg = c_avggg + Yi_avg(i)*Phi_avg(i)
+        end do
+    #:endif
+
 #:enddef roe_avg
 
 #:def compute_average_state()

src/simulation/m_global_parameters.fpp

@@ -571,6 +571,7 @@ contains
 
         chem_params%diffusion = .false.
         chem_params%reactions = .false.
+        chem_params%gamma_method = 1
 
         bc_x%beg = dflt_int; bc_x%end = dflt_int
         bc_y%beg = dflt_int; bc_y%end = dflt_int

src/simulation/m_mpi_proxy.fpp

@@ -211,6 +211,10 @@ contains
             #:for VAR in [ 'diffusion', 'reactions' ]
                 call MPI_BCAST(chem_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
             #:endfor
+
+            #:for VAR in [ 'gamma_method' ]
+                call MPI_BCAST(chem_params%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            #:endfor
         end if
 
         #:for VAR in [ 'dt','weno_eps','teno_CT','pref','rhoref','R0ref','Web','Ca', 'sigma', &