Skip to content

Lodi multi component #838

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 15 commits into from
May 22, 2025
Merged

Lodi multi component #838

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
bcaaa35
Chemistry Bug Fixes
DimAdam-01 May 6, 2025
4833d7b
Chemistry Bug Fixes
DimAdam-01 May 6, 2025
c8e7f9f
Format
DimAdam-01 May 6, 2025
652d0fa
Testing-Changes-In-CodeGenerator-For-Pyrometheus
DimAdam-01 May 8, 2025
e63683a
Merge branch 'MFlowCode:master' into my_new_branch
DimAdam-01 May 8, 2025
cc20ee9
Extend Frontier Run Time
DimAdam-01 May 8, 2025
b720ecc
Debugggg
DimAdam-01 May 17, 2025
1182c2c
Small Changes
DimAdam-01 May 17, 2025
24a324c
Lodi for Multicomponent
DimAdam-01 May 19, 2025
b02547e
Merge
DimAdam-01 May 19, 2025
308d919
Merge_2
DimAdam-01 May 19, 2025
043b8a7
Upstream
DimAdam-01 May 21, 2025
f847598
Merge branch 'master' into Lodi_MultiComponent
DimAdam-01 May 21, 2025
a5a0198
Images
DimAdam-01 May 21, 2025
55c9549
Merge remote-tracking branch 'origin/Lodi_MultiComponent' into Lodi_M…
DimAdam-01 May 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file removed examples/1D_reactive_shocktube/initial.png
View file
Open in desktop
Binary file not shown.
Binary file removed examples/1D_reactive_shocktube/result.png
View file
Open in desktop
Binary file not shown.
31 changes: 28 additions & 3 deletions src/common/m_variables_conversion.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1467,6 +1467,8 @@ contains
real(wp) :: qv_K
real(wp), dimension(2) :: Re_K
real(wp) :: G_K
real(wp), dimension(num_species) :: Y_K
real(wp) :: T_K, mix_mol_weight, R_gas

integer :: i, j, k, l !< Generic loop iterators

Expand All @@ -1479,7 +1481,7 @@ contains
! Computing the flux variables from the primitive variables, without
! accounting for the contribution of either viscosity or capillarity
#ifdef MFC_SIMULATION
!$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_K, vel_K, alpha_K, Re_K)
!$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_K, vel_K, alpha_K, Re_K, Y_K)
do l = is3b, is3e
do k = is2b, is2e
do j = is1b, is1e
Expand Down Expand Up @@ -1518,8 +1520,23 @@ contains
end if

! Computing the energy from the pressure
E_K = gamma_K*pres_K + pi_inf_K &
+ 5e-1_wp*rho_K*vel_K_sum + qv_K

if (chemistry) then
!$acc loop seq
do i = chemxb, chemxe
Y_K(i - chemxb + 1) = qK_prim_vf(j, k, l, i)
end do
!Computing the energy from the internal energy of the mixture
call get_mixture_molecular_weight(Y_k, mix_mol_weight)
R_gas = gas_constant/mix_mol_weight
T_K = pres_K/rho_K/R_gas
call get_mixture_energy_mass(T_K, Y_K, E_K)
E_K = rho_K*E_K + 5e-1_wp*rho_K*vel_K_sum
else
! Computing the energy from the pressure
E_K = gamma_K*pres_K + pi_inf_K &
+ 5e-1_wp*rho_K*vel_K_sum + qv_K
end if
Comment on lines +1524 to +1539
Copy link
Collaborator

@wilfonba wilfonba May 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you make sure this doesn't have any performance impact on non-chemistry cases? This variables conversion subroutine can be really touchy on AMD hardware...


! mass flux, this should be \alpha_i \rho_i u_i
!$acc loop seq
Expand All @@ -1538,6 +1555,14 @@ contains
! energy flux, u(E+p)
FK_vf(j, k, l, E_idx) = vel_K(dir_idx(1))*(E_K + pres_K)

! Species advection Flux, \rho*u*Y
if (chemistry) then
!$acc loop seq
do i = 1, num_species
FK_vf(j, k, l, i - 1 + chemxb) = vel_K(dir_idx(1))*(rho_K*Y_K(i))
end do
end if

