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