q_T_sf: Chemistry Temperature Optimization (#758)

Author: henryleberre

src/common/m_chemistry.fpp

@@ -0,0 +1,106 @@
+!!>
+!! @file   m_chemistry.f90
+!! @brief  Contains module m_chemistry
+!! @author Henry Le Berre <[email protected]>
+
+#:include 'macros.fpp'
+#:include 'case.fpp'
+
+module m_chemistry
+
+    use m_thermochem, only: &
+        num_species, mol_weights, get_temperature, get_net_production_rates
+
+    use m_global_parameters
+
+    implicit none
+
+contains
+
+    subroutine s_compute_q_T_sf(q_T_sf, q_cons_vf, bounds)
+
+        ! Initialize the temperature field at the start of the simulation to
+        ! reasonable values. Temperature is computed the regular way using the
+        ! conservative variables.
+
+        type(scalar_field), intent(inout) :: q_T_sf
+        type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
+        type(int_bounds_info), dimension(1:3), intent(in) :: bounds
+
+        integer :: x, y, z, eqn
+        real(wp) :: energy, mean_molecular_weight
+        real(wp), dimension(num_species) :: Ys
+
+        do z = bounds(3)%beg, bounds(3)%end
+            do y = bounds(2)%beg, bounds(2)%end
+                do x = bounds(1)%beg, bounds(1)%end
+                    !$acc loop seq
+                    do eqn = chemxb, chemxe
+                        Ys(eqn - chemxb + 1) = &
+                            q_cons_vf(eqn)%sf(x, y, z)/q_cons_vf(contxb)%sf(x, y, z)
+                    end do
+
+                    ! e = E - 1/2*|u|^2
+                    ! cons. E_idx     = \rho E
+                    ! cons. contxb    = \rho         (1-fluid model)
+                    ! cons. momxb + i = \rho u_i
+                    energy = q_cons_vf(E_idx)%sf(x, y, z)/q_cons_vf(contxb)%sf(x, y, z)
+                    !$acc loop seq
+                    do eqn = momxb, momxe
+                        energy = energy - &
+                                 0.5_wp*(q_cons_vf(eqn)%sf(x, y, z)/q_cons_vf(contxb)%sf(x, y, z))**2._wp
+                    end do
+
+                    call get_temperature(energy, dflt_T_guess, Ys, .true., q_T_sf%sf(x, y, z))
+                end do
+            end do
+        end do
+
+    end subroutine s_compute_q_T_sf
+
+    subroutine s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp, q_T_sf, q_prim_qp, bounds)
+
+        type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
+        type(scalar_field), intent(inout) :: q_T_sf
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_qp, q_prim_qp
+        type(int_bounds_info), dimension(1:3), intent(in) :: bounds
+
+        integer :: x, y, z
+        integer :: eqn
+        real(wp) :: T
+        real(wp) :: rho, omega_m
+        real(wp), dimension(num_species) :: Ys
+        real(wp), dimension(num_species) :: omega
+
+        !$acc parallel loop collapse(3) gang vector default(present) &
+        !$acc private(Ys, omega)
+        do z = bounds(3)%beg, bounds(3)%end
+            do y = bounds(2)%beg, bounds(2)%end
+                do x = bounds(1)%beg, bounds(1)%end
+
+                    !$acc loop seq
+                    do eqn = chemxb, chemxe
+                        Ys(eqn - chemxb + 1) = q_prim_qp(eqn)%sf(x, y, z)
+                    end do
+
+                    rho = q_cons_qp(contxe)%sf(x, y, z)
+                    T = q_T_sf%sf(x, y, z)
+
+                    call get_net_production_rates(rho, T, Ys, omega)
+
+                    !$acc loop seq
+                    do eqn = chemxb, chemxe
+
+                        omega_m = mol_weights(eqn - chemxb + 1)*omega(eqn - chemxb + 1)
+
+                        rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m
+
+                    end do
+
+                end do
+            end do
+        end do
+
+    end subroutine s_compute_chemistry_reaction_flux
+
+end module m_chemistry

src/common/m_constants.fpp

@@ -30,6 +30,9 @@ module m_constants
     real(wp), parameter :: broadband_spectral_level_constant = 20._wp !< The constant to scale the spectral level at the lower frequency bound
     real(wp), parameter :: broadband_spectral_level_growth_rate = 10._wp !< The spectral level constant to correct the magnitude at each frqeuency to ensure the source is overall broadband
 
