Skip to content

Commit 55a4cb8

Browse files
sbryngelsonsbryngelson
authored
Merge branch 'master' into MPI_FFTW
2 parents 179e3b7 + eb152c5 commit 55a4cb8

File tree

26 files changed

+778
-2000
lines changed

26 files changed

+778
-2000
lines changed

‎examples/1D_inert_shocktube/case.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -70,8 +70,8 @@
7070
"mapped_weno": "T",
7171
"mp_weno": "T",
7272
"riemann_solver": 2,
73-
"wave_speeds": 2,
74-
"avg_state": 1,
73+
"wave_speeds": 1,
74+
"avg_state": 2,
7575
"bc_x%beg": -2,
7676
"bc_x%end": -3,
7777
# Chemistry

‎src/common/m_chemistry.fpp‎

Lines changed: 166 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,14 +10,37 @@ module m_chemistry
1010

1111
use m_thermochem, only: &
1212
num_species, molecular_weights, get_temperature, get_net_production_rates, &
13-
gas_constant, get_mixture_molecular_weight
13+
get_mole_fractions, get_species_binary_mass_diffusivities, &
14+
get_species_mass_diffusivities_mixavg, gas_constant, get_mixture_molecular_weight, &
15+
get_mixture_energy_mass, get_mixture_thermal_conductivity_mixavg, get_species_enthalpies_rt, &
16+
get_mixture_viscosity_mixavg
1417

1518
use m_global_parameters
1619

1720
implicit none
1821

22+
type(int_bounds_info) :: isc1, isc2, isc3
23+
$:GPU_DECLARE(create='[isc1, isc2, isc3]')
24+
integer, dimension(3) :: offsets
25+
$:GPU_DECLARE(create='[offsets]')
26+
1927
contains
2028

29+
subroutine compute_viscosity_and_inversion(T_L, Ys_L, T_R, Ys_R, Re_L, Re_R)
30+
31+
$:GPU_ROUTINE(function_name='compute_viscosity_and_inversion',parallelism='[seq]', &
32+
& cray_inline=True)
33+
34+
real(wp), intent(inout) :: T_L, T_R, Re_L, Re_R
35+
real(wp), dimension(num_species), intent(inout) :: Ys_R, Ys_L
36+
37+
call get_mixture_viscosity_mixavg(T_L, Ys_L, Re_L)
38+
call get_mixture_viscosity_mixavg(T_R, Ys_R, Re_R)
39+
Re_L = 1.0_wp/Re_L
40+
Re_R = 1.0_wp/Re_R
41+
42+
end subroutine compute_viscosity_and_inversion
43+
2144
subroutine s_compute_q_T_sf(q_T_sf, q_cons_vf, bounds)
2245

2346
! Initialize the temperature field at the start of the simulation to
@@ -129,4 +152,146 @@ contains
129152

130153
end subroutine s_compute_chemistry_reaction_flux
131154

