Multicomponent diffusion fluxes, thermal conduction, and mixture viscosity

src/common/m_chemistry.fpp

@@ -10,14 +10,34 @@ module m_chemistry
 
     use m_thermochem, only: &
         num_species, molecular_weights, get_temperature, get_net_production_rates, &
-        gas_constant, get_mixture_molecular_weight
+        get_mole_fractions, get_species_binary_mass_diffusivities, &
+        get_species_mass_diffusivities_mixavg, gas_constant, get_mixture_molecular_weight, &
+        get_mixture_energy_mass, get_mixture_thermal_conductivity_mixavg, get_species_enthalpies_rt, &
+        get_mixture_viscosity_mixavg
 
     use m_global_parameters
 
     implicit none
 
+    type(int_bounds_info) :: isc1, isc2, isc3
+    !$acc declare create(isc1, isc2, isc3)
+
+    integer, dimension(3) :: offsets
+    !$acc declare create(offsets)
 contains
 
+    subroutine compute_viscosity_and_inversion(T_L, Ys_L, T_R, Ys_R, Re_L, Re_R)
+
+        real(wp), intent(inout) :: T_L, T_R, Re_L, Re_R
+        real(wp), dimension(num_species), intent(inout) :: Ys_R, Ys_L
+
+        call get_mixture_viscosity_mixavg(T_L, Ys_L, Re_L)
+        call get_mixture_viscosity_mixavg(T_R, Ys_R, Re_R)
+        Re_L = 1.0_wp/Re_L
+        Re_R = 1.0_wp/Re_R
+
+    end subroutine compute_viscosity_and_inversion
+
     subroutine s_compute_q_T_sf(q_T_sf, q_cons_vf, bounds)
 
         ! Initialize the temperature field at the start of the simulation to
@@ -130,4 +150,145 @@ contains
 
     end subroutine s_compute_chemistry_reaction_flux
 
