q_T_sf: Chemistry Temperature Optimization

Author: henryleberre

src/common/m_chemistry.fpp

@@ -0,0 +1,104 @@
+!!>
+!! @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)
+
+        ! 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
+
+        integer :: x, y, z, eqn
+        real(kind(0d0)) :: energy, mean_molecular_weight
+        real(kind(0d0)), dimension(num_species) :: Ys
+
+        !$acc parallel loop collapse(3) gang vector default(present) &
+        !$acc private(Ys)
+        do z = idwbuff(3)%beg, idwbuff(3)%end
+            do y = idwbuff(2)%beg, idwbuff(2)%end
+                do x = idwbuff(1)%beg, idwbuff(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)
+                    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 - &
+                                 5d-1*(q_cons_vf(eqn)%sf(x, y, z)/q_cons_vf(contxb)%sf(x, y, z))**2d0
+                    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)
+
+        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
+        integer :: x, y, z
+        integer :: eqn
+        real(kind(0d0)) :: T
+        real(kind(0d0)) :: rho, omega_m
+        real(kind(0d0)), dimension(num_species) :: Ys
+        real(kind(0d0)), dimension(num_species) :: omega
+
+        !$acc parallel loop collapse(3) gang vector default(present) &
+        !$acc private(Ys, omega)
+        do z = idwint(3)%beg, idwint(3)%end
+            do y = idwint(2)%beg, idwint(2)%end
+                do x = idwint(1)%beg, idwint(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

@@ -28,6 +28,9 @@ module m_constants
     real(kind(0d0)), parameter :: broadband_spectral_level_constant = 20d0 !< The constant to scale the spectral level at the lower frequency bound
     real(kind(0d0)), parameter :: broadband_spectral_level_growth_rate = 10d0 !< The spectral level constant to correct the magnitude at each frqeuency to ensure the source is overall broadband
 
+    ! Chemistry
+    real(kind(0d0)), parameter :: dflt_T_guess = 1200d0 ! 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,7 +128,8 @@ contains
         real(kind(0d0)), intent(in) :: energy, alf
         real(kind(0d0)), intent(in) :: dyn_p
         real(kind(0d0)), intent(in) :: pi_inf, gamma, rho, qv
-        real(kind(0d0)), intent(out) :: pres, T
+        real(kind(0d0)), intent(out) :: pres
+        real(kind(0d0)), intent(inout) :: T
         real(kind(0d0)), intent(in), optional :: stress, mom, G
 
         ! Chemistry
@@ -183,7 +184,7 @@ contains
             e_Per_Kg = energy/rho
             Pdyn_Per_Kg = dyn_p/rho
 
-            call get_temperature(e_Per_Kg - Pdyn_Per_Kg, 1200d0, Y_rs, .true., T)
+            call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T, Y_rs, .true., T)
             call get_pressure(rho, T, Y_rs, pres)
 
         #:endif
@@ -814,11 +815,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), &
@@ -845,11 +848,11 @@ contains
 
         real(kind(0d0)) :: G_K
 
-        real(kind(0d0)) :: pres, Yksum, T
+        real(kind(0d0)) :: pres, Yksum
 
         integer :: i, j, k, l, q !< Generic loop iterators
 
-        real(kind(0.d0)) :: ntmp
+        real(kind(0.d0)) :: ntmp, T
 
         #:if MFC_CASE_OPTIMIZATION
 #ifndef MFC_SIMULATION
@@ -922,8 +925,6 @@ contains
                         do i = chemxb, chemxe
                             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(T_idx)%sf(j, k, l) = qK_cons_vf(T_idx)%sf(j, k, l)
                     else
                         !$acc loop seq
                         do i = 1, contxe
@@ -953,15 +954,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
@@ -1119,8 +1124,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

@@ -129,7 +129,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).
@@ -633,13 +632,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)
+
         ! 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

src/pre_process/m_data_output.fpp

@@ -128,13 +128,15 @@ contains
         real(kind(0d0)) :: nbub                         !< Temporary bubble number density
         real(kind(0d0)) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params
         real(kind(0d0)) :: rho                          !< Temporary density
-        real(kind(0d0)) :: pres, Temp                         !< Temporary pressure
+        real(kind(0d0)) :: pres, T                      !< Temporary pressure and temperature
 
         real(kind(0d0)) :: nR3
         real(kind(0d0)) :: ntmp
 
         real(kind(0d0)) :: rhoYks(1:num_species) !< Temporary species mass fractions
 
+        T = dflt_T_guess
+
         t_step = 0
 
         ! Outputting the Locations of the Cell-boundaries ==================
@@ -297,8 +299,6 @@ contains
                                  ((i >= adv_idx%beg) .and. (i <= adv_idx%end)) &
                                  .or. &
                                  ((i >= chemxb) .and. (i <= chemxe)) &
-                                 .or. &
-                                 ((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
@@ -310,7 +310,7 @@ contains
                                 q_cons_vf(E_idx)%sf(j, 0, 0), &
                                 q_cons_vf(alf_idx)%sf(j, 0, 0), &
                                 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho, &
-                                pi_inf, gamma, rho, qv, rhoYks, pres, Temp)
+                                pi_inf, gamma, rho, qv, rhoYks, pres, T)
                             write (2, FMT) x_cb(j), pres
                         else if ((i >= bub_idx%beg) .and. (i <= bub_idx%end) .and. bubbles) then
 
@@ -910,8 +910,6 @@ contains
             do i = 1, num_species
                 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)') T_idx, "T", "T"
         end if
 
         write (1, '(A)') ""
@@ -923,7 +921,6 @@ 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 (T_idx /= 0) write (1, '("[",I2,",",I2,"]",A)') T_idx, T_idx, " Temperature"
 
         close (1)