MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_helper.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
2# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
3! This file exists so that Fypp can be run without generating case.fpp files for
4! each target. This is useful when generating documentation, for example. This
5! should also let MFC be built with CMake directly, without invoking mfc.sh.
6
7! For pre-process.
8# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
9
10! For moving immersed boundaries in simulation
11# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
12# 2 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp" 2
13# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
14# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
15# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
16# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
19# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21
22# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39
40# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45! New line at end of file is required for FYPP
46# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
47# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
48# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
49# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54
55# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64
65# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66
67# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68
69# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
70
71# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
72
73# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
74
75# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
76
77# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
78! New line at end of file is required for FYPP
79# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
80
81# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
82# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
83# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
84# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
85# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86
87# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88
89# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
94
95# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
96
97# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
98
99# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
100
101# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102
103# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104
105# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
106
107# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145
146# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147
148# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149
150# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
151
152# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
153
154# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
155
156# 403 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
157! New line at end of file is required for FYPP
158# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
159# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
160# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
161# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
164# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
165# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166
167# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170
171# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
173# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
174
175# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
176
177# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
178
179# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
180
181# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
182
183# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
184
185# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
186
187# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
188
189# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
190! New line at end of file is required for FYPP
191# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
192
193# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234
235# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
236
237# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
238
239# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
240
241# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
242
243# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
244
245# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
246
247# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
248! New line at end of file is required for FYPP
249# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
250
251! GPU parallel region (scalar reductions, maxval/minval)
252# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
253
254! GPU parallel loop over threads (most common GPU macro)
255# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257! Required closing for GPU_PARALLEL_LOOP
258# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
259
260! Mark routine for device compilation
261# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263! Declare device-resident data
264# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
265
266! Inner loop within a GPU parallel region
267# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
268
269! Scoped GPU data region
270# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
271
272! Host code with device pointers (for MPI with GPU buffers)
273# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
274
275! Allocate device memory (unscoped)
276# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
277
278! Free device memory
279# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
280
281! Atomic operation on device
282# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
283
284! End atomic capture block
285# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
286
287! Copy data between host and device
288# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
289
290! Synchronization barrier
291# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
292
293! Import GPU library module (openacc or omp_lib)
294# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
295
296! Emit code only for AMD compiler
297# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
298
299! Emit code for non-Cray compilers
300# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
301
302! Emit code only for Cray compiler
303# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
304
305! Emit code for non-NVIDIA compilers
306# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
307
308# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
309# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
310! New line at end of file is required for FYPP
311# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
312
313# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
314
315! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
316! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
317! example see misc/nvidia_uvm/bind.sh.
318# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319
320! Allocate and create GPU device memory
321# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
322
323! Free GPU device memory and deallocate
324# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
325
326! Cray-specific GPU pointer setup for vector fields
327# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
328
329! Cray-specific GPU pointer setup for scalar fields
330# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
331
332! Cray-specific GPU pointer setup for acoustic source spatials
333# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
334
335# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
336
337# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
338! New line at end of file is required for FYPP
339# 3 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp" 2
340
341!>
342!! @file
343!! @brief Contains module m_helper
344
345!> @brief Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions
347
350 use m_constants, only: bc_periodic
351 use ieee_arithmetic !< for checking nan
352
353 implicit none
354
355 private
361
362contains
363
364 !> Computes the bubble number density n from the primitive variables
365 subroutine s_comp_n_from_prim(vftmp, Rtmp, ntmp, weights)
366
367
368# 30 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
369#if MFC_OpenACC
370# 30 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
371!$acc routine seq
372# 30 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
373#elif MFC_OpenMP
374# 30 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
375
376# 30 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
377
378# 30 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
379!$omp declare target device_type(any)
380# 30 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
381#endif
382 real(wp), intent(in) :: vftmp
383 real(wp), dimension(nb), intent(in) :: rtmp
384 real(wp), intent(out) :: ntmp
385 real(wp), dimension(nb), intent(in) :: weights
386 real(wp) :: r3
387
388 r3 = dot_product(weights, rtmp**3._wp)
389 ntmp = (3._wp/(4._wp*pi))*vftmp/r3
390
391 end subroutine s_comp_n_from_prim
392
393 !> Compute the bubble number density from the conservative void fraction and weighted bubble radii.
394 subroutine s_comp_n_from_cons(vftmp, nRtmp, ntmp, weights)
395
396
397# 45 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
398#if MFC_OpenACC
399# 45 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
400!$acc routine seq
401# 45 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
402#elif MFC_OpenMP
403# 45 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
404
405# 45 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
406
407# 45 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
408!$omp declare target device_type(any)
409# 45 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
410#endif
411 real(wp), intent(in) :: vftmp
412 real(wp), dimension(nb), intent(in) :: nrtmp
413 real(wp), intent(out) :: ntmp
414 real(wp), dimension(nb), intent(in) :: weights
415 real(wp) :: nr3
416
417 nr3 = dot_product(weights, nrtmp**3._wp)
418 ntmp = sqrt((4._wp*pi/3._wp)*nr3/vftmp)
419
420 end subroutine s_comp_n_from_cons
421
422 !> Print a 2D real array to standard output, optionally dividing each element by a given scalar.
423 impure subroutine s_print_2d_array(A, div)
424
425 real(wp), dimension(:,:), intent(in) :: a
426 real(wp), optional, intent(in) :: div
427 integer :: i, j
428 integer :: local_m, local_n
429 real(wp) :: c
430
431 local_m = size(a, 1)
432 local_n = size(a, 2)
433
434 if (present(div)) then
435 c = div
436 else
437 c = 1._wp
438 end if
439
440 print *, local_m, local_n
441
442 do i = 1, local_m
443 do j = 1, local_n
444 write (*, fmt="(F12.4)", advance="no") a(i, j)/c
445 end do
446 write (*, fmt="(A1)") " "
447 end do
448 write (*, fmt="(A1)") " "
449
450 end subroutine s_print_2d_array
451
452 !> Initialize bubble model arrays for Euler or Lagrangian bubbles with polytropic or non-polytropic gas.
453 impure subroutine s_initialize_bubbles_model()
454
455 ! Allocate memory
456 if (bubbles_euler) then
457#ifdef MFC_DEBUG
458# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
459 block
460# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
461 use iso_fortran_env, only: output_unit
462# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
463
464# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
465 print *, 'm_helper.fpp:92: ', '@:ALLOCATE(weight(nb), R0(nb))'
466# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
467
468# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
469 call flush (output_unit)
470# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
471 end block
472# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
473#endif
474# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
475 allocate (weight(nb), r0(nb))
476# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
477
478# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
479
480# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
481
482# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
483#if defined(MFC_OpenACC)
484# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
485!$acc enter data create(weight, R0)
486# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
487#elif defined(MFC_OpenMP)
488# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
489!$omp target enter data map(always,alloc:weight, R0)
490# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
491#endif
492 if (.not. polytropic) then
493#ifdef MFC_DEBUG
494# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
495 block
496# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
497 use iso_fortran_env, only: output_unit
498# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
499
500# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
501 print *, 'm_helper.fpp:94: ', '@:ALLOCATE(pb0(nb), Pe_T(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))'
502# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
503
504# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
505 call flush (output_unit)
506# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
507 end block
508# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
509#endif
510# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
511 allocate (pb0(nb), pe_t(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))
512# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
513
514# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
515
516# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
517
518# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
519
520# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
521
522# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
523
524# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
525
526# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
527#if defined(MFC_OpenACC)
528# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
529!$acc enter data create(pb0, Pe_T, k_g, k_v, mass_g0, mass_v0)
530# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
531#elif defined(MFC_OpenMP)
532# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
533!$omp target enter data map(always,alloc:pb0, Pe_T, k_g, k_v, mass_g0, mass_v0)
534# 94 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
535#endif
536#ifdef MFC_DEBUG
537# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
538 block
539# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
540 use iso_fortran_env, only: output_unit
541# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
542
543# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
544 print *, 'm_helper.fpp:95: ', '@:ALLOCATE(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb))'
545# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
546
547# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
548 call flush (output_unit)
549# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
550 end block
551# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
552#endif
553# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
554 allocate (re_trans_t(nb), re_trans_c(nb), im_trans_t(nb), im_trans_c(nb))
555# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
556
557# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
558
559# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
560
561# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
562
563# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
564
565# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
566#if defined(MFC_OpenACC)
567# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
568!$acc enter data create(Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c)
569# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
570#elif defined(MFC_OpenMP)
571# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
572!$omp target enter data map(always,alloc:Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c)
573# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
574#endif
575 else if (qbmm) then
576#ifdef MFC_DEBUG
577# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
578 block
579# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
580 use iso_fortran_env, only: output_unit
581# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
582
583# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
584 print *, 'm_helper.fpp:97: ', '@:ALLOCATE(pb0(nb))'
585# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
586
587# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
588 call flush (output_unit)
589# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
590 end block
591# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
592#endif
593# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
594 allocate (pb0(nb))
595# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
596
597# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
598
599# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
600#if defined(MFC_OpenACC)
601# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
602!$acc enter data create(pb0)
603# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
604#elif defined(MFC_OpenMP)
605# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
606!$omp target enter data map(always,alloc:pb0)
607# 97 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
608#endif
609 end if
610
611 ! Compute quadrature weights and nodes for polydisperse simulations
612 if (nb > 1) then
613 call s_simpson(weight, r0)
614 else if (nb == 1) then
615 r0 = 1._wp
616 weight = 1._wp
617 else
618 stop 'Invalid value of nb'
619 end if
620 r0 = r0*bub_pp%R0ref
621 end if
622
623 ! Initialize bubble variables
625
626 end subroutine s_initialize_bubbles_model
627
628 !> Set bubble physical parameters and nondimensional numbers from the input configuration.
629 impure subroutine s_initialize_bubble_vars()
630
631 r0ref = bub_pp%R0ref; p0ref = bub_pp%p0ref
632 rho0ref = bub_pp%rho0ref
633 ss = bub_pp%ss; pv = bub_pp%pv; vd = bub_pp%vd
634 mu_l = bub_pp%mu_l; mu_v = bub_pp%mu_v; mu_g = bub_pp%mu_g
635 gam_v = bub_pp%gam_v; gam_g = bub_pp%gam_g
636 if (.not. polytropic) then
637 if (bubbles_euler) then
638 m_v = bub_pp%M_v; m_g = bub_pp%M_g
639 k_v = bub_pp%k_v; k_g = bub_pp%k_g
640 end if
641 r_v = bub_pp%R_v; r_g = bub_pp%R_g
642 tw = bub_pp%T0ref
643 end if
644 if (bubbles_lagrange) then
645 cp_v = bub_pp%cp_v; cp_g = bub_pp%cp_g
646 k_vl = bub_pp%k_v; k_gl = bub_pp%k_g
647 end if
648
649 ! Input quantities
650 if (bubbles_euler .and. (.not. polytropic)) then
651 if (thermal == 2) then
652 gam_m = 1._wp
653 else
654 gam_m = gam_g
655 end if
656 end if
657
658 ! Nondimensional numbers
659 eu = p0ref
660 ca = eu - pv
661 if (.not. f_is_default(bub_pp%ss)) web = 1._wp/ss
662 if (.not. f_is_default(bub_pp%mu_l)) re_inv = mu_l
663 if (.not. polytropic) pe_c = 1._wp/vd
664
665 if (bubbles_euler) then
666 ! Initialize variables for non-polytropic (Preston) model
667 if (.not. polytropic) then
669 end if
670 ! Initialize pb based on surface tension for qbmm (polytropic)
671 if (qbmm .and. polytropic) then
672 pb0 = eu
673 if (.not. f_is_default(web)) then
674 pb0 = pb0 + 2._wp/web/r0
675 end if
676 end if
677 end if
678
679 end subroutine s_initialize_bubble_vars
680
681 !> Initializes non-polydisperse bubble modeling
682 impure subroutine s_initialize_nonpoly()
683
684 integer :: ir
685 real(wp), dimension(nb) :: chi_vw0, cp_m0, k_m0, rho_m0, x_vw, omegan
686 real(wp), parameter :: k_poly = 1._wp !< polytropic index used to compute isothermal natural frequency
687 ! Chapman-Enskog transport coefficients for vapor-gas mixture, Ando JAS (2010)
688
689 phi_vg = (1._wp + sqrt(mu_v/mu_g)*(m_g/m_v)**(0.25_wp))**2/(sqrt(8._wp)*sqrt(1._wp + m_v/m_g))
690 phi_gv = (1._wp + sqrt(mu_g/mu_v)*(m_v/m_g)**(0.25_wp))**2/(sqrt(8._wp)*sqrt(1._wp + m_g/m_v))
691
692 ! Initial internal bubble pressure (Euler number + Laplace pressure)
693 pb0 = eu + 2._wp/web/r0
694
695 ! Vapor mass fraction at bubble wall, Ando JAS (2010)
696 chi_vw0 = 1._wp/(1._wp + r_v/r_g*(pb0/pv - 1._wp))
697
698 ! Mixture specific heat from mass-weighted vapor/gas contributions
699 cp_m0 = chi_vw0*r_v*gam_v/(gam_v - 1._wp) + (1._wp - chi_vw0)*r_g*gam_g/(gam_g - 1._wp)
700
701 ! mole fraction of vapor (Eq. 2.23 in Ando 2010)
702 x_vw = m_g*chi_vw0/(m_v + (m_g - m_v)*chi_vw0)
703
704 ! thermal conductivity for gas/vapor mixture (Eq. 2.21 in Ando 2010)
705 k_m0 = x_vw*k_v/(x_vw + (1._wp - x_vw)*phi_vg) + (1._wp - x_vw)*k_g/(x_vw*phi_gv + 1._wp - x_vw)
706 k_g(:) = k_g(:)/k_m0(:)
707 k_v(:) = k_v(:)/k_m0(:)
708
709 ! mixture density (Eq. 2.20 in Ando 2010)
710 rho_m0 = pv/(chi_vw0*r_v*tw)
711
712 ! mass of gas/vapor
713 mass_g0(:) = (4._wp*pi/3._wp)*(pb0(:) - pv)/(r_g*tw)*r0(:)**3
714 mass_v0(:) = (4._wp*pi/3._wp)*pv/(r_v*tw)*r0(:)**3
715
716 ! Peclet numbers
717 pe_t(:) = rho_m0*cp_m0(:)/k_m0(:)
718
719 ! Bubble natural frequency, Ando JAS (2010)
720 omegan(:) = sqrt(3._wp*k_poly*ca + 2._wp*(3._wp*k_poly - 1._wp)/(web*r0))/r0/sqrt(rho0ref)
721 do ir = 1, nb
722 call s_transcoeff(omegan(ir)*r0(ir), pe_t(ir)*r0(ir), re_trans_t(ir), im_trans_t(ir))
723 call s_transcoeff(omegan(ir)*r0(ir), pe_c*r0(ir), re_trans_c(ir), im_trans_c(ir))
724 end do
725 im_trans_t = 0._wp
726
727 end subroutine s_initialize_nonpoly
728
729 !> Computes the transfer coefficient for the non-polytropic bubble compression process
730 elemental subroutine s_transcoeff(omega, peclet, Re_trans, Im_trans)
731
732 real(wp), intent(in) :: omega, peclet
733 real(wp), intent(out) :: re_trans, im_trans
734 complex(wp) :: imag, trans, c1, c2, c3
735
736 imag = (0._wp, 1._wp)
737
738 c1 = imag*omega*peclet
739 c2 = sqrt(c1)
740 c3 = (exp(c2) - exp(-c2))/(exp(c2) + exp(-c2)) ! TANH(c2)
741 trans = ((c2/c3 - 1._wp)**(-1) - 3._wp/c1)**(-1) ! transfer function
742
743 re_trans = trans
744 im_trans = aimag(trans)
745
746 end subroutine s_transcoeff
747
748 !> Convert an integer to its trimmed string representation.
749 elemental subroutine s_int_to_str(i, res)
750
751 integer, intent(in) :: i
752 character(len=*), intent(inout) :: res
753
754 write (res, '(I0)') i
755 res = trim(res)
756
757 end subroutine s_int_to_str
758
759 !> Compute the Simpson weights for quadrature
760 subroutine s_simpson(local_weight, local_R0)
761
762 real(wp), dimension(:), intent(inout) :: local_weight
763 real(wp), dimension(:), intent(inout) :: local_r0
764 integer :: ir
765 real(wp) :: r0mn, r0mx, dphi, tmp, sd
766 real(wp), dimension(nb) :: phi
767
768 sd = poly_sigma
769 r0mn = 0.8_wp*exp(-2.8_wp*sd)
770 r0mx = 0.2_wp*exp(9.5_wp*sd) + 1._wp
771
772 ! phi = ln( R0 ) & return R0
773 do ir = 1, nb
774 phi(ir) = log(r0mn) + (ir - 1._wp)*log(r0mx/r0mn)/(nb - 1._wp)
775 local_r0(ir) = exp(phi(ir))
776 end do
777
778# 268 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
779 dphi = phi(2) - phi(1)
780# 270 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
781
782 ! weights for quadrature using Simpson's rule
783 do ir = 2, nb - 1
784 ! Gaussian
785 tmp = exp(-0.5_wp*(phi(ir)/sd)**2)/sqrt(2._wp*pi)/sd
786 if (mod(ir, 2) == 0) then
787 local_weight(ir) = tmp*4._wp*dphi/3._wp
788 else
789 local_weight(ir) = tmp*2._wp*dphi/3._wp
790 end if
791 end do
792 tmp = exp(-0.5_wp*(phi(1)/sd)**2)/sqrt(2._wp*pi)/sd
793 local_weight(1) = tmp*dphi/3._wp
794 tmp = exp(-0.5_wp*(phi(nb)/sd)**2)/sqrt(2._wp*pi)/sd
795 local_weight(nb) = tmp*dphi/3._wp
796
797 end subroutine s_simpson
798
799 !> Compute the cross product of two vectors.
800 pure function f_cross(a, b) result(c)
801
802
803# 291 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
804#if MFC_OpenACC
805# 291 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
806!$acc routine seq
807# 291 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
808#elif MFC_OpenMP
809# 291 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
810
811# 291 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
812
813# 291 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
814!$omp declare target device_type(any)
815# 291 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
816#endif
817
818 real(wp), dimension(3), intent(in) :: a, b
819 real(wp), dimension(3) :: c
820
821 c(1) = a(2)*b(3) - a(3)*b(2)
822 c(2) = a(3)*b(1) - a(1)*b(3)
823 c(3) = a(1)*b(2) - a(2)*b(1)
824
825 end function f_cross
826
827 !> Generate a unit vector uniformly distributed on the sphere from two random parameters.
828 function f_unit_vector(theta, eta) result(vec)
829
830
831# 305 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
832#if MFC_OpenACC
833# 305 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
834!$acc routine seq
835# 305 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
836#elif MFC_OpenMP
837# 305 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
838
839# 305 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
840
841# 305 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
842!$omp declare target device_type(any)
843# 305 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
844#endif
845 real(wp), intent(in) :: theta, eta
846 real(wp) :: zeta, xi
847 real(wp), dimension(3) :: vec
848
849 xi = 2._wp*pi*theta
850 zeta = acos(2._wp*eta - 1._wp)
851 vec(1) = sin(zeta)*cos(xi)
852 vec(2) = sin(zeta)*sin(xi)
853 vec(3) = cos(zeta)
854
855 end function f_unit_vector
856
857 !> Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
858 subroutine s_prng(var, seed)
859
860
861# 321 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
862#if MFC_OpenACC
863# 321 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
864!$acc routine seq
865# 321 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
866#elif MFC_OpenMP
867# 321 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
868
869# 321 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
870
871# 321 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
872!$omp declare target device_type(any)
873# 321 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
874#endif
875 integer, intent(inout) :: seed
876 real(wp), intent(out) :: var
877
878 seed = mod(modmul(seed), modulus)
879 var = seed/real(modulus, wp)
880
881 end subroutine s_prng
882
883 !> Compute a modular multiplication step for the linear congruential pseudo-random number generator.
884 function modmul(a) result(val)
885
886
887# 333 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
888#if MFC_OpenACC
889# 333 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
890!$acc routine seq
891# 333 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
892#elif MFC_OpenMP
893# 333 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
894
895# 333 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
896
897# 333 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
898!$omp declare target device_type(any)
899# 333 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
900#endif
901 integer, intent(in) :: a
902 integer :: val
903 real(wp) :: x, y
904
905 x = (multiplier/real(modulus, wp))*a + (increment/real(modulus, wp))
906 y = nint((x - floor(x))*decimal_trim)/decimal_trim
907 val = nint(y*modulus)
908
909 end function modmul
910
911 !> Compute the cross product c = a x b of two 3D vectors.
912 subroutine s_cross_product(a, b, c)
913
914
915# 347 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
916#if MFC_OpenACC
917# 347 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
918!$acc routine seq
919# 347 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
920#elif MFC_OpenMP
921# 347 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
922
923# 347 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
924
925# 347 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
926!$omp declare target device_type(any)
927# 347 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
928#endif
929 real(wp), intent(in) :: a(3), b(3)
930 real(wp), intent(out) :: c(3)
931
932 c(1) = a(2)*b(3) - a(3)*b(2)
933 c(2) = a(3)*b(1) - a(1)*b(3)
934 c(3) = a(1)*b(2) - a(2)*b(1)
935
936 end subroutine s_cross_product
937
938 !> Swap two real numbers.
939 elemental subroutine s_swap(lhs, rhs)
940
941 real(wp), intent(inout) :: lhs, rhs
942 real(wp) :: ltemp
943
944 ltemp = lhs
945 lhs = rhs
946 rhs = ltemp
947
948 end subroutine s_swap
949
950 !> Create a transformation matrix.
951 function f_create_transform_matrix(param, center) result(out_matrix)
952
953 type(ic_model_parameters), intent(in) :: param
954 real(wp), dimension(1:3), optional, intent(in) :: center
955 real(wp), dimension(1:4,1:4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
956
957 sc = transpose(reshape([param%scale(1), 0._wp, 0._wp, 0._wp, 0._wp, param%scale(2), 0._wp, 0._wp, 0._wp, 0._wp, &
958 & param%scale(3), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(sc)))
959
960 rz = transpose(reshape([cos(param%rotate(3)), -sin(param%rotate(3)), 0._wp, 0._wp, sin(param%rotate(3)), &
961 & cos(param%rotate(3)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp, 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(rz)))
962
963 rx = transpose(reshape([1._wp, 0._wp, 0._wp, 0._wp, 0._wp, cos(param%rotate(1)), -sin(param%rotate(1)), 0._wp, 0._wp, &
964 & sin(param%rotate(1)), cos(param%rotate(1)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(rx)))
965
966 ry = transpose(reshape([cos(param%rotate(2)), 0._wp, sin(param%rotate(2)), 0._wp, 0._wp, 1._wp, 0._wp, 0._wp, &
967 & -sin(param%rotate(2)), 0._wp, cos(param%rotate(2)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(ry)))
968
969 tr = transpose(reshape([1._wp, 0._wp, 0._wp, param%translate(1), 0._wp, 1._wp, 0._wp, param%translate(2), 0._wp, 0._wp, &
970 & 1._wp, param%translate(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
971
972 if (present(center)) then
973 ! Translation matrix to move center to the origin
974 t_to_origin = transpose(reshape([1._wp, 0._wp, 0._wp, -center(1), 0._wp, 1._wp, 0._wp, -center(2), 0._wp, 0._wp, &
975 & 1._wp, -center(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
976
977 ! Translation matrix to move center back to original position
978 t_back = transpose(reshape([1._wp, 0._wp, 0._wp, center(1), 0._wp, 1._wp, 0._wp, center(2), 0._wp, 0._wp, 1._wp, &
979 & center(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
980
981 out_matrix = matmul(tr, matmul(t_back, matmul(ry, matmul(rx, matmul(rz, matmul(sc, t_to_origin))))))
982 else
983 out_matrix = matmul(ry, matmul(rx, rz))
984 end if
985
986 end function f_create_transform_matrix
987
988 !> Transform a vector by a matrix.
989 subroutine s_transform_vec(vec, matrix)
990
991 real(wp), dimension(1:3), intent(inout) :: vec
992 real(wp), dimension(1:4,1:4), intent(in) :: matrix
993 real(wp), dimension(1:4) :: tmp
994
995 tmp = matmul(matrix, [vec(1), vec(2), vec(3), 1._wp])
996 vec = tmp(1:3)
997
998 end subroutine s_transform_vec
999
1000 !> Transform a triangle by a matrix, one vertex at a time.
1001 subroutine s_transform_triangle(triangle, matrix, matrix_n)
1002
1003 type(t_triangle), intent(inout) :: triangle
1004 real(wp), dimension(1:4,1:4), intent(in) :: matrix, matrix_n
1005 integer :: i
1006
1007 do i = 1, 3
1008 call s_transform_vec(triangle%v(i,:), matrix)
1009 end do
1010
1011 call s_transform_vec(triangle%n(1:3), matrix_n)
1012
1013 end subroutine s_transform_triangle
1014
1015 !> Transform a model by a matrix, one triangle at a time.
1016 subroutine s_transform_model(model, matrix, matrix_n)
1017
1018 type(t_model), intent(inout) :: model
1019 real(wp), dimension(1:4,1:4), intent(in) :: matrix, matrix_n
1020 integer :: i
1021
1022 do i = 1, size(model%trs)
1023 call s_transform_triangle(model%trs(i), matrix, matrix_n)
1024 end do
1025
1026 end subroutine s_transform_model
1027
1028 !> Create a bounding box for a model.
1029 function f_create_bbox(model) result(bbox)
1030
1031 type(t_model), intent(in) :: model
1032 type(t_bbox) :: bbox
1033 integer :: i, j
1034
1035 if (size(model%trs) == 0) then
1036 bbox%min = 0._wp
1037 bbox%max = 0._wp
1038 return
1039 end if
1040
1041 bbox%min = model%trs(1)%v(1,:)
1042 bbox%max = model%trs(1)%v(1,:)
1043
1044 do i = 1, size(model%trs)
1045 do j = 1, 3
1046 bbox%min = min(bbox%min, model%trs(i)%v(j,:))
1047 bbox%max = max(bbox%max, model%trs(i)%v(j,:))
1048 end do
1049 end do
1050
1051 end function f_create_bbox
1052
1053 !> Perform XOR on lhs and rhs.
1054 elemental function f_xor(lhs, rhs) result(res)
1055
1056 logical, intent(in) :: lhs, rhs
1057 logical :: res
1058
1059 res = (lhs .and. .not. rhs) .or. (.not. lhs .and. rhs)
1060
1061 end function f_xor
1062
1063 !> Convert a logical to 1 or 0.
1064 elemental function f_logical_to_int(predicate) result(int)
1065
1066 logical, intent(in) :: predicate
1067 integer :: int
1068
1069 if (predicate) then
1070 int = 1
1071 else
1072 int = 0
1073 end if
1074
1075 end function f_logical_to_int
1076
1077 !> Real spherical harmonic Y_lm(theta, phi). theta = polar angle from +z (acos(z/r)), phi = atan2(y,x). Uses associated Legendre
1078 !! P_l^|m|(cos theta). Standard normalisation.
1079 function real_ylm(theta, phi, l, m) result(Y)
1080
1081 integer, intent(in) :: l, m
1082 real(wp), intent(in) :: theta, phi
1083 real(wp) :: y, x, prefac
1084 integer :: m_abs
1085
1086 m_abs = abs(m)
1087 if (m_abs > l) then
1088 y = 0._wp
1089 return
1090 end if
1091 x = cos(theta)
1092 prefac = sqrt((2*l + 1)*real(factorial(l - m_abs), wp)/real(factorial(l + m_abs), wp)/(4._wp*pi))
1093 if (m == 0) then
1094 y = prefac*associated_legendre(x, l, 0)
1095 else if (m > 0) then
1096 y = prefac*sqrt(2._wp)*associated_legendre(x, l, m_abs)*cos(m*phi)
1097 else
1098 y = prefac*sqrt(2._wp)*associated_legendre(x, l, m_abs)*sin(m_abs*phi)
1099 end if
1100
1101 end function real_ylm
1102
1103 !> Associated Legendre polynomial P_l^m(x) (Ferrers function, Condon-Shortley phase). Valid for integer l >= 0, 0 <= m <= l, and
1104 !! x in [-1,1]. Returns 0 for |m| > l or l < 0. Formulas: DLMF 14.10.3 (recurrence in degree), Wikipedia "Associated Legendre
1105 !! polynomials" (P_l^l and P_l^{l-1} identities). Recurrence: (l-m)P_l^m = (2l-1)x P_{l-1}^m - (l+m-1)P_{l-2}^m.
1106 !! @param x argument (typically cos(theta)), should be in [-1,1]
1107 !! @param l degree (>= 0)
1108 !! @param m_order order (0 <= m_order <= l)
1109 recursive function associated_legendre(x, l, m_order) result(result_P)
1110
1111 integer, intent(in) :: l, m_order
1112 real(wp), intent(in) :: x
1113 real(wp) :: result_p
1114 real(wp) :: one_minus_x2
1115
1116 ! Out-of-domain: P_l^m = 0 for |m| > l or l < 0 (standard convention)
1117
1118 if (l < 0 .or. m_order < 0 .or. m_order > l) then
1119 result_p = 0._wp
1120 return
1121 end if
1122
1123 if (m_order <= 0 .and. l <= 0) then
1124 result_p = 1._wp
1125 else if (l == 1 .and. m_order <= 0) then
1126 result_p = x
1127 else if (l == 1 .and. m_order == 1) then
1128 one_minus_x2 = max(0._wp, 1._wp - x**2)
1129 result_p = -sqrt(one_minus_x2)
1130 else if (m_order == l) then
1131 ! P_l^l(x) = (-1)^l (2l-1)!! (1-x^2)^(l/2). Use real exponent for odd l
1132 one_minus_x2 = max(0._wp, 1._wp - x**2)
1133 result_p = (-1)**l*real(double_factorial(2*l - 1), wp)*one_minus_x2**(0.5_wp*real(l, wp))
1134 else if (m_order == l - 1) then
1135 result_p = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1)
1136 else
1137 result_p = ((2*l - 1)*x*associated_legendre(x, l - 1, m_order) - (l + m_order - 1)*associated_legendre(x, l - 2, &
1138 & m_order))/(l - m_order)
1139 end if
1140
1141 end function associated_legendre
1142
1143 !> Calculate the double factorial of an integer
1144 elemental function double_factorial(n_in) result(R_result)
1145
1146 integer, intent(in) :: n_in
1147 integer, parameter :: int64_kind = selected_int_kind(18) !< 18 bytes for 64-bit integer
1148 integer(kind=int64_kind) :: r_result
1149 integer :: i
1150
1151 r_result = product((/(i, i=n_in, 1, -2)/))
1152
1153 end function double_factorial
1154
1155 !> Calculate the factorial of an integer
1156 elemental function factorial(n_in) result(R_result)
1157
1158 integer, intent(in) :: n_in
1159 integer, parameter :: int64_kind = selected_int_kind(18) !< 18 bytes for 64-bit integer
1160 integer(kind=int64_kind) :: r_result
1161 integer :: i
1162
1163 r_result = product((/(i, i=n_in, 1, -1)/))
1164
1165 end function factorial
1166
1167 !> Calculate a smooth cut-on function that is zero for x values smaller than zero and goes to one, for generating smooth initial
1168 !! conditions
1169 function f_cut_on(x, eps) result(fx)
1170
1171 real(wp), intent(in) :: x, eps
1172 real(wp) :: fx
1173
1174 fx = 1 - f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
1175
1176 end function f_cut_on
1177
1178 !> Calculate a smooth cut-off function that is one for x values smaller than zero and goes to zero, for generating smooth
1179 !! initial conditions
1180 function f_cut_off(x, eps) result(fx)
1181
1182 real(wp), intent(in) :: x, eps
1183 real(wp) :: fx
1184
1185 fx = f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
1186
1187 end function f_cut_off
1188
1189 !> Helper function for f_cut_on and f_cut_off
1190 function f_gx(x) result(gx)
1191
1192 real(wp), intent(in) :: x
1193 real(wp) :: gx
1194
1195 if (x > 0) then
1196 gx = exp(-1._wp/x)
1197 else
1198 gx = 0._wp
1199 end if
1200
1201 end function f_gx
1202
1203 !> Downsample conservative variable fields by a factor of 3 in each direction using volume averaging.
1204 subroutine s_downsample_data(q_cons_vf, q_cons_temp, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
1205
1206 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf, q_cons_temp
1207
1208 ! Down sampling variables
1209 integer :: i, j, k, l
1210 integer :: ix, iy, iz, x_id, y_id, z_id
1211 integer, intent(inout) :: m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds
1212
1213 m_ds = int((m + 1)/3) - 1
1214 n_ds = int((n + 1)/3) - 1
1215 p_ds = int((p + 1)/3) - 1
1216
1217 m_glb_ds = int((m_glb + 1)/3) - 1
1218 n_glb_ds = int((n_glb + 1)/3) - 1
1219 p_glb_ds = int((p_glb + 1)/3) - 1
1220
1221 do l = -1, p_ds + 1
1222 do k = -1, n_ds + 1
1223 do j = -1, m_ds + 1
1224 x_id = 3*j + 1
1225 y_id = 3*k + 1
1226 z_id = 3*l + 1
1227 do i = 1, sys_size
1228 q_cons_temp(i)%sf(j, k, l) = 0
1229
1230 do iz = -1, 1
1231 do iy = -1, 1
1232 do ix = -1, 1
1233 q_cons_temp(i)%sf(j, k, l) = q_cons_temp(i)%sf(j, k, &
1234 & l) + (1._wp/27._wp)*q_cons_vf(i)%sf(x_id + ix, y_id + iy, z_id + iz)
1235 end do
1236 end do
1237 end do
1238 end do
1239 end do
1240 end do
1241 end do
1242
1243 end subroutine s_downsample_data
1244
1245 !> Upsample conservative variable fields from a coarsened grid back to the original resolution using interpolation.
1246 subroutine s_upsample_data(q_cons_vf, q_cons_temp)
1247
1248 type(scalar_field), intent(inout), dimension(sys_size) :: q_cons_vf, q_cons_temp
1249 integer :: i, j, k, l
1250 integer :: ix, iy, iz
1251 integer :: x_id, y_id, z_id
1252 real(wp), dimension(4) :: temp
1253
1254 do l = 0, p
1255 do k = 0, n
1256 do j = 0, m
1257 do i = 1, sys_size
1258 ix = int(j/3._wp)
1259 iy = int(k/3._wp)
1260 iz = int(l/3._wp)
1261
1262 x_id = j - int(3*ix) - 1
1263 y_id = k - int(3*iy) - 1
1264 z_id = l - int(3*iz) - 1
1265
1266 temp(1) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy, iz) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, iy, iz)
1267 temp(2) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy + y_id, iz) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, &
1268 & iy + y_id, iz)
1269 temp(3) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
1270
1271 temp(1) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy, iz + z_id) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, &
1272 & iy, iz + z_id)
1273 temp(2) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy + y_id, &
1274 & iz + z_id) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, iy + y_id, iz + z_id)
1275 temp(4) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
1276
1277 q_cons_vf(i)%sf(j, k, l) = (2._wp/3._wp)*temp(3) + (1._wp/3._wp)*temp(4)
1278 end do
1279 end do
1280 end do
1281 end do
1282
1283 end subroutine s_upsample_data
1284
1285 !> @brief True if `location` falls within this rank's own subdomain (a strict partition - each location is owned by exactly one
1286 !! rank, unlike an overlapping multi-rank ghost-stencil neighborhood). Used by both pre_process (to decide which
1287 !! generated/namelist IBs a rank writes to its own restart_data/ib_state_0.dat chunk) and simulation (to decide which IBs a rank
1288 !! owns for its own periodic IB-state writes) so the two stay consistent. glb_bounds_in is the global domain extent, used only
1289 !! to project a location that falls just outside the domain (floating-point edge case) onto the domain so some rank still claims
1290 !! it.
1291 function f_local_rank_owns_location(location, glb_bounds_in) result(owns_location)
1292
1293
1294# 712 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1295#if MFC_OpenACC
1296# 712 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1297!$acc routine seq
1298# 712 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1299#elif MFC_OpenMP
1300# 712 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1301
1302# 712 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1303
1304# 712 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1305!$omp declare target device_type(any)
1306# 712 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1307#endif
1308
1309 real(wp), dimension(3), intent(in) :: location
1310 type(bounds_info), dimension(3), intent(in) :: glb_bounds_in
1311 logical :: owns_location
1312 real(wp), dimension(3) :: projected_location
1313
1314 owns_location = .true.
1315
1316#ifdef MFC_MPI
1317 if (num_procs > 1) then
1318 projected_location(:) = location(:)
1319
1320 ! catch the edge case where the location lies just outside the computational domain
1321# 727 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1322 if (num_dims >= 1) then
1323 if (ib_bc_x%beg /= bc_periodic) then
1324 ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
1325 if (location(1) < glb_bounds_in(1)%beg) then
1326 projected_location(1) = glb_bounds_in(1)%beg
1327 else if (glb_bounds_in(1)%end < location(1)) then
1328 projected_location(1) = glb_bounds_in(1)%end - 1.0e-10_wp
1329 end if
1330 end if
1331 owns_location = owns_location .and. x_cb(-1) <= projected_location(1) &
1332 & .and. projected_location(1) < x_cb(m)
1333 end if
1334# 727 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1335 if (num_dims >= 2) then
1336 if (ib_bc_y%beg /= bc_periodic) then
1337 ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
1338 if (location(2) < glb_bounds_in(2)%beg) then
1339 projected_location(2) = glb_bounds_in(2)%beg
1340 else if (glb_bounds_in(2)%end < location(2)) then
1341 projected_location(2) = glb_bounds_in(2)%end - 1.0e-10_wp
1342 end if
1343 end if
1344 owns_location = owns_location .and. y_cb(-1) <= projected_location(2) &
1345 & .and. projected_location(2) < y_cb(n)
1346 end if
1347# 727 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1348 if (num_dims >= 3) then
1349 if (ib_bc_z%beg /= bc_periodic) then
1350 ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
1351 if (location(3) < glb_bounds_in(3)%beg) then
1352 projected_location(3) = glb_bounds_in(3)%beg
1353 else if (glb_bounds_in(3)%end < location(3)) then
1354 projected_location(3) = glb_bounds_in(3)%end - 1.0e-10_wp
1355 end if
1356 end if
1357 owns_location = owns_location .and. z_cb(-1) <= projected_location(3) &
1358 & .and. projected_location(3) < z_cb(p)
1359 end if
1360# 740 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
1361 end if
1362#endif
1363
1364 end function f_local_rank_owns_location
1365
1366end module m_helper
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
integer, intent(in) k
integer, intent(in) j
integer, intent(in) l
Compile-time constant parameters: default values, tolerances, and physical constants.
real(wp), parameter pi
Pi.
integer, parameter bc_periodic
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Defines global parameters for the computational domain, simulation algorithm, and initial conditions.
real(wp), dimension(:), allocatable weight
real(wp), dimension(:), allocatable r0
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
subroutine, public s_comp_n_from_prim(vftmp, rtmp, ntmp, weights)
Computes the bubble number density n from the primitive variables.
subroutine, public s_cross_product(a, b, c)
Compute the cross product c = a x b of two 3D vectors.
integer function, public modmul(a)
Compute a modular multiplication step for the linear congruential pseudo-random number generator.
impure subroutine, public s_initialize_nonpoly()
Initializes non-polydisperse bubble modeling.
recursive real(wp) function, public associated_legendre(x, l, m_order)
Associated Legendre polynomial P_l^m(x) (Ferrers function, Condon-Shortley phase)....
subroutine, public s_transform_triangle(triangle, matrix, matrix_n)
Transform a triangle by a matrix, one vertex at a time.
real(wp) function, dimension(3), public f_unit_vector(theta, eta)
Generate a unit vector uniformly distributed on the sphere from two random parameters.
logical function, public f_local_rank_owns_location(location, glb_bounds_in)
True if location falls within this rank's own subdomain (a strict partition - each location is owned ...
real(wp) function, public f_cut_on(x, eps)
Calculate a smooth cut-on function that is zero for x values smaller than zero and goes to one,...
subroutine, public s_transform_model(model, matrix, matrix_n)
Transform a model by a matrix, one triangle at a time.
type(t_bbox) function, public f_create_bbox(model)
Create a bounding box for a model.
real(wp) function f_gx(x)
Helper function for f_cut_on and f_cut_off.
subroutine, public s_simpson(local_weight, local_r0)
Compute the Simpson weights for quadrature.
elemental integer(kind=int64_kind) function, public double_factorial(n_in)
Calculate the double factorial of an integer.
impure subroutine s_initialize_bubble_vars()
Set bubble physical parameters and nondimensional numbers from the input configuration.
elemental integer(kind=int64_kind) function, public factorial(n_in)
Calculate the factorial of an integer.
subroutine, public s_upsample_data(q_cons_vf, q_cons_temp)
Upsample conservative variable fields from a coarsened grid back to the original resolution using int...
real(wp) function, dimension(1:4, 1:4), public f_create_transform_matrix(param, center)
Create a transformation matrix.
real(wp) function, public f_cut_off(x, eps)
Calculate a smooth cut-off function that is one for x values smaller than zero and goes to zero,...
subroutine, public s_transform_vec(vec, matrix)
Transform a vector by a matrix.
impure subroutine, public s_initialize_bubbles_model()
Initialize bubble model arrays for Euler or Lagrangian bubbles with polytropic or non-polytropic gas.
impure subroutine, public s_print_2d_array(a, div)
Print a 2D real array to standard output, optionally dividing each element by a given scalar.
elemental subroutine, public s_transcoeff(omega, peclet, re_trans, im_trans)
Computes the transfer coefficient for the non-polytropic bubble compression process.
subroutine, public s_prng(var, seed)
Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
pure real(wp) function, dimension(3), public f_cross(a, b)
Compute the cross product of two vectors.
subroutine, public s_downsample_data(q_cons_vf, q_cons_temp, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
Downsample conservative variable fields by a factor of 3 in each direction using volume averaging.
elemental logical function, public f_xor(lhs, rhs)
Perform XOR on lhs and rhs.
subroutine, public s_comp_n_from_cons(vftmp, nrtmp, ntmp, weights)
Compute the bubble number density from the conservative void fraction and weighted bubble radii.
real(wp) function, public real_ylm(theta, phi, l, m)
Real spherical harmonic Y_lm(theta, phi). theta = polar angle from +z (acos(z/r)),...
elemental integer function, public f_logical_to_int(predicate)
Convert a logical to 1 or 0.
elemental subroutine, public s_int_to_str(i, res)
Convert an integer to its trimmed string representation.
elemental subroutine, public s_swap(lhs, rhs)
Swap two real numbers.