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