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