Skip to content

Commit 96bb43d

Browse files
henryleberreDimAdam-01
henryleberre
and
DimAdam-01
committed
CheMFC fixes for HLL+Advection & Pressure/Energy
Co-authored-by: Dimitrios Adam <[email protected]>
1 parent b0d9461 commit 96bb43d

File tree

12 files changed

+133
-129
lines changed

12 files changed

+133
-129
lines changed

‎src/common/m_thermochem.fpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,8 @@ module m_thermochem
77
num_species, species_names, gas_constant, mol_weights, &
88
get_temperature, get_net_production_rates, get_pressure, &
99
get_mixture_molecular_weight, get_mixture_energy_mass, &
10-
get_temperature
10+
get_mixture_specific_heat_cv_mass, get_mixture_specific_heat_cp_mass, &
11+
get_species_enthalpies_rt
1112
#:endif
1213

1314
implicit none

‎src/common/m_variables_conversion.fpp‎

Lines changed: 18 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -127,20 +127,21 @@ contains
127127
!! @param pres Pressure to calculate
128128
!! @param stress Shear Stress
129129
!! @param mom Momentum
130-
subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, stress, mom, G)
130+
subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G)
131+
131132
!$acc routine seq
132133

133134
real(kind(0d0)), intent(in) :: energy, alf
134135
real(kind(0d0)), intent(in) :: dyn_p
135136
real(kind(0d0)), intent(in) :: pi_inf, gamma, rho, qv
136-
real(kind(0d0)), intent(out) :: pres
137+
real(kind(0d0)), intent(out) :: pres, T
137138
real(kind(0d0)), intent(in), optional :: stress, mom, G
138139

139140
! Chemistry
140141
integer :: i
141142
real(kind(0d0)), dimension(1:num_species), intent(in) :: rhoYks
142143
real(kind(0d0)) :: E_e
143-
real(kind(0d0)) :: T
144+
real(kind(0d0)) :: e_Per_Kg, Pdyn_Per_Kg
144145
real(kind(0d0)), dimension(1:num_species) :: Y_rs
145146

146147
integer :: s !< Generic loop iterator
@@ -189,12 +190,11 @@ contains
189190
Y_rs(i) = rhoYks(i)/rho
190191
end do
191192

192-
if (sum(Y_rs) > 1d-16) then
193-
call get_temperature(.true., energy - dyn_p, 1200d0, Y_rs, T)
194-
call get_pressure(rho, T, Y_rs, pres)
195-
else
196-
pres = 0d0
197-
end if
193+
e_Per_Kg = energy/rho
194+
Pdyn_Per_Kg = dyn_p/rho
195+
196+
call get_temperature(e_Per_Kg - Pdyn_Per_Kg, 1200d0, Y_rs, .true., T)
197+
call get_pressure(rho, T, Y_rs, pres)
198198

199199
#:endif
200200

@@ -886,7 +886,7 @@ contains
886886

887887
real(kind(0d0)) :: G_K
888888

889-
real(kind(0d0)) :: pres, Yksum
889+
real(kind(0d0)) :: pres, Yksum, T
890890

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

@@ -1004,9 +1004,12 @@ contains
10041004

10051005
call s_compute_pressure(qK_cons_vf(E_idx)%sf(j, k, l), &
10061006
qK_cons_vf(alf_idx)%sf(j, k, l), &
1007-
dyn_pres_K, pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres)
1007+
dyn_pres_K, pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres, T)
10081008

10091009
qK_prim_vf(E_idx)%sf(j, k, l) = pres
1010+
if (chemistry) then
1011+
qK_prim_vf(tempxb)%sf(j, k, l) = T
1012+
end if
10101013

10111014
if (bubbles) then
10121015
!$acc loop seq
@@ -1117,7 +1120,7 @@ contains
11171120
integer :: spec
11181121

11191122
real(kind(0d0)), dimension(num_species) :: Ys
1120-
real(kind(0d0)) :: temperature, e_mix, mix_mol_weight, T
1123+
real(kind(0d0)) :: e_mix, mix_mol_weight, T
11211124

11221125
#ifndef MFC_SIMULATION
11231126
! Converting the primitive variables to the conservative variables
@@ -1158,13 +1161,13 @@ contains
11581161
end do
11591162

11601163
call get_mixture_molecular_weight(Ys, mix_mol_weight)
1161-
T = q_prim_vf(E_idx)%sf(j, k, l)*mix_mol_weight/(gas_constant*rho)
1164+
T = q_prim_vf(tempxb)%sf(j, k, l)
11621165
call get_mixture_energy_mass(T, Ys, e_mix)
11631166

