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