if (riemann_solver == 1 .or. riemann_solver == 4) then
!$acc loop seq
do i = advxb, advxe
Expand Down
139 changes: 112 additions & 27 deletions src/simulation/m_cbc.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,13 @@ module m_cbc

use m_compute_cbc

use m_thermochem, only: &
get_mixture_energy_mass, get_mixture_specific_heat_cv_mass, &
get_mixture_specific_heat_cp_mass, gas_constant, &
get_mixture_molecular_weight, get_species_enthalpies_rt, &
molecular_weights, get_species_specific_heats_r, &
get_mole_fractions, get_species_specific_heats_r

implicit none

private; public :: s_initialize_cbc_module, s_cbc, s_finalize_cbc_module
Expand Down Expand Up @@ -96,7 +103,8 @@ module m_cbc
integer :: dj
integer :: bcxb, bcxe, bcyb, bcye, bczb, bcze
integer :: cbc_dir, cbc_loc
!$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc)
integer :: flux_cbc_index
!$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc, flux_cbc_index)

!! GRCBC inputs for subsonic inflow and outflow conditions consisting of
!! inflow velocities, pressure, density and void fraction as well as
Expand Down Expand Up @@ -124,6 +132,13 @@ contains
integer :: i
logical :: is_cbc

if (chemistry) then
flux_cbc_index = sys_size
else
flux_cbc_index = adv_idx%end
end if
!$acc update device(flux_cbc_index)

call s_any_cbc_boundaries(is_cbc)

if (is_cbc .eqv. .false.) return
Expand Down Expand Up @@ -153,7 +168,7 @@ contains

@:ALLOCATE(F_rsx_vf(0:buff_size, &
is2%beg:is2%end, &
is3%beg:is3%end, 1:adv_idx%end))
is3%beg:is3%end, 1:flux_cbc_index))

@:ALLOCATE(F_src_rsx_vf(0:buff_size, &
is2%beg:is2%end, &
Expand All @@ -163,7 +178,7 @@ contains

@:ALLOCATE(flux_rsx_vf_l(-1:buff_size, &
is2%beg:is2%end, &
is3%beg:is3%end, 1:adv_idx%end))
is3%beg:is3%end, 1:flux_cbc_index))

@:ALLOCATE(flux_src_rsx_vf_l(-1:buff_size, &
is2%beg:is2%end, &
Expand Down Expand Up @@ -196,7 +211,7 @@ contains

@:ALLOCATE(F_rsy_vf(0:buff_size, &
is2%beg:is2%end, &
is3%beg:is3%end, 1:adv_idx%end))
is3%beg:is3%end, 1:flux_cbc_index))

@:ALLOCATE(F_src_rsy_vf(0:buff_size, &
is2%beg:is2%end, &
Expand All @@ -206,7 +221,7 @@ contains

@:ALLOCATE(flux_rsy_vf_l(-1:buff_size, &
is2%beg:is2%end, &
is3%beg:is3%end, 1:adv_idx%end))
is3%beg:is3%end, 1:flux_cbc_index))

@:ALLOCATE(flux_src_rsy_vf_l(-1:buff_size, &
is2%beg:is2%end, &
Expand Down Expand Up @@ -241,7 +256,7 @@ contains

@:ALLOCATE(F_rsz_vf(0:buff_size, &
is2%beg:is2%end, &
is3%beg:is3%end, 1:adv_idx%end))
is3%beg:is3%end, 1:flux_cbc_index))

@:ALLOCATE(F_src_rsz_vf(0:buff_size, &
is2%beg:is2%end, &
Expand All @@ -251,7 +266,7 @@ contains

@:ALLOCATE(flux_rsz_vf_l(-1:buff_size, &
is2%beg:is2%end, &
is3%beg:is3%end, 1:adv_idx%end))
is3%beg:is3%end, 1:flux_cbc_index))

