MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_qbmm.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
2!>
3!! @file
4!! @brief Contains module m_qbmm
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
7! This file exists so that Fypp can be run without generating case.fpp files for
8! each target. This is useful when generating documentation, for example. This
9! should also let MFC be built with CMake directly, without invoking mfc.sh.
10
11! For pre-process.
12# 9 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
13
14! For moving immersed boundaries in simulation
15# 14 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
16# 6 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp" 2
17# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
18# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
19# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
20# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39
40# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43! New line at end of file is required for FYPP
44# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
45# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
46# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
47# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52
53# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56
57# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60
61# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64
65# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66
67# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68
69# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
70! New line at end of file is required for FYPP
71# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
72
73# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
74# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
76# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
78
79# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
80
81# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
82
83# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
84
85# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86
87# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88
89# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 76 "/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# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145
146# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147
148# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149! New line at end of file is required for FYPP
150# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
151# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
152# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
153# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
156# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
158
159# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
160# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162
163# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
164
165# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166
167# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168
169# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170
171# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172
173# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
174
175# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
176! New line at end of file is required for FYPP
177# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
178
179# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
180
181# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
182
183# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
184
185# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
186
187# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
188
189# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
190
191# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
192
193# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234! New line at end of file is required for FYPP
235# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
236
237! GPU parallel region (scalar reductions, maxval/minval)
238# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
239
240! GPU parallel loop over threads (most common GPU macro)
241# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
242
243! Required closing for GPU_PARALLEL_LOOP
244# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
245
246! Mark routine for device compilation
247# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
248
249! Declare device-resident data
250# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
251
252! Inner loop within a GPU parallel region
253# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
254
255! Scoped GPU data region
256# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
257
258! Host code with device pointers (for MPI with GPU buffers)
259# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
260
261! Allocate device memory (unscoped)
262# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
263
264! Free device memory
265# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
266
267! Atomic operation on device
268# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
269
270! End atomic capture block
271# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
272
273! Copy data between host and device
274# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
275
276! Synchronization barrier
277# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
278
279! Import GPU library module (openacc or omp_lib)
280# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
281
282! Emit code only for AMD compiler
283# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
284
285! Emit code for non-Cray compilers
286# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
287
288! Emit code only for Cray compiler
289# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
290
291! Emit code for non-NVIDIA compilers
292# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
293
294# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
295# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
296! New line at end of file is required for FYPP
297# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
298
299# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
302! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
303! example see misc/nvidia_uvm/bind.sh. NVIDIA unified memory page placement hint
304# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
305
306! Allocate and create GPU device memory
307# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
308
309! Free GPU device memory and deallocate
310# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
311
312! Cray-specific GPU pointer setup for vector fields
313# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
314
315! Cray-specific GPU pointer setup for scalar fields
316# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318! Cray-specific GPU pointer setup for acoustic source spatials
319# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
320
321# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
322
323# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
324! New line at end of file is required for FYPP
325# 7 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp" 2
326
327!> @brief Quadrature-based moment methods (QBMM) for polydisperse bubble moment inversion and transport
328module m_qbmm
329
332 use m_mpi_proxy
335 use m_helper
336
337 implicit none
338
340
341 real(wp), allocatable, dimension(:,:,:,:,:) :: momrhs
342
343# 23 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
344#if defined(MFC_OpenACC)
345# 23 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
346!$acc declare create(momrhs)
347# 23 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
348#elif defined(MFC_OpenMP)
349# 23 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
350!$omp declare target (momrhs)
351# 23 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
352#endif
353
354# 28 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
355 integer :: nterms
356
357# 29 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
358#if defined(MFC_OpenACC)
359# 29 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
360!$acc declare create(nterms)
361# 29 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
362#elif defined(MFC_OpenMP)
363# 29 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
364!$omp declare target (nterms)
365# 29 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
366#endif
367# 31 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
368
370
371# 33 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
372#if defined(MFC_OpenACC)
373# 33 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
374!$acc declare create(is1_qbmm, is2_qbmm, is3_qbmm)
375# 33 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
376#elif defined(MFC_OpenMP)
377# 33 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
378!$omp declare target (is1_qbmm, is2_qbmm, is3_qbmm)
379# 33 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
380#endif
381
382 integer, allocatable, dimension(:) :: bubrs_qbmm
383 integer, allocatable, dimension(:,:) :: bubmoms
384
385# 37 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
386#if defined(MFC_OpenACC)
387# 37 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
388!$acc declare create(bubrs_qbmm, bubmoms)
389# 37 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
390#elif defined(MFC_OpenMP)
391# 37 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
392!$omp declare target (bubrs_qbmm, bubmoms)
393# 37 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
394#endif
395
396contains
397
398 !> Initialize the QBMM module
399 impure subroutine s_initialize_qbmm_module
400
401 integer :: i1, i2, q, i, j
402
403# 47 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
404 if (bubble_model == 2) then
405 ! Keller-Miksis without viscosity/surface tension
406 nterms = 32
407 else if (bubble_model == 3) then
408 ! Rayleigh-Plesset with viscosity/surface tension
409 nterms = 7
410 end if
411
412
413# 55 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
414#if defined(MFC_OpenACC)
415# 55 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
416!$acc enter data copyin(nterms)
417# 55 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
418#elif defined(MFC_OpenMP)
419# 55 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
420!$omp target enter data map(to:nterms)
421# 55 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
422#endif
423
424# 56 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
425#if defined(MFC_OpenACC)
426# 56 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
427!$acc update device(nterms)
428# 56 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
429#elif defined(MFC_OpenMP)
430# 56 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
431!$omp target update to(nterms)
432# 56 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
433#endif
434# 58 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
435
436#ifdef MFC_DEBUG
437# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
438 block
439# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
440 use iso_fortran_env, only: output_unit
441# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
442
443# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
444 print *, 'm_qbmm.fpp:59: ', '@:ALLOCATE(momrhs(1:3, 0:2, 0:2, 1:nterms, 1:nb))'
445# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
446
447# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
448 call flush (output_unit)
449# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
450 end block
451# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
452#endif
453# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
454 allocate (momrhs(1:3, 0:2, 0:2, 1:nterms, 1:nb))
455# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
456
457# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
458
459# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
460#if defined(MFC_OpenACC)
461# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
462!$acc enter data create(momrhs)
463# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
464#elif defined(MFC_OpenMP)
465# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
466!$omp target enter data map(always,alloc:momrhs)
467# 59 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
468#endif
469 momrhs = 0._wp
470
471 ! Assigns the required RHS moments for moment transport equations The rhs%(:,3) is only to be used for R0 quadrature, not
472 ! for computing X/Y indices Accounts for different governing equations in polytropic and non-polytropic models
473 if (.not. polytropic) then
474 do q = 1, nb
475 do i1 = 0, 2; do i2 = 0, 2
476 if ((i1 + i2) <= 2) then
477 if (bubble_model == 3) then
478 momrhs(1, i1, i2, 1, q) = -1._wp + i1
479 momrhs(2, i1, i2, 1, q) = -1._wp + i2
480 momrhs(3, i1, i2, 1, q) = 0._wp
481
482 momrhs(1, i1, i2, 2, q) = -1._wp + i1
483 momrhs(2, i1, i2, 2, q) = 1._wp + i2
484 momrhs(3, i1, i2, 2, q) = 0._wp
485
486 momrhs(1, i1, i2, 3, q) = -1._wp + i1
487 momrhs(2, i1, i2, 3, q) = -1._wp + i2
488 momrhs(3, i1, i2, 3, q) = 0._wp
489
490 momrhs(1, i1, i2, 4, q) = -1._wp + i1
491 momrhs(2, i1, i2, 4, q) = 1._wp + i2
492 momrhs(3, i1, i2, 4, q) = 0._wp
493
494 if (.not. f_is_default(re_inv)) then
495 ! add viscosity
496 momrhs(1, i1, i2, 5, q) = -2._wp + i1
497 momrhs(2, i1, i2, 5, q) = i2
498 momrhs(3, i1, i2, 5, q) = 0._wp
499 end if
500
501 if (.not. f_is_default(web)) then
502 ! add surface tension
503 momrhs(1, i1, i2, 6, q) = -2._wp + i1
504 momrhs(2, i1, i2, 6, q) = -1._wp + i2
505 momrhs(3, i1, i2, 6, q) = 0._wp
506 end if
507
508 momrhs(1, i1, i2, 7, q) = -1._wp + i1
509 momrhs(2, i1, i2, 7, q) = -1._wp + i2
510 momrhs(3, i1, i2, 7, q) = 0._wp
511 else if (bubble_model == 2) then
512 ! KM with approximation of 1/(1-V/C) = 1+V/C
513 momrhs(1, i1, i2, 1, q) = -1._wp + i1
514 momrhs(2, i1, i2, 1, q) = 1._wp + i2
515 momrhs(3, i1, i2, 1, q) = 0._wp
516
517 momrhs(1, i1, i2, 2, q) = -1._wp + i1
518 momrhs(2, i1, i2, 2, q) = 2._wp + i2
519 momrhs(3, i1, i2, 2, q) = 0._wp
520
521 momrhs(1, i1, i2, 3, q) = -1._wp + i1
522 momrhs(2, i1, i2, 3, q) = 3._wp + i2
523 momrhs(3, i1, i2, 3, q) = 0._wp
524
525 momrhs(1, i1, i2, 4, q) = -1._wp + i1
526 momrhs(2, i1, i2, 4, q) = -1._wp + i2
527 momrhs(3, i1, i2, 4, q) = 0._wp
528
529 momrhs(1, i1, i2, 5, q) = -1._wp + i1
530 momrhs(2, i1, i2, 5, q) = i2
531 momrhs(3, i1, i2, 5, q) = 0._wp
532
533 momrhs(1, i1, i2, 6, q) = -1._wp + i1
534 momrhs(2, i1, i2, 6, q) = 1._wp + i2
535 momrhs(3, i1, i2, 6, q) = 0._wp
536
537 momrhs(1, i1, i2, 7, q) = -1._wp + i1
538 momrhs(2, i1, i2, 7, q) = -1._wp + i2
539 momrhs(3, i1, i2, 7, q) = 0._wp
540
541 momrhs(1, i1, i2, 8, q) = -1._wp + i1
542 momrhs(2, i1, i2, 8, q) = i2
543 momrhs(3, i1, i2, 8, q) = 0._wp
544
545 momrhs(1, i1, i2, 9, q) = -1._wp + i1
546 momrhs(2, i1, i2, 9, q) = 1._wp + i2
547 momrhs(3, i1, i2, 9, q) = 0._wp
548
549 momrhs(1, i1, i2, 10, q) = -1._wp + i1
550 momrhs(2, i1, i2, 10, q) = i2
551 momrhs(3, i1, i2, 10, q) = 0._wp
552
553 momrhs(1, i1, i2, 11, q) = -1._wp + i1
554 momrhs(2, i1, i2, 11, q) = 1._wp + i2
555 momrhs(3, i1, i2, 11, q) = 0._wp
556
557 momrhs(1, i1, i2, 12, q) = -1._wp + i1
558 momrhs(2, i1, i2, 12, q) = 1._wp + i2
559 momrhs(3, i1, i2, 12, q) = 0._wp
560
561 momrhs(1, i1, i2, 13, q) = -1._wp + i1
562 momrhs(2, i1, i2, 13, q) = -1._wp + i2
563 momrhs(3, i1, i2, 13, q) = 0._wp
564
565 momrhs(1, i1, i2, 14, q) = -1._wp + i1
566 momrhs(2, i1, i2, 14, q) = i2
567 momrhs(3, i1, i2, 14, q) = 0._wp
568
569 momrhs(1, i1, i2, 15, q) = -1._wp + i1
570 momrhs(2, i1, i2, 15, q) = 1._wp + i2
571 momrhs(3, i1, i2, 15, q) = 0._wp
572
573 momrhs(1, i1, i2, 16, q) = -2._wp + i1
574 momrhs(2, i1, i2, 16, q) = i2
575 momrhs(3, i1, i2, 16, q) = 0._wp
576
577 momrhs(1, i1, i2, 17, q) = -2._wp + i1
578 momrhs(2, i1, i2, 17, q) = -1._wp + i2
579 momrhs(3, i1, i2, 17, q) = 0._wp
580
581 momrhs(1, i1, i2, 18, q) = -2._wp + i1
582 momrhs(2, i1, i2, 18, q) = 1._wp + i2
583 momrhs(3, i1, i2, 18, q) = 0._wp
584
585 momrhs(1, i1, i2, 19, q) = -2._wp + i1
586 momrhs(2, i1, i2, 19, q) = 2._wp + i2
587 momrhs(3, i1, i2, 19, q) = 0._wp
588
589 momrhs(1, i1, i2, 20, q) = -2._wp + i1
590 momrhs(2, i1, i2, 20, q) = -1._wp + i2
591 momrhs(3, i1, i2, 20, q) = 0._wp
592
593 momrhs(1, i1, i2, 21, q) = -2._wp + i1
594 momrhs(2, i1, i2, 21, q) = i2
595 momrhs(3, i1, i2, 21, q) = 0._wp
596
597 momrhs(1, i1, i2, 22, q) = -2._wp + i1
598 momrhs(2, i1, i2, 22, q) = -1._wp + i2
599 momrhs(3, i1, i2, 22, q) = 0._wp
600
601 momrhs(1, i1, i2, 23, q) = -2._wp + i1
602 momrhs(2, i1, i2, 23, q) = i2
603 momrhs(3, i1, i2, 23, q) = 0._wp
604
605 momrhs(1, i1, i2, 24, q) = -3._wp + i1
606 momrhs(2, i1, i2, 24, q) = i2
607 momrhs(3, i1, i2, 24, q) = 0._wp
608
609 momrhs(1, i1, i2, 25, q) = -3._wp + i1
610 momrhs(2, i1, i2, 25, q) = -1._wp + i2
611 momrhs(3, i1, i2, 25, q) = 0._wp
612
613 momrhs(1, i1, i2, 26, q) = -2._wp + i1
614 momrhs(2, i1, i2, 26, q) = i2
615 momrhs(3, i1, i2, 26, q) = 0._wp
616
617 momrhs(1, i1, i2, 27, q) = -1._wp + i1
618 momrhs(2, i1, i2, 27, q) = -1._wp + i2
619 momrhs(3, i1, i2, 27, q) = 0._wp
620
621 momrhs(1, i1, i2, 28, q) = -1._wp + i1
622 momrhs(2, i1, i2, 28, q) = i2
623 momrhs(3, i1, i2, 28, q) = 0._wp
624
625 momrhs(1, i1, i2, 29, q) = -2._wp + i1
626 momrhs(2, i1, i2, 29, q) = i2
627 momrhs(3, i1, i2, 29, q) = 0._wp
628
629 momrhs(1, i1, i2, 30, q) = -1._wp + i1
630 momrhs(2, i1, i2, 30, q) = -1._wp + i2
631 momrhs(3, i1, i2, 30, q) = 0._wp
632
633 momrhs(1, i1, i2, 31, q) = -1._wp + i1
634 momrhs(2, i1, i2, 31, q) = i2
635 momrhs(3, i1, i2, 31, q) = 0._wp
636
637 momrhs(1, i1, i2, 32, q) = -2._wp + i1
638 momrhs(2, i1, i2, 32, q) = i2
639 momrhs(3, i1, i2, 32, q) = 0._wp
640 end if
641 end if
642 end do; end do
643 end do
644 else
645 do q = 1, nb
646 do i1 = 0, 2; do i2 = 0, 2
647 if ((i1 + i2) <= 2) then
648 if (bubble_model == 3) then
649 momrhs(1, i1, i2, 1, q) = -1._wp + i1
650 momrhs(2, i1, i2, 1, q) = -1._wp + i2
651 momrhs(3, i1, i2, 1, q) = 0._wp
652
653 momrhs(1, i1, i2, 2, q) = -1._wp + i1
654 momrhs(2, i1, i2, 2, q) = 1._wp + i2
655 momrhs(3, i1, i2, 2, q) = 0._wp
656
657 momrhs(1, i1, i2, 3, q) = -1._wp + i1 - 3._wp*gam
658 momrhs(2, i1, i2, 3, q) = -1._wp + i2
659 momrhs(3, i1, i2, 3, q) = 3._wp*gam
660
661 momrhs(1, i1, i2, 4, q) = -1._wp + i1
662 momrhs(2, i1, i2, 4, q) = 1._wp + i2
663 momrhs(3, i1, i2, 4, q) = 0._wp
664
665 if (.not. f_is_default(re_inv)) then
666 ! add viscosity
667 momrhs(1, i1, i2, 5, q) = -2._wp + i1
668 momrhs(2, i1, i2, 5, q) = i2
669 momrhs(3, i1, i2, 5, q) = 0._wp
670 end if
671
672 if (.not. f_is_default(web)) then
673 ! add surface tension
674 momrhs(1, i1, i2, 6, q) = -2._wp + i1
675 momrhs(2, i1, i2, 6, q) = -1._wp + i2
676 momrhs(3, i1, i2, 6, q) = 0._wp
677 end if
678
679 momrhs(1, i1, i2, 7, q) = -1._wp + i1
680 momrhs(2, i1, i2, 7, q) = -1._wp + i2
681 momrhs(3, i1, i2, 7, q) = 0._wp
682 else if (bubble_model == 2) then
683 ! KM with approximation of 1/(1-V/C) = 1+V/C
684 momrhs(1, i1, i2, 1, q) = -1._wp + i1
685 momrhs(2, i1, i2, 1, q) = 1._wp + i2
686 momrhs(3, i1, i2, 1, q) = 0._wp
687
688 momrhs(1, i1, i2, 2, q) = -1._wp + i1
689 momrhs(2, i1, i2, 2, q) = 2._wp + i2
690 momrhs(3, i1, i2, 2, q) = 0._wp
691
692 momrhs(1, i1, i2, 3, q) = -1._wp + i1
693 momrhs(2, i1, i2, 3, q) = 3._wp + i2
694 momrhs(3, i1, i2, 3, q) = 0._wp
695
696 momrhs(1, i1, i2, 4, q) = -1._wp + i1
697 momrhs(2, i1, i2, 4, q) = -1._wp + i2
698 momrhs(3, i1, i2, 4, q) = 0._wp
699
700 momrhs(1, i1, i2, 5, q) = -1._wp + i1
701 momrhs(2, i1, i2, 5, q) = i2
702 momrhs(3, i1, i2, 5, q) = 0._wp
703
704 momrhs(1, i1, i2, 6, q) = -1._wp + i1
705 momrhs(2, i1, i2, 6, q) = 1._wp + i2
706 momrhs(3, i1, i2, 6, q) = 0._wp
707
708 momrhs(1, i1, i2, 7, q) = -1._wp + i1 - 3._wp*gam
709 momrhs(2, i1, i2, 7, q) = -1._wp + i2
710 momrhs(3, i1, i2, 7, q) = 3._wp*gam
711
712 momrhs(1, i1, i2, 8, q) = -1._wp + i1 - 3._wp*gam
713 momrhs(2, i1, i2, 8, q) = i2
714 momrhs(3, i1, i2, 8, q) = 3._wp*gam
715
716 momrhs(1, i1, i2, 9, q) = -1._wp + i1 - 3._wp*gam
717 momrhs(2, i1, i2, 9, q) = 1._wp + i2
718 momrhs(3, i1, i2, 9, q) = 3._wp*gam
719
720 momrhs(1, i1, i2, 10, q) = -1._wp + i1 - 3._wp*gam
721 momrhs(2, i1, i2, 10, q) = i2
722 momrhs(3, i1, i2, 10, q) = 3._wp*gam
723
724 momrhs(1, i1, i2, 11, q) = -1._wp + i1 - 3._wp*gam
725 momrhs(2, i1, i2, 11, q) = 1._wp + i2
726 momrhs(3, i1, i2, 11, q) = 3._wp*gam
727
728 momrhs(1, i1, i2, 12, q) = -1._wp + i1
729 momrhs(2, i1, i2, 12, q) = 1._wp + i2
730 momrhs(3, i1, i2, 12, q) = 0._wp
731
732 momrhs(1, i1, i2, 13, q) = -1._wp + i1
733 momrhs(2, i1, i2, 13, q) = -1._wp + i2
734 momrhs(3, i1, i2, 13, q) = 0._wp
735
736 momrhs(1, i1, i2, 14, q) = -1._wp + i1
737 momrhs(2, i1, i2, 14, q) = i2
738 momrhs(3, i1, i2, 14, q) = 0._wp
739
740 momrhs(1, i1, i2, 15, q) = -1._wp + i1
741 momrhs(2, i1, i2, 15, q) = 1._wp + i2
742 momrhs(3, i1, i2, 15, q) = 0._wp
743
744 momrhs(1, i1, i2, 16, q) = -2._wp + i1
745 momrhs(2, i1, i2, 16, q) = i2
746 momrhs(3, i1, i2, 16, q) = 0._wp
747
748 momrhs(1, i1, i2, 17, q) = -2._wp + i1
749 momrhs(2, i1, i2, 17, q) = -1._wp + i2
750 momrhs(3, i1, i2, 17, q) = 0._wp
751
752 momrhs(1, i1, i2, 18, q) = -2._wp + i1
753 momrhs(2, i1, i2, 18, q) = 1._wp + i2
754 momrhs(3, i1, i2, 18, q) = 0._wp
755
756 momrhs(1, i1, i2, 19, q) = -2._wp + i1
757 momrhs(2, i1, i2, 19, q) = 2._wp + i2
758 momrhs(3, i1, i2, 19, q) = 0._wp
759
760 momrhs(1, i1, i2, 20, q) = -2._wp + i1
761 momrhs(2, i1, i2, 20, q) = -1._wp + i2
762 momrhs(3, i1, i2, 20, q) = 0._wp
763
764 momrhs(1, i1, i2, 21, q) = -2._wp + i1
765 momrhs(2, i1, i2, 21, q) = i2
766 momrhs(3, i1, i2, 21, q) = 0._wp
767
768 momrhs(1, i1, i2, 22, q) = -2._wp + i1 - 3._wp*gam
769 momrhs(2, i1, i2, 22, q) = -1._wp + i2
770 momrhs(3, i1, i2, 22, q) = 3._wp*gam
771
772 momrhs(1, i1, i2, 23, q) = -2._wp + i1 - 3._wp*gam
773 momrhs(2, i1, i2, 23, q) = i2
774 momrhs(3, i1, i2, 23, q) = 3._wp*gam
775
776 momrhs(1, i1, i2, 24, q) = -3._wp + i1
777 momrhs(2, i1, i2, 24, q) = i2
778 momrhs(3, i1, i2, 24, q) = 0._wp
779
780 momrhs(1, i1, i2, 25, q) = -3._wp + i1
781 momrhs(2, i1, i2, 25, q) = -1._wp + i2
782 momrhs(3, i1, i2, 25, q) = 0._wp
783
784 momrhs(1, i1, i2, 26, q) = -2._wp + i1 - 3._wp*gam
785 momrhs(2, i1, i2, 26, q) = i2
786 momrhs(3, i1, i2, 26, q) = 3._wp*gam
787 end if
788 end if
789 end do; end do
790 end do
791 end if
792
793
794# 384 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
795#if defined(MFC_OpenACC)
796# 384 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
797!$acc update device(momrhs)
798# 384 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
799#elif defined(MFC_OpenMP)
800# 384 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
801!$omp target update to(momrhs)
802# 384 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
803#endif
804
805#ifdef MFC_DEBUG
806# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
807 block
808# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
809 use iso_fortran_env, only: output_unit
810# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
811
812# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
813 print *, 'm_qbmm.fpp:386: ', '@:ALLOCATE(bubrs_qbmm(1:nb))'
814# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
815
816# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
817 call flush (output_unit)
818# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
819 end block
820# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
821#endif
822# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
823 allocate (bubrs_qbmm(1:nb))
824# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
825
826# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
827
828# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
829#if defined(MFC_OpenACC)
830# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
831!$acc enter data create(bubrs_qbmm)
832# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
833#elif defined(MFC_OpenMP)
834# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
835!$omp target enter data map(always,alloc:bubrs_qbmm)
836# 386 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
837#endif
838#ifdef MFC_DEBUG
839# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
840 block
841# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
842 use iso_fortran_env, only: output_unit
843# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
844
845# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
846 print *, 'm_qbmm.fpp:387: ', '@:ALLOCATE(bubmoms(1:nb, 1:nmom))'
847# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
848
849# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
850 call flush (output_unit)
851# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
852 end block
853# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
854#endif
855# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
856 allocate (bubmoms(1:nb, 1:nmom))
857# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
858
859# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
860
861# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
862#if defined(MFC_OpenACC)
863# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
864!$acc enter data create(bubmoms)
865# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
866#elif defined(MFC_OpenMP)
867# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
868!$omp target enter data map(always,alloc:bubmoms)
869# 387 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
870#endif
871
872 do i = 1, nb
873 bubrs_qbmm(i) = bub_idx%rs(i)
874 end do
875
876# 392 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
877#if defined(MFC_OpenACC)
878# 392 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
879!$acc update device(bubrs_qbmm)
880# 392 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
881#elif defined(MFC_OpenMP)
882# 392 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
883!$omp target update to(bubrs_qbmm)
884# 392 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
885#endif
886
887 do j = 1, nmom
888 do i = 1, nb
889 bubmoms(i, j) = bub_idx%moms(i, j)
890 end do
891 end do
892
893# 399 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
894#if defined(MFC_OpenACC)
895# 399 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
896!$acc update device(bubmoms)
897# 399 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
898#elif defined(MFC_OpenMP)
899# 399 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
900!$omp target update to(bubmoms)
901# 399 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
902#endif
903
904 end subroutine s_initialize_qbmm_module
905
906 !> Compute the QBMM right-hand side source terms for bubble moment transport equations
907 subroutine s_compute_qbmm_rhs(idir, q_cons_vf, q_prim_vf, rhs_vf, flux_n_vf, pb, rhs_pb)
908
909 integer, intent(in) :: idir
910 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf, q_prim_vf
911 type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
912 type(scalar_field), dimension(sys_size), intent(in) :: flux_n_vf
913 real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb
914
915 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), &
916 & intent(inout) :: rhs_pb ! TODO :: I think that this should be stp as well.
917
918 integer :: i, j, k, l, q
919 real(wp) :: nb_q, nb_dot, r, r2, nr, nr2, nr_dot, nr2_dot, var, ax
920 logical :: is_axisym
921
922 select case (idir)
923 case (1)
924 is_axisym = .false.
925 case (2)
926 is_axisym = .false.
927 case (3)
928 is_axisym = (grid_geometry == 3)
929 end select
930
931 if (.not. polytropic) then
932
933# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
934
935# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
936#if defined(MFC_OpenACC)
937# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
938!$acc parallel loop collapse(5) gang vector default(present) private(i, j, k, l, q, nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX)
939# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
940#elif defined(MFC_OpenMP)
941# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
942
943# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
944
945# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
946
947# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
948!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(5) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, l, q, nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var, AX)
949# 429 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
950#endif
951 do i = 1, nb
952 do q = 1, nnode
953 do l = 0, p
954 do k = 0, n
955 do j = 0, m
956 nb_q = q_cons_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
957 nr = q_cons_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
958 nr2 = q_cons_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)
959 r = q_prim_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
960 r2 = q_prim_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)
961 var = max(r2 - r**2._wp, sgm_eps)
962 if (q <= 2) then
963 ax = r - sqrt(var)
964 else
965 ax = r + sqrt(var)
966 end if
967
968 select case (idir)
969 case (1)
970 nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j - 1, k, &
971 & l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
972 nr_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j - 1, k, &
973 & l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
974 nr2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j - 1, k, &
975 & l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)
976 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
977 & i) - 3._wp*gam/(dx(j)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, q, i))
978 case (2)
979 nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k - 1, &
980 & l) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
981 nr_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k - 1, &
982 & l) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
983 nr2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k - 1, &
984 & l) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)
985 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
986 & i) - 3._wp*gam/(dy(k)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, q, i))
987 case (3)
988 if (is_axisym) then
989 nb_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, &
990 & l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l))
991 nr_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, &
992 & k, l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l))
993 nr2_dot = q_prim_vf(contxe + idir)%sf(j, k, l)*(flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, &
994 & k, l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l))
995 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
996 & i) - 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, &
997 & q, i))
998 else
999 nb_dot = flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, &
1000 & l - 1) - flux_n_vf(bubxb + (i - 1)*nmom)%sf(j, k, l)
1001 nr_dot = flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, &
1002 & l - 1) - flux_n_vf(bubxb + 1 + (i - 1)*nmom)%sf(j, k, l)
1003 nr2_dot = flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, &
1004 & l - 1) - flux_n_vf(bubxb + 3 + (i - 1)*nmom)%sf(j, k, l)
1005 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1006 & i) - 3._wp*gam/(dz(l)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, q, i))
1007 end if
1008 end select
1009 if (q <= 2) then
1010 select case (idir)
1011 case (1)
1012 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1013 & i) + 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1014 & *(pb(j, k, l, q, i))
1015 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1016 & i) + 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
1017 & - nr*nb_dot))*(pb(j, k, l, q, i))
1018 case (2)
1019 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1020 & i) + 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1021 & *(pb(j, k, l, q, i))
1022 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1023 & i) + 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
1024 & - nr*nb_dot))*(pb(j, k, l, q, i))
1025 case (3)
1026 if (is_axisym) then
1027 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1028 & i) + 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q &
1029 & - nr2*nb_dot)*(pb(j, k, l, q, i))
1030 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1031 & i) + 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
1032 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
1033 else
1034 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1035 & i) + 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1036 & *(pb(j, k, l, q, i))
1037 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1038 & i) + 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
1039 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
1040 end if
1041 end select
1042 else
1043 select case (idir)
1044 case (1)
1045 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1046 & i) - 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1047 & *(pb(j, k, l, q, i))
1048 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1049 & i) - 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
1050 & - nr*nb_dot))*(pb(j, k, l, q, i))
1051 case (2)
1052 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1053 & i) - 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1054 & *(pb(j, k, l, q, i))
1055 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1056 & i) - 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
1057 & - nr*nb_dot))*(pb(j, k, l, q, i))
1058 case (3)
1059 if (is_axisym) then
1060 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1061 & i) - 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q &
1062 & - nr2*nb_dot)*(pb(j, k, l, q, i))
1063 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1064 & i) - 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
1065 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
1066 else
1067 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1068 & i) - 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1069 & *(pb(j, k, l, q, i))
1070 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1071 & i) - 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
1072 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
1073 end if
1074 end select
1075 end if
1076 end do
1077 end do
1078 end do
1079 end do
1080 end do
1081
1082# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1083#if defined(MFC_OpenACC)
1084# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1085!$acc end parallel loop
1086# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1087#elif defined(MFC_OpenMP)
1088# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1089
1090# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1091!$omp end target teams loop
1092# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1093#endif
1094 end if
1095
1096 ! The following block is not repeated and is left as is
1097 if (idir == 1) then
1098
1099# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1100
1101# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1102#if defined(MFC_OpenACC)
1103# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1104!$acc parallel loop collapse(3) gang vector default(present) private(i, l, q)
1105# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1106#elif defined(MFC_OpenMP)
1107# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1108
1109# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1110
1111# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1112
1113# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1114!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, l, q)
1115# 565 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1116#endif
1117 do l = 0, p
1118 do q = 0, n
1119 do i = 0, m
1120 rhs_vf(alf_idx)%sf(i, q, l) = rhs_vf(alf_idx)%sf(i, q, l) + mom_sp(2)%sf(i, q, l)
1121 j = bubxb
1122
1123# 571 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1124#if defined(MFC_OpenACC)
1125# 571 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1126!$acc loop seq
1127# 571 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1128#elif defined(MFC_OpenMP)
1129# 571 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1130
1131# 571 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1132#endif
1133 do k = 1, nb
1134 rhs_vf(j)%sf(i, q, l) = rhs_vf(j)%sf(i, q, l) + mom_3d(0, 0, k)%sf(i, q, l)
1135 rhs_vf(j + 1)%sf(i, q, l) = rhs_vf(j + 1)%sf(i, q, l) + mom_3d(1, 0, k)%sf(i, q, l)
1136 rhs_vf(j + 2)%sf(i, q, l) = rhs_vf(j + 2)%sf(i, q, l) + mom_3d(0, 1, k)%sf(i, q, l)
1137 rhs_vf(j + 3)%sf(i, q, l) = rhs_vf(j + 3)%sf(i, q, l) + mom_3d(2, 0, k)%sf(i, q, l)
1138 rhs_vf(j + 4)%sf(i, q, l) = rhs_vf(j + 4)%sf(i, q, l) + mom_3d(1, 1, k)%sf(i, q, l)
1139 rhs_vf(j + 5)%sf(i, q, l) = rhs_vf(j + 5)%sf(i, q, l) + mom_3d(0, 2, k)%sf(i, q, l)
1140 j = j + 6
1141 end do
1142 end do
1143 end do
1144 end do
1145
1146# 584 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1147#if defined(MFC_OpenACC)
1148# 584 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1149!$acc end parallel loop
1150# 584 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1151#elif defined(MFC_OpenMP)
1152# 584 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1153
1154# 584 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1155!$omp end target teams loop
1156# 584 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1157#endif
1158 end if
1159
1160 end subroutine s_compute_qbmm_rhs
1161
1162 !> Build the coefficient array for the non-polytropic bubble model
1163 subroutine s_coeff_nonpoly(pres, rho, c, coeffs)
1164
1165
1166# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1167#ifdef _CRAYFTN
1168# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1169#if MFC_OpenACC
1170# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1171!$acc routine seq
1172# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1173#elif MFC_OpenMP
1174# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1175
1176# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1177
1178# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1179!$omp declare target device_type(any)
1180# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1181#else
1182# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1183!DIR$ INLINEALWAYS s_coeff_nonpoly
1184# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1185#endif
1186# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1187#elif MFC_OpenACC
1188# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1189!$acc routine seq
1190# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1191#elif MFC_OpenMP
1192# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1193
1194# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1195
1196# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1197!$omp declare target device_type(any)
1198# 592 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1199#endif
1200
1201 real(wp), intent(in) :: pres, rho, c
1202# 598 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1203 real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeffs
1204# 600 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1205
1206 integer :: i1, i2
1207
1208 coeffs(:,:,:) = 0._wp
1209
1210 do i2 = 0, 2; do i1 = 0, 2
1211 if ((i1 + i2) <= 2) then
1212 if (bubble_model == 3) then
1213# 609 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1214 ! RPE
1215 coeffs(1, i1, i2) = -1._wp*i2*pres/rho
1216 coeffs(2, i1, i2) = -3._wp*i2/2._wp
1217 coeffs(3, i1, i2) = i2/rho
1218 coeffs(4, i1, i2) = i1
1219 if (.not. f_is_default(re_inv)) coeffs(5, i1, i2) = -4._wp*i2*re_inv/rho
1220 if (.not. f_is_default(web)) coeffs(6, i1, i2) = -2._wp*i2/web/rho
1221 coeffs(7, i1, i2) = 0._wp
1222# 618 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1223 else if (bubble_model == 2) then
1224 ! KM with approximation of 1/(1-V/C) = 1+V/C
1225# 621 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1226 coeffs(1, i1, i2) = -3._wp*i2/2._wp
1227 coeffs(2, i1, i2) = -i2/c
1228 coeffs(3, i1, i2) = i2/(2._wp*c*c)
1229 coeffs(4, i1, i2) = -i2*pres/rho
1230 coeffs(5, i1, i2) = -2._wp*i2*pres/(c*rho)
1231 coeffs(6, i1, i2) = -i2*pres/(c*c*rho)
1232 coeffs(7, i1, i2) = i2/rho
1233 coeffs(8, i1, i2) = 2._wp*i2/(c*rho)
1234 coeffs(9, i1, i2) = i2/(c*c*rho)
1235 coeffs(10, i1, i2) = -3._wp*i2*gam/(c*rho)
1236 coeffs(11, i1, i2) = -3._wp*i2*gam/(c*c*rho)
1237 coeffs(12, i1, i2) = i1
1238 coeffs(13, i1, i2) = 0._wp
1239 coeffs(14, i1, i2) = 0._wp
1240 coeffs(15, i1, i2) = 0._wp
1241 if (.not. f_is_default(re_inv)) coeffs(16, i1, i2) = -i2*4._wp*re_inv/rho
1242 if (.not. f_is_default(web)) coeffs(17, i1, i2) = -i2*2._wp/web/rho
1243 if (.not. f_is_default(re_inv)) then
1244 coeffs(18, i1, i2) = i2*6._wp*re_inv/(rho*c)
1245 coeffs(19, i1, i2) = -i2*2._wp*re_inv/(rho*c*c)
1246 coeffs(20, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c)
1247 coeffs(21, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c*c)
1248 coeffs(22, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c)
1249 coeffs(23, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c*c)
1250 coeffs(24, i1, i2) = i2*16._wp*re_inv*re_inv/(rho*rho*c)
1251 if (.not. f_is_default(web)) then
1252 coeffs(25, i1, i2) = i2*8._wp*re_inv/web/(rho*rho*c)
1253 end if
1254 coeffs(26, i1, i2) = -12._wp*i2*gam*re_inv/(rho*rho*c*c)
1255 end if
1256 coeffs(27, i1, i2) = 3._wp*i2*gam*r_v*tw/(c*rho)
1257 coeffs(28, i1, i2) = 3._wp*i2*gam*r_v*tw/(c*c*rho)
1258 if (.not. f_is_default(re_inv)) then
1259 coeffs(29, i1, i2) = 12._wp*i2*gam*r_v*tw*re_inv/(rho*rho*c*c)
1260 end if
1261 coeffs(30, i1, i2) = 3._wp*i2*gam/(c*rho)
1262 coeffs(31, i1, i2) = 3._wp*i2*gam/(c*c*rho)
1263 if (.not. f_is_default(re_inv)) then
1264 coeffs(32, i1, i2) = 12._wp*i2*gam*re_inv/(rho*rho*c*c)
1265 end if
1266# 662 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1267 end if
1268 end if
1269 end do; end do
1270
1271 end subroutine s_coeff_nonpoly
1272
1273 !> Build the coefficient array for the polytropic bubble model
1274 subroutine s_coeff(pres, rho, c, coeffs)
1275
1276
1277# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1278#ifdef _CRAYFTN
1279# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1280#if MFC_OpenACC
1281# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1282!$acc routine seq
1283# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1284#elif MFC_OpenMP
1285# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1286
1287# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1288
1289# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1290!$omp declare target device_type(any)
1291# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1292#else
1293# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1294!DIR$ INLINEALWAYS s_coeff
1295# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1296#endif
1297# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1298#elif MFC_OpenACC
1299# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1300!$acc routine seq
1301# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1302#elif MFC_OpenMP
1303# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1304
1305# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1306
1307# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1308!$omp declare target device_type(any)
1309# 671 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1310#endif
1311
1312 real(wp), intent(in) :: pres, rho, c
1313# 677 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1314 real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeffs
1315# 679 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1316
1317 integer :: i1, i2
1318
1319 coeffs(:,:,:) = 0._wp
1320
1321 do i2 = 0, 2; do i1 = 0, 2
1322 if ((i1 + i2) <= 2) then
1323 if (bubble_model == 3) then
1324 ! RPE
1325# 689 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1326 coeffs(1, i1, i2) = -1._wp*i2*pres/rho
1327 coeffs(2, i1, i2) = -3._wp*i2/2._wp
1328 coeffs(3, i1, i2) = i2/rho
1329 coeffs(4, i1, i2) = i1
1330 if (.not. f_is_default(re_inv)) coeffs(5, i1, i2) = -4._wp*i2*re_inv/rho
1331 if (.not. f_is_default(web)) coeffs(6, i1, i2) = -2._wp*i2/web/rho
1332 coeffs(7, i1, i2) = i2*pv/rho
1333# 697 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1334 else if (bubble_model == 2) then
1335 ! KM with approximation of 1/(1-V/C) = 1+V/C
1336# 700 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1337 coeffs(1, i1, i2) = -3._wp*i2/2._wp
1338 coeffs(2, i1, i2) = -i2/c
1339 coeffs(3, i1, i2) = i2/(2._wp*c*c)
1340 coeffs(4, i1, i2) = -i2*pres/rho
1341 coeffs(5, i1, i2) = -2._wp*i2*pres/(c*rho)
1342 coeffs(6, i1, i2) = -i2*pres/(c*c*rho)
1343 coeffs(7, i1, i2) = i2/rho
1344 coeffs(8, i1, i2) = 2._wp*i2/(c*rho)
1345 coeffs(9, i1, i2) = i2/(c*c*rho)
1346 coeffs(10, i1, i2) = -3._wp*i2*gam/(c*rho)
1347 coeffs(11, i1, i2) = -3._wp*i2*gam/(c*c*rho)
1348 coeffs(12, i1, i2) = i1
1349 coeffs(13, i1, i2) = i2*(pv)/rho
1350 coeffs(14, i1, i2) = 2._wp*i2*(pv)/(c*rho)
1351 coeffs(15, i1, i2) = i2*(pv)/(c*c*rho)
1352 if (.not. f_is_default(re_inv)) coeffs(16, i1, i2) = -i2*4._wp*re_inv/rho
1353 if (.not. f_is_default(web)) coeffs(17, i1, i2) = -i2*2._wp/web/rho
1354 if (.not. f_is_default(re_inv)) then
1355 coeffs(18, i1, i2) = i2*6._wp*re_inv/(rho*c)
1356 coeffs(19, i1, i2) = -i2*2._wp*re_inv/(rho*c*c)
1357 coeffs(20, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c)
1358 coeffs(21, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c*c)
1359 coeffs(22, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c)
1360 coeffs(23, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c*c)
1361 coeffs(24, i1, i2) = i2*16._wp*re_inv*re_inv/(rho*rho*c)
1362 if (.not. f_is_default(web)) then
1363 coeffs(25, i1, i2) = i2*8._wp*re_inv/web/(rho*rho*c)
1364 end if
1365 coeffs(26, i1, i2) = -12._wp*i2*gam*re_inv/(rho*rho*c*c)
1366 end if
1367# 731 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1368 end if
1369 end if
1370 end do; end do
1371
1372 end subroutine s_coeff
1373
1374 !> Perform moment inversion to recover quadrature weights and abscissas and evaluate bubble source terms
1375 subroutine s_mom_inv(q_cons_vf, q_prim_vf, momsp, moms3d, pb, rhs_pb, mv, rhs_mv, ix, iy, iz)
1376
1377 type(scalar_field), dimension(:), intent(inout) :: q_cons_vf, q_prim_vf
1378 type(scalar_field), dimension(:), intent(inout) :: momsp
1379 type(scalar_field), dimension(0:,0:,:), intent(inout) :: moms3d
1380 real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb
1381 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_pb
1382 real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: mv
1383 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_mv
1384 type(int_bounds_info), intent(in) :: ix, iy, iz
1385
1386# 753 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1387 real(wp), dimension(nmom) :: moms, msum
1388 real(wp), dimension(nnode, nb) :: wght, abscx, abscy, wght_pb, wght_mv, wght_ht, ht
1389# 756 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1390# 759 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1391 real(wp), dimension(nterms,0:2,0:2) :: coeff
1392# 761 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1393 real(wp) :: pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, grad_t
1394 real(wp) :: n_tait, b_tait
1395 integer :: id1, id2, id3, i1, i2, j, q, r
1396
1397 is1_qbmm = ix; is2_qbmm = iy; is3_qbmm = iz
1398
1399# 766 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1400#if defined(MFC_OpenACC)
1401# 766 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1402!$acc update device(is1_qbmm, is2_qbmm, is3_qbmm)
1403# 766 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1404#elif defined(MFC_OpenMP)
1405# 766 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1406!$omp target update to(is1_qbmm, is2_qbmm, is3_qbmm)
1407# 766 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1408#endif
1409
1410
1411# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1412
1413# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1414#if defined(MFC_OpenACC)
1415# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1416!$acc parallel loop collapse(3) gang vector default(present) private(id1, id2, id3, moms, msum, wght, abscX, abscY, wght_pb, wght_mv, wght_ht, coeff, ht, r, q, n_tait, B_tait, pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, grad_T, i1, i2, j)
1417# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1418#elif defined(MFC_OpenMP)
1419# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1420
1421# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1422
1423# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1424
1425# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1426!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(id1, id2, id3, moms, msum, wght, abscX, abscY, wght_pb, wght_mv, wght_ht, coeff, ht, r, q, n_tait, B_tait, pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, grad_T, i1, i2, j)
1427# 768 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1428#endif
1429# 771 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1430 do id3 = is3_qbmm%beg, is3_qbmm%end
1431 do id2 = is2_qbmm%beg, is2_qbmm%end
1432 do id1 = is1_qbmm%beg, is1_qbmm%end
1433 alf = q_prim_vf(alf_idx)%sf(id1, id2, id3)
1434 pres = q_prim_vf(e_idx)%sf(id1, id2, id3)
1435 rho = q_prim_vf(contxb)%sf(id1, id2, id3)
1436
1437 if (bubble_model == 2) then
1438 n_tait = 1._wp/gammas(1) + 1._wp
1439 b_tait = pi_infs(1)*(n_tait - 1)/n_tait
1440 c = n_tait*(pres + b_tait)*(1._wp - alf)/(rho)
1441 c = merge(sqrt(c), sgm_eps, c > 0._wp)
1442 end if
1443
1444 call s_coeff_selector(pres, rho, c, coeff, polytropic)
1445
1446 if (alf > small_alf) then
1447 nbub = q_cons_vf(bubxb)%sf(id1, id2, id3)
1448
1449# 789 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1450#if defined(MFC_OpenACC)
1451# 789 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1452!$acc loop seq
1453# 789 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1454#elif defined(MFC_OpenMP)
1455# 789 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1456
1457# 789 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1458#endif
1459 do q = 1, nb
1460 ! Gather moments for this bubble bin
1461
1462# 792 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1463#if defined(MFC_OpenACC)
1464# 792 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1465!$acc loop seq
1466# 792 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1467#elif defined(MFC_OpenMP)
1468# 792 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1469
1470# 792 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1471#endif
1472 do r = 2, nmom
1473 moms(r) = q_prim_vf(bubmoms(q, r))%sf(id1, id2, id3)
1474 end do
1475 moms(1) = 1._wp
1476 call s_chyqmom(moms, wght(:,q), abscx(:,q), abscy(:,q))
1477
1478 if (polytropic) then
1479
1480# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1481#if defined(MFC_OpenACC)
1482# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1483!$acc loop seq
1484# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1485#elif defined(MFC_OpenMP)
1486# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1487
1488# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1489#endif
1490 do j = 1, nnode
1491 wght_pb(j, q) = wght(j, q)*(pb0(q) - pv)
1492 end do
1493 else
1494
1495# 805 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1496#if defined(MFC_OpenACC)
1497# 805 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1498!$acc loop seq
1499# 805 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1500#elif defined(MFC_OpenMP)
1501# 805 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1502
1503# 805 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1504#endif
1505 do j = 1, nnode
1506 chi_vw = 1._wp/(1._wp + r_v/r_g*(pb(id1, id2, id3, j, q)/pv - 1._wp))
1507 x_vw = m_g*chi_vw/(m_v + (m_g - m_v)*chi_vw)
1508 k_mw = x_vw*k_v(q)/(x_vw + (1._wp - x_vw)*phi_vg) + (1._wp - x_vw)*k_g(q)/(x_vw*phi_gv &
1509 & + 1._wp - x_vw)
1510 rho_mw = pv/(chi_vw*r_v*tw)
1511 rhs_mv(id1, id2, id3, j, q) = -re_trans_c(q)*((mv(id1, id2, id3, j, q)/(mv(id1, id2, id3, j, &
1512 & q) + mass_g0(q))) - chi_vw)
1513 rhs_mv(id1, id2, id3, j, q) = rho_mw*rhs_mv(id1, id2, id3, j, &
1514 & q)/pe_c/(1._wp - chi_vw)/abscx(j, q)
1515 grad_t = -re_trans_t(q)*((pb(id1, id2, id3, j, q)/pb0(q))*(abscx(j, &
1516 & q)/r0(q))**3*(mass_g0(q) + mass_v0(q))/(mass_g0(q) + mv(id1, id2, id3, &
1517 & j, q)) - 1._wp)
1518 ht(j, q) = pb0(q)*k_mw*grad_t/pe_t(q)/abscx(j, q)
1519 wght_pb(j, q) = wght(j, q)*(pb(id1, id2, id3, j, q))
1520 wght_mv(j, q) = wght(j, q)*(rhs_mv(id1, id2, id3, j, q))
1521 wght_ht(j, q) = wght(j, q)*ht(j, q)
1522 end do
1523 end if
1524
1525 ! Compute change in moments due to bubble dynamics
1526 r = 1
1527
1528# 828 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1529#if defined(MFC_OpenACC)
1530# 828 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1531!$acc loop seq
1532# 828 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1533#elif defined(MFC_OpenMP)
1534# 828 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1535
1536# 828 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1537#endif
1538 do i2 = 0, 2
1539
1540# 830 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1541#if defined(MFC_OpenACC)
1542# 830 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1543!$acc loop seq
1544# 830 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1545#elif defined(MFC_OpenMP)
1546# 830 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1547
1548# 830 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1549#endif
1550 do i1 = 0, 2
1551 if ((i1 + i2) <= 2) then
1552 momsum = 0._wp
1553
1554# 834 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1555#if defined(MFC_OpenACC)
1556# 834 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1557!$acc loop seq
1558# 834 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1559#elif defined(MFC_OpenMP)
1560# 834 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1561
1562# 834 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1563#endif
1564 do j = 1, nterms
1565 select case (bubble_model)
1566 case (3)
1567 if (j == 3) then
1568 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1569 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_pb(:,q), &
1570 & momrhs(:,i1, i2, j, q))
1571 else
1572 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1573 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght(:,q), &
1574 & momrhs(:,i1, i2, j, q))
1575 end if
1576 case (2)
1577 if ((j >= 7 .and. j <= 9) .or. (j >= 22 .and. j <= 23) &
1578 & .or. (j >= 10 .and. j <= 11) .or. (j == 26)) then
1579 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1580 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_pb(:,q), &
1581 & momrhs(:,i1, i2, j, q))
1582 else if ((j >= 27 .and. j <= 29) .and. (.not. polytropic)) then
1583 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1584 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_mv(:,q), &
1585 & momrhs(:,i1, i2, j, q))
1586 else if ((j >= 30 .and. j <= 32) .and. (.not. polytropic)) then
1587 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1588 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_ht(:,q), &
1589 & momrhs(:,i1, i2, j, q))
1590 else
1591 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1592 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght(:,q), &
1593 & momrhs(:,i1, i2, j, q))
1594 end if
1595 end select
1596 end do
1597 moms3d(i1, i2, q)%sf(id1, id2, id3) = nbub*momsum
1598 msum(r) = momsum
1599 r = r + 1
1600 end if
1601 end do
1602 end do
1603
1604 ! Compute change in pb and mv for non-polytropic model
1605 if (.not. polytropic) then
1606
1607# 877 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1608#if defined(MFC_OpenACC)
1609# 877 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1610!$acc loop seq
1611# 877 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1612#elif defined(MFC_OpenMP)
1613# 877 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1614
1615# 877 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1616#endif
1617 do j = 1, nnode
1618 drdt = msum(2)
1619 drdt2 = merge(-1._wp, 1._wp, j == 1 .or. j == 2)/(2._wp*sqrt(merge(moms(4) - moms(2)**2._wp, &
1620 & sgm_eps, moms(4) - moms(2)**2._wp > 0._wp)))
1621 drdt2 = drdt2*(msum(3) - 2._wp*moms(2)*msum(2))
1622 drdt = drdt + drdt2
1623 rhs_pb(id1, id2, id3, j, q) = (-3._wp*gam*drdt/abscx(j, q))*(pb(id1, id2, id3, j, q))
1624 rhs_pb(id1, id2, id3, j, q) = rhs_pb(id1, id2, id3, j, q) + (3._wp*gam/abscx(j, &
1625 & q))*rhs_mv(id1, id2, id3, j, q)*r_v*tw
1626 rhs_pb(id1, id2, id3, j, q) = rhs_pb(id1, id2, id3, j, q) + (3._wp*gam/abscx(j, q))*ht(j, q)
1627 rhs_mv(id1, id2, id3, j, q) = rhs_mv(id1, id2, id3, j, q)*(4._wp*pi*abscx(j, q)**2._wp)
1628 end do
1629 end if
1630 end do
1631
1632 ! Compute special high-order moments
1633 momsp(1)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght, 3._wp, 0._wp, 0._wp)
1634 momsp(2)%sf(id1, id2, id3) = 4._wp*pi*nbub*f_quad(abscx, abscy, wght, 2._wp, 1._wp, 0._wp)
1635 momsp(3)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght, 3._wp, 2._wp, 0._wp)
1636 if (abs(gam - 1._wp) <= 1.e-4_wp) then
1637 momsp(4)%sf(id1, id2, id3) = 1._wp
1638 else
1639 if (polytropic) then
1640 momsp(4)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght_pb, 3._wp*(1._wp - gam), 0._wp, &
1641 & 3._wp*gam) + pv*f_quad(abscx, abscy, wght, 3._wp, 0._wp, &
1642 & 0._wp) - 4._wp*re_inv*f_quad(abscx, abscy, wght, 2._wp, 1._wp, &
1643 & 0._wp) - (2._wp/web)*f_quad(abscx, abscy, wght, 2._wp, 0._wp, 0._wp)
1644 else
1645 momsp(4)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght_pb, 3._wp, 0._wp, &
1646 & 0._wp) - 4._wp*re_inv*f_quad(abscx, abscy, wght, 2._wp, 1._wp, &
1647 & 0._wp) - (2._wp/web)*f_quad(abscx, abscy, wght, 2._wp, 0._wp, 0._wp)
1648 end if
1649 end if
1650 else
1651
1652# 912 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1653#if defined(MFC_OpenACC)
1654# 912 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1655!$acc loop seq
1656# 912 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1657#elif defined(MFC_OpenMP)
1658# 912 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1659
1660# 912 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1661#endif
1662 do q = 1, nb
1663
1664# 914 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1665#if defined(MFC_OpenACC)
1666# 914 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1667!$acc loop seq
1668# 914 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1669#elif defined(MFC_OpenMP)
1670# 914 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1671
1672# 914 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1673#endif
1674 do i1 = 0, 2
1675
1676# 916 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1677#if defined(MFC_OpenACC)
1678# 916 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1679!$acc loop seq
1680# 916 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1681#elif defined(MFC_OpenMP)
1682# 916 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1683
1684# 916 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1685#endif
1686 do i2 = 0, 2
1687 moms3d(i1, i2, q)%sf(id1, id2, id3) = 0._wp
1688 end do
1689 end do
1690 end do
1691 momsp(1)%sf(id1, id2, id3) = 0._wp
1692 momsp(2)%sf(id1, id2, id3) = 0._wp
1693 momsp(3)%sf(id1, id2, id3) = 0._wp
1694 momsp(4)%sf(id1, id2, id3) = 0._wp
1695 end if
1696 end do
1697 end do
1698 end do
1699
1700# 930 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1701#if defined(MFC_OpenACC)
1702# 930 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1703!$acc end parallel loop
1704# 930 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1705#elif defined(MFC_OpenMP)
1706# 930 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1707
1708# 930 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1709!$omp end target teams loop
1710# 930 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1711#endif
1712
1713 contains
1714 !> Select the polytropic or non-polytropic coefficient routine
1715 subroutine s_coeff_selector(pres, rho, c, coeff, polytropic)
1716
1717
1718# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1719#ifdef _CRAYFTN
1720# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1721#if MFC_OpenACC
1722# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1723!$acc routine seq
1724# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1725#elif MFC_OpenMP
1726# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1727
1728# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1729
1730# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1731!$omp declare target device_type(any)
1732# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1733#else
1734# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1735!DIR$ INLINEALWAYS s_coeff_selector
1736# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1737#endif
1738# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1739#elif MFC_OpenACC
1740# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1741!$acc routine seq
1742# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1743#elif MFC_OpenMP
1744# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1745
1746# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1747
1748# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1749!$omp declare target device_type(any)
1750# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1751#endif
1752 real(wp), intent(in) :: pres, rho, c
1753# 941 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1754 real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeff
1755# 943 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1756 logical, intent(in) :: polytropic
1757 if (polytropic) then
1758 call s_coeff(pres, rho, c, coeff)
1759 else
1760 call s_coeff_nonpoly(pres, rho, c, coeff)
1761 end if
1762
1763 end subroutine s_coeff_selector
1764
1765 !> Perform CHyQMOM inversion for bivariate moments
1766 subroutine s_chyqmom(momin, wght, abscX, abscY)
1767
1768
1769# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1770#ifdef _CRAYFTN
1771# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1772#if MFC_OpenACC
1773# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1774!$acc routine seq
1775# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1776#elif MFC_OpenMP
1777# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1778
1779# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1780
1781# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1782!$omp declare target device_type(any)
1783# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1784#else
1785# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1786!DIR$ INLINEALWAYS s_chyqmom
1787# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1788#endif
1789# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1790#elif MFC_OpenACC
1791# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1792!$acc routine seq
1793# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1794#elif MFC_OpenMP
1795# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1796
1797# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1798
1799# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1800!$omp declare target device_type(any)
1801# 955 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1802#endif
1803
1804 real(wp), dimension(nmom), intent(in) :: momin
1805 real(wp), dimension(nnode), intent(inout) :: wght, abscX, abscY
1806
1807 ! Local variables
1808 real(wp), dimension(0:2,0:2) :: moms
1809 real(wp), dimension(3) :: M1, M3
1810 real(wp), dimension(2) :: myrho, myrho3, up, up3, Vf
1811 real(wp) :: bu, bv, d20, d11, d_02, c20, c11, c02
1812 real(wp) :: mu2, vp21, vp22, rho21, rho22
1813
1814 ! Assign moments to 2D array for clarity
1815 moms(0, 0) = momin(1)
1816 moms(1, 0) = momin(2)
1817 moms(0, 1) = momin(3)
1818 moms(2, 0) = momin(4)
1819 moms(1, 1) = momin(5)
1820 moms(0, 2) = momin(6)
1821
1822 ! Compute means and central moments
1823 bu = moms(1, 0)/moms(0, 0)
1824 bv = moms(0, 1)/moms(0, 0)
1825 d20 = moms(2, 0)/moms(0, 0)
1826 d11 = moms(1, 1)/moms(0, 0)
1827 d_02 = moms(0, 2)/moms(0, 0)
1828
1829 c20 = d20 - bu**2._wp
1830 c11 = d11 - bu*bv
1831 c02 = d_02 - bv**2._wp
1832
1833 ! First 1D quadrature (X direction)
1834 m1 = (/1._wp, 0._wp, c20/)
1835 call s_hyqmom(myrho, up, m1)
1836 vf = c11*up/c20
1837
1838 ! Second 1D quadrature (Y direction, conditional on X)
1839 mu2 = max(0._wp, c02 - sum(myrho*(vf**2._wp)))
1840 m3 = (/1._wp, 0._wp, mu2/)
1841 call s_hyqmom(myrho3, up3, m3)
1842
1843 ! Assign roots and weights for 2D quadrature
1844 vp21 = up3(1)
1845 vp22 = up3(2)
1846 rho21 = myrho3(1)
1847 rho22 = myrho3(2)
1848
1849 ! Compute weights (vectorized)
1850 wght = moms(0, 0)*[myrho(1)*rho21, myrho(1)*rho22, myrho(2)*rho21, myrho(2)*rho22]
1851
1852 ! Compute abscissas (vectorized)
1853 abscx = bu + [up(1), up(1), up(2), up(2)]
1854 abscy = bv + [vf(1) + vp21, vf(1) + vp22, vf(2) + vp21, vf(2) + vp22]
1855
1856 end subroutine s_chyqmom
1857
1858 !> Perform HyQMOM inversion for univariate moments
1859 subroutine s_hyqmom(frho, fup, fmom)
1860
1861
1862# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1863#ifdef _CRAYFTN
1864# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1865#if MFC_OpenACC
1866# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1867!$acc routine seq
1868# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1869#elif MFC_OpenMP
1870# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1871
1872# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1873
1874# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1875!$omp declare target device_type(any)
1876# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1877#else
1878# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1879!DIR$ INLINEALWAYS s_hyqmom
1880# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1881#endif
1882# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1883#elif MFC_OpenACC
1884# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1885!$acc routine seq
1886# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1887#elif MFC_OpenMP
1888# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1889
1890# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1891
1892# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1893!$omp declare target device_type(any)
1894# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1895#endif
1896
1897 real(wp), dimension(2), intent(inout) :: frho, fup
1898 real(wp), dimension(3), intent(in) :: fmom
1899 real(wp) :: bu, d2, c2
1900
1901 bu = fmom(2)/fmom(1)
1902 d2 = fmom(3)/fmom(1)
1903 c2 = d2 - bu**2._wp
1904 frho(1) = fmom(1)/2._wp
1905 frho(2) = fmom(1)/2._wp
1906 c2 = maxval((/c2, sgm_eps/))
1907 fup(1) = bu - sqrt(c2)
1908 fup(2) = bu + sqrt(c2)
1909
1910 end subroutine s_hyqmom
1911
1912 !> Evaluate a weighted quadrature sum over all bubble size bins and nodes
1913 function f_quad(abscX, abscY, wght_in, q, r, s)
1914
1915
1916# 1034 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1917#if MFC_OpenACC
1918# 1034 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1919!$acc routine seq
1920# 1034 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1921#elif MFC_OpenMP
1922# 1034 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1923
1924# 1034 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1925
1926# 1034 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1927!$omp declare target device_type(any)
1928# 1034 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1929#endif
1930# 1038 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1931 real(wp), dimension(nnode, nb), intent(in) :: abscx, abscy, wght_in
1932# 1040 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1933 real(wp), intent(in) :: q, r, s
1934 real(wp) :: f_quad_rv, f_quad
1935 integer :: i, i1
1936
1937 f_quad = 0._wp
1938
1939# 1045 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1940#if defined(MFC_OpenACC)
1941# 1045 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1942!$acc loop seq
1943# 1045 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1944#elif defined(MFC_OpenMP)
1945# 1045 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1946
1947# 1045 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1948#endif
1949 do i = 1, nb
1950 f_quad_rv = 0._wp
1951
1952# 1048 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1953#if defined(MFC_OpenACC)
1954# 1048 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1955!$acc loop seq
1956# 1048 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1957#elif defined(MFC_OpenMP)
1958# 1048 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1959
1960# 1048 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1961#endif
1962 do i1 = 1, nnode
1963 f_quad_rv = f_quad_rv + wght_in(i1, i)*(abscx(i1, i)**q)*(abscy(i1, i)**r)
1964 end do
1965 f_quad = f_quad + weight(i)*(r0(i)**s)*f_quad_rv
1966 end do
1967
1968 end function f_quad
1969
1970 !> Evaluate a weighted 2D quadrature sum over quadrature nodes for a single size bin
1971 function f_quad2d(abscX, abscY, wght_in, pow)
1972
1973
1974# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1975#if MFC_OpenACC
1976# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1977!$acc routine seq
1978# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1979#elif MFC_OpenMP
1980# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1981
1982# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1983
1984# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1985!$omp declare target device_type(any)
1986# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1987#endif
1988# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1989 real(wp), dimension(nnode), intent(in) :: abscx, abscy, wght_in
1990# 1066 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1991 real(wp), dimension(3), intent(in) :: pow
1992 real(wp) :: f_quad2d
1993 integer :: i
1994
1995 f_quad2d = 0._wp
1996
1997# 1071 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1998#if defined(MFC_OpenACC)
1999# 1071 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
2000!$acc loop seq
2001# 1071 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
2002#elif defined(MFC_OpenMP)
2003# 1071 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
2004
2005# 1071 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
2006#endif
2007 do i = 1, nnode
2008 f_quad2d = f_quad2d + wght_in(i)*(abscx(i)**pow(1))*(abscy(i)**pow(2))
2009 end do
2010
2011 end function f_quad2d
2012
2013 end subroutine s_mom_inv
2014
2015end module m_qbmm
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
integer, intent(in) k
integer, intent(in) j
integer, intent(in) l
subroutine s_hyqmom(frho, fup, fmom)
Perform HyQMOM inversion for univariate moments.
subroutine s_coeff_selector(pres, rho, c, coeff, polytropic)
Select the polytropic or non-polytropic coefficient routine.
subroutine s_chyqmom(momin, wght, abscx, abscy)
Perform CHyQMOM inversion for bivariate moments.
real(wp) function f_quad(abscx, abscy, wght_in, q, r, s)
Evaluate a weighted quadrature sum over all bubble size bins and nodes.
real(wp) function f_quad2d(abscx, abscy, wght_in, pow)
Evaluate a weighted 2D quadrature sum over quadrature nodes for a single size bin.
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
real(wp) re_inv
Inverse Reynolds number.
logical polytropic
Polytropic switch.
integer nb
Number of eq. bubble sizes.
integer bubble_model
Gilmore or Keller–Miksis bubble model.
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
logical elemental function, public f_is_default(var)
Checks if a real(wp) variable is of default value.
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.
Quadrature-based moment methods (QBMM) for polydisperse bubble moment inversion and transport.
integer, dimension(:,:), allocatable bubmoms
integer, dimension(:), allocatable bubrs_qbmm
type(int_bounds_info) is1_qbmm
subroutine, public s_compute_qbmm_rhs(idir, q_cons_vf, q_prim_vf, rhs_vf, flux_n_vf, pb, rhs_pb)
Compute the QBMM right-hand side source terms for bubble moment transport equations.
integer nterms
subroutine s_coeff_nonpoly(pres, rho, c, coeffs)
Build the coefficient array for the non-polytropic bubble model.
real(wp), dimension(:,:,:,:,:), allocatable momrhs
subroutine, public s_mom_inv(q_cons_vf, q_prim_vf, momsp, moms3d, pb, rhs_pb, mv, rhs_mv, ix, iy, iz)
Perform moment inversion to recover quadrature weights and abscissas and evaluate bubble source terms...
type(int_bounds_info) is3_qbmm
type(int_bounds_info) is2_qbmm
subroutine, public s_coeff(pres, rho, c, coeffs)
Build the coefficient array for the polytropic bubble model.
impure subroutine, public s_initialize_qbmm_module
Initialize the QBMM module.
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Integer bounds for variables.