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 ! TODO :: I think that this should be stp as well.
867 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_pb
868 integer :: i, j, k, l, q
869 real(wp) :: nb_q, nb_dot, r, r2, nr, nr2, nr_dot, nr2_dot, var, ax
870 logical :: is_axisym
871
872 select case (idir)
873 case (1)
874 is_axisym = .false.
875 case (2)
876 is_axisym = .false.
877 case (3)
878 is_axisym = (grid_geometry == 3)
879 end select
880
881 if (.not. polytropic) then
882
883# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
884
885# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
886#if defined(MFC_OpenACC)
887# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
888!$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)
889# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
890#elif defined(MFC_OpenMP)
891# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
892
893# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
894
895# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
896
897# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
898!$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)
899# 421 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
900#endif
901 do i = 1, nb
902 do q = 1, nnode
903 do l = 0, p
904 do k = 0, n
905 do j = 0, m
906 nb_q = q_cons_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k, l)
907 nr = q_cons_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l)
908 nr2 = q_cons_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l)
909 r = q_prim_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l)
910 r2 = q_prim_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l)
911 var = max(r2 - r**2._wp, sgm_eps)
912 if (q <= 2) then
913 ax = r - sqrt(var)
914 else
915 ax = r + sqrt(var)
916 end if
917
918 select case (idir)
919 case (1)
920 nb_dot = flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j - 1, k, &
921 & l) - flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k, l)
922 nr_dot = flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j - 1, k, &
923 & l) - flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l)
924 nr2_dot = flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j - 1, k, &
925 & l) - flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l)
926 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
927 & i) - 3._wp*gam/(dx(j)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, q, i))
928 case (2)
929 nb_dot = flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k - 1, &
930 & l) - flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k, l)
931 nr_dot = flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k - 1, &
932 & l) - flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l)
933 nr2_dot = flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k - 1, &
934 & l) - flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l)
935 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
936 & i) - 3._wp*gam/(dy(k)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, q, i))
937 case (3)
938 if (is_axisym) then
939 nb_dot = q_prim_vf(eqn_idx%cont%end + idir)%sf(j, k, &
940 & l)*(flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k, &
941 & l - 1) - flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k, l))
942 nr_dot = q_prim_vf(eqn_idx%cont%end + idir)%sf(j, k, &
943 & l)*(flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, &
944 & l - 1) - flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l))
945 nr2_dot = q_prim_vf(eqn_idx%cont%end + idir)%sf(j, k, &
946 & l)*(flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, &
947 & l - 1) - flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l))
948 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
949 & i) - 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, &
950 & q, i))
951 else
952 nb_dot = flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k, &
953 & l - 1) - flux_n_vf(eqn_idx%bub%beg + (i - 1)*nmom)%sf(j, k, l)
954 nr_dot = flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, &
955 & l - 1) - flux_n_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l)
956 nr2_dot = flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, &
957 & l - 1) - flux_n_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l)
958 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
959 & i) - 3._wp*gam/(dz(l)*ax*nb_q**2)*(nr_dot*nb_q - nr*nb_dot)*(pb(j, k, l, q, i))
960 end if
961 end select
962 if (q <= 2) then
963 select case (idir)
964 case (1)
965 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
966 & i) + 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
967 & *(pb(j, k, l, q, i))
968 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
969 & i) + 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
970 & - nr*nb_dot))*(pb(j, k, l, q, i))
971 case (2)
972 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
973 & i) + 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
974 & *(pb(j, k, l, q, i))
975 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
976 & i) + 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
977 & - nr*nb_dot))*(pb(j, k, l, q, i))
978 case (3)
979 if (is_axisym) then
980 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
981 & i) + 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q &
982 & - nr2*nb_dot)*(pb(j, k, l, q, i))
983 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
984 & i) + 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
985 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
986 else
987 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
988 & i) + 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
989 & *(pb(j, k, l, q, i))
990 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
991 & i) + 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
992 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
993 end if
994 end select
995 else
996 select case (idir)
997 case (1)
998 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
999 & i) - 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1000 & *(pb(j, k, l, q, i))
1001 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1002 & i) - 3._wp*gam/(dx(j)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
1003 & - nr*nb_dot))*(pb(j, k, l, q, i))
1004 case (2)
1005 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1006 & i) - 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1007 & *(pb(j, k, l, q, i))
1008 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1009 & i) - 3._wp*gam/(dy(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q)*(nr_dot*nb_q &
1010 & - nr*nb_dot))*(pb(j, k, l, q, i))
1011 case (3)
1012 if (is_axisym) then
1013 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1014 & i) - 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q &
1015 & - nr2*nb_dot)*(pb(j, k, l, q, i))
1016 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1017 & i) - 3._wp*gam/(dz(l)*y_cc(k)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
1018 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
1019 else
1020 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1021 & i) - 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(nr2_dot*nb_q - nr2*nb_dot) &
1022 & *(pb(j, k, l, q, i))
1023 rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, &
1024 & i) - 3._wp*gam/(dz(l)*ax*nb_q**2*sqrt(var)*2._wp)*(-2._wp*(nr/nb_q) &
1025 & *(nr_dot*nb_q - nr*nb_dot))*(pb(j, k, l, q, i))
1026 end if
1027 end select
1028 end if
1029 end do
1030 end do
1031 end do
1032 end do
1033 end do
1034
1035# 555 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1036#if defined(MFC_OpenACC)
1037# 555 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1038!$acc end parallel loop
1039# 555 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1040#elif defined(MFC_OpenMP)
1041# 555 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1042
1043# 555 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1044!$omp end target teams loop
1045# 555 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1046#endif
1047 end if
1048
1049 ! The following block is not repeated and is left as is
1050 if (idir == 1) then
1051
1052# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1053
1054# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1055#if defined(MFC_OpenACC)
1056# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1057!$acc parallel loop collapse(3) gang vector default(present) private(i, l, q)
1058# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1059#elif defined(MFC_OpenMP)
1060# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1061
1062# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1063
1064# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1065
1066# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1067!$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)
1068# 560 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1069#endif
1070 do l = 0, p
1071 do q = 0, n
1072 do i = 0, m
1073 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)
1074 j = eqn_idx%bub%beg
1075
1076# 566 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1077#if defined(MFC_OpenACC)
1078# 566 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1079!$acc loop seq
1080# 566 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1081#elif defined(MFC_OpenMP)
1082# 566 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1083
1084# 566 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1085#endif
1086 do k = 1, nb
1087 rhs_vf(j)%sf(i, q, l) = rhs_vf(j)%sf(i, q, l) + mom_3d(0, 0, k)%sf(i, q, l)
1088 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)
1089 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)
1090 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)
1091 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)
1092 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)
1093 j = j + 6
1094 end do
1095 end do
1096 end do
1097 end do
1098
1099# 579 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1100#if defined(MFC_OpenACC)
1101# 579 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1102!$acc end parallel loop
1103# 579 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1104#elif defined(MFC_OpenMP)
1105# 579 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1106
1107# 579 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1108!$omp end target teams loop
1109# 579 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1110#endif
1111 end if
1112
1113 end subroutine s_compute_qbmm_rhs
1114
1115 !> Build the coefficient array for the non-polytropic bubble model
1116 subroutine s_coeff_nonpoly(pres, rho, c, coeffs)
1117
1118
1119# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1120#ifdef _CRAYFTN
1121# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1122#if MFC_OpenACC
1123# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1124!$acc routine seq
1125# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1126#elif MFC_OpenMP
1127# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1128
1129# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1130
1131# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1132!$omp declare target device_type(any)
1133# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1134#else
1135# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1136!DIR$ INLINEALWAYS s_coeff_nonpoly
1137# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1138#endif
1139# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1140#elif MFC_OpenACC
1141# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1142!$acc routine seq
1143# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1144#elif MFC_OpenMP
1145# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1146
1147# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1148
1149# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1150!$omp declare target device_type(any)
1151# 587 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1152#endif
1153
1154 real(wp), intent(in) :: pres, rho, c
1155# 593 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1156 real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeffs
1157# 595 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1158
1159 integer :: i1, i2
1160
1161 coeffs(:,:,:) = 0._wp
1162
1163 do i2 = 0, 2; do i1 = 0, 2
1164 if ((i1 + i2) <= 2) then
1165 if (bubble_model == 3) then
1166# 604 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1167 ! RPE
1168 coeffs(1, i1, i2) = -1._wp*i2*pres/rho
1169 coeffs(2, i1, i2) = -3._wp*i2/2._wp
1170 coeffs(3, i1, i2) = i2/rho
1171 coeffs(4, i1, i2) = i1
1172 if (.not. f_is_default(re_inv)) coeffs(5, i1, i2) = -4._wp*i2*re_inv/rho
1173 if (.not. f_is_default(web)) coeffs(6, i1, i2) = -2._wp*i2/web/rho
1174 coeffs(7, i1, i2) = 0._wp
1175# 613 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1176 else if (bubble_model == 2) then
1177 ! KM with approximation of 1/(1-V/C) = 1+V/C
1178# 616 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1179 coeffs(1, i1, i2) = -3._wp*i2/2._wp
1180 coeffs(2, i1, i2) = -i2/c
1181 coeffs(3, i1, i2) = i2/(2._wp*c*c)
1182 coeffs(4, i1, i2) = -i2*pres/rho
1183 coeffs(5, i1, i2) = -2._wp*i2*pres/(c*rho)
1184 coeffs(6, i1, i2) = -i2*pres/(c*c*rho)
1185 coeffs(7, i1, i2) = i2/rho
1186 coeffs(8, i1, i2) = 2._wp*i2/(c*rho)
1187 coeffs(9, i1, i2) = i2/(c*c*rho)
1188 coeffs(10, i1, i2) = -3._wp*i2*gam/(c*rho)
1189 coeffs(11, i1, i2) = -3._wp*i2*gam/(c*c*rho)
1190 coeffs(12, i1, i2) = i1
1191 coeffs(13, i1, i2) = 0._wp
1192 coeffs(14, i1, i2) = 0._wp
1193 coeffs(15, i1, i2) = 0._wp
1194 if (.not. f_is_default(re_inv)) coeffs(16, i1, i2) = -i2*4._wp*re_inv/rho
1195 if (.not. f_is_default(web)) coeffs(17, i1, i2) = -i2*2._wp/web/rho
1196 if (.not. f_is_default(re_inv)) then
1197 coeffs(18, i1, i2) = i2*6._wp*re_inv/(rho*c)
1198 coeffs(19, i1, i2) = -i2*2._wp*re_inv/(rho*c*c)
1199 coeffs(20, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c)
1200 coeffs(21, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c*c)
1201 coeffs(22, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c)
1202 coeffs(23, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c*c)
1203 coeffs(24, i1, i2) = i2*16._wp*re_inv*re_inv/(rho*rho*c)
1204 if (.not. f_is_default(web)) then
1205 coeffs(25, i1, i2) = i2*8._wp*re_inv/web/(rho*rho*c)
1206 end if
1207 coeffs(26, i1, i2) = -12._wp*i2*gam*re_inv/(rho*rho*c*c)
1208 end if
1209 coeffs(27, i1, i2) = 3._wp*i2*gam*r_v*tw/(c*rho)
1210 coeffs(28, i1, i2) = 3._wp*i2*gam*r_v*tw/(c*c*rho)
1211 if (.not. f_is_default(re_inv)) then
1212 coeffs(29, i1, i2) = 12._wp*i2*gam*r_v*tw*re_inv/(rho*rho*c*c)
1213 end if
1214 coeffs(30, i1, i2) = 3._wp*i2*gam/(c*rho)
1215 coeffs(31, i1, i2) = 3._wp*i2*gam/(c*c*rho)
1216 if (.not. f_is_default(re_inv)) then
1217 coeffs(32, i1, i2) = 12._wp*i2*gam*re_inv/(rho*rho*c*c)
1218 end if
1219# 657 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1220 end if
1221 end if
1222 end do; end do
1223
1224 end subroutine s_coeff_nonpoly
1225
1226 !> Build the coefficient array for the polytropic bubble model
1227 subroutine s_coeff(pres, rho, c, coeffs)
1228
1229
1230# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1231#ifdef _CRAYFTN
1232# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1233#if MFC_OpenACC
1234# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1235!$acc routine seq
1236# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1237#elif MFC_OpenMP
1238# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1239
1240# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1241
1242# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1243!$omp declare target device_type(any)
1244# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1245#else
1246# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1247!DIR$ INLINEALWAYS s_coeff
1248# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1249#endif
1250# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1251#elif MFC_OpenACC
1252# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1253!$acc routine seq
1254# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1255#elif MFC_OpenMP
1256# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1257
1258# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1259
1260# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1261!$omp declare target device_type(any)
1262# 666 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1263#endif
1264
1265 real(wp), intent(in) :: pres, rho, c
1266# 672 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1267 real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeffs
1268# 674 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1269
1270 integer :: i1, i2
1271
1272 coeffs(:,:,:) = 0._wp
1273
1274 do i2 = 0, 2; do i1 = 0, 2
1275 if ((i1 + i2) <= 2) then
1276 if (bubble_model == 3) then
1277 ! RPE
1278# 684 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1279 coeffs(1, i1, i2) = -1._wp*i2*pres/rho
1280 coeffs(2, i1, i2) = -3._wp*i2/2._wp
1281 coeffs(3, i1, i2) = i2/rho
1282 coeffs(4, i1, i2) = i1
1283 if (.not. f_is_default(re_inv)) coeffs(5, i1, i2) = -4._wp*i2*re_inv/rho
1284 if (.not. f_is_default(web)) coeffs(6, i1, i2) = -2._wp*i2/web/rho
1285 coeffs(7, i1, i2) = i2*pv/rho
1286# 692 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1287 else if (bubble_model == 2) then
1288 ! KM with approximation of 1/(1-V/C) = 1+V/C
1289# 695 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1290 coeffs(1, i1, i2) = -3._wp*i2/2._wp
1291 coeffs(2, i1, i2) = -i2/c
1292 coeffs(3, i1, i2) = i2/(2._wp*c*c)
1293 coeffs(4, i1, i2) = -i2*pres/rho
1294 coeffs(5, i1, i2) = -2._wp*i2*pres/(c*rho)
1295 coeffs(6, i1, i2) = -i2*pres/(c*c*rho)
1296 coeffs(7, i1, i2) = i2/rho
1297 coeffs(8, i1, i2) = 2._wp*i2/(c*rho)
1298 coeffs(9, i1, i2) = i2/(c*c*rho)
1299 coeffs(10, i1, i2) = -3._wp*i2*gam/(c*rho)
1300 coeffs(11, i1, i2) = -3._wp*i2*gam/(c*c*rho)
1301 coeffs(12, i1, i2) = i1
1302 coeffs(13, i1, i2) = i2*(pv)/rho
1303 coeffs(14, i1, i2) = 2._wp*i2*(pv)/(c*rho)
1304 coeffs(15, i1, i2) = i2*(pv)/(c*c*rho)
1305 if (.not. f_is_default(re_inv)) coeffs(16, i1, i2) = -i2*4._wp*re_inv/rho
1306 if (.not. f_is_default(web)) coeffs(17, i1, i2) = -i2*2._wp/web/rho
1307 if (.not. f_is_default(re_inv)) then
1308 coeffs(18, i1, i2) = i2*6._wp*re_inv/(rho*c)
1309 coeffs(19, i1, i2) = -i2*2._wp*re_inv/(rho*c*c)
1310 coeffs(20, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c)
1311 coeffs(21, i1, i2) = i2*4._wp*pres*re_inv/(rho*rho*c*c)
1312 coeffs(22, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c)
1313 coeffs(23, i1, i2) = -i2*4._wp*re_inv/(rho*rho*c*c)
1314 coeffs(24, i1, i2) = i2*16._wp*re_inv*re_inv/(rho*rho*c)
1315 if (.not. f_is_default(web)) then
1316 coeffs(25, i1, i2) = i2*8._wp*re_inv/web/(rho*rho*c)
1317 end if
1318 coeffs(26, i1, i2) = -12._wp*i2*gam*re_inv/(rho*rho*c*c)
1319 end if
1320# 726 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1321 end if
1322 end if
1323 end do; end do
1324
1325 end subroutine s_coeff
1326
1327 !> Perform moment inversion to recover quadrature weights and abscissas and evaluate bubble source terms
1328 subroutine s_mom_inv(q_cons_vf, q_prim_vf, momsp, moms3d, pb, rhs_pb, mv, rhs_mv, ix, iy, iz)
1329
1330 type(scalar_field), dimension(:), intent(inout) :: q_cons_vf, q_prim_vf
1331 type(scalar_field), dimension(:), intent(inout) :: momsp
1332 type(scalar_field), dimension(0:,0:,:), intent(inout) :: moms3d
1333 real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb
1334 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_pb
1335 real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: mv
1336 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_mv
1337 type(int_bounds_info), intent(in) :: ix, iy, iz
1338
1339# 748 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1340 real(wp), dimension(nmom) :: moms, msum
1341 real(wp), dimension(nnode, nb) :: wght, abscx, abscy, wght_pb, wght_mv, wght_ht, ht
1342# 751 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1343# 754 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1344 real(wp), dimension(nterms,0:2,0:2) :: coeff
1345# 756 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1346 real(wp) :: pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, grad_t
1347 real(wp) :: n_tait, b_tait
1348 integer :: id1, id2, id3, i1, i2, j, q, r
1349
1350 is1_qbmm = ix; is2_qbmm = iy; is3_qbmm = iz
1351
1352# 761 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1353#if defined(MFC_OpenACC)
1354# 761 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1355!$acc update device(is1_qbmm, is2_qbmm, is3_qbmm)
1356# 761 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1357#elif defined(MFC_OpenMP)
1358# 761 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1359!$omp target update to(is1_qbmm, is2_qbmm, is3_qbmm)
1360# 761 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1361#endif
1362
1363
1364# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1365
1366# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1367#if defined(MFC_OpenACC)
1368# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1369!$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)
1370# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1371#elif defined(MFC_OpenMP)
1372# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1373
1374# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1375
1376# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1377
1378# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1379!$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)
1380# 763 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1381#endif
1382# 766 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1383 do id3 = is3_qbmm%beg, is3_qbmm%end
1384 do id2 = is2_qbmm%beg, is2_qbmm%end
1385 do id1 = is1_qbmm%beg, is1_qbmm%end
1386 alf = q_prim_vf(eqn_idx%alf)%sf(id1, id2, id3)
1387 pres = q_prim_vf(eqn_idx%E)%sf(id1, id2, id3)
1388 rho = q_prim_vf(eqn_idx%cont%beg)%sf(id1, id2, id3)
1389
1390 if (bubble_model == 2) then
1391 n_tait = 1._wp/gammas(1) + 1._wp
1392 b_tait = pi_infs(1)*(n_tait - 1)/n_tait
1393 c = n_tait*(pres + b_tait)*(1._wp - alf)/(rho)
1394 c = merge(sqrt(c), sgm_eps, c > 0._wp)
1395 end if
1396
1397 call s_coeff_selector(pres, rho, c, coeff, polytropic)
1398
1399 if (alf > small_alf) then
1400 nbub = q_cons_vf(eqn_idx%bub%beg)%sf(id1, id2, id3)
1401
1402# 784 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1403#if defined(MFC_OpenACC)
1404# 784 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1405!$acc loop seq
1406# 784 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1407#elif defined(MFC_OpenMP)
1408# 784 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1409
1410# 784 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1411#endif
1412 do q = 1, nb
1413 ! Gather moments for this bubble bin
1414
1415# 787 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1416#if defined(MFC_OpenACC)
1417# 787 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1418!$acc loop seq
1419# 787 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1420#elif defined(MFC_OpenMP)
1421# 787 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1422
1423# 787 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1424#endif
1425 do r = 2, nmom
1426 moms(r) = q_prim_vf(bubmoms(q, r))%sf(id1, id2, id3)
1427 end do
1428 moms(1) = 1._wp
1429 call s_chyqmom(moms, wght(:,q), abscx(:,q), abscy(:,q))
1430
1431 if (polytropic) then
1432
1433# 795 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1434#if defined(MFC_OpenACC)
1435# 795 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1436!$acc loop seq
1437# 795 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1438#elif defined(MFC_OpenMP)
1439# 795 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1440
1441# 795 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1442#endif
1443 do j = 1, nnode
1444 wght_pb(j, q) = wght(j, q)*(pb0(q) - pv)
1445 end do
1446 else
1447
1448# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1449#if defined(MFC_OpenACC)
1450# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1451!$acc loop seq
1452# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1453#elif defined(MFC_OpenMP)
1454# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1455
1456# 800 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1457#endif
1458 do j = 1, nnode
1459 chi_vw = 1._wp/(1._wp + r_v/r_g*(pb(id1, id2, id3, j, q)/pv - 1._wp))
1460 x_vw = m_g*chi_vw/(m_v + (m_g - m_v)*chi_vw)
1461 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 &
1462 & + 1._wp - x_vw)
1463 rho_mw = pv/(chi_vw*r_v*tw)
1464 rhs_mv(id1, id2, id3, j, q) = -re_trans_c(q)*((mv(id1, id2, id3, j, q)/(mv(id1, id2, id3, j, &
1465 & q) + mass_g0(q))) - chi_vw)
1466 rhs_mv(id1, id2, id3, j, q) = rho_mw*rhs_mv(id1, id2, id3, j, &
1467 & q)/pe_c/(1._wp - chi_vw)/abscx(j, q)
1468 grad_t = -re_trans_t(q)*((pb(id1, id2, id3, j, q)/pb0(q))*(abscx(j, &
1469 & q)/r0(q))**3*(mass_g0(q) + mass_v0(q))/(mass_g0(q) + mv(id1, id2, id3, &
1470 & j, q)) - 1._wp)
1471 ht(j, q) = pb0(q)*k_mw*grad_t/pe_t(q)/abscx(j, q)
1472 wght_pb(j, q) = wght(j, q)*(pb(id1, id2, id3, j, q))
1473 wght_mv(j, q) = wght(j, q)*(rhs_mv(id1, id2, id3, j, q))
1474 wght_ht(j, q) = wght(j, q)*ht(j, q)
1475 end do
1476 end if
1477
1478 ! Compute change in moments due to bubble dynamics
1479 r = 1
1480
1481# 823 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1482#if defined(MFC_OpenACC)
1483# 823 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1484!$acc loop seq
1485# 823 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1486#elif defined(MFC_OpenMP)
1487# 823 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1488
1489# 823 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1490#endif
1491 do i2 = 0, 2
1492
1493# 825 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1494#if defined(MFC_OpenACC)
1495# 825 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1496!$acc loop seq
1497# 825 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1498#elif defined(MFC_OpenMP)
1499# 825 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1500
1501# 825 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1502#endif
1503 do i1 = 0, 2
1504 if ((i1 + i2) <= 2) then
1505 momsum = 0._wp
1506
1507# 829 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1508#if defined(MFC_OpenACC)
1509# 829 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1510!$acc loop seq
1511# 829 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1512#elif defined(MFC_OpenMP)
1513# 829 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1514
1515# 829 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1516#endif
1517 do j = 1, nterms
1518 select case (bubble_model)
1519 case (3)
1520 if (j == 3) then
1521 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1522 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_pb(:,q), &
1523 & momrhs(:,i1, i2, j, q))
1524 else
1525 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1526 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght(:,q), &
1527 & momrhs(:,i1, i2, j, q))
1528 end if
1529 case (2)
1530 if ((j >= 7 .and. j <= 9) .or. (j >= 22 .and. j <= 23) .or. (j >= 10 &
1531 & .and. j <= 11) .or. (j == 26)) then
1532 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1533 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_pb(:,q), &
1534 & momrhs(:,i1, i2, j, q))
1535 else if ((j >= 27 .and. j <= 29) .and. (.not. polytropic)) then
1536 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1537 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_mv(:,q), &
1538 & momrhs(:,i1, i2, j, q))
1539 else if ((j >= 30 .and. j <= 32) .and. (.not. polytropic)) then
1540 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1541 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght_ht(:,q), &
1542 & momrhs(:,i1, i2, j, q))
1543 else
1544 momsum = momsum + coeff(j, i1, i2)*(r0(q)**momrhs(3, i1, i2, j, &
1545 & q))*f_quad2d(abscx(:,q), abscy(:,q), wght(:,q), &
1546 & momrhs(:,i1, i2, j, q))
1547 end if
1548 end select
1549 end do
1550 moms3d(i1, i2, q)%sf(id1, id2, id3) = nbub*momsum
1551 msum(r) = momsum
1552 r = r + 1
1553 end if
1554 end do
1555 end do
1556
1557 ! Compute change in pb and mv for non-polytropic model
1558 if (.not. polytropic) then
1559
1560# 872 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1561#if defined(MFC_OpenACC)
1562# 872 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1563!$acc loop seq
1564# 872 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1565#elif defined(MFC_OpenMP)
1566# 872 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1567
1568# 872 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1569#endif
1570 do j = 1, nnode
1571 drdt = msum(2)
1572 drdt2 = merge(-1._wp, 1._wp, j == 1 .or. j == 2)/(2._wp*sqrt(merge(moms(4) - moms(2)**2._wp, &
1573 & sgm_eps, moms(4) - moms(2)**2._wp > 0._wp)))
1574 drdt2 = drdt2*(msum(3) - 2._wp*moms(2)*msum(2))
1575 drdt = drdt + drdt2
1576 rhs_pb(id1, id2, id3, j, q) = (-3._wp*gam*drdt/abscx(j, q))*(pb(id1, id2, id3, j, q))
1577 rhs_pb(id1, id2, id3, j, q) = rhs_pb(id1, id2, id3, j, q) + (3._wp*gam/abscx(j, &
1578 & q))*rhs_mv(id1, id2, id3, j, q)*r_v*tw
1579 rhs_pb(id1, id2, id3, j, q) = rhs_pb(id1, id2, id3, j, q) + (3._wp*gam/abscx(j, q))*ht(j, q)
1580 rhs_mv(id1, id2, id3, j, q) = rhs_mv(id1, id2, id3, j, q)*(4._wp*pi*abscx(j, q)**2._wp)
1581 end do
1582 end if
1583 end do
1584
1585 ! Compute special high-order moments
1586 momsp(1)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght, 3._wp, 0._wp, 0._wp)
1587 momsp(2)%sf(id1, id2, id3) = 4._wp*pi*nbub*f_quad(abscx, abscy, wght, 2._wp, 1._wp, 0._wp)
1588 momsp(3)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght, 3._wp, 2._wp, 0._wp)
1589 if (abs(gam - 1._wp) <= 1.e-4_wp) then
1590 momsp(4)%sf(id1, id2, id3) = 1._wp
1591 else
1592 if (polytropic) then
1593 momsp(4)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght_pb, 3._wp*(1._wp - gam), 0._wp, &
1594 & 3._wp*gam) + pv*f_quad(abscx, abscy, wght, 3._wp, 0._wp, &
1595 & 0._wp) - 4._wp*re_inv*f_quad(abscx, abscy, wght, 2._wp, 1._wp, &
1596 & 0._wp) - (2._wp/web)*f_quad(abscx, abscy, wght, 2._wp, 0._wp, 0._wp)
1597 else
1598 momsp(4)%sf(id1, id2, id3) = f_quad(abscx, abscy, wght_pb, 3._wp, 0._wp, &
1599 & 0._wp) - 4._wp*re_inv*f_quad(abscx, abscy, wght, 2._wp, 1._wp, &
1600 & 0._wp) - (2._wp/web)*f_quad(abscx, abscy, wght, 2._wp, 0._wp, 0._wp)
1601 end if
1602 end if
1603 else
1604
1605# 907 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1606#if defined(MFC_OpenACC)
1607# 907 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1608!$acc loop seq
1609# 907 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1610#elif defined(MFC_OpenMP)
1611# 907 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1612
1613# 907 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1614#endif
1615 do q = 1, nb
1616
1617# 909 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1618#if defined(MFC_OpenACC)
1619# 909 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1620!$acc loop seq
1621# 909 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1622#elif defined(MFC_OpenMP)
1623# 909 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1624
1625# 909 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1626#endif
1627 do i1 = 0, 2
1628
1629# 911 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1630#if defined(MFC_OpenACC)
1631# 911 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1632!$acc loop seq
1633# 911 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1634#elif defined(MFC_OpenMP)
1635# 911 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1636
1637# 911 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1638#endif
1639 do i2 = 0, 2
1640 moms3d(i1, i2, q)%sf(id1, id2, id3) = 0._wp
1641 end do
1642 end do
1643 end do
1644 momsp(1)%sf(id1, id2, id3) = 0._wp
1645 momsp(2)%sf(id1, id2, id3) = 0._wp
1646 momsp(3)%sf(id1, id2, id3) = 0._wp
1647 momsp(4)%sf(id1, id2, id3) = 0._wp
1648 end if
1649 end do
1650 end do
1651 end do
1652
1653# 925 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1654#if defined(MFC_OpenACC)
1655# 925 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1656!$acc end parallel loop
1657# 925 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1658#elif defined(MFC_OpenMP)
1659# 925 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1660
1661# 925 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1662!$omp end target teams loop
1663# 925 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1664#endif
1665
1666 contains
1667 !> Select the polytropic or non-polytropic coefficient routine
1668 subroutine s_coeff_selector(pres, rho, c, coeff, polytropic)
1669
1670
1671# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1672#ifdef _CRAYFTN
1673# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1674#if MFC_OpenACC
1675# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1676!$acc routine seq
1677# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1678#elif MFC_OpenMP
1679# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1680
1681# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1682
1683# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1684!$omp declare target device_type(any)
1685# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1686#else
1687# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1688!DIR$ INLINEALWAYS s_coeff_selector
1689# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1690#endif
1691# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1692#elif MFC_OpenACC
1693# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1694!$acc routine seq
1695# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1696#elif MFC_OpenMP
1697# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1698
1699# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1700
1701# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1702!$omp declare target device_type(any)
1703# 931 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1704#endif
1705 real(wp), intent(in) :: pres, rho, c
1706# 936 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1707 real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeff
1708# 938 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1709 logical, intent(in) :: polytropic
1710 if (polytropic) then
1711 call s_coeff(pres, rho, c, coeff)
1712 else
1713 call s_coeff_nonpoly(pres, rho, c, coeff)
1714 end if
1715
1716 end subroutine s_coeff_selector
1717
1718 !> Perform CHyQMOM inversion for bivariate moments
1719 subroutine s_chyqmom(momin, wght, abscX, abscY)
1720
1721
1722# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1723#ifdef _CRAYFTN
1724# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1725#if MFC_OpenACC
1726# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1727!$acc routine seq
1728# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1729#elif MFC_OpenMP
1730# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1731
1732# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1733
1734# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1735!$omp declare target device_type(any)
1736# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1737#else
1738# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1739!DIR$ INLINEALWAYS s_chyqmom
1740# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1741#endif
1742# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1743#elif MFC_OpenACC
1744# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1745!$acc routine seq
1746# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1747#elif MFC_OpenMP
1748# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1749
1750# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1751
1752# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1753!$omp declare target device_type(any)
1754# 950 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1755#endif
1756
1757 real(wp), dimension(nmom), intent(in) :: momin
1758 real(wp), dimension(nnode), intent(inout) :: wght, abscX, abscY
1759
1760 ! Local variables
1761 real(wp), dimension(0:2,0:2) :: moms
1762 real(wp), dimension(3) :: M1, M3
1763 real(wp), dimension(2) :: myrho, myrho3, up, up3, Vf
1764 real(wp) :: bu, bv, d20, d11, d_02, c20, c11, c02
1765 real(wp) :: mu2, vp21, vp22, rho21, rho22
1766
1767 ! Assign moments to 2D array for clarity
1768 moms(0, 0) = momin(1)
1769 moms(1, 0) = momin(2)
1770 moms(0, 1) = momin(3)
1771 moms(2, 0) = momin(4)
1772 moms(1, 1) = momin(5)
1773 moms(0, 2) = momin(6)
1774
1775 ! Compute means and central moments
1776 bu = moms(1, 0)/moms(0, 0)
1777 bv = moms(0, 1)/moms(0, 0)
1778 d20 = moms(2, 0)/moms(0, 0)
1779 d11 = moms(1, 1)/moms(0, 0)
1780 d_02 = moms(0, 2)/moms(0, 0)
1781
1782 c20 = d20 - bu**2._wp
1783 c11 = d11 - bu*bv
1784 c02 = d_02 - bv**2._wp
1785
1786 ! First 1D quadrature (X direction)
1787 m1 = (/1._wp, 0._wp, c20/)
1788 call s_hyqmom(myrho, up, m1)
1789 vf = c11*up/c20
1790
1791 ! Second 1D quadrature (Y direction, conditional on X)
1792 mu2 = max(0._wp, c02 - sum(myrho*(vf**2._wp)))
1793 m3 = (/1._wp, 0._wp, mu2/)
1794 call s_hyqmom(myrho3, up3, m3)
1795
1796 ! Assign roots and weights for 2D quadrature
1797 vp21 = up3(1)
1798 vp22 = up3(2)
1799 rho21 = myrho3(1)
1800 rho22 = myrho3(2)
1801
1802 ! Compute weights (vectorized)
1803 wght = moms(0, 0)*[myrho(1)*rho21, myrho(1)*rho22, myrho(2)*rho21, myrho(2)*rho22]
1804
1805 ! Compute abscissas (vectorized)
1806 abscx = bu + [up(1), up(1), up(2), up(2)]
1807 abscy = bv + [vf(1) + vp21, vf(1) + vp22, vf(2) + vp21, vf(2) + vp22]
1808
1809 end subroutine s_chyqmom
1810
1811 !> Perform HyQMOM inversion for univariate moments
1812 subroutine s_hyqmom(frho, fup, fmom)
1813
1814
1815# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1816#ifdef _CRAYFTN
1817# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1818#if MFC_OpenACC
1819# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1820!$acc routine seq
1821# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1822#elif MFC_OpenMP
1823# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1824
1825# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1826
1827# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1828!$omp declare target device_type(any)
1829# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1830#else
1831# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1832!DIR$ INLINEALWAYS s_hyqmom
1833# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1834#endif
1835# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1836#elif MFC_OpenACC
1837# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1838!$acc routine seq
1839# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1840#elif MFC_OpenMP
1841# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1842
1843# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1844
1845# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1846!$omp declare target device_type(any)
1847# 1009 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1848#endif
1849
1850 real(wp), dimension(2), intent(inout) :: frho, fup
1851 real(wp), dimension(3), intent(in) :: fmom
1852 real(wp) :: bu, d2, c2
1853
1854 bu = fmom(2)/fmom(1)
1855 d2 = fmom(3)/fmom(1)
1856 c2 = d2 - bu**2._wp
1857 frho(1) = fmom(1)/2._wp
1858 frho(2) = fmom(1)/2._wp
1859 c2 = maxval((/c2, sgm_eps/))
1860 fup(1) = bu - sqrt(c2)
1861 fup(2) = bu + sqrt(c2)
1862
1863 end subroutine s_hyqmom
1864
1865 !> Evaluate a weighted quadrature sum over all bubble size bins and nodes
1866 function f_quad(abscX, abscY, wght_in, q, r, s)
1867
1868
1869# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1870#if MFC_OpenACC
1871# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1872!$acc routine seq
1873# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1874#elif MFC_OpenMP
1875# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1876
1877# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1878
1879# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1880!$omp declare target device_type(any)
1881# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1882#endif
1883# 1033 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1884 real(wp), dimension(nnode, nb), intent(in) :: abscx, abscy, wght_in
1885# 1035 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1886 real(wp), intent(in) :: q, r, s
1887 real(wp) :: f_quad_rv, f_quad
1888 integer :: i, i1
1889
1890 f_quad = 0._wp
1891
1892# 1040 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1893#if defined(MFC_OpenACC)
1894# 1040 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1895!$acc loop seq
1896# 1040 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1897#elif defined(MFC_OpenMP)
1898# 1040 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1899
1900# 1040 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1901#endif
1902 do i = 1, nb
1903 f_quad_rv = 0._wp
1904
1905# 1043 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1906#if defined(MFC_OpenACC)
1907# 1043 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1908!$acc loop seq
1909# 1043 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1910#elif defined(MFC_OpenMP)
1911# 1043 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1912
1913# 1043 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1914#endif
1915 do i1 = 1, nnode
1916 f_quad_rv = f_quad_rv + wght_in(i1, i)*(abscx(i1, i)**q)*(abscy(i1, i)**r)
1917 end do
1918 f_quad = f_quad + weight(i)*(r0(i)**s)*f_quad_rv
1919 end do
1920
1921 end function f_quad
1922
1923 !> Evaluate a weighted 2D quadrature sum over quadrature nodes for a single size bin
1924 function f_quad2d(abscX, abscY, wght_in, pow)
1925
1926
1927# 1055 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1928#if MFC_OpenACC
1929# 1055 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1930!$acc routine seq
1931# 1055 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1932#elif MFC_OpenMP
1933# 1055 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1934
1935# 1055 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1936
1937# 1055 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1938!$omp declare target device_type(any)
1939# 1055 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1940#endif
1941# 1059 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1942 real(wp), dimension(nnode), intent(in) :: abscx, abscy, wght_in
1943# 1061 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1944 real(wp), dimension(3), intent(in) :: pow
1945 real(wp) :: f_quad2d
1946 integer :: i
1947
1948 f_quad2d = 0._wp
1949
1950# 1066 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1951#if defined(MFC_OpenACC)
1952# 1066 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1953!$acc loop seq
1954# 1066 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1955#elif defined(MFC_OpenMP)
1956# 1066 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1957
1958# 1066 "/home/runner/work/MFC/MFC/src/simulation/m_qbmm.fpp"
1959#endif
1960 do i = 1, nnode
1961 f_quad2d = f_quad2d + wght_in(i)*(abscx(i)**pow(1))*(abscy(i)**pow(2))
1962 end do
1963
1964 end function f_quad2d
1965
1966 end subroutine s_mom_inv
1967
1968end 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:1813
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:1669
s_chyqmom
subroutine s_chyqmom(momin, wght, abscx, abscy)
Perform CHyQMOM inversion for bivariate moments.
Definition m_qbmm.fpp.f90:1720
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:1867
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:1925
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:1045
m_global_parameters::re_inv
real(wp) re_inv
Inverse Reynolds number.
Definition m_global_parameters.fpp.f90:1046
m_global_parameters::polytropic
logical polytropic
Polytropic switch.
Definition m_global_parameters.fpp.f90:1074
m_global_parameters::nb
integer nb
Number of eq. bubble sizes.
Definition m_global_parameters.fpp.f90:1040
m_global_parameters::bubble_model
integer bubble_model
Gilmore or Keller–Miksis bubble model.
Definition m_global_parameters.fpp.f90:1104
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:1117
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:1329
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:1228
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