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