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