@:ALLOCATE(flux_src_rsz_vf_l(-1:buff_size, &
is2%beg:is2%end, &
Expand Down Expand Up @@ -651,6 +666,9 @@ contains
real(wp) :: qv !< Cell averaged fluid reference energy
real(wp) :: c
real(wp) :: Ma
real(wp) :: T, sum_Enthalpies
real(wp) :: Cv, Cp, e_mix, Mw, R_gas
real(wp), dimension(num_species) :: Ys, h_k, dYs_dt, dYs_ds, Xs, Gamma_i, Cp_i

real(wp) :: vel_K_sum, vel_dv_dt_sum

Expand Down Expand Up @@ -684,7 +702,7 @@ contains
is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg)

!$acc parallel loop collapse(3) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do r = is3%beg, is3%end
do k = is2%beg, is2%end
flux_rs${XYZ}$_vf_l(0, k, r, i) = F_rs${XYZ}$_vf(0, k, r, i) &
Expand Down Expand Up @@ -715,7 +733,7 @@ contains
is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg)

!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do j = 0, 1
do r = is3%beg, is3%end
do k = is2%beg, is2%end
Expand Down Expand Up @@ -757,7 +775,7 @@ contains
end if

! FD2 or FD4 of RHS at j = 0
!$acc parallel loop collapse(2) gang vector default(present) private(alpha_rho, vel, adv, mf, dvel_ds, dadv_ds, Re_cbc, dalpha_rho_ds,dvel_dt, dadv_dt, dalpha_rho_dt,L, lambda)
!$acc parallel loop collapse(2) gang vector default(present) private(alpha_rho, vel, adv, mf, dvel_ds, dadv_ds, Re_cbc, dalpha_rho_ds,dvel_dt, dadv_dt, dalpha_rho_dt,L, lambda,Ys,dYs_dt,dYs_ds,h_k,Cp_i,Gamma_i,Xs)
do r = is3%beg, is3%end
do k = is2%beg, is2%end

Expand Down Expand Up @@ -796,7 +814,33 @@ contains
mf(i) = alpha_rho(i)/rho
end do

E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_K_sum
if (chemistry) then
!$acc loop seq
do i = chemxb, chemxe
Ys(i - chemxb + 1) = q_prim_rs${XYZ}$_vf(0, k, r, i)
end do

call get_mixture_molecular_weight(Ys, Mw)
R_gas = gas_constant/Mw
T = pres/rho/R_gas
call get_mixture_specific_heat_cp_mass(T, Ys, Cp)
call get_mixture_energy_mass(T, Ys, e_mix)
E = rho*e_mix + 5e-1_wp*rho*vel_K_sum
if (chem_params%gamma_method == 1) then
!> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
call get_mole_fractions(Mw, Ys, Xs)
call get_species_specific_heats_r(T, Cp_i)
Gamma_i = Cp_i/(Cp_i - 1.0_wp)
gamma = sum(Xs(:)/(Gamma_i(:) - 1.0_wp))
else if (chem_params%gamma_method == 2) then
!> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
call get_mixture_specific_heat_cv_mass(T, Ys, Cv)
gamma = 1.0_wp/(Cp/Cv - 1.0_wp)
end if
else
E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_K_sum
end if

H = (E + pres)/rho

! Compute mixture sound speed
Expand All @@ -820,6 +864,13 @@ contains
dadv_ds(i) = 0._wp
end do

if (chemistry) then
!$acc loop seq
do i = 1, num_species
dYs_ds(i) = 0._wp
end do
end if

!$acc loop seq
do j = 0, buff_size

Expand All @@ -845,6 +896,15 @@ contains
fd_coef_${XYZ}$ (j, cbc_loc) + &
dadv_ds(i)
end do

if (chemistry) then
!$acc loop seq
do i = 1, num_species
dYs_ds(i) = q_prim_rs${XYZ}$_vf(j, k, r, chemxb - 1 + i)* &
fd_coef_${XYZ}$ (j, cbc_loc) + &
dYs_ds(i)
end do
end if
end do