+    ! Chemistry
+    real(wp), parameter :: dflt_T_guess = 1200._wp ! Default guess for temperature (when a previous value is not available)
+
     ! IBM+STL interpolation constants
     integer, parameter :: Ifactor_2D = 50 !< Multiple factor of the ratio (edge to cell width) for interpolation along edges for 2D models
     integer, parameter :: Ifactor_3D = 5 !< Multiple factor of the ratio (edge to cell width) for interpolation along edges for 3D models

src/common/m_variables_conversion.fpp

@@ -128,13 +128,15 @@ contains
         real(wp), intent(in) :: energy, alf
         real(wp), intent(in) :: dyn_p
         real(wp), intent(in) :: pi_inf, gamma, rho, qv
-        real(wp), intent(out) :: pres, T
+        real(wp), intent(out) :: pres
+        real(wp), intent(inout) :: T
         real(wp), intent(in), optional :: stress, mom, G
 
         ! Chemistry
         real(wp), dimension(1:num_species), intent(in) :: rhoYks
         real(wp) :: E_e
         real(wp) :: e_Per_Kg, Pdyn_Per_Kg
+        real(wp) :: T_guess
         real(wp), dimension(1:num_species) :: Y_rs
 
         integer :: s !< Generic loop iterator
@@ -183,7 +185,9 @@ contains
             e_Per_Kg = energy/rho
             Pdyn_Per_Kg = dyn_p/rho
 
-            call get_temperature(e_Per_Kg - Pdyn_Per_Kg, 1200._wp, Y_rs, .true., T)
+            T_guess = T
+
+            call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T_guess, Y_rs, .true., T)
             call get_pressure(rho, T, Y_rs, pres)
 
         #:endif
@@ -816,11 +820,13 @@ contains
         !! @param iy Index bounds in second coordinate direction
         !! @param iz Index bounds in third coordinate direction
     subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, &
+                                                             q_T_sf, &
                                                              qK_prim_vf, &
                                                              ibounds, &
                                                              gm_alphaK_vf)
 
         type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf
+        type(scalar_field), intent(inout) :: q_T_sf
         type(scalar_field), dimension(sys_size), intent(inout) :: qK_prim_vf
         type(int_bounds_info), dimension(1:3), intent(in) :: ibounds
         type(scalar_field), &
@@ -847,11 +853,11 @@ contains
 
         real(wp) :: G_K
 
-        real(wp) :: pres, Yksum, T
+        real(wp) :: pres, Yksum
 
         integer :: i, j, k, l, q !< Generic loop iterators
 
-        real(wp) :: ntmp
+        real(wp) :: ntmp, T
 
         #:if MFC_CASE_OPTIMIZATION
 #ifndef MFC_SIMULATION
@@ -924,8 +930,6 @@ contains
                         do i = chemxb, chemxe
                             qK_prim_vf(i)%sf(j, k, l) = max(0._wp, qK_cons_vf(i)%sf(j, k, l)/rho_K)
                         end do
-
-                        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
@@ -955,15 +959,19 @@ contains
                         do i = 1, num_species
                             rhoYks(i) = qK_cons_vf(chemxb + i - 1)%sf(j, k, l)
                         end do
+
+                        T = q_T_sf%sf(j, k, l)
                     end if
 
                     call s_compute_pressure(qK_cons_vf(E_idx)%sf(j, k, l), &
                                             qK_cons_vf(alf_idx)%sf(j, k, l), &
-                                            dyn_pres_K, pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres, T)
+                                            dyn_pres_K, pi_inf_K, gamma_K, rho_K, &
+                                            qv_K, rhoYks, pres, T)
 
                     qK_prim_vf(E_idx)%sf(j, k, l) = pres
+
                     if (chemistry) then
-                        qK_prim_vf(T_idx)%sf(j, k, l) = T
+                        q_T_sf%sf(j, k, l) = T
                     end if
 
                     if (bubbles) then
@@ -1121,8 +1129,6 @@ contains
 
                         q_cons_vf(E_idx)%sf(j, k, l) = &
                             dyn_pres + rho*e_mix
-
-                        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

src/post_process/m_data_input.f90

@@ -51,6 +51,9 @@ end subroutine s_read_abstract_data_files
     type(scalar_field), allocatable, dimension(:), public :: q_prim_vf !<
     !! Primitive variables
 