11641167
q_cons_vf(E_idx)%sf(j, k, l) = &
1165-
dyn_pres + e_mix
1168+
dyn_pres + rho*e_mix
11661169

1167-
q_cons_vf(tempxb)%sf(j, k, l) = T
1170+
q_cons_vf(tempxb)%sf(j, k, l) = q_prim_vf(tempxb)%sf(j, k, l)
11681171
#:else
11691172
! Computing the energy from the pressure
11701173
if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then

‎src/post_process/m_start_up.f90‎

Lines changed: 21 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -296,19 +296,32 @@ subroutine s_save_data(t_step, varname, pres, c, H)
296296
! ----------------------------------------------------------------------
297297

298298
! Adding the species' concentrations to the formatted database file ----
299-
do i = 1, num_species
300-
if (chem_wrt_Y(i) .or. prim_vars_wrt) then
301-
q_sf = q_prim_vf(chemxb + i - 1)%sf(-offset_x%beg:m + offset_x%end, &
302-
-offset_y%beg:n + offset_y%end, &
303-
-offset_z%beg:p + offset_z%end)
299+
if (chemistry) then
300+
do i = 1, num_species
301+
if (chem_wrt_Y(i) .or. prim_vars_wrt) then
302+
q_sf = q_prim_vf(chemxb + i - 1)%sf(-offset_x%beg:m + offset_x%end, &
303+
-offset_y%beg:n + offset_y%end, &
304+
-offset_z%beg:p + offset_z%end)
305+
306+
write (varname, '(A,A)') 'Y_', trim(species_names(i))
307+
call s_write_variable_to_formatted_database_file(varname, t_step)
308+
309+
varname(:) = ' '
310+
311+
end if
312+
end do
313+
314+
if (chem_wrt_T) then
315+
q_sf = q_prim_vf(tempxb)%sf(-offset_x%beg:m + offset_x%end, &
316+
-offset_y%beg:n + offset_y%end, &
317+
-offset_z%beg:p + offset_z%end)
304318

305-
write (varname, '(A,A)') 'Y_', trim(species_names(i))
319+
write (varname, '(A)') 'T'
306320
call s_write_variable_to_formatted_database_file(varname, t_step)
307321

308322
varname(:) = ' '
309-
310323
end if
311-
end do
324+
end if
312325

313326
! Adding the flux limiter function to the formatted database file
314327
do i = 1, E_idx - mom_idx%beg

‎src/pre_process/m_data_output.fpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -104,7 +104,7 @@ contains
104104
real(kind(0d0)) :: nbub !< Temporary bubble number density
105105
real(kind(0d0)) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params
106106
real(kind(0d0)) :: rho !< Temporary density
107-
real(kind(0d0)) :: pres !< Temporary pressure
107+
real(kind(0d0)) :: pres, Temp !< Temporary pressure
108108

109109
real(kind(0d0)) :: nR3
110110
real(kind(0d0)) :: ntmp
@@ -273,7 +273,7 @@ contains
273273
q_cons_vf(E_idx)%sf(j, 0, 0), &
274274
q_cons_vf(alf_idx)%sf(j, 0, 0), &
275275
0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho, &
276-
pi_inf, gamma, rho, qv, rhoYks, pres)
276+
pi_inf, gamma, rho, qv, rhoYks, pres, Temp)
277277
write (2, FMT) x_cb(j), pres
278278
else if ((i >= bub_idx%beg) .and. (i <= bub_idx%end) .and. bubbles) then
279279

‎src/simulation/m_cbc.fpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -646,6 +646,7 @@ contains
646646
real(kind(0d0)), dimension(num_fluids) :: adv, dadv_ds
647647
real(kind(0d0)), dimension(sys_size) :: L
648648
real(kind(0d0)), dimension(3) :: lambda
649+
real(kind(0d0)), dimension(num_species) :: Y_s
649650

650651
real(kind(0d0)) :: rho !< Cell averaged density
651652
real(kind(0d0)) :: pres !< Cell averaged pressure
@@ -803,7 +804,6 @@ contains
803804
end do
804805

805806
E = gamma*pres + pi_inf + 5d-1*rho*vel_K_sum
806-
807807
H = (E + pres)/rho
808808

809809
! Compute mixture sound speed

‎src/simulation/m_chemistry.fpp‎

Lines changed: 17 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -55,41 +55,6 @@ contains
5555

5656
end subroutine s_finalize_chemistry_module
5757

