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