+    type(scalar_field), public :: q_T_sf !<
+    !! Temperature field
+
     ! type(scalar_field), public :: ib_markers !<
     type(integer_field), public :: ib_markers
 
@@ -1173,6 +1176,12 @@ subroutine s_initialize_data_input_module
                                             -buff_size:p + buff_size))
                 end if
 
+                if (chemistry) then
+                    allocate (q_T_sf%sf(-buff_size:m + buff_size, &
+                                        -buff_size:n + buff_size, &
+                                        -buff_size:p + buff_size))
+                end if
+
                 ! Simulation is 2D
             else
 
@@ -1190,6 +1199,12 @@ subroutine s_initialize_data_input_module
                                             -buff_size:n + buff_size, &
                                             0:0))
                 end if
+
+                if (chemistry) then
+                    allocate (q_T_sf%sf(-buff_size:m + buff_size, &
+                                        -buff_size:n + buff_size, &
+                                        0:0))
+                end if
             end if
 
             ! Simulation is 1D
@@ -1208,6 +1223,10 @@ subroutine s_initialize_data_input_module
                 allocate (ib_markers%sf(-buff_size:m + buff_size, 0:0, 0:0))
             end if
 
+            if (chemistry) then
+                allocate (q_T_sf%sf(-buff_size:m + buff_size, 0:0, 0:0))
+            end if
+
         end if
 
         if (parallel_io .neqv. .true.) then
@@ -1236,6 +1255,10 @@ subroutine s_finalize_data_input_module
             deallocate (ib_markers%sf)
         end if
 
+        if (chemistry) then
+            deallocate (q_T_sf%sf)
+        end if
+
         s_read_data_files => null()
 
     end subroutine s_finalize_data_input_module

src/post_process/m_global_parameters.fpp

@@ -130,7 +130,6 @@ 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.
-    integer :: T_idx                               !< Index of temperature eqn.
     !> @}
 
     ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
@@ -635,13 +634,9 @@ contains
             species_idx%beg = sys_size + 1
             species_idx%end = sys_size + num_species
             sys_size = species_idx%end
-
-            T_idx = sys_size + 1
-            sys_size = T_idx
         else
             species_idx%beg = 1
             species_idx%end = 1
-            T_idx = 1
         end if
 
         momxb = mom_idx%beg

src/post_process/m_start_up.f90

@@ -39,6 +39,8 @@ module m_start_up
 
     use m_finite_differences
 
+    use m_chemistry
+
     ! ==========================================================================
 
     implicit none
@@ -168,6 +170,7 @@ subroutine s_perform_time_step(t_step)
 
         ! Populating the grid and conservative variables
         call s_read_data_files(t_step)
+
         ! Populating the buffer regions of the grid variables
         if (buff_size > 0) then
             call s_populate_grid_variables_buffer_regions()
@@ -178,8 +181,11 @@ subroutine s_perform_time_step(t_step)
             call s_populate_conservative_variables_buffer_regions()
         end if
 
+        ! Initialize the Temperature cache.
+        if (chemistry) call s_compute_q_T_sf(q_T_sf, q_cons_vf, idwbuff)
+
         ! Converting the conservative variables to the primitive ones
-        call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf, idwbuff)
+        call s_convert_conservative_to_primitive_variables(q_cons_vf, q_T_sf, q_prim_vf, idwbuff)
 
     end subroutine s_perform_time_step
 
@@ -311,9 +317,9 @@ subroutine s_save_data(t_step, varname, pres, c, H)
             end do
 
             if (chem_wrt_T) then
-                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)
+                q_sf = q_T_sf%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)

src/pre_process/m_assign_variables.fpp

@@ -189,11 +189,6 @@ contains
                     Ys(i) = q_prim_vf(chemxb + i - 1)%sf(j, k, l)
                 end do
             end block
-
-            call get_mixture_molecular_weight(Ys, mean_molecular_weight)
-            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))
         end if
 
         ! Updating the patch identities bookkeeping variable
@@ -568,10 +563,6 @@ contains
                     Ys(i) = q_prim_vf(chemxb + i - 1)%sf(j, k, l)
                 end do
             end block
-
-            call get_mixture_molecular_weight(Ys, mean_molecular_weight)
-            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))
         end if
 
         ! Set streamwise velocity to hyperbolic tangent function of y