58-
#:for NORM_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
59-
subroutine s_compute_chemistry_rhs_${XYZ}$ (flux_n, rhs_vf, flux_src_vf, q_prim_vf)
60-
61-
type(vector_field), dimension(:), intent(IN) :: flux_n
62-
type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf, flux_src_vf, q_prim_vf
63-
type(int_bounds_info) :: ix, iy, iz
64-
65-
integer :: x, y, z
66-
integer :: eqn
67-
68-
integer, parameter :: mx = ${1 if NORM_DIR == 1 else 0}$
69-
integer, parameter :: my = ${1 if NORM_DIR == 2 else 0}$
70-
integer, parameter :: mz = ${1 if NORM_DIR == 3 else 0}$
71-
72-
!$acc parallel loop collapse(4) present(rhs_vf, flux_n)
73-
do x = 0, m
74-
do y = 0, n
75-
do z = 0, p
76-
77-
do eqn = chemxb, chemxe
78-
79-
! \nabla \cdot (F)
80-
rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + &
81-
(flux_n(${NORM_DIR}$)%vf(eqn)%sf(x - mx, y - my, z - mz) - &
82-
flux_n(${NORM_DIR}$)%vf(eqn)%sf(x, y, z))/d${XYZ}$ (${XYZ}$)
83-
84-
end do
85-
86-
end do
87-
end do
88-
end do
89-
90-
end subroutine s_compute_chemistry_rhs_${XYZ}$
91-
#:endfor
92-
9358
subroutine s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp, q_prim_qp)
9459

9560
type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf, q_cons_qp, q_prim_qp
@@ -104,50 +69,47 @@ contains
10469
real(kind(0d0)) :: dyn_pres
10570
real(kind(0d0)) :: E
10671

107-
real(kind(0d0)) :: rho
108-
real(kind(1.d0)), dimension(num_species) :: Ys
109-
real(kind(1.d0)), dimension(num_species) :: omega
72+
real(kind(0d0)) :: rho, omega_m, omega_T
73+
real(kind(0d0)), dimension(num_species) :: Ys
74+
real(kind(0d0)), dimension(num_species) :: omega
11075
real(kind(0d0)), dimension(num_species) :: enthalpies
11176
real(kind(0d0)) :: cp_mix
11277

11378
#:if chemistry
11479

115-
!$acc parallel loop collapse(4) private(rho)
80+
!$acc parallel loop collapse(3) private(rho)
11681
do x = 0, m
11782
do y = 0, n
11883
do z = 0, p
11984

120-
! Maybe use q_prim_vf instead?
12185
rho = 0d0
12286
do eqn = chemxb, chemxe
12387
rho = rho + q_cons_qp(eqn)%sf(x, y, z)
12488
end do
12589

12690
do eqn = chemxb, chemxe
127-
Ys(eqn - chemxb + 1) = q_cons_qp(eqn)%sf(x, y, z)/rho
91+
Ys(eqn - chemxb + 1) = q_prim_qp(eqn)%sf(x, y, z)
12892
end do
12993

13094
dyn_pres = 0d0
131-
13295
do i = momxb, momxe
133-
dyn_pres = dyn_pres + rho*q_cons_qp(i)%sf(x, y, z)* &
134-
q_cons_qp(i)%sf(x, y, z)/2d0
96+
dyn_pres = dyn_pres + q_cons_qp(i)%sf(x, y, z)* &
97+
q_prim_qp(i)%sf(x, y, z)/2d0
13598
end do
13699

137-
call get_temperature(.true., q_cons_qp(E_idx)%sf(x, y, z) - dyn_pres, &
138-
& q_prim_qp(tempxb)%sf(x, y, z), Ys, T)
100+
call get_temperature(q_cons_qp(E_idx)%sf(x, y, z)/rho - dyn_pres/rho, &
101+
& 1200d0, Ys, .true., T)
139102

140103
call get_net_production_rates(rho, T, Ys, omega)
141104

142-
q_cons_qp(tempxb)%sf(x, y, z) = T
143-
q_prim_qp(tempxb)%sf(x, y, z) = T
144-
145-
!print*, x, y, z, T, rho, Ys, omega, q_cons_qp(E_idx)%sf(x, y, z), dyn_pres
105+
call get_species_enthalpies_rt(T, enthalpies)
146106

147107
do eqn = chemxb, chemxe
148108

149-
rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + &
150-
mol_weights(eqn - chemxb + 1)*omega(eqn - chemxb + 1)
109+
omega_m = mol_weights(eqn - chemxb + 1)*omega(eqn - chemxb + 1)
110+
omega_T = omega_m*enthalpies(eqn - chemxb + 1)*gas_constant*T
111+
112+
rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m
151113

152114
end do
153115

@@ -171,9 +133,9 @@ contains
171133
integer :: eqn
172134

173135
!$acc parallel loop collapse(4)
174-
do x = ix%beg, ix%end
175-
do y = iy%beg, iy%end
176-
do z = iz%beg, iz%end
136+
do x = 0, m
137+
do y = 0, n
138+
do z = 0, p
177139

