@@ -130,35 +130,35 @@ contains
130130 real (wp), dimension (num_species) :: omega
131131
132132 $:GPU_PARALLEL_LOOP(collapse= 3 , private= ' [x,y,z,Ys, omega, T]' )
133- do z = bounds(3 )%beg, bounds(3 )%end
134- do y = bounds(2 )%beg, bounds(2 )%end
135- do x = bounds(1 )%beg, bounds(1 )%end
136-
137- $:GPU_LOOP(parallelism= ' [seq]' )
138- do eqn = chemxb, chemxe
139- Ys(eqn - chemxb + 1 ) = q_prim_qp(eqn)%sf(x, y, z)
140- end do
133+ do z = bounds(3 )%beg, bounds(3 )%end
134+ do y = bounds(2 )%beg, bounds(2 )%end
135+ do x = bounds(1 )%beg, bounds(1 )%end
141136
142- rho = q_cons_qp(contxe)%sf(x, y, z)
143- T = q_T_sf%sf(x, y, z)
137+ $:GPU_LOOP(parallelism= ' [seq]' )
138+ do eqn = chemxb, chemxe
139+ Ys(eqn - chemxb + 1 ) = q_prim_qp(eqn)%sf(x, y, z)
140+ end do
144141
145- call get_net_production_rates(rho, T, Ys, omega)
142+ rho = q_cons_qp(contxe)%sf(x, y, z)
143+ T = q_T_sf%sf(x, y, z)
146144
147- $:GPU_LOOP(parallelism= ' [seq]' )
148- do eqn = chemxb, chemxe
149- #:block UNDEF_AMD
150- omega_m = molecular_weights(eqn - chemxb + 1 )* omega(eqn - chemxb + 1 )
151- #:endblock UNDEF_AMD
152- #:block DEF_AMD
153- omega_m = molecular_weights_nonparameter(eqn - chemxb + 1 )* omega(eqn - chemxb + 1 )
154- #:endblock DEF_AMD
155- rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m
145+ call get_net_production_rates(rho, T, Ys, omega)
156146
157- end do
147+ $:GPU_LOOP(parallelism= ' [seq]' )
148+ do eqn = chemxb, chemxe
149+ #:block UNDEF_AMD
150+ omega_m = molecular_weights(eqn - chemxb + 1 )* omega(eqn - chemxb + 1 )
151+ #:endblock UNDEF_AMD
152+ #:block DEF_AMD
153+ omega_m = molecular_weights_nonparameter(eqn - chemxb + 1 )* omega(eqn - chemxb + 1 )
154+ #:endblock DEF_AMD
155+ rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m
158156
159157 end do
158+
160159 end do
161160 end do
161+ end do
162162 $:END_GPU_PARALLEL_LOOP()
163163
164164 end subroutine s_compute_chemistry_reaction_flux
@@ -192,112 +192,112 @@ contains
192192 offsets(idir) = 1
193193
194194 $:GPU_PARALLEL_LOOP(collapse= 3 , private= ' [x,y,z,Ys_L, Ys_R, Ys_cell, Xs_L, Xs_R, mass_diffusivities_mixavg1, mass_diffusivities_mixavg2, mass_diffusivities_mixavg_Cell, h_l, h_r, Xs_cell, h_k, dXk_dxi,Mass_Diffu_Flux]' , copyin= ' [offsets]' )
195- do z = isc3%beg, isc3%end
196- do y = isc2%beg, isc2%end
197- do x = isc1%beg, isc1%end
198- ! Calculate grid spacing using direction- based indexing
199- select case (idir)
200- case (1 )
201- grid_spacing = x_cc(x + 1 ) - x_cc(x)
202- case (2 )
203- grid_spacing = y_cc(y + 1 ) - y_cc(y)
204- case (3 )
205- grid_spacing = z_cc(z + 1 ) - z_cc(z)
206- end select
207-
208- ! Extract species mass fractions
209- $:GPU_LOOP(parallelism= ' [seq]' )
210- do i = chemxb, chemxe
211- Ys_L(i - chemxb + 1 ) = q_prim_qp(i)%sf(x, y, z)
212- Ys_R(i - chemxb + 1 ) = q_prim_qp(i)%sf(x + offsets(1 ), y + offsets(2 ), z + offsets(3 ))
213- Ys_cell(i - chemxb + 1 ) = 0.5_wp * (Ys_L(i - chemxb + 1 ) + Ys_R(i - chemxb + 1 ))
214- end do
215-
216- ! Calculate molecular weights and mole fractions
217- call get_mixture_molecular_weight(Ys_L, MW_L)
218- call get_mixture_molecular_weight(Ys_R, MW_R)
219- MW_cell = 0.5_wp * (MW_L + MW_R)
220-
221- call get_mole_fractions(MW_L, Ys_L, Xs_L)
222- call get_mole_fractions(MW_R, Ys_R, Xs_R)
223-
224- ! Calculate gas constants and thermodynamic properties
225- Rgas_L = gas_constant/ MW_L
226- Rgas_R = gas_constant/ MW_R
227-
228- P_L = q_prim_qp(E_idx)%sf(x, y, z)
229- P_R = q_prim_qp(E_idx)%sf(x + offsets(1 ), y + offsets(2 ), z + offsets(3 ))
230-
231- rho_L = q_prim_qp(1 )%sf(x, y, z)
232- rho_R = q_prim_qp(1 )%sf(x + offsets(1 ), y + offsets(2 ), z + offsets(3 ))
233-
234- T_L = P_L/ rho_L/ Rgas_L
235- T_R = P_R/ rho_R/ Rgas_R
236-
237- rho_cell = 0.5_wp * (rho_L + rho_R)
238- dT_dxi = (T_R - T_L)/ grid_spacing
239-
240- ! Get transport properties
241- call get_species_mass_diffusivities_mixavg(P_L, T_L, Ys_L, mass_diffusivities_mixavg1)
242- call get_species_mass_diffusivities_mixavg(P_R, T_R, Ys_R, mass_diffusivities_mixavg2)
243-
244- call get_mixture_thermal_conductivity_mixavg(T_L, Ys_L, lambda_L)
245- call get_mixture_thermal_conductivity_mixavg(T_R, Ys_R, lambda_R)
246-
247- call get_species_enthalpies_rt(T_L, h_l)
248- call get_species_enthalpies_rt(T_R, h_r)
249-
250- ! Calculate species properties and gradients
251- $:GPU_LOOP(parallelism= ' [seq]' )
252- do i = chemxb, chemxe
253- h_l(i - chemxb + 1 ) = h_l(i - chemxb + 1 )* gas_constant* T_L/ molecular_weights(i - chemxb + 1 )
254- h_r(i - chemxb + 1 ) = h_r(i - chemxb + 1 )* gas_constant* T_R/ molecular_weights(i - chemxb + 1 )
255- Xs_cell(i - chemxb + 1 ) = 0.5_wp * (Xs_L(i - chemxb + 1 ) + Xs_R(i - chemxb + 1 ))
256- h_k(i - chemxb + 1 ) = 0.5_wp * (h_l(i - chemxb + 1 ) + h_r(i - chemxb + 1 ))
257- dXk_dxi(i - chemxb + 1 ) = (Xs_R(i - chemxb + 1 ) - Xs_L(i - chemxb + 1 ))/ grid_spacing
258- end do
259-
260- ! Calculate mixture- averaged diffusivities
261- $:GPU_LOOP(parallelism= ' [seq]' )
262- do i = chemxb, chemxe
263- mass_diffusivities_mixavg_Cell(i - chemxb + 1 ) = &
264- (mass_diffusivities_mixavg2(i - chemxb + 1 ) + mass_diffusivities_mixavg1(i - chemxb + 1 ))/ 2.0_wp
265- end do
266-
267- lambda_Cell = 0.5_wp * (lambda_R + lambda_L)
268-
269- ! Calculate mass diffusion fluxes
270- rho_Vic = 0.0_wp
271- Mass_Diffu_Energy = 0.0_wp
272-
273- $:GPU_LOOP(parallelism= ' [seq]' )
274- do eqn = chemxb, chemxe
275- Mass_Diffu_Flux(eqn - chemxb + 1 ) = rho_cell* mass_diffusivities_mixavg_Cell(eqn - chemxb + 1 )* &
276- molecular_weights(eqn - chemxb + 1 )/ MW_cell* dXk_dxi(eqn - chemxb + 1 )
277- rho_Vic = rho_Vic + Mass_Diffu_Flux(eqn - chemxb + 1 )
278- Mass_Diffu_Energy = Mass_Diffu_Energy + h_k(eqn - chemxb + 1 )* Mass_Diffu_Flux(eqn - chemxb + 1 )
279- end do
280-
281- ! Apply corrections for mass conservation
282- $:GPU_LOOP(parallelism= ' [seq]' )
283- do eqn = chemxb, chemxe
284- Mass_Diffu_Energy = Mass_Diffu_Energy - h_k(eqn - chemxb + 1 )* Ys_cell(eqn - chemxb + 1 )* rho_Vic
285- Mass_Diffu_Flux(eqn - chemxb + 1 ) = Mass_Diffu_Flux(eqn - chemxb + 1 ) - rho_Vic* Ys_cell(eqn - chemxb + 1 )
286- end do
287-
288- ! Add thermal conduction contribution
289- Mass_Diffu_Energy = lambda_Cell* dT_dxi + Mass_Diffu_Energy
290-
291- ! Update flux arrays
292- flux_src_vf(E_idx)%sf(x, y, z) = flux_src_vf(E_idx)%sf(x, y, z) - Mass_Diffu_Energy
293-
294- $:GPU_LOOP(parallelism= ' [seq]' )
295- do eqn = chemxb, chemxe
296- flux_src_vf(eqn)%sf(x, y, z) = flux_src_vf(eqn)%sf(x, y, z) - Mass_diffu_Flux(eqn - chemxb + 1 )
297- end do
195+ do z = isc3%beg, isc3%end
196+ do y = isc2%beg, isc2%end
197+ do x = isc1%beg, isc1%end
198+ ! Calculate grid spacing using direction- based indexing
199+ select case (idir)
200+ case (1 )
201+ grid_spacing = x_cc(x + 1 ) - x_cc(x)
202+ case (2 )
203+ grid_spacing = y_cc(y + 1 ) - y_cc(y)
204+ case (3 )
205+ grid_spacing = z_cc(z + 1 ) - z_cc(z)
206+ end select
207+
208+ ! Extract species mass fractions
209+ $:GPU_LOOP(parallelism= ' [seq]' )
210+ do i = chemxb, chemxe
211+ Ys_L(i - chemxb + 1 ) = q_prim_qp(i)%sf(x, y, z)
212+ Ys_R(i - chemxb + 1 ) = q_prim_qp(i)%sf(x + offsets(1 ), y + offsets(2 ), z + offsets(3 ))
213+ Ys_cell(i - chemxb + 1 ) = 0.5_wp * (Ys_L(i - chemxb + 1 ) + Ys_R(i - chemxb + 1 ))
214+ end do
215+
216+ ! Calculate molecular weights and mole fractions
217+ call get_mixture_molecular_weight(Ys_L, MW_L)
218+ call get_mixture_molecular_weight(Ys_R, MW_R)
219+ MW_cell = 0.5_wp * (MW_L + MW_R)
220+
221+ call get_mole_fractions(MW_L, Ys_L, Xs_L)
222+ call get_mole_fractions(MW_R, Ys_R, Xs_R)
223+
224+ ! Calculate gas constants and thermodynamic properties
225+ Rgas_L = gas_constant/ MW_L
226+ Rgas_R = gas_constant/ MW_R
227+
228+ P_L = q_prim_qp(E_idx)%sf(x, y, z)
229+ P_R = q_prim_qp(E_idx)%sf(x + offsets(1 ), y + offsets(2 ), z + offsets(3 ))
230+
231+ rho_L = q_prim_qp(1 )%sf(x, y, z)
232+ rho_R = q_prim_qp(1 )%sf(x + offsets(1 ), y + offsets(2 ), z + offsets(3 ))
233+
234+ T_L = P_L/ rho_L/ Rgas_L
235+ T_R = P_R/ rho_R/ Rgas_R
236+
237+ rho_cell = 0.5_wp * (rho_L + rho_R)
238+ dT_dxi = (T_R - T_L)/ grid_spacing
239+
240+ ! Get transport properties
241+ call get_species_mass_diffusivities_mixavg(P_L, T_L, Ys_L, mass_diffusivities_mixavg1)
242+ call get_species_mass_diffusivities_mixavg(P_R, T_R, Ys_R, mass_diffusivities_mixavg2)
243+
244+ call get_mixture_thermal_conductivity_mixavg(T_L, Ys_L, lambda_L)
245+ call get_mixture_thermal_conductivity_mixavg(T_R, Ys_R, lambda_R)
246+
247+ call get_species_enthalpies_rt(T_L, h_l)
248+ call get_species_enthalpies_rt(T_R, h_r)
249+
250+ ! Calculate species properties and gradients
251+ $:GPU_LOOP(parallelism= ' [seq]' )
252+ do i = chemxb, chemxe
253+ h_l(i - chemxb + 1 ) = h_l(i - chemxb + 1 )* gas_constant* T_L/ molecular_weights(i - chemxb + 1 )
254+ h_r(i - chemxb + 1 ) = h_r(i - chemxb + 1 )* gas_constant* T_R/ molecular_weights(i - chemxb + 1 )
255+ Xs_cell(i - chemxb + 1 ) = 0.5_wp * (Xs_L(i - chemxb + 1 ) + Xs_R(i - chemxb + 1 ))
256+ h_k(i - chemxb + 1 ) = 0.5_wp * (h_l(i - chemxb + 1 ) + h_r(i - chemxb + 1 ))
257+ dXk_dxi(i - chemxb + 1 ) = (Xs_R(i - chemxb + 1 ) - Xs_L(i - chemxb + 1 ))/ grid_spacing
258+ end do
259+
260+ ! Calculate mixture- averaged diffusivities
261+ $:GPU_LOOP(parallelism= ' [seq]' )
262+ do i = chemxb, chemxe
263+ mass_diffusivities_mixavg_Cell(i - chemxb + 1 ) = &
264+ (mass_diffusivities_mixavg2(i - chemxb + 1 ) + mass_diffusivities_mixavg1(i - chemxb + 1 ))/ 2.0_wp
265+ end do
266+
267+ lambda_Cell = 0.5_wp * (lambda_R + lambda_L)
268+
269+ ! Calculate mass diffusion fluxes
270+ rho_Vic = 0.0_wp
271+ Mass_Diffu_Energy = 0.0_wp
272+
273+ $:GPU_LOOP(parallelism= ' [seq]' )
274+ do eqn = chemxb, chemxe
275+ Mass_Diffu_Flux(eqn - chemxb + 1 ) = rho_cell* mass_diffusivities_mixavg_Cell(eqn - chemxb + 1 )* &
276+ molecular_weights(eqn - chemxb + 1 )/ MW_cell* dXk_dxi(eqn - chemxb + 1 )
277+ rho_Vic = rho_Vic + Mass_Diffu_Flux(eqn - chemxb + 1 )
278+ Mass_Diffu_Energy = Mass_Diffu_Energy + h_k(eqn - chemxb + 1 )* Mass_Diffu_Flux(eqn - chemxb + 1 )
279+ end do
280+
281+ ! Apply corrections for mass conservation
282+ $:GPU_LOOP(parallelism= ' [seq]' )
283+ do eqn = chemxb, chemxe
284+ Mass_Diffu_Energy = Mass_Diffu_Energy - h_k(eqn - chemxb + 1 )* Ys_cell(eqn - chemxb + 1 )* rho_Vic
285+ Mass_Diffu_Flux(eqn - chemxb + 1 ) = Mass_Diffu_Flux(eqn - chemxb + 1 ) - rho_Vic* Ys_cell(eqn - chemxb + 1 )
286+ end do
287+
288+ ! Add thermal conduction contribution
289+ Mass_Diffu_Energy = lambda_Cell* dT_dxi + Mass_Diffu_Energy
290+
291+ ! Update flux arrays
292+ flux_src_vf(E_idx)%sf(x, y, z) = flux_src_vf(E_idx)%sf(x, y, z) - Mass_Diffu_Energy
293+
294+ $:GPU_LOOP(parallelism= ' [seq]' )
295+ do eqn = chemxb, chemxe
296+ flux_src_vf(eqn)%sf(x, y, z) = flux_src_vf(eqn)%sf(x, y, z) - Mass_diffu_Flux(eqn - chemxb + 1 )
298297 end do
299298 end do
300299 end do
300+ end do
301301 $:END_GPU_PARALLEL_LOOP()
302302 end if
303303
0 commit comments