! First-Order Temporal Derivatives of Primitive Variables
Expand All @@ -859,7 +919,7 @@ contains
call s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_BUFFER) .or. &
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_BUFFER)) then
call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_INFLOW) .or. &
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_INFLOW)) then
call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
Expand All @@ -883,7 +943,7 @@ contains
end if
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_OUTFLOW) .or. &
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_OUTFLOW)) then
call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
! Add GRCBC for Subsonic Outflow (Pressure)
if (bc_${XYZ}$%grcbc_out) then
L(advxe) = c*(1._wp - Ma)*(pres - pres_out(${CBC_DIR}$))/Del_out(${CBC_DIR}$)
Expand All @@ -904,7 +964,7 @@ contains
call s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_SUP_OUTFLOW) .or. &
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_SUP_OUTFLOW)) then
call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
end if

! Be careful about the cylindrical coordinate!
Expand Down Expand Up @@ -935,6 +995,13 @@ contains
vel_dv_dt_sum = vel_dv_dt_sum + vel(i)*dvel_dt(i)
end do

if (chemistry) then
!$acc loop seq
do i = 1, num_species
dYs_dt(i) = -1._wp*L(chemxb + i - 1)
end do
end if

! The treatment of void fraction source is unclear
if (cyl_coord .and. cbc_dir == 2 .and. cbc_loc == 1) then
!$acc loop seq
Expand Down Expand Up @@ -978,13 +1045,31 @@ contains
+ rho*dvel_dt(i - contxe))
end do

flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) &
+ ds(0)*(pres*dgamma_dt &
+ gamma*dpres_dt &
+ dpi_inf_dt &
+ dqv_dt &
+ rho*vel_dv_dt_sum &
+ 5e-1_wp*drho_dt*vel_K_sum)
if (chemistry) then
! Evolution of LODI equation of energy for real gases adjusted to perfect gas, doi:10.1006/jcph.2002.6990
call get_species_enthalpies_rt(T, h_k)
sum_Enthalpies = 0._wp
!$acc loop seq
do i = 1, num_species
h_k(i) = h_k(i)*gas_constant/molecular_weights(i)*T
sum_Enthalpies = sum_Enthalpies + (rho*h_k(i) - pres*Mw/molecular_weights(i)*Cp/R_gas)*dYs_dt(i)
end do
flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) &
+ ds(0)*((E/rho + pres/rho)*drho_dt + rho*vel_dv_dt_sum + Cp*T*L(2)/(c*c) + sum_Enthalpies)
!$acc loop seq
do i = 1, num_species
flux_rs${XYZ}$_vf_l(-1, k, r, i - 1 + chemxb) = flux_rs${XYZ}$_vf_l(0, k, r, chemxb + i - 1) &
+ ds(0)*(drho_dt*Ys(i) + rho*dYs_dt(i))
end do
else
flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) &
+ ds(0)*(pres*dgamma_dt &
+ gamma*dpres_dt &
+ dpi_inf_dt &
+ dqv_dt &
+ rho*vel_dv_dt_sum &
+ 5e-1_wp*drho_dt*vel_K_sum)
end if

if (riemann_solver == 1) then
!$acc loop seq
Expand Down Expand Up @@ -1106,7 +1191,7 @@ contains
end do

!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do r = is3%beg, is3%end
do k = is2%beg, is2%end
do j = -1, buff_size
Expand Down Expand Up @@ -1182,7 +1267,7 @@ contains
end do

!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do r = is3%beg, is3%end
do k = is2%beg, is2%end
do j = -1, buff_size
Expand Down Expand Up @@ -1258,7 +1343,7 @@ contains
end do

!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do r = is3%beg, is3%end
do k = is2%beg, is2%end
do j = -1, buff_size
Expand Down Expand Up @@ -1339,7 +1424,7 @@ contains
if (cbc_dir == 1) then

!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do r = is3%beg, is3%end
do k = is2%beg, is2%end
do j = -1, buff_size
Expand Down Expand Up @@ -1390,7 +1475,7 @@ contains
elseif (cbc_dir == 2) then

!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do r = is3%beg, is3%end
do k = is2%beg, is2%end
do j = -1, buff_size
Expand Down Expand Up @@ -1443,7 +1528,7 @@ contains
else

!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, advxe
do i = 1, flux_cbc_index
do r = is3%beg, is3%end
do k = is2%beg, is2%end
do j = -1, buff_size
Expand Down
Loading
Loading