MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_chemistry.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
2!>
3!! @file
4!! @brief Contains module m_chemistry
5!! @author Henry Le Berre <hberre3@gatech.edu>
6
7# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
9# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
10# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
15
16# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
19
20# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21
22# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23
24# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27
28# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35! New line at end of file is required for FYPP
36# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
37# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
38# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
39# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
44
45# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48
49# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50
51# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52
53# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54
55# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56
57# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60
61# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64! New line at end of file is required for FYPP
65# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
66
67# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
68# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
69# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
70# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
71# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
72
73# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
74
75# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
76
77# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
78
79# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
80
81# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
82
83# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
84
85# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86
87# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88
89# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
94
95# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
96
97# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
98
99# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
100
101# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102
103# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104
105# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
106
107# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143! New line at end of file is required for FYPP
144# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
145# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
146# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
147# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
148# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
152
153# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
156
157# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
158
159# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
160
161# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162
163# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
164
165# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166
167# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168
169# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170
171# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172! New line at end of file is required for FYPP
173# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
174
175# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
176
177# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
178
179# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
180
181# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
182
183# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
184
185# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
186
187# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
188
189# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
190
191# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
192
193# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230! New line at end of file is required for FYPP
231# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
232
233! GPU parallel region (scalar reductions, maxval/minval)
234# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
235
236! GPU parallel loop over threads (most common GPU macro)
237# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
238
239! Required closing for GPU_PARALLEL_LOOP
240# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
241
242! Mark routine for device compilation
243# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
244
245! Declare device-resident data
246# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
247
248! Inner loop within a GPU parallel region
249# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
250
251! Scoped GPU data region
252# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
253
254! Host code with device pointers (for MPI with GPU buffers)
255# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257! Allocate device memory (unscoped)
258# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
259
260! Free device memory
261# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263! Atomic operation on device
264# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
265
266! End atomic capture block
267# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
268
269! Copy data between host and device
270# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
271
272! Synchronization barrier
273# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
274
275! Import GPU library module (openacc or omp_lib)
276# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
277
278! Emit code only for AMD compiler
279# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
280
281! Emit code for non-Cray compilers
282# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
283
284! Emit code only for Cray compiler
285# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
286
287! Emit code for non-NVIDIA compilers
288# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
289
290# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
292! New line at end of file is required for FYPP
293# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
294
295# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
296
297! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
298! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
299! example see misc/nvidia_uvm/bind.sh.
300# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
301
302! Allocate and create GPU device memory
303# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
304
305! Free GPU device memory and deallocate
306# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
307
308! Cray-specific GPU pointer setup for vector fields
309# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
310
311! Cray-specific GPU pointer setup for scalar fields
312# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
313
314! Cray-specific GPU pointer setup for acoustic source spatials
315# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
316
317# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
318
319# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
320! New line at end of file is required for FYPP
321# 7 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp" 2
322# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
323! This file exists so that Fypp can be run without generating case.fpp files for
324! each target. This is useful when generating documentation, for example. This
325! should also let MFC be built with CMake directly, without invoking mfc.sh.
326
327! For pre-process.
328# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
329
330! For moving immersed boundaries in simulation
331# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
332# 8 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp" 2
333
334!> @brief Multi-species chemistry interface for thermodynamic properties, reaction rates, and transport coefficients
335module m_chemistry
336
337 use m_thermochem, only: num_species, molecular_weights, get_temperature, get_net_production_rates, &
338 & get_creation_destruction_rates, get_mole_fractions, get_species_binary_mass_diffusivities, &
339 & get_species_mass_diffusivities_mixavg, gas_constant, get_mixture_molecular_weight, get_mixture_energy_mass, &
340 & get_mixture_thermal_conductivity_mixavg, get_species_enthalpies_rt, get_mixture_viscosity_mixavg, &
341 & get_mixture_specific_heat_cp_mass, get_mixture_enthalpy_mass
342
343 use m_global_parameters
344
345 implicit none
346
347 type(int_bounds_info) :: isc1, isc2, isc3
348
349# 23 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
350#if defined(MFC_OpenACC)
351# 23 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
352!$acc declare create(isc1, isc2, isc3)
353# 23 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
354#elif defined(MFC_OpenMP)
355# 23 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
356!$omp declare target (isc1, isc2, isc3)
357# 23 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
358#endif
359 integer, dimension(3) :: offsets
360
361# 25 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
362#if defined(MFC_OpenACC)
363# 25 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
364!$acc declare create(offsets)
365# 25 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
366#elif defined(MFC_OpenMP)
367# 25 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
368!$omp declare target (offsets)
369# 25 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
370#endif
371
372contains
373
374 !> Compute mixture viscosities for left and right states and invert them for use as reciprocal Reynolds numbers.
375 subroutine compute_viscosity_and_inversion(T_L, Ys_L, T_R, Ys_R, Re_L, Re_R)
376
377
378# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
379#ifdef _CRAYFTN
380# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
381#if MFC_OpenACC
382# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
383!$acc routine seq
384# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
385#elif MFC_OpenMP
386# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
387
388# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
389
390# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
391!$omp declare target device_type(any)
392# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
393#else
394# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
395!DIR$ INLINEALWAYS compute_viscosity_and_inversion
396# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
397#endif
398# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
399#elif MFC_OpenACC
400# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
401!$acc routine seq
402# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
403#elif MFC_OpenMP
404# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
405
406# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
407
408# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
409!$omp declare target device_type(any)
410# 32 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
411#endif
412
413 real(wp), intent(inout) :: T_L, T_R, Re_L, Re_R
414 real(wp), dimension(num_species), intent(inout) :: Ys_R, Ys_L
415
416 call get_mixture_viscosity_mixavg(t_l, ys_l, re_l)
417 call get_mixture_viscosity_mixavg(t_r, ys_r, re_r)
418 ! Convert dynamic viscosity to inverse (MFC stores 1/mu for Reynolds number convention)
419 re_l = 1.0_wp/re_l
420 re_r = 1.0_wp/re_r
421
422 end subroutine compute_viscosity_and_inversion
423
424 !> Initialize the temperature field from conservative variables by inverting the energy equation.
425 subroutine s_compute_q_t_sf(q_T_sf, q_cons_vf, bounds)
426
427 ! Initialize the temperature field at the start of the simulation to reasonable values. Temperature is computed the regular
428 ! way using the conservative variables.
429
430 type(scalar_field), intent(inout) :: q_T_sf
431 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
432 type(int_bounds_info), dimension(1:3), intent(in) :: bounds
433 integer :: x, y, z, eqn
434 real(wp) :: energy, T_in
435 real(wp), dimension(num_species) :: Ys
436
437 do z = bounds(3)%beg, bounds(3)%end
438 do y = bounds(2)%beg, bounds(2)%end
439 do x = bounds(1)%beg, bounds(1)%end
440 do eqn = eqn_idx%species%beg, eqn_idx%species%end
441 ys(eqn - eqn_idx%species%beg + 1) = q_cons_vf(eqn)%sf(x, y, z)/q_cons_vf(eqn_idx%cont%beg)%sf(x, y, z)
442 end do
443
444 ! e = E - 1/2*|u|^2 cons. eqn_idx%E = \rho E cons. eqn_idx%cont%beg = \rho (1-fluid model) cons. eqn_idx%mom%beg
445 ! + i = \rho u_i
446 energy = q_cons_vf(eqn_idx%E)%sf(x, y, z)/q_cons_vf(eqn_idx%cont%beg)%sf(x, y, z)
447 do eqn = eqn_idx%mom%beg, eqn_idx%mom%end
448 energy = energy - 0.5_wp*(q_cons_vf(eqn)%sf(x, y, z)/q_cons_vf(eqn_idx%cont%beg)%sf(x, y, z))**2._wp
449 end do
450
451 t_in = real(q_t_sf%sf(x, y, z), kind=wp)
452 call get_temperature(energy, dflt_t_guess, ys, .true., t_in)
453 q_t_sf%sf(x, y, z) = t_in
454 end do
455 end do
456 end do
457
458 end subroutine s_compute_q_t_sf
459
460 !> Compute the temperature field from primitive variables using the ideal gas law and mixture molecular weight.
461 subroutine s_compute_t_from_primitives(q_T_sf, q_prim_vf, bounds)
462
463 type(scalar_field), intent(inout) :: q_T_sf
464 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
465 type(int_bounds_info), dimension(1:3), intent(in) :: bounds
466 integer :: x, y, z, i
467 real(wp), dimension(num_species) :: Ys
468 real(wp) :: mix_mol_weight
469
470 do z = bounds(3)%beg, bounds(3)%end
471 do y = bounds(2)%beg, bounds(2)%end
472 do x = bounds(1)%beg, bounds(1)%end
473 do i = eqn_idx%species%beg, eqn_idx%species%end
474 ys(i - eqn_idx%species%beg + 1) = q_prim_vf(i)%sf(x, y, z)
475 end do
476
477 call get_mixture_molecular_weight(ys, mix_mol_weight)
478 q_t_sf%sf(x, y, z) = q_prim_vf(eqn_idx%E)%sf(x, y, z)*mix_mol_weight/(gas_constant*q_prim_vf(1)%sf(x, y, z))
479 end do
480 end do
481 end do
482
483 end subroutine s_compute_t_from_primitives
484
485 !> Add chemical reaction source terms to the species transport RHS using net production rates.
486 subroutine s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp, q_T_sf, q_prim_qp, bounds)
487
488 type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
489 type(scalar_field), intent(inout) :: q_T_sf
490 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_qp, q_prim_qp
491 type(int_bounds_info), dimension(1:3), intent(in) :: bounds
492 integer :: x, y, z
493 integer :: eqn
494 real(wp) :: T
495 real(wp) :: rho, omega_m
496
497# 122 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
498 real(wp), dimension(num_species) :: Ys
499 real(wp), dimension(num_species) :: omega
500# 125 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
501
502
503# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
504
505# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
506#if defined(MFC_OpenACC)
507# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
508!$acc parallel loop collapse(3) gang vector default(present) private(Ys, omega, eqn, T, rho, omega_m) copyin(bounds)
509# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
510#elif defined(MFC_OpenMP)
511# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
512
513# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
514
515# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
516
517# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
518!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
519# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
520!$omp& private(Ys, omega, eqn, T, rho, omega_m) map(to:bounds)
521# 126 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
522#endif
523 do z = bounds(3)%beg, bounds(3)%end
524 do y = bounds(2)%beg, bounds(2)%end
525 do x = bounds(1)%beg, bounds(1)%end
526
527# 130 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
528#if defined(MFC_OpenACC)
529# 130 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
530!$acc loop seq
531# 130 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
532#elif defined(MFC_OpenMP)
533# 130 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
534
535# 130 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
536#endif
537 do eqn = eqn_idx%species%beg, eqn_idx%species%end
538 ys(eqn - eqn_idx%species%beg + 1) = q_prim_qp(eqn)%sf(x, y, z)
539 end do
540
541 rho = q_cons_qp(eqn_idx%cont%end)%sf(x, y, z)
542 t = q_t_sf%sf(x, y, z)
543
544 call get_net_production_rates(rho, t, ys, omega)
545
546
547# 140 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
548#if defined(MFC_OpenACC)
549# 140 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
550!$acc loop seq
551# 140 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
552#elif defined(MFC_OpenMP)
553# 140 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
554
555# 140 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
556#endif
557 do eqn = eqn_idx%species%beg, eqn_idx%species%end
558 omega_m = molecular_weights(eqn - eqn_idx%species%beg + 1)*omega(eqn - eqn_idx%species%beg + 1)
559 rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m
560 end do
561 end do
562 end do
563 end do
564
565# 148 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
566#if defined(MFC_OpenACC)
567# 148 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
568!$acc end parallel loop
569# 148 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
570#elif defined(MFC_OpenMP)
571# 148 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
572
573# 148 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
574!$omp end target teams loop
575# 148 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
576#endif
577
578 end subroutine s_compute_chemistry_reaction_flux
579
580 !> Operator-split integration of the reaction source with an alpha-QSS (Mott quasi-steady-state) reactor. Called after the flow
581 !! update: each cell's constant-(rho, e) reactor is advanced over dtime with chem_params%reaction_substeps predictor-corrector
582 !! sub-steps -- or, when adap_substeps = T, an adaptive count nsub clamped to [reaction_substeps, reaction_substeps_max] and
583 !! sized once per rank from the rank's stiffest cell -- updating the species partial densities and temperature in place. Mixture
584 !! density, momentum, and total energy are unchanged (reactions convert chemical to thermal energy at fixed internal energy).
585 !! The alpha-QSS update treats each species' destruction as a pseudo-first-order loss, so it is stable for stiff ignition where
586 !! an explicit source would overshoot and diverge, and relaxes to the correct chemical equilibrium rather than over-heating.
587 !! Reaction is split from the flow update at first order (Lie-Trotter), and each sub-step renormalizes the mass fractions to sum
588 !! to one (which does not strictly conserve elemental composition).
589 subroutine s_chemistry_reaction_substep(q_cons_vf, q_T_sf, dtime, bounds)
590
591 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
592 type(scalar_field), intent(inout) :: q_T_sf
593 real(wp), intent(in) :: dtime
594 type(int_bounds_info), dimension(1:3), intent(in) :: bounds
595 integer :: x, y, z, eqn, s, nsub
596 real(wp) :: rho, energy, T, T_new, dt_sub, Ysum
597 real(wp) :: r, r2, wr, loss_i, prod_p, loss_p, Lbar, pbar
598 real(wp) :: stiff_max, cell_stiff
599 real(wp), parameter :: y_floor = 1.e-16_wp
600 ! stiff_target: fractional net composition change per sub-step targeted when sizing the adaptive
601 ! nsub. An uncalibrated engineering default -- alpha-QSS is unconditionally stable, so it trades
602 ! accuracy for cost (never stability), and the cost is bounded by reaction_substeps_max.
603 real(wp), parameter :: stiff_target = 0.5_wp
604
605# 180 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
606 real(wp), dimension(num_species) :: Ys, cdot, ddot, y0, prod0, Lloss, alp
607# 182 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
608
609 if (chem_params%adap_substeps) then
610 ! Pass 1: per-rank local stiffness probe -> adapt nsub for this step, no MPI. Each rank
611 ! sizes its own work from the largest fractional net species change any of its cells sees.
612 stiff_max = 0._wp
613
614# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
615
616# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
617#if defined(MFC_OpenACC)
618# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
619!$acc parallel loop collapse(3) gang vector default(present) private(Ys, cdot, eqn, rho, T, T_new, energy, wr, cell_stiff) reduction(MAX:stiff_max) copyin(bounds, dtime)
620# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
621#elif defined(MFC_OpenMP)
622# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
623
624# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
625
626# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
627
628# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
629!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
630# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
631!$omp& private(Ys, cdot, eqn, rho, T, T_new, energy, wr, cell_stiff) reduction(MAX:stiff_max) map(to:bounds, dtime)
632# 187 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
633#endif
634# 189 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
635 do z = bounds(3)%beg, bounds(3)%end
636 do y = bounds(2)%beg, bounds(2)%end
637 do x = bounds(1)%beg, bounds(1)%end
638 rho = q_cons_vf(eqn_idx%cont%beg)%sf(x, y, z)
639
640# 193 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
641#if defined(MFC_OpenACC)
642# 193 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
643!$acc loop seq
644# 193 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
645#elif defined(MFC_OpenMP)
646# 193 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
647
648# 193 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
649#endif
650 do eqn = eqn_idx%species%beg, eqn_idx%species%end
651 ys(eqn - eqn_idx%species%beg + 1) = q_cons_vf(eqn)%sf(x, y, z)/rho
652 end do
653 ! q_T_sf still holds the pre-update RK-stage temperature; re-solve T from the fresh
654 ! post-advection internal energy so a just-shock-heated cell is probed at its true
655 ! (hot) temperature. Otherwise the stiffness is under-read and nsub under-sizes
656 ! exactly at an ignition front.
657 energy = q_cons_vf(eqn_idx%E)%sf(x, y, z)/rho
658
659# 202 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
660#if defined(MFC_OpenACC)
661# 202 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
662!$acc loop seq
663# 202 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
664#elif defined(MFC_OpenMP)
665# 202 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
666
667# 202 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
668#endif
669 do eqn = eqn_idx%mom%beg, eqn_idx%mom%end
670 energy = energy - 0.5_wp*(q_cons_vf(eqn)%sf(x, y, z)/rho)**2
671 end do
672 t = q_t_sf%sf(x, y, z)
673 call get_temperature(energy, t, ys, .true., t_new)
674 t = t_new
675 ! Net rate (creation - destruction) on purpose: nsub sizes the accuracy of the
676 ! composition trajectory, which is set by how fast Ys actually moves, not by the
677 ! raw forward/reverse magnitudes. In fast partial equilibrium the net is ~0 and the
678 ! composition is static, so the floor is adequate; the alpha-QSS update is itself
679 ! unconditionally stable, so an under-sized nsub loses accuracy, never stability.
680 call get_net_production_rates(rho, t, ys, cdot) ! net omega in cdot
681 cell_stiff = 0._wp
682
683# 216 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
684#if defined(MFC_OpenACC)
685# 216 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
686!$acc loop seq
687# 216 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
688#elif defined(MFC_OpenMP)
689# 216 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
690
691# 216 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
692#endif
693 do eqn = 1, num_species
694 wr = molecular_weights(eqn)/rho
695 cell_stiff = max(cell_stiff, dtime*abs(wr*cdot(eqn))/max(ys(eqn), y_floor))
696 end do
697 stiff_max = max(stiff_max, cell_stiff)
698 end do
699 end do
700 end do
701
702# 225 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
703#if defined(MFC_OpenACC)
704# 225 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
705!$acc end parallel loop
706# 225 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
707#elif defined(MFC_OpenMP)
708# 225 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
709
710# 225 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
711!$omp end target teams loop
712# 225 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
713#endif
714 nsub = ceiling(max(real(chem_params%reaction_substeps, wp), min(real(chem_params%reaction_substeps_max, wp), &
715 & stiff_max/stiff_target)))
716 else
717 nsub = chem_params%reaction_substeps
718 end if
719 dt_sub = dtime/real(nsub, wp)
720
721
722# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
723
724# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
725#if defined(MFC_OpenACC)
726# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
727!$acc parallel loop collapse(3) gang vector default(present) private(Ys, cdot, ddot, y0, prod0, Lloss, alp, eqn, s, rho, energy, T, T_new, Ysum, r, r2, wr, loss_i, prod_p, loss_p, Lbar, pbar) &
728# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
729!$acc& copyin(bounds, dt_sub, nsub)
730# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
731#elif defined(MFC_OpenMP)
732# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
733
734# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
735
736# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
737
738# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
739!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
740# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
741!$omp& private(Ys, cdot, ddot, y0, prod0, Lloss, alp, eqn, s, rho, energy, T, T_new, Ysum, r, r2, wr, loss_i, prod_p, loss_p, Lbar, pbar) map(to:bounds, dt_sub, nsub)
742# 233 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
743#endif
744# 235 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
745 do z = bounds(3)%beg, bounds(3)%end
746 do y = bounds(2)%beg, bounds(2)%end
747 do x = bounds(1)%beg, bounds(1)%end
748 rho = q_cons_vf(eqn_idx%cont%beg)%sf(x, y, z)
749
750
751# 240 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
752#if defined(MFC_OpenACC)
753# 240 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
754!$acc loop seq
755# 240 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
756#elif defined(MFC_OpenMP)
757# 240 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
758
759# 240 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
760#endif
761 do eqn = eqn_idx%species%beg, eqn_idx%species%end
762 ys(eqn - eqn_idx%species%beg + 1) = q_cons_vf(eqn)%sf(x, y, z)/rho
763 end do
764
765 ! internal energy per mass, held fixed through the reactor sub-steps
766 energy = q_cons_vf(eqn_idx%E)%sf(x, y, z)/rho
767
768# 247 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
769#if defined(MFC_OpenACC)
770# 247 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
771!$acc loop seq
772# 247 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
773#elif defined(MFC_OpenMP)
774# 247 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
775
776# 247 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
777#endif
778 do eqn = eqn_idx%mom%beg, eqn_idx%mom%end
779 energy = energy - 0.5_wp*(q_cons_vf(eqn)%sf(x, y, z)/rho)**2
780 end do
781
782 ! re-solve T from the fresh internal energy so the first predictor sub-step starts
783 ! from the post-advection state (q_T_sf holds the pre-update RK-stage temperature).
784 t = q_t_sf%sf(x, y, z)
785 call get_temperature(energy, t, ys, .true., t_new)
786 t = t_new
787
788 do s = 1, nsub
789 ! predictor: rates at the start of the sub-step (one fused pass fills both)
790 call get_creation_destruction_rates(rho, t, ys, cdot, ddot)
791
792# 261 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
793#if defined(MFC_OpenACC)
794# 261 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
795!$acc loop seq
796# 261 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
797#elif defined(MFC_OpenMP)
798# 261 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
799
800# 261 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
801#endif
802 do eqn = 1, num_species
803 y0(eqn) = ys(eqn)
804 wr = molecular_weights(eqn)/rho
805 prod0(eqn) = wr*cdot(eqn) ! mass-fraction production
806 loss_i = wr*ddot(eqn) ! mass-fraction loss
807 lloss(eqn) = loss_i/max(ys(eqn), y_floor) ! pseudo-first-order loss rate
808 r = dt_sub*lloss(eqn); r2 = r*r
809 alp(eqn) = (180._wp + 60._wp*r + 11._wp*r2 + r2*r)/(360._wp + 60._wp*r + 12._wp*r2 + r2*r)
810 ys(eqn) = y0(eqn) + dt_sub*(prod0(eqn) - loss_i)/(1._wp + alp(eqn)*dt_sub*lloss(eqn))
811 if (ys(eqn) < 0._wp) ys(eqn) = 0._wp
812 end do
813 ! corrector: re-evaluate rates at the predicted state (T is the Newton guess; T_new is intent(out))
814 call get_temperature(energy, t, ys, .true., t_new)
815 call get_creation_destruction_rates(rho, t_new, ys, cdot, ddot)
816 ysum = 0._wp
817
818# 277 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
819#if defined(MFC_OpenACC)
820# 277 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
821!$acc loop seq
822# 277 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
823#elif defined(MFC_OpenMP)
824# 277 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
825
826# 277 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
827#endif
828 do eqn = 1, num_species
829 wr = molecular_weights(eqn)/rho
830 prod_p = wr*cdot(eqn)
831 loss_p = wr*ddot(eqn)
832 lbar = 0.5_wp*(lloss(eqn) + loss_p/max(ys(eqn), y_floor))
833 ! reuse the predictor's alp(eqn) here (and in the denominator) rather than
834 ! recomputing alpha from the averaged loss Lbar -- a deliberate CHEMEQ2
835 ! simplification, stability-neutral (alpha in [0.5, 1]) and mitigated by sub-stepping.
836 pbar = alp(eqn)*prod_p + (1._wp - alp(eqn))*prod0(eqn)
837 ys(eqn) = y0(eqn) + dt_sub*(pbar - lbar*y0(eqn))/(1._wp + alp(eqn)*dt_sub*lbar)
838 if (ys(eqn) < 0._wp) ys(eqn) = 0._wp
839 ysum = ysum + ys(eqn)
840 end do
841 if (ysum > y_floor) then
842
843# 292 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
844#if defined(MFC_OpenACC)
845# 292 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
846!$acc loop seq
847# 292 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
848#elif defined(MFC_OpenMP)
849# 292 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
850
851# 292 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
852#endif
853 do eqn = 1, num_species
854 ys(eqn) = ys(eqn)/ysum
855 end do
856 else
857 ! Degenerate corrector (every species clipped to zero): fall back to the
858 ! sub-step's starting composition rather than dividing by a vanishing sum.
859
860# 299 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
861#if defined(MFC_OpenACC)
862# 299 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
863!$acc loop seq
864# 299 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
865#elif defined(MFC_OpenMP)
866# 299 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
867
868# 299 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
869#endif
870 do eqn = 1, num_species
871 ys(eqn) = y0(eqn)
872 end do
873 end if
874 call get_temperature(energy, t, ys, .true., t_new)
875 t = t_new
876 end do
877
878
879# 308 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
880#if defined(MFC_OpenACC)
881# 308 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
882!$acc loop seq
883# 308 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
884#elif defined(MFC_OpenMP)
885# 308 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
886
887# 308 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
888#endif
889 do eqn = eqn_idx%species%beg, eqn_idx%species%end
890 q_cons_vf(eqn)%sf(x, y, z) = rho*ys(eqn - eqn_idx%species%beg + 1)
891 end do
892 q_t_sf%sf(x, y, z) = t
893 end do
894 end do
895 end do
896
897# 316 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
898#if defined(MFC_OpenACC)
899# 316 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
900!$acc end parallel loop
901# 316 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
902#elif defined(MFC_OpenMP)
903# 316 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
904
905# 316 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
906!$omp end target teams loop
907# 316 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
908#endif
909
910 end subroutine s_chemistry_reaction_substep
911
912 !> Compute species mass diffusion fluxes at cell interfaces using mixture-averaged diffusivities.
913 subroutine s_compute_chemistry_diffusion_flux(idir, q_prim_qp, flux_src_vf, irx, iry, irz, q_T_sf)
914
915 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_qp
916 type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf
917 type(int_bounds_info), intent(in) :: irx, iry, irz
918 integer, intent(in) :: idir
919 type(scalar_field), intent(in) :: q_T_sf
920
921# 335 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
922 real(wp), dimension(num_species) :: Xs_L, Xs_R, Xs_cell, Ys_L, Ys_R, Ys_cell
923 real(wp), dimension(num_species) :: mass_diffusivities_mixavg1, mass_diffusivities_mixavg2
924 real(wp), dimension(num_species) :: mass_diffusivities_mixavg_Cell, dXk_dxi, h_l, h_r, h_k
925 real(wp), dimension(num_species) :: Mass_Diffu_Flux, dYk_dxi
926# 340 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
927
928 real(wp) :: Mass_Diffu_Energy
929 real(wp) :: MW_L, MW_R, MW_cell, T_L, T_R, P_L, P_R, rho_L, rho_R, rho_cell, rho_Vic
930 real(wp) :: lambda_L, lambda_R, lambda_Cell, dT_dxi, grid_spacing
931 real(wp) :: Cp_L, Cp_R
932 real(wp) :: diffusivity_L, diffusivity_R, diffusivity_cell
933 real(wp) :: hmix_L, hmix_R, dh_dxi
934 integer :: x, y, z, i, n, eqn
935 integer, dimension(3) :: offsets
936
937 isc1 = irx; isc2 = iry; isc3 = irz
938
939
940# 352 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
941#if defined(MFC_OpenACC)
942# 352 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
943!$acc update device(isc1, isc2, isc3)
944# 352 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
945#elif defined(MFC_OpenMP)
946# 352 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
947!$omp target update to(isc1, isc2, isc3)
948# 352 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
949#endif
950
951 if (chemistry) then
952 ! Set offsets based on direction using array indexing
953 offsets = 0
954 offsets(idir) = 1
955 ! Model 1: Mixture-Average Transport
956 if (chem_params%transport_model == 1) then
957 ! Note: Added 'i' and 'eqn' to private list.
958
959# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
960
961# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
962#if defined(MFC_OpenACC)
963# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
964!$acc parallel loop collapse(3) gang vector default(present) &
965# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
966!$acc& private(x, y, z, i, eqn, 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, Mass_Diffu_Energy, MW_L, MW_R, MW_cell, T_L, T_R, P_L, P_R, rho_L, rho_R, rho_cell, rho_Vic, lambda_L, lambda_R, lambda_Cell, dT_dxi, grid_spacing) &
967# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
968!$acc& copyin(offsets)
969# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
970#elif defined(MFC_OpenMP)
971# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
972
973# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
974
975# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
976
977# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
978!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
979# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
980!$omp& private(x, y, z, i, eqn, 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, Mass_Diffu_Energy, MW_L, MW_R, MW_cell, T_L, T_R, P_L, P_R, rho_L, rho_R, rho_cell, rho_Vic, lambda_L, lambda_R, lambda_Cell, dT_dxi, grid_spacing) &
981# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
982!$omp& map(to:offsets)
983# 361 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
984#endif
985# 366 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
986 do z = isc3%beg, isc3%end
987 do y = isc2%beg, isc2%end
988 do x = isc1%beg, isc1%end
989 ! Calculate grid spacing using direction-based indexing
990 select case (idir)
991 case (1)
992 grid_spacing = x_cc(x + 1) - x_cc(x)
993 case (2)
994 grid_spacing = y_cc(y + 1) - y_cc(y)
995 case (3)
996 grid_spacing = z_cc(z + 1) - z_cc(z)
997 end select
998
999 ! Extract species mass fractions
1000
1001# 380 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1002#if defined(MFC_OpenACC)
1003# 380 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1004!$acc loop seq
1005# 380 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1006#elif defined(MFC_OpenMP)
1007# 380 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1008
1009# 380 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1010#endif
1011 do i = eqn_idx%species%beg, eqn_idx%species%end
1012 ys_l(i - eqn_idx%species%beg + 1) = q_prim_qp(i)%sf(x, y, z)
1013 ys_r(i - eqn_idx%species%beg + 1) = q_prim_qp(i)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1014 ys_cell(i - eqn_idx%species%beg + 1) = 0.5_wp*(ys_l(i - eqn_idx%species%beg + 1) + ys_r(i &
1015 & - eqn_idx%species%beg + 1))
1016 end do
1017
1018 ! Calculate molecular weights and mole fractions
1019 call get_mixture_molecular_weight(ys_l, mw_l)
1020 call get_mixture_molecular_weight(ys_r, mw_r)
1021 mw_cell = 0.5_wp*(mw_l + mw_r)
1022
1023 call get_mole_fractions(mw_l, ys_l, xs_l)
1024 call get_mole_fractions(mw_r, ys_r, xs_r)
1025
1026 p_l = q_prim_qp(eqn_idx%E)%sf(x, y, z)
1027 p_r = q_prim_qp(eqn_idx%E)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1028
1029 rho_l = q_prim_qp(1)%sf(x, y, z)
1030 rho_r = q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1031
1032 t_l = q_t_sf%sf(x, y, z)
1033 t_r = q_t_sf%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1034
1035 rho_cell = 0.5_wp*(rho_l + rho_r)
1036 dt_dxi = (t_r - t_l)/grid_spacing
1037
1038 ! Get transport properties
1039 call get_species_mass_diffusivities_mixavg(p_l, t_l, ys_l, mass_diffusivities_mixavg1)
1040 call get_species_mass_diffusivities_mixavg(p_r, t_r, ys_r, mass_diffusivities_mixavg2)
1041
1042 call get_mixture_thermal_conductivity_mixavg(t_l, ys_l, lambda_l)
1043 call get_mixture_thermal_conductivity_mixavg(t_r, ys_r, lambda_r)
1044
1045 call get_species_enthalpies_rt(t_l, h_l)
1046 call get_species_enthalpies_rt(t_r, h_r)
1047
1048 ! Calculate species properties and gradients
1049
1050# 419 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1051#if defined(MFC_OpenACC)
1052# 419 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1053!$acc loop seq
1054# 419 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1055#elif defined(MFC_OpenMP)
1056# 419 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1057
1058# 419 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1059#endif
1060 do i = eqn_idx%species%beg, eqn_idx%species%end
1061 h_l(i - eqn_idx%species%beg + 1) = h_l(i - eqn_idx%species%beg + 1) &
1062 & *gas_constant*t_l/molecular_weights(i - eqn_idx%species%beg + 1)
1063 h_r(i - eqn_idx%species%beg + 1) = h_r(i - eqn_idx%species%beg + 1) &
1064 & *gas_constant*t_r/molecular_weights(i - eqn_idx%species%beg + 1)
1065 xs_cell(i - eqn_idx%species%beg + 1) = 0.5_wp*(xs_l(i - eqn_idx%species%beg + 1) + xs_r(i &
1066 & - eqn_idx%species%beg + 1))
1067 h_k(i - eqn_idx%species%beg + 1) = 0.5_wp*(h_l(i - eqn_idx%species%beg + 1) + h_r(i &
1068 & - eqn_idx%species%beg + 1))
1069 dxk_dxi(i - eqn_idx%species%beg + 1) = (xs_r(i - eqn_idx%species%beg + 1) - xs_l(i &
1070 & - eqn_idx%species%beg + 1))/grid_spacing
1071 end do
1072
1073 ! Calculate mixture-averaged diffusivities
1074
1075# 434 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1076#if defined(MFC_OpenACC)
1077# 434 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1078!$acc loop seq
1079# 434 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1080#elif defined(MFC_OpenMP)
1081# 434 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1082
1083# 434 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1084#endif
1085 do i = eqn_idx%species%beg, eqn_idx%species%end
1086 mass_diffusivities_mixavg_cell(i - eqn_idx%species%beg + 1) = (mass_diffusivities_mixavg2(i &
1087 & - eqn_idx%species%beg + 1) + mass_diffusivities_mixavg1(i &
1088 & - eqn_idx%species%beg + 1))/2.0_wp
1089 end do
1090
1091 lambda_cell = 0.5_wp*(lambda_r + lambda_l)
1092
1093 ! Calculate mass diffusion fluxes
1094 rho_vic = 0.0_wp
1095 mass_diffu_energy = 0.0_wp
1096
1097
1098# 447 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1099#if defined(MFC_OpenACC)
1100# 447 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1101!$acc loop seq
1102# 447 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1103#elif defined(MFC_OpenMP)
1104# 447 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1105
1106# 447 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1107#endif
1108 do eqn = eqn_idx%species%beg, eqn_idx%species%end
1109 mass_diffu_flux(eqn - eqn_idx%species%beg + 1) = rho_cell*mass_diffusivities_mixavg_cell(eqn &
1110 & - eqn_idx%species%beg + 1)*molecular_weights(eqn - eqn_idx%species%beg + 1) &
1111 & /mw_cell*dxk_dxi(eqn - eqn_idx%species%beg + 1)
1112 rho_vic = rho_vic + mass_diffu_flux(eqn - eqn_idx%species%beg + 1)
1113 mass_diffu_energy = mass_diffu_energy + h_k(eqn - eqn_idx%species%beg + 1)*mass_diffu_flux(eqn &
1114 & - eqn_idx%species%beg + 1)
1115 end do
1116
1117 ! Apply corrections for mass conservation
1118
1119# 458 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1120#if defined(MFC_OpenACC)
1121# 458 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1122!$acc loop seq
1123# 458 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1124#elif defined(MFC_OpenMP)
1125# 458 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1126
1127# 458 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1128#endif
1129 do eqn = eqn_idx%species%beg, eqn_idx%species%end
1130 mass_diffu_energy = mass_diffu_energy - h_k(eqn - eqn_idx%species%beg + 1)*ys_cell(eqn &
1131 & - eqn_idx%species%beg + 1)*rho_vic
1132 mass_diffu_flux(eqn - eqn_idx%species%beg + 1) = mass_diffu_flux(eqn - eqn_idx%species%beg + 1) &
1133 & - rho_vic*ys_cell(eqn - eqn_idx%species%beg + 1)
1134 end do
1135
1136 ! Add thermal conduction contribution
1137 mass_diffu_energy = lambda_cell*dt_dxi + mass_diffu_energy
1138
1139 ! Update flux arrays
1140 flux_src_vf(eqn_idx%E)%sf(x, y, z) = flux_src_vf(eqn_idx%E)%sf(x, y, z) - mass_diffu_energy
1141
1142
1143# 472 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1144#if defined(MFC_OpenACC)
1145# 472 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1146!$acc loop seq
1147# 472 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1148#elif defined(MFC_OpenMP)
1149# 472 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1150
1151# 472 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1152#endif
1153 do eqn = eqn_idx%species%beg, eqn_idx%species%end
1154 flux_src_vf(eqn)%sf(x, y, z) = flux_src_vf(eqn)%sf(x, y, &
1155 & z) - mass_diffu_flux(eqn - eqn_idx%species%beg + 1)
1156 end do
1157 end do
1158 end do
1159 end do
1160
1161# 480 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1162#if defined(MFC_OpenACC)
1163# 480 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1164!$acc end parallel loop
1165# 480 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1166#elif defined(MFC_OpenMP)
1167# 480 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1168
1169# 480 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1170!$omp end target teams loop
1171# 480 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1172#endif
1173
1174 ! Model 2: Unity Lewis Number
1175 else if (chem_params%transport_model == 2) then
1176 ! Note: Added ALL scalars and 'i'/'eqn' to private list to prevent race conditions.
1177
1178# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1179
1180# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1181#if defined(MFC_OpenACC)
1182# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1183!$acc parallel loop collapse(3) gang vector default(present) &
1184# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1185!$acc& private(x, y, z, i, eqn, Ys_L, Ys_R, Ys_cell, dYk_dxi, Mass_Diffu_Flux, grid_spacing, MW_L, MW_R, MW_cell, P_L, P_R, rho_L, rho_R, rho_cell, T_L, T_R, Cp_L, Cp_R, hmix_L, hmix_R, dh_dxi, lambda_L, lambda_R, lambda_Cell, diffusivity_L, diffusivity_R, diffusivity_cell, Mass_Diffu_Energy) &
1186# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1187!$acc& copyin(offsets)
1188# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1189#elif defined(MFC_OpenMP)
1190# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1191
1192# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1193
1194# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1195
1196# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1197!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
1198# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1199!$omp& private(x, y, z, i, eqn, Ys_L, Ys_R, Ys_cell, dYk_dxi, Mass_Diffu_Flux, grid_spacing, MW_L, MW_R, MW_cell, P_L, P_R, rho_L, rho_R, rho_cell, T_L, T_R, Cp_L, Cp_R, hmix_L, hmix_R, dh_dxi, lambda_L, lambda_R, lambda_Cell, diffusivity_L, diffusivity_R, diffusivity_cell, Mass_Diffu_Energy) &
1200# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1201!$omp& map(to:offsets)
1202# 485 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1203#endif
1204# 489 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1205 do z = isc3%beg, isc3%end
1206 do y = isc2%beg, isc2%end
1207 do x = isc1%beg, isc1%end
1208 ! Calculate grid spacing using direction-based indexing
1209 select case (idir)
1210 case (1)
1211 grid_spacing = x_cc(x + 1) - x_cc(x)
1212 case (2)
1213 grid_spacing = y_cc(y + 1) - y_cc(y)
1214 case (3)
1215 grid_spacing = z_cc(z + 1) - z_cc(z)
1216 end select
1217
1218 ! Extract species mass fractions
1219
1220# 503 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1221#if defined(MFC_OpenACC)
1222# 503 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1223!$acc loop seq
1224# 503 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1225#elif defined(MFC_OpenMP)
1226# 503 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1227
1228# 503 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1229#endif
1230 do i = eqn_idx%species%beg, eqn_idx%species%end
1231 ys_l(i - eqn_idx%species%beg + 1) = q_prim_qp(i)%sf(x, y, z)
1232 ys_r(i - eqn_idx%species%beg + 1) = q_prim_qp(i)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1233 ys_cell(i - eqn_idx%species%beg + 1) = 0.5_wp*(ys_l(i - eqn_idx%species%beg + 1) + ys_r(i &
1234 & - eqn_idx%species%beg + 1))
1235 end do
1236
1237 ! Calculate molecular weights and mole fractions
1238 call get_mixture_molecular_weight(ys_l, mw_l)
1239 call get_mixture_molecular_weight(ys_r, mw_r)
1240 mw_cell = 0.5_wp*(mw_l + mw_r)
1241
1242 p_l = q_prim_qp(eqn_idx%E)%sf(x, y, z)
1243 p_r = q_prim_qp(eqn_idx%E)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1244
1245 rho_l = q_prim_qp(1)%sf(x, y, z)
1246 rho_r = q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1247
1248 t_l = q_t_sf%sf(x, y, z)
1249 t_r = q_t_sf%sf(x + offsets(1), y + offsets(2), z + offsets(3))
1250
1251 rho_cell = 0.5_wp*(rho_l + rho_r)
1252
1253 call get_mixture_specific_heat_cp_mass(t_l, ys_l, cp_l)
1254 call get_mixture_specific_heat_cp_mass(t_r, ys_r, cp_r)
1255 call get_mixture_enthalpy_mass(t_l, ys_l, hmix_l)
1256 call get_mixture_enthalpy_mass(t_r, ys_r, hmix_r)
1257 dh_dxi = (hmix_r - hmix_l)/grid_spacing
1258
1259 ! Get transport properties
1260 call get_mixture_thermal_conductivity_mixavg(t_l, ys_l, lambda_l)
1261 call get_mixture_thermal_conductivity_mixavg(t_r, ys_r, lambda_r)
1262
1263 ! Calculate species properties and gradients
1264
1265# 538 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1266#if defined(MFC_OpenACC)
1267# 538 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1268!$acc loop seq
1269# 538 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1270#elif defined(MFC_OpenMP)
1271# 538 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1272
1273# 538 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1274#endif
1275 do i = eqn_idx%species%beg, eqn_idx%species%end
1276 dyk_dxi(i - eqn_idx%species%beg + 1) = (ys_r(i - eqn_idx%species%beg + 1) - ys_l(i &
1277 & - eqn_idx%species%beg + 1))/grid_spacing
1278 end do
1279
1280 ! Calculate mixture-averaged diffusivities
1281 diffusivity_l = lambda_l/rho_l/cp_l
1282 diffusivity_r = lambda_r/rho_r/cp_r
1283
1284 lambda_cell = 0.5_wp*(lambda_r + lambda_l)
1285 diffusivity_cell = 0.5_wp*(diffusivity_r + diffusivity_l)
1286
1287 ! Calculate mass diffusion fluxes
1288 mass_diffu_energy = 0.0_wp
1289
1290
1291# 554 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1292#if defined(MFC_OpenACC)
1293# 554 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1294!$acc loop seq
1295# 554 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1296#elif defined(MFC_OpenMP)
1297# 554 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1298
1299# 554 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1300#endif
1301 do eqn = eqn_idx%species%beg, eqn_idx%species%end
1302 mass_diffu_flux(eqn - eqn_idx%species%beg + 1) = rho_cell*diffusivity_cell*dyk_dxi(eqn &
1303 & - eqn_idx%species%beg + 1)
1304 end do
1305 mass_diffu_energy = rho_cell*diffusivity_cell*dh_dxi
1306
1307 ! Update flux arrays
1308 flux_src_vf(eqn_idx%E)%sf(x, y, z) = flux_src_vf(eqn_idx%E)%sf(x, y, z) - mass_diffu_energy
1309
1310
1311# 564 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1312#if defined(MFC_OpenACC)
1313# 564 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1314!$acc loop seq
1315# 564 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1316#elif defined(MFC_OpenMP)
1317# 564 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1318
1319# 564 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1320#endif
1321 do eqn = eqn_idx%species%beg, eqn_idx%species%end
1322 flux_src_vf(eqn)%sf(x, y, z) = flux_src_vf(eqn)%sf(x, y, &
1323 & z) - mass_diffu_flux(eqn - eqn_idx%species%beg + 1)
1324 end do
1325 end do
1326 end do
1327 end do
1328
1329# 572 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1330#if defined(MFC_OpenACC)
1331# 572 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1332!$acc end parallel loop
1333# 572 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1334#elif defined(MFC_OpenMP)
1335# 572 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1336
1337# 572 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1338!$omp end target teams loop
1339# 572 "/home/runner/work/MFC/MFC/src/common/m_chemistry.fpp"
1340#endif
1341 end if
1342 end if
1343
1344 end subroutine s_compute_chemistry_diffusion_flux
1345
1346end module m_chemistry
m_chemistry
Multi-species chemistry interface for thermodynamic properties, reaction rates, and transport coeffic...
Definition m_chemistry.fpp.f90:335
m_chemistry::s_compute_chemistry_diffusion_flux
subroutine s_compute_chemistry_diffusion_flux(idir, q_prim_qp, flux_src_vf, irx, iry, irz, q_t_sf)
Compute species mass diffusion fluxes at cell interfaces using mixture-averaged diffusivities.
Definition m_chemistry.fpp.f90:914
m_chemistry::s_compute_q_t_sf
subroutine s_compute_q_t_sf(q_t_sf, q_cons_vf, bounds)
Initialize the temperature field from conservative variables by inverting the energy equation.
Definition m_chemistry.fpp.f90:426
m_chemistry::isc2
type(int_bounds_info) isc2
Definition m_chemistry.fpp.f90:347
m_chemistry::s_compute_chemistry_reaction_flux
subroutine s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp, q_t_sf, q_prim_qp, bounds)
Add chemical reaction source terms to the species transport RHS using net production rates.
Definition m_chemistry.fpp.f90:487
m_chemistry::offsets
integer, dimension(3) offsets
Definition m_chemistry.fpp.f90:359
m_chemistry::isc1
type(int_bounds_info) isc1
Definition m_chemistry.fpp.f90:347
m_chemistry::s_chemistry_reaction_substep
subroutine s_chemistry_reaction_substep(q_cons_vf, q_t_sf, dtime, bounds)
Operator-split integration of the reaction source with an alpha-QSS (Mott quasi-steady-state) reactor...
Definition m_chemistry.fpp.f90:590
m_chemistry::isc3
type(int_bounds_info) isc3
Definition m_chemistry.fpp.f90:347
m_chemistry::compute_viscosity_and_inversion
subroutine compute_viscosity_and_inversion(t_l, ys_l, t_r, ys_r, re_l, re_r)
Compute mixture viscosities for left and right states and invert them for use as reciprocal Reynolds ...
Definition m_chemistry.fpp.f90:376
m_chemistry::s_compute_t_from_primitives
subroutine s_compute_t_from_primitives(q_t_sf, q_prim_vf, bounds)
Compute the temperature field from primitive variables using the ideal gas law and mixture molecular ...
Definition m_chemistry.fpp.f90:462
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:334
m_global_parameters::y_cc
real(wp), dimension(:), allocatable, target y_cc
Definition m_global_parameters.fpp.f90:400
m_global_parameters::z_cc
real(wp), dimension(:), allocatable, target z_cc
Definition m_global_parameters.fpp.f90:400
m_global_parameters::x_cc
real(wp), dimension(:), allocatable, target x_cc
Definition m_global_parameters.fpp.f90:400