178140
do eqn = chemxb, chemxe
179141
q_cons_qp(eqn)%sf(x, y, z) = max(0d0, q_cons_qp(eqn)%sf(x, y, z))

‎src/simulation/m_data_output.fpp‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -960,7 +960,7 @@ contains
960960
real(kind(0d0)) :: E_e
961961
real(kind(0d0)), dimension(6) :: tau_e
962962
real(kind(0d0)) :: G
963-
real(kind(0d0)) :: dyn_p
963+
real(kind(0d0)) :: dyn_p, Temp
964964

965965
integer :: i, j, k, l, s, q, d !< Generic loop iterator
966966

@@ -1054,14 +1054,14 @@ contains
10541054
call s_compute_pressure( &
10551055
q_cons_vf(1)%sf(j - 2, k, l), &
10561056
q_cons_vf(alf_idx)%sf(j - 2, k, l), &
1057-
dyn_p, pi_inf, gamma, rho, qv, rhoYks(:), pres, &
1057+
dyn_p, pi_inf, gamma, rho, qv, rhoYks(:), pres, Temp, &
10581058
q_cons_vf(stress_idx%beg)%sf(j - 2, k, l), &
10591059
q_cons_vf(mom_idx%beg)%sf(j - 2, k, l), G)
10601060
else
10611061
call s_compute_pressure( &
10621062
q_cons_vf(1)%sf(j - 2, k, l), &
10631063
q_cons_vf(alf_idx)%sf(j - 2, k, l), &
1064-
dyn_p, pi_inf, gamma, rho, qv, rhoYks(:), pres)
1064+
dyn_p, pi_inf, gamma, rho, qv, rhoYks(:), pres, Temp)
10651065
end if
10661066

10671067
if (model_eqns == 4) then
@@ -1163,13 +1163,14 @@ contains
11631163
dyn_p, pi_inf, gamma, rho, qv, &
11641164
rhoYks, &
11651165
pres, &
1166+
Temp, &
11661167
q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l), &
11671168
q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l), G)
11681169
else
11691170
call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l), &
11701171
q_cons_vf(alf_idx)%sf(j - 2, k - 2, l), &
11711172
dyn_p, pi_inf, gamma, rho, qv, &
1172-
rhoYks, pres)
1173+
rhoYks, pres, Temp)
11731174
end if
11741175

11751176
if (model_eqns == 4) then
@@ -1253,14 +1254,14 @@ contains
12531254
q_cons_vf(1)%sf(j - 2, k - 2, l - 2), &
12541255
q_cons_vf(alf_idx)%sf(j - 2, k - 2, l - 2), &
12551256
dyn_p, pi_inf, gamma, rho, qv, &
1256-
rhoYks, pres, &
1257+
rhoYks, pres, Temp, &
12571258
q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l - 2), &
12581259
q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l - 2), G)
12591260
else
12601261
call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l - 2), &
12611262
q_cons_vf(alf_idx)%sf(j - 2, k - 2, l - 2), &
12621263
dyn_p, pi_inf, gamma, rho, qv, &
1263-
rhoYks, pres)
1264+
rhoYks, pres, Temp)
12641265
end if
12651266

12661267
! Compute mixture sound speed

‎src/simulation/m_global_parameters.fpp‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -979,6 +979,16 @@ contains
979979
end if
980980
! END: Volume Fraction Model =======================================
981981

982+
if (chemistry) then
983+
species_idx%beg = sys_size + 1
984+
species_idx%end = sys_size + num_species
985+
sys_size = species_idx%end
986+
987+
temperature_idx%beg = sys_size + 1
988+
temperature_idx%end = sys_size + 1
989+
sys_size = temperature_idx%end
990+
end if
991+
982992
if (qbmm .and. .not. polytropic) then
983993
allocate (MPI_IO_DATA%view(1:sys_size + 2*nb*4))
984994
allocate (MPI_IO_DATA%var(1:sys_size + 2*nb*4))
@@ -1067,16 +1077,6 @@ contains
10671077
grid_geometry = 3
10681078
end if
10691079

1070-
if (chemistry) then
1071-
species_idx%beg = sys_size + 1
1072-
species_idx%end = sys_size + num_species
1073-
sys_size = species_idx%end
1074-
1075-
temperature_idx%beg = sys_size + 1
1076-
temperature_idx%end = sys_size + 1
1077-
sys_size = temperature_idx%end
1078-
end if
1079-
10801080
momxb = mom_idx%beg
10811081
momxe = mom_idx%end
10821082
advxb = adv_idx%beg

0 commit comments

Comments
 (0)