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 all 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
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 @@
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 @@
! 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 @@
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)

Check warning on line 1527 in src/common/m_variables_conversion.fpp

View check run for this annotation

Codecov / codecov/patch

src/common/m_variables_conversion.fpp#L1526-L1527 

Added lines #L1526 - L1527 were not covered by tests
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

Check warning on line 1534 in src/common/m_variables_conversion.fpp

View check run for this annotation

Codecov / codecov/patch

src/common/m_variables_conversion.fpp#L1530-L1534 

Added lines #L1530 - L1534 were not covered by tests
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 @@
! 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))

Check warning on line 1562 in src/common/m_variables_conversion.fpp

View check run for this annotation

Codecov / codecov/patch

src/common/m_variables_conversion.fpp#L1561-L1562 

Added lines #L1561 - L1562 were not covered by tests
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 @@

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 @@
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 @@
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 @@

@: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 @@

@: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 @@

@: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 @@

@: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 @@

@: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 @@

@: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 @@
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 @@
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 @@
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 @@
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 @@
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)

Check warning on line 820 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L819-L820 

Added lines #L819 - L820 were not covered by tests
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

Check warning on line 828 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L823-L828 

Added lines #L823 - L828 were not covered by tests
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

Check warning on line 835 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L831-L835 

Added lines #L831 - L835 were not covered by tests
!> 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)

Check warning on line 838 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L837-L838 

Added lines #L837 - L838 were not covered by tests
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 @@
dadv_ds(i) = 0._wp
end do

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

Check warning on line 870 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L869-L870 

Added lines #L869 - L870 were not covered by tests
end do
end if

!$acc loop seq
do j = 0, buff_size

Expand All @@ -845,6 +896,15 @@
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)* &

Check warning on line 903 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L902-L903 

Added lines #L902 - L903 were not covered by tests
fd_coef_${XYZ}$ (j, cbc_loc) + &
dYs_ds(i)

Check warning on line 905 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L905 

Added line #L905 was not covered by tests
end do
end if
end do

! First-Order Temporal Derivatives of Primitive Variables
Expand All @@ -859,7 +919,7 @@
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 @@
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 @@
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 @@
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)

Check warning on line 1001 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L1000-L1001 

Added lines #L1000 - L1001 were not covered by tests
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 @@
+ 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)

Check warning on line 1050 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L1050 

Added line #L1050 was not covered by tests
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)

Check warning on line 1055 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L1053-L1055 

Added lines #L1053 - L1055 were not covered by tests
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)

Check warning on line 1058 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L1057-L1058 

Added lines #L1057 - L1058 were not covered by tests
!$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))

Check warning on line 1062 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L1060-L1062 

Added lines #L1060 - L1062 were not covered by tests
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 &

Check warning on line 1066 in src/simulation/m_cbc.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_cbc.fpp#L1066 

Added line #L1066 was not covered by tests
+ 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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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