+    subroutine s_compute_chemistry_diffusion_flux(idir, q_prim_qp, flux_src_vf, irx, iry, irz)
+
+        type(scalar_field), dimension(sys_size), intent(in) :: q_prim_qp
+        type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf
+        type(int_bounds_info), intent(in) :: irx, iry, irz
+
+        integer, intent(in) :: idir
+
+        real(wp), dimension(num_species) :: Xs_L, Xs_R, Xs_cell, Ys_L, Ys_R, Ys_cell
+        real(wp), dimension(num_species) :: mass_diffusivities_mixavg1, mass_diffusivities_mixavg2
+        real(wp), dimension(num_species) :: mass_diffusivities_mixavg_Cell, dXk_dxi, h_l, h_r, h_k
+        real(wp), dimension(num_species) :: Mass_Diffu_Flux
+        real(wp) :: Mass_Diffu_Energy
+        real(wp) :: MW_L, MW_R, MW_cell, Rgas_L, Rgas_R, T_L, T_R, P_L, P_R, rho_L, rho_R, rho_cell, rho_Vic
+        real(wp) :: lambda_L, lambda_R, lambda_Cell, dT_dxi, grid_spacing
+
+        integer :: x, y, z, i, n, eqn
+        integer, dimension(3) :: offsets
+
+        isc1 = irx; isc2 = iry; isc3 = irz
+
+        !$acc update device(isc1, isc2, isc3)
+
+        if (chemistry) then
+            ! Set offsets based on direction using array indexing
+            offsets = 0
+            offsets(idir) = 1
+
+            !$acc parallel loop collapse(3) gang vector default(present) copyin(offsets) &
+            !$acc private(Ys_L,Ys_R,Ys_cell,Xs_L,Xs_R,mass_diffusivities_mixavg1,mass_diffusivities_mixavg2,mass_diffusivities_mixavg_Cell,h_l,h_r,Xs_cell,h_k,dXk_dxi,Mass_Diffu_Flux)
+            do z = isc3%beg, isc3%end
+                do y = isc2%beg, isc2%end
+                    do x = isc1%beg, isc1%end
+                        ! Calculate grid spacing using direction-based indexing
+                        select case (idir)
+                        case (1)
+                            grid_spacing = x_cc(x + 1) - x_cc(x)
+                        case (2)
+                            grid_spacing = y_cc(y + 1) - y_cc(y)
+                        case (3)
+                            grid_spacing = z_cc(z + 1) - z_cc(z)
+                        end select
+
+                        ! Extract species mass fractions
+                        !$acc loop seq
+                        do i = chemxb, chemxe
+                            Ys_L(i - chemxb + 1) = q_prim_qp(i)%sf(x, y, z)
+                            Ys_R(i - chemxb + 1) = q_prim_qp(i)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
+                            Ys_cell(i - chemxb + 1) = 0.5_wp*(Ys_L(i - chemxb + 1) + Ys_R(i - chemxb + 1))
+                        end do
+
+                        ! Calculate molecular weights and mole fractions
+                        call get_mixture_molecular_weight(Ys_L, MW_L)
+                        call get_mixture_molecular_weight(Ys_R, MW_R)
+                        MW_cell = 0.5_wp*(MW_L + MW_R)
+
+                        call get_mole_fractions(MW_L, Ys_L, Xs_L)
+                        call get_mole_fractions(MW_R, Ys_R, Xs_R)
+
+                        ! Calculate gas constants and thermodynamic properties
+                        Rgas_L = gas_constant/MW_L
+                        Rgas_R = gas_constant/MW_R
+
+                        P_L = q_prim_qp(E_idx)%sf(x, y, z)
+                        P_R = q_prim_qp(E_idx)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
+
+                        rho_L = q_prim_qp(1)%sf(x, y, z)
+                        rho_R = q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
+
+                        T_L = P_L/rho_L/Rgas_L
+                        T_R = P_R/rho_R/Rgas_R
+
+                        rho_cell = 0.5_wp*(rho_L + rho_R)
+                        dT_dxi = (T_R - T_L)/grid_spacing
+
+                        ! Get transport properties
+                        call get_species_mass_diffusivities_mixavg(P_L, T_L, Ys_L, mass_diffusivities_mixavg1)
+                        call get_species_mass_diffusivities_mixavg(P_R, T_R, Ys_R, mass_diffusivities_mixavg2)
+
+                        call get_mixture_thermal_conductivity_mixavg(T_L, Ys_L, lambda_L)
+                        call get_mixture_thermal_conductivity_mixavg(T_R, Ys_R, lambda_R)
+
+                        call get_species_enthalpies_rt(T_L, h_l)
+                        call get_species_enthalpies_rt(T_R, h_r)
+
+                        ! Calculate species properties and gradients
+                        !$acc loop seq
+                        do i = chemxb, chemxe
+                            h_l(i - chemxb + 1) = h_l(i - chemxb + 1)*gas_constant*T_L/molecular_weights(i - chemxb + 1)
+                            h_r(i - chemxb + 1) = h_r(i - chemxb + 1)*gas_constant*T_R/molecular_weights(i - chemxb + 1)
+                            Xs_cell(i - chemxb + 1) = 0.5_wp*(Xs_L(i - chemxb + 1) + Xs_R(i - chemxb + 1))
+                            h_k(i - chemxb + 1) = 0.5_wp*(h_l(i - chemxb + 1) + h_r(i - chemxb + 1))
+                            dXk_dxi(i - chemxb + 1) = (Xs_R(i - chemxb + 1) - Xs_L(i - chemxb + 1))/grid_spacing
+                        end do
+
+                        ! Calculate mixture-averaged diffusivities
+                        !$acc loop seq
+                        do i = chemxb, chemxe
+                            mass_diffusivities_mixavg_Cell(i - chemxb + 1) = &
+                                (mass_diffusivities_mixavg2(i - chemxb + 1) + mass_diffusivities_mixavg1(i - chemxb + 1))/ &
+                                2.0_wp*(1.0_wp - Xs_cell(i - chemxb + 1))/(1.0_wp - Ys_cell(i - chemxb + 1))
+                        end do
+
+                        lambda_Cell = 0.5_wp*(lambda_R + lambda_L)
+
+                        ! Calculate mass diffusion fluxes
+                        rho_Vic = 0.0_wp
+                        Mass_Diffu_Energy = 0.0_wp
+
+                        !$acc loop seq
+                        do eqn = chemxb, chemxe
+                            Mass_Diffu_Flux(eqn - chemxb + 1) = rho_cell*mass_diffusivities_mixavg_Cell(eqn - chemxb + 1)* &
+                                                                molecular_weights(eqn - chemxb + 1)/MW_cell*dXk_dxi(eqn - chemxb + 1)
+                            rho_Vic = rho_Vic + Mass_Diffu_Flux(eqn - chemxb + 1)
+                            Mass_Diffu_Energy = Mass_Diffu_Energy + h_k(eqn - chemxb + 1)*Mass_Diffu_Flux(eqn - chemxb + 1)
+                        end do
+
+                        ! Apply corrections for mass conservation
+                        !$acc loop seq
+                        do eqn = chemxb, chemxe
+                            Mass_Diffu_Energy = Mass_Diffu_Energy - h_k(eqn - chemxb + 1)*Ys_cell(eqn - chemxb + 1)*rho_Vic
+                            Mass_Diffu_Flux(eqn - chemxb + 1) = Mass_Diffu_Flux(eqn - chemxb + 1) - rho_Vic*Ys_cell(eqn - chemxb + 1)
+                        end do
+
+                        ! Add thermal conduction contribution
+                        Mass_Diffu_Energy = lambda_Cell*dT_dxi + Mass_Diffu_Energy
+
+                        ! Update flux arrays
+                        flux_src_vf(E_idx)%sf(x, y, z) = flux_src_vf(E_idx)%sf(x, y, z) - Mass_Diffu_Energy
+
+                        !$acc loop seq
+                        do eqn = chemxb, chemxe
+                            flux_src_vf(eqn)%sf(x, y, z) = flux_src_vf(eqn)%sf(x, y, z) - Mass_diffu_Flux(eqn - chemxb + 1)
+                        end do
+                    end do
+                end do
+            end do
+        end if
+
+    end subroutine s_compute_chemistry_diffusion_flux
+
 end module m_chemistry

