Skip to content

Commit 59d9431

Browse files
DimAdam-01DimAdam-01
authored
Lodi multi component (#838)
1 parent 4949c07 commit 59d9431

File tree

3 files changed

+177
-34
lines changed

3 files changed

+177
-34
lines changed

‎src/common/m_variables_conversion.fpp

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1467,6 +1467,8 @@ contains
14671467
real(wp) :: qv_K
14681468
real(wp), dimension(2) :: Re_K
14691469
real(wp) :: G_K
1470+
real(wp), dimension(num_species) :: Y_K
1471+
real(wp) :: T_K, mix_mol_weight, R_gas
14701472

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

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

15201522
! Computing the energy from the pressure
1521-
E_K = gamma_K*pres_K + pi_inf_K &
1522-
+ 5e-1_wp*rho_K*vel_K_sum + qv_K
1523+
1524+
if (chemistry) then
1525+
!$acc loop seq
1526+
do i = chemxb, chemxe
1527+
Y_K(i - chemxb + 1) = qK_prim_vf(j, k, l, i)
1528+
end do
1529+
!Computing the energy from the internal energy of the mixture
1530+
call get_mixture_molecular_weight(Y_k, mix_mol_weight)
1531+
R_gas = gas_constant/mix_mol_weight
1532+
T_K = pres_K/rho_K/R_gas
1533+
call get_mixture_energy_mass(T_K, Y_K, E_K)
1534+
E_K = rho_K*E_K + 5e-1_wp*rho_K*vel_K_sum
1535+
else
1536+
! Computing the energy from the pressure
1537+
E_K = gamma_K*pres_K + pi_inf_K &
1538+
+ 5e-1_wp*rho_K*vel_K_sum + qv_K
1539+
end if
15231540

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

1558+
! Species advection Flux, \rho*u*Y
1559+
if (chemistry) then
1560+
!$acc loop seq
1561+
do i = 1, num_species
1562+
FK_vf(j, k, l, i - 1 + chemxb) = vel_K(dir_idx(1))*(rho_K*Y_K(i))
1563+
end do
1564+
end if
1565+
15411566
if (riemann_solver == 1 .or. riemann_solver == 4) then
15421567
!$acc loop seq
15431568
do i = advxb, advxe

‎src/simulation/m_cbc.fpp

Lines changed: 112 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,13 @@ module m_cbc
3030

3131
use m_compute_cbc
3232

33+
use m_thermochem, only: &
34+
get_mixture_energy_mass, get_mixture_specific_heat_cv_mass, &
35+
get_mixture_specific_heat_cp_mass, gas_constant, &
36+
get_mixture_molecular_weight, get_species_enthalpies_rt, &
37+
molecular_weights, get_species_specific_heats_r, &
38+
get_mole_fractions, get_species_specific_heats_r
39+
3340
implicit none
3441

3542
private; public :: s_initialize_cbc_module, s_cbc, s_finalize_cbc_module
@@ -96,7 +103,8 @@ module m_cbc
96103
integer :: dj
97104
integer :: bcxb, bcxe, bcyb, bcye, bczb, bcze
98105
integer :: cbc_dir, cbc_loc
99-
!$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc)
106+
integer :: flux_cbc_index
107+
!$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc, flux_cbc_index)
100108
101109
!! GRCBC inputs for subsonic inflow and outflow conditions consisting of
102110
!! inflow velocities, pressure, density and void fraction as well as
@@ -124,6 +132,13 @@ contains
124132
integer :: i
125133
logical :: is_cbc
126134
135+
if (chemistry) then
136+
flux_cbc_index = sys_size
137+
else
138+
flux_cbc_index = adv_idx%end
139+
end if
140+
!$acc update device(flux_cbc_index)
141+
127142
call s_any_cbc_boundaries(is_cbc)
128143
129144
if (is_cbc .eqv. .false.) return
@@ -153,7 +168,7 @@ contains
153168
154169
@:ALLOCATE(F_rsx_vf(0:buff_size, &
155170
is2%beg:is2%end, &
156-
is3%beg:is3%end, 1:adv_idx%end))
171+
is3%beg:is3%end, 1:flux_cbc_index))
157172
158173
@:ALLOCATE(F_src_rsx_vf(0:buff_size, &
159174
is2%beg:is2%end, &
@@ -163,7 +178,7 @@ contains
163178
164179
@:ALLOCATE(flux_rsx_vf_l(-1:buff_size, &
165180
is2%beg:is2%end, &
166-
is3%beg:is3%end, 1:adv_idx%end))
181+
is3%beg:is3%end, 1:flux_cbc_index))
167182
168183
@:ALLOCATE(flux_src_rsx_vf_l(-1:buff_size, &
169184
is2%beg:is2%end, &
@@ -196,7 +211,7 @@ contains
196211
197212
@:ALLOCATE(F_rsy_vf(0:buff_size, &
198213
is2%beg:is2%end, &
199-
is3%beg:is3%end, 1:adv_idx%end))
214+
is3%beg:is3%end, 1:flux_cbc_index))
200215
201216
@:ALLOCATE(F_src_rsy_vf(0:buff_size, &
202217
is2%beg:is2%end, &
@@ -206,7 +221,7 @@ contains
206221
207222
@:ALLOCATE(flux_rsy_vf_l(-1:buff_size, &
208223
is2%beg:is2%end, &
209-
is3%beg:is3%end, 1:adv_idx%end))
224+
is3%beg:is3%end, 1:flux_cbc_index))
210225
211226
@:ALLOCATE(flux_src_rsy_vf_l(-1:buff_size, &
212227
is2%beg:is2%end, &
@@ -241,7 +256,7 @@ contains
241256
242257
@:ALLOCATE(F_rsz_vf(0:buff_size, &
243258
is2%beg:is2%end, &
244-
is3%beg:is3%end, 1:adv_idx%end))
259+
is3%beg:is3%end, 1:flux_cbc_index))
245260
246261
@:ALLOCATE(F_src_rsz_vf(0:buff_size, &
247262
is2%beg:is2%end, &
@@ -251,7 +266,7 @@ contains
251266
252267
@:ALLOCATE(flux_rsz_vf_l(-1:buff_size, &
253268
is2%beg:is2%end, &
254-
is3%beg:is3%end, 1:adv_idx%end))
269+
is3%beg:is3%end, 1:flux_cbc_index))
255270
256271
@:ALLOCATE(flux_src_rsz_vf_l(-1:buff_size, &
257272
is2%beg:is2%end, &
@@ -651,6 +666,9 @@ contains
651666
real(wp) :: qv !< Cell averaged fluid reference energy
652667
real(wp) :: c
653668
real(wp) :: Ma
669+
real(wp) :: T, sum_Enthalpies
670+
real(wp) :: Cv, Cp, e_mix, Mw, R_gas
671+
real(wp), dimension(num_species) :: Ys, h_k, dYs_dt, dYs_ds, Xs, Gamma_i, Cp_i
654672
655673
real(wp) :: vel_K_sum, vel_dv_dt_sum
656674
@@ -684,7 +702,7 @@ contains
684702
is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg)
685703
686704
!$acc parallel loop collapse(3) gang vector default(present)
687-
do i = 1, advxe
705+
do i = 1, flux_cbc_index
688706
do r = is3%beg, is3%end
689707
do k = is2%beg, is2%end
690708
flux_rs${XYZ}$_vf_l(0, k, r, i) = F_rs${XYZ}$_vf(0, k, r, i) &
@@ -715,7 +733,7 @@ contains
715733
is1, is2, is3, idwbuff(2)%beg, idwbuff(3)%beg)
716734
717735
!$acc parallel loop collapse(4) gang vector default(present)
718-
do i = 1, advxe
736+
do i = 1, flux_cbc_index
719737
do j = 0, 1
720738
do r = is3%beg, is3%end
721739
do k = is2%beg, is2%end
@@ -757,7 +775,7 @@ contains
757775
end if
758776
759777
! FD2 or FD4 of RHS at j = 0
760-
!$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)
778+
!$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)
761779
do r = is3%beg, is3%end
762780
do k = is2%beg, is2%end
763781
@@ -796,7 +814,33 @@ contains
796814
mf(i) = alpha_rho(i)/rho
797815
end do
798816
799-
E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_K_sum
817+
if (chemistry) then
818+
!$acc loop seq
819+
do i = chemxb, chemxe
820+
Ys(i - chemxb + 1) = q_prim_rs${XYZ}$_vf(0, k, r, i)
821+
end do
822+
823+
call get_mixture_molecular_weight(Ys, Mw)
824+
R_gas = gas_constant/Mw
825+
T = pres/rho/R_gas
826+
call get_mixture_specific_heat_cp_mass(T, Ys, Cp)
827+
call get_mixture_energy_mass(T, Ys, e_mix)
828+
E = rho*e_mix + 5e-1_wp*rho*vel_K_sum
829+
if (chem_params%gamma_method == 1) then
830+
!> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
831+
call get_mole_fractions(Mw, Ys, Xs)
832+
call get_species_specific_heats_r(T, Cp_i)
833+
Gamma_i = Cp_i/(Cp_i - 1.0_wp)
834+
gamma = sum(Xs(:)/(Gamma_i(:) - 1.0_wp))
835+
else if (chem_params%gamma_method == 2) then
836+
!> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
837+
call get_mixture_specific_heat_cv_mass(T, Ys, Cv)
838+
gamma = 1.0_wp/(Cp/Cv - 1.0_wp)
839+
end if
840+
else
841+
E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_K_sum
842+
end if
843+
800844
H = (E + pres)/rho
801845
802846
! Compute mixture sound speed
@@ -820,6 +864,13 @@ contains
820864
dadv_ds(i) = 0._wp
821865
end do
822866
867+
if (chemistry) then
868+
!$acc loop seq
869+
do i = 1, num_species
870+
dYs_ds(i) = 0._wp
871+
end do
872+
end if
873+
823874
!$acc loop seq
824875
do j = 0, buff_size
825876
@@ -845,6 +896,15 @@ contains
845896
fd_coef_${XYZ}$ (j, cbc_loc) + &
846897
dadv_ds(i)
847898
end do
899+
900+
if (chemistry) then
901+
!$acc loop seq
902+
do i = 1, num_species
903+
dYs_ds(i) = q_prim_rs${XYZ}$_vf(j, k, r, chemxb - 1 + i)* &
904+
fd_coef_${XYZ}$ (j, cbc_loc) + &
905+
dYs_ds(i)
906+
end do
907+
end if
848908
end do
849909
850910
! First-Order Temporal Derivatives of Primitive Variables
@@ -859,7 +919,7 @@ contains
859919
call s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
860920
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_BUFFER) .or. &
861921
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_BUFFER)) then
862-
call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
922+
call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
863923
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_INFLOW) .or. &
864924
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_INFLOW)) then
865925
call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
@@ -883,7 +943,7 @@ contains
883943
end if
884944
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_OUTFLOW) .or. &
885945
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_OUTFLOW)) then
886-
call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
946+
call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
887947
! Add GRCBC for Subsonic Outflow (Pressure)
888948
if (bc_${XYZ}$%grcbc_out) then
889949
L(advxe) = c*(1._wp - Ma)*(pres - pres_out(${CBC_DIR}$))/Del_out(${CBC_DIR}$)
@@ -904,7 +964,7 @@ contains
904964
call s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
905965
else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_SUP_OUTFLOW) .or. &
906966
(cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_SUP_OUTFLOW)) then
907-
call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
967+
call s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds, dYs_ds)
908968
end if
909969
910970
! Be careful about the cylindrical coordinate!
@@ -935,6 +995,13 @@ contains
935995
vel_dv_dt_sum = vel_dv_dt_sum + vel(i)*dvel_dt(i)
936996
end do
937997
998+
if (chemistry) then
999+
!$acc loop seq
1000+
do i = 1, num_species
1001+
dYs_dt(i) = -1._wp*L(chemxb + i - 1)
1002+
end do
1003+
end if
1004+
9381005
! The treatment of void fraction source is unclear
9391006
if (cyl_coord .and. cbc_dir == 2 .and. cbc_loc == 1) then
9401007
!$acc loop seq
@@ -978,13 +1045,31 @@ contains
9781045
+ rho*dvel_dt(i - contxe))
9791046
end do
9801047
981-
flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) &
982-
+ ds(0)*(pres*dgamma_dt &
983-
+ gamma*dpres_dt &
984-
+ dpi_inf_dt &
985-
+ dqv_dt &
986-
+ rho*vel_dv_dt_sum &
987-
+ 5e-1_wp*drho_dt*vel_K_sum)
1048+
if (chemistry) then
1049+
! Evolution of LODI equation of energy for real gases adjusted to perfect gas, doi:10.1006/jcph.2002.6990
1050+
call get_species_enthalpies_rt(T, h_k)
1051+
sum_Enthalpies = 0._wp
1052+
!$acc loop seq
1053+
do i = 1, num_species
1054+
h_k(i) = h_k(i)*gas_constant/molecular_weights(i)*T
1055+
sum_Enthalpies = sum_Enthalpies + (rho*h_k(i) - pres*Mw/molecular_weights(i)*Cp/R_gas)*dYs_dt(i)
1056+
end do
1057+
flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) &
1058+
+ ds(0)*((E/rho + pres/rho)*drho_dt + rho*vel_dv_dt_sum + Cp*T*L(2)/(c*c) + sum_Enthalpies)
1059+
!$acc loop seq
1060+
do i = 1, num_species
1061+
flux_rs${XYZ}$_vf_l(-1, k, r, i - 1 + chemxb) = flux_rs${XYZ}$_vf_l(0, k, r, chemxb + i - 1) &
1062+
+ ds(0)*(drho_dt*Ys(i) + rho*dYs_dt(i))
1063+
end do
1064+
else
1065+
flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) &
1066+
+ ds(0)*(pres*dgamma_dt &
1067+
+ gamma*dpres_dt &
1068+
+ dpi_inf_dt &
1069+
+ dqv_dt &
1070+
+ rho*vel_dv_dt_sum &
1071+
+ 5e-1_wp*drho_dt*vel_K_sum)
1072+
end if
9881073
9891074
if (riemann_solver == 1) then
9901075
!$acc loop seq
@@ -1106,7 +1191,7 @@ contains
11061191
end do
11071192
11081193
!$acc parallel loop collapse(4) gang vector default(present)
1109-
do i = 1, advxe
1194+
do i = 1, flux_cbc_index
11101195
do r = is3%beg, is3%end
11111196
do k = is2%beg, is2%end
11121197
do j = -1, buff_size
@@ -1182,7 +1267,7 @@ contains
11821267
end do
11831268
11841269
!$acc parallel loop collapse(4) gang vector default(present)
1185-
do i = 1, advxe
1270+
do i = 1, flux_cbc_index
11861271
do r = is3%beg, is3%end
11871272
do k = is2%beg, is2%end
11881273
do j = -1, buff_size
@@ -1258,7 +1343,7 @@ contains
12581343
end do
12591344
12601345
!$acc parallel loop collapse(4) gang vector default(present)
1261-
do i = 1, advxe
1346+
do i = 1, flux_cbc_index
12621347
do r = is3%beg, is3%end
12631348
do k = is2%beg, is2%end
12641349
do j = -1, buff_size
@@ -1339,7 +1424,7 @@ contains
13391424
if (cbc_dir == 1) then
13401425
13411426
!$acc parallel loop collapse(4) gang vector default(present)
1342-
do i = 1, advxe
1427+
do i = 1, flux_cbc_index
13431428
do r = is3%beg, is3%end
13441429
do k = is2%beg, is2%end
13451430
do j = -1, buff_size
@@ -1390,7 +1475,7 @@ contains
13901475
elseif (cbc_dir == 2) then
13911476
13921477
!$acc parallel loop collapse(4) gang vector default(present)
1393-
do i = 1, advxe
1478+
do i = 1, flux_cbc_index
13941479
do r = is3%beg, is3%end
13951480
do k = is2%beg, is2%end
13961481
do j = -1, buff_size
@@ -1443,7 +1528,7 @@ contains
14431528
else
14441529
14451530
!$acc parallel loop collapse(4) gang vector default(present)
1446-
do i = 1, advxe
1531+
do i = 1, flux_cbc_index
14471532
do r = is3%beg, is3%end
14481533
do k = is2%beg, is2%end
14491534
do j = -1, buff_size

0 commit comments

Comments
 (0)