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