src/simulation/m_rhs.fpp

@@ -531,6 +531,12 @@ contains
                                  & idwbuff(2)%beg:idwbuff(2)%end, &
                                  & idwbuff(3)%beg:idwbuff(3)%end))
                     end do
+                    if (chem_params%diffusion .and. .not. viscous) then
+                        @:ALLOCATE(flux_src_n(i)%vf(E_idx)%sf( &
+                                 & idwbuff(1)%beg:idwbuff(1)%end, &
+                                 & idwbuff(2)%beg:idwbuff(2)%end, &
+                                 & idwbuff(3)%beg:idwbuff(3)%end))
+                    end if
                 end if
 
             else
@@ -795,9 +801,15 @@ contains
                                                                q_prim_qp%vf, &
                                                                rhs_vf)
             call nvtxEndRange
+            ! RHS for diffusion
+            if (chemistry .and. chem_params%diffusion) then
+                call nvtxStartRange("RHS-CHEM-DIFFUSION")
+                call s_compute_chemistry_diffusion_flux(id, q_prim_qp%vf, flux_src_n(id)%vf, irx, iry, irz)
+                call nvtxEndRange
+            end if
 
             ! RHS additions for viscosity
-            if (viscous .or. surface_tension) then
+            if (viscous .or. surface_tension .or. chem_params%diffusion) then
                 call nvtxStartRange("RHS-ADD-PHYSICS")
                 call s_compute_additional_physics_rhs(id, &
                                                       q_prim_qp%vf, &
@@ -1414,20 +1426,46 @@ contains
                 end do
             end if
 
-            !$acc parallel loop collapse(3) gang vector default(present)
-            do l = 0, p
-                do k = 0, n
-                    do j = 0, m
-                        !$acc loop seq
-                        do i = momxb, E_idx
-                            rhs_vf(i)%sf(j, k, l) = &
-                                rhs_vf(i)%sf(j, k, l) + 1._wp/dx(j)* &
-                                (flux_src_n(i)%sf(j - 1, k, l) &
-                                 - flux_src_n(i)%sf(j, k, l))
+            if (surface_tension .or. viscous) then
+                !$acc parallel loop collapse(3) gang vector default(present)
+                do l = 0, p
+                    do k = 0, n
+                        do j = 0, m
+                            !$acc loop seq
+                            do i = momxb, E_idx
+                                rhs_vf(i)%sf(j, k, l) = &
+                                    rhs_vf(i)%sf(j, k, l) + 1._wp/dx(j)* &
+                                    (flux_src_n(i)%sf(j - 1, k, l) &
+                                     - flux_src_n(i)%sf(j, k, l))
+                            end do
                         end do
                     end do
                 end do
-            end do
+            end if
+
+            if (chem_params%diffusion) then
+                !$acc parallel loop collapse(3) gang vector default(present)
+                do l = 0, p
+                    do k = 0, n
+                        do j = 0, m
+                            !$acc loop seq
+                            do i = chemxb, chemxe
+                                rhs_vf(i)%sf(j, k, l) = &
+                                    rhs_vf(i)%sf(j, k, l) + 1._wp/dx(j)* &
+                                    (flux_src_n(i)%sf(j - 1, k, l) &
+                                     - flux_src_n(i)%sf(j, k, l))
+                            end do
+
+                            if (.not. viscous) then
+                                rhs_vf(E_idx)%sf(j, k, l) = &
+                                    rhs_vf(E_idx)%sf(j, k, l) + 1._wp/dx(j)* &
+                                    (flux_src_n(E_idx)%sf(j - 1, k, l) &
+                                     - flux_src_n(E_idx)%sf(j, k, l))
+                            end if
+                        end do
+                    end do
+                end do
+            end if
 
         elseif (idir == 2) then ! y-direction
 
@@ -1495,20 +1533,46 @@ contains
                 end do
 
             else
-                !$acc parallel loop collapse(3) gang vector default(present)
-                do l = 0, p
-                    do k = 0, n
-                        do j = 0, m
-                            !$acc loop seq
-                            do i = momxb, E_idx
-                                rhs_vf(i)%sf(j, k, l) = &
-                                    rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* &
-                                    (flux_src_n(i)%sf(j, k - 1, l) &
-                                     - flux_src_n(i)%sf(j, k, l))
+
+                if (viscous .or. surface_tension) then
+                    !$acc parallel loop collapse(3) gang vector default(present)
+                    do l = 0, p
+                        do k = 0, n
+                            do j = 0, m
+                                !$acc loop seq
+                                do i = momxb, E_idx
+                                    rhs_vf(i)%sf(j, k, l) = &
+                                        rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* &
+                                        (flux_src_n(i)%sf(j, k - 1, l) &
+                                         - flux_src_n(i)%sf(j, k, l))
+                                end do
                             end do
                         end do
                     end do
-                end do
+                end if
+
+                if (chem_params%diffusion) then
+                    !$acc parallel loop collapse(3) gang vector default(present)
+                    do l = 0, p
+                        do k = 0, n
+                            do j = 0, m
+                                !$acc loop seq
+                                do i = chemxb, chemxe
+                                    rhs_vf(i)%sf(j, k, l) = &
+                                        rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* &
+                                        (flux_src_n(i)%sf(j, k - 1, l) &
+                                         - flux_src_n(i)%sf(j, k, l))
+                                end do
+                                if (.not. viscous) then
+                                    rhs_vf(E_idx)%sf(j, k, l) = &
+                                        rhs_vf(E_idx)%sf(j, k, l) + 1._wp/dy(k)* &
+                                        (flux_src_n(E_idx)%sf(j, k - 1, l) &
+                                         - flux_src_n(E_idx)%sf(j, k, l))
+                                end if
+                            end do
+                        end do
+                    end do
+                end if
             end if
 
             ! Applying the geometrical viscous Riemann source fluxes calculated as average
@@ -1581,20 +1645,45 @@ contains
                 end do
             end if
 
-            !$acc parallel loop collapse(3) gang vector default(present)
-            do l = 0, p
-                do k = 0, n
-                    do j = 0, m
-                        !$acc loop seq
-                        do i = momxb, E_idx
-                            rhs_vf(i)%sf(j, k, l) = &
-                                rhs_vf(i)%sf(j, k, l) + 1._wp/dz(l)* &
-                                (flux_src_n(i)%sf(j, k, l - 1) &
-                                 - flux_src_n(i)%sf(j, k, l))
+            if (viscous .or. surface_tension) then
+                !$acc parallel loop collapse(3) gang vector default(present)
+                do l = 0, p
+                    do k = 0, n
+                        do j = 0, m
+                            !$acc loop seq
+                            do i = momxb, E_idx
+                                rhs_vf(i)%sf(j, k, l) = &
+                                    rhs_vf(i)%sf(j, k, l) + 1._wp/dz(l)* &
+                                    (flux_src_n(i)%sf(j, k, l - 1) &
+                                     - flux_src_n(i)%sf(j, k, l))
+                            end do
                         end do
                     end do
                 end do
-            end do
+            end if
+
+            if (chem_params%diffusion) then
+                !$acc parallel loop collapse(3) gang vector default(present)
+                do l = 0, p
+                    do k = 0, n
+                        do j = 0, m
+                            !$acc loop seq
+                            do i = chemxb, chemxe
+                                rhs_vf(i)%sf(j, k, l) = &
+                                    rhs_vf(i)%sf(j, k, l) + 1._wp/dz(l)* &
+                                    (flux_src_n(i)%sf(j, k, l - 1) &
+                                     - flux_src_n(i)%sf(j, k, l))
+                            end do
+                            if (.not. viscous) then
+                                rhs_vf(E_idx)%sf(j, k, l) = &
+                                    rhs_vf(E_idx)%sf(j, k, l) + 1._wp/dz(l)* &
+                                    (flux_src_n(E_idx)%sf(j, k, l - 1) &
+                                     - flux_src_n(E_idx)%sf(j, k, l))
+                            end if
+                        end do
+                    end do
+                end do
+            end if
 
             if (grid_geometry == 3) then
                 !$acc parallel loop collapse(3) gang vector default(present)
@@ -1895,6 +1984,10 @@ contains
                     end do
                 end if
 
+                if (chem_params%diffusion .and. .not. viscous) then
+                    @:DEALLOCATE(flux_src_n(i)%vf(E_idx)%sf)
+                end if
+
                 if (riemann_solver == 1 .or. riemann_solver == 4) then
                     do l = adv_idx%beg + 1, adv_idx%end
                         @:DEALLOCATE(flux_src_n(i)%vf(l)%sf)