Fix Chemistry advection fluxes (#752)

Author: henryleberre

examples/nD_perfect_reactor/case.py

@@ -26,7 +26,7 @@
 
 sol.TPX = 1_600, ct.one_atm, 'H2:0.04, O2:0.02, AR:0.94'
 
-Nx   = 25 * args.scale
+Nx   = int(25 * args.scale)
 Tend = 1e-4
 s    = 1e-2
 dt   = 1e-7

src/common/m_variables_conversion.fpp

@@ -23,7 +23,7 @@ module m_variables_conversion
     use m_helper
 
     use m_thermochem, only: &
-        num_species, get_temperature, get_pressure, &
+        num_species, get_temperature, get_pressure, gas_constant, &
         get_mixture_molecular_weight, get_mixture_energy_mass
 
     ! ==========================================================================
@@ -1114,7 +1114,7 @@ contains
                         end do
 
                         call get_mixture_molecular_weight(Ys, mix_mol_weight)
-                        T = q_prim_vf(T_idx)%sf(j, k, l)
+                        T = q_prim_vf(E_idx)%sf(j, k, l)*mix_mol_weight/(gas_constant*rho)
                         call get_mixture_energy_mass(T, Ys, e_mix)
 
                         q_cons_vf(E_idx)%sf(j, k, l) = &

src/post_process/p_main.fpp

@@ -44,6 +44,13 @@ program p_main
 
     ! Time-Marching Loop =======================================================
     do
+        ! If all time-steps are not ready to be post-processed and one rank is
+        ! faster than another, the slower rank processing the last available
+        ! step might be killed when the faster rank attempts to process the
+        ! first missing step, before the slower rank finishes writing the last
+        ! available step. To avoid this, we force synchronization here.
+        call s_mpi_barrier()
+
         call s_perform_time_step(t_step)
 
         call s_save_data(t_step, varname, pres, c, H)

src/simulation/m_chemistry.fpp

@@ -51,53 +51,6 @@ contains
 
     end subroutine s_finalize_chemistry_module
 
-    subroutine s_compute_chemistry_advection_flux(flux_n, rhs_vf)
-
-        type(vector_field), dimension(:), intent(IN) :: flux_n
-        type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf
-        integer :: x, y, z
-        integer :: eqn
-
-        real(kind(0d0)) :: flux_x, flux_y, flux_z
-
-        #:for num_dims in range(1, 4)
-            if (num_dims == ${num_dims}$) then
-                !$acc parallel loop collapse(4) gang vector default(present) &
-                !$acc private(flux_x, flux_y, flux_z)
-                do z = idwint(3)%beg, idwint(3)%end
-                    do y = idwint(2)%beg, idwint(2)%end
-                        do x = idwint(1)%beg, idwint(1)%end
-                            do eqn = chemxb, chemxe
-                                ! \nabla \cdot (F)
-                                flux_x = (flux_n(1)%vf(eqn)%sf(x - 1, y, z) - &
-                                          flux_n(1)%vf(eqn)%sf(x, y, z))/dx(x)
-
-                                #:if num_dims >= 2
-                                    flux_y = (flux_n(2)%vf(eqn)%sf(x, y - 1, z) - &
-                                              flux_n(2)%vf(eqn)%sf(x, y, z))/dy(y)
-                                #:else
-                                    flux_y = 0d0
-                                #:endif
-
-                                #:if num_dims == 3
-                                    flux_z = (flux_n(3)%vf(eqn)%sf(x, y, z - 1) - &
-                                              flux_n(3)%vf(eqn)%sf(x, y, z))/dz(z)
-                                #:else
-                                    flux_z = 0d0
-                                #:endif
-
-                                rhs_vf(eqn)%sf(x, y, z) = flux_x + flux_y + flux_z
-                            end do
-
-                            rhs_vf(T_idx)%sf(x, y, z) = 0d0
-                        end do
-                    end do
-                end do
-            end if
-        #:endfor
-
-    end subroutine s_compute_chemistry_advection_flux
-
     subroutine s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp, q_prim_qp)
 
         type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf, q_cons_qp, q_prim_qp

src/simulation/m_rhs.fpp

@@ -860,12 +860,6 @@ contains
         end do
         ! END: Dimensional Splitting Loop =================================
 
-        if (chemistry) then
-            call nvtxStartRange("RHS-CHEM-ADVECTION")
-            call s_compute_chemistry_advection_flux(flux_n, rhs_vf)
-            call nvtxEndRange
-        end if
-
         if (ib) then
             !$acc parallel loop collapse(3) gang vector default(present)
             do l = 0, p
@@ -904,6 +898,10 @@ contains
         end if
 
         if (chemistry) then
+            !$acc kernels
+            rhs_vf(T_idx)%sf(:, :, :) = 0d0
+            !$acc end kernels
+
             if (chem_params%reactions) then
                 call nvtxStartRange("RHS-CHEM-REACTIONS")
                 call s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp%vf, q_prim_qp%vf)

src/simulation/m_riemann_solvers.fpp

@@ -797,8 +797,8 @@ contains
                                     Y_L = qL_prim_rs${XYZ}$_vf(j, k, l, i)
                                     Y_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)
 
-                                    flux_rs${XYZ}$_vf(j, k, l, i) = (s_M*Y_R*rho_R*vel_R(dir_idx(norm_dir)) &
-                                                                     - s_P*Y_L*rho_L*vel_L(dir_idx(norm_dir)) &
+                                    flux_rs${XYZ}$_vf(j, k, l, i) = (s_M*Y_R*rho_R*vel_R(dir_idx(1)) &
+                                                                     - s_P*Y_L*rho_L*vel_L(dir_idx(1)) &
                                                                      + s_M*s_P*(Y_L*rho_L - Y_R*rho_R)) &
                                                                     /(s_M - s_P)
                                     flux_src_rs${XYZ}$_vf(j, k, l, i) = 0d0