155+
subroutine s_compute_chemistry_diffusion_flux(idir, q_prim_qp, flux_src_vf, irx, iry, irz)
156+
157+
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_qp
158+
type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf
159+
type(int_bounds_info), intent(in) :: irx, iry, irz
160+
161+
integer, intent(in) :: idir
162+
163+
real(wp), dimension(num_species) :: Xs_L, Xs_R, Xs_cell, Ys_L, Ys_R, Ys_cell
164+
real(wp), dimension(num_species) :: mass_diffusivities_mixavg1, mass_diffusivities_mixavg2
165+
real(wp), dimension(num_species) :: mass_diffusivities_mixavg_Cell, dXk_dxi, h_l, h_r, h_k
166+
real(wp), dimension(num_species) :: Mass_Diffu_Flux
167+
real(wp) :: Mass_Diffu_Energy
168+
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
169+
real(wp) :: lambda_L, lambda_R, lambda_Cell, dT_dxi, grid_spacing
170+
171+
integer :: x, y, z, i, n, eqn
172+
integer, dimension(3) :: offsets
173+
174+
isc1 = irx; isc2 = iry; isc3 = irz
175+
176+
$:GPU_UPDATE(device='[isc1,isc2,isc3]')
177+
178+
if (chemistry) then
179+
! Set offsets based on direction using array indexing
180+
offsets = 0
181+
offsets(idir) = 1
182+
183+
$:GPU_PARALLEL_LOOP(collapse=3, private='[Ys_L, Ys_R, Ys_cell, &
184+
& Xs_L, Xs_R, mass_diffusivities_mixavg1, mass_diffusivities_mixavg2, &
185+
& mass_diffusivities_mixavg_Cell, h_l, h_r, Xs_cell, h_k, &
186+
& dXk_dxi,Mass_Diffu_Flux]', copyin='[offsets]')
187+
do z = isc3%beg, isc3%end
188+
do y = isc2%beg, isc2%end
189+
do x = isc1%beg, isc1%end
190+
! Calculate grid spacing using direction-based indexing
191+
select case (idir)
192+
case (1)
193+
grid_spacing = x_cc(x + 1) - x_cc(x)
194+
case (2)
195+
grid_spacing = y_cc(y + 1) - y_cc(y)
196+
case (3)
197+
grid_spacing = z_cc(z + 1) - z_cc(z)
198+
end select
199+
200+
! Extract species mass fractions
201+
$:GPU_LOOP(parallelism='[seq]')
202+
do i = chemxb, chemxe
203+
Ys_L(i - chemxb + 1) = q_prim_qp(i)%sf(x, y, z)
204+
Ys_R(i - chemxb + 1) = q_prim_qp(i)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
205+
Ys_cell(i - chemxb + 1) = 0.5_wp*(Ys_L(i - chemxb + 1) + Ys_R(i - chemxb + 1))
206+
end do
207+
208+
! Calculate molecular weights and mole fractions
209+
call get_mixture_molecular_weight(Ys_L, MW_L)
210+
call get_mixture_molecular_weight(Ys_R, MW_R)
211+
MW_cell = 0.5_wp*(MW_L + MW_R)
212+
213+
call get_mole_fractions(MW_L, Ys_L, Xs_L)
214+
call get_mole_fractions(MW_R, Ys_R, Xs_R)
215+
216+
! Calculate gas constants and thermodynamic properties
217+
Rgas_L = gas_constant/MW_L
218+
Rgas_R = gas_constant/MW_R
219+
220+
P_L = q_prim_qp(E_idx)%sf(x, y, z)
221+
P_R = q_prim_qp(E_idx)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
222+
223+
rho_L = q_prim_qp(1)%sf(x, y, z)
224+
rho_R = q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
225+
226+
T_L = P_L/rho_L/Rgas_L
227+
T_R = P_R/rho_R/Rgas_R
228+
229+
rho_cell = 0.5_wp*(rho_L + rho_R)
230+
dT_dxi = (T_R - T_L)/grid_spacing
231+
232+
! Get transport properties
233+
call get_species_mass_diffusivities_mixavg(P_L, T_L, Ys_L, mass_diffusivities_mixavg1)
234+
call get_species_mass_diffusivities_mixavg(P_R, T_R, Ys_R, mass_diffusivities_mixavg2)
235+
236+
call get_mixture_thermal_conductivity_mixavg(T_L, Ys_L, lambda_L)
237+
call get_mixture_thermal_conductivity_mixavg(T_R, Ys_R, lambda_R)
238+
239+
call get_species_enthalpies_rt(T_L, h_l)
240+
call get_species_enthalpies_rt(T_R, h_r)
241+
242+
! Calculate species properties and gradients
243+
$:GPU_LOOP(parallelism='[seq]')
244+
do i = chemxb, chemxe
245+
h_l(i - chemxb + 1) = h_l(i - chemxb + 1)*gas_constant*T_L/molecular_weights(i - chemxb + 1)
246+
h_r(i - chemxb + 1) = h_r(i - chemxb + 1)*gas_constant*T_R/molecular_weights(i - chemxb + 1)
247+
Xs_cell(i - chemxb + 1) = 0.5_wp*(Xs_L(i - chemxb + 1) + Xs_R(i - chemxb + 1))
248+
h_k(i - chemxb + 1) = 0.5_wp*(h_l(i - chemxb + 1) + h_r(i - chemxb + 1))
249+
dXk_dxi(i - chemxb + 1) = (Xs_R(i - chemxb + 1) - Xs_L(i - chemxb + 1))/grid_spacing
250+
end do
251+
252+
! Calculate mixture-averaged diffusivities
253+
$:GPU_LOOP(parallelism='[seq]')
254+
do i = chemxb, chemxe
255+
mass_diffusivities_mixavg_Cell(i - chemxb + 1) = &
256+
(mass_diffusivities_mixavg2(i - chemxb + 1) + mass_diffusivities_mixavg1(i - chemxb + 1))/2.0_wp
257+
end do
258+
259+
lambda_Cell = 0.5_wp*(lambda_R + lambda_L)
260+
261+
! Calculate mass diffusion fluxes
262+
rho_Vic = 0.0_wp
263+
Mass_Diffu_Energy = 0.0_wp
264+
265+
$:GPU_LOOP(parallelism='[seq]')
266+
do eqn = chemxb, chemxe
267+
Mass_Diffu_Flux(eqn - chemxb + 1) = rho_cell*mass_diffusivities_mixavg_Cell(eqn - chemxb + 1)* &
268+
molecular_weights(eqn - chemxb + 1)/MW_cell*dXk_dxi(eqn - chemxb + 1)
269+
rho_Vic = rho_Vic + Mass_Diffu_Flux(eqn - chemxb + 1)
270+
Mass_Diffu_Energy = Mass_Diffu_Energy + h_k(eqn - chemxb + 1)*Mass_Diffu_Flux(eqn - chemxb + 1)
271+
end do
272+
273+
! Apply corrections for mass conservation
274+
$:GPU_LOOP(parallelism='[seq]')
275+
do eqn = chemxb, chemxe
276+
Mass_Diffu_Energy = Mass_Diffu_Energy - h_k(eqn - chemxb + 1)*Ys_cell(eqn - chemxb + 1)*rho_Vic
277+
Mass_Diffu_Flux(eqn - chemxb + 1) = Mass_Diffu_Flux(eqn - chemxb + 1) - rho_Vic*Ys_cell(eqn - chemxb + 1)
278+
end do
279+
280+
! Add thermal conduction contribution
281+
Mass_Diffu_Energy = lambda_Cell*dT_dxi + Mass_Diffu_Energy
282+
283+
! Update flux arrays
284+
flux_src_vf(E_idx)%sf(x, y, z) = flux_src_vf(E_idx)%sf(x, y, z) - Mass_Diffu_Energy
285+
286+
$:GPU_LOOP(parallelism='[seq]')
287+
do eqn = chemxb, chemxe
288+
flux_src_vf(eqn)%sf(x, y, z) = flux_src_vf(eqn)%sf(x, y, z) - Mass_diffu_Flux(eqn - chemxb + 1)
289+
end do
290+
end do
291+
end do
292+
end do
293+
end if
294+
295+
end subroutine s_compute_chemistry_diffusion_flux
296+
132297
end module m_chemistry

0 commit comments

Comments
 (0)