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