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# 207 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 232 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 243 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 245 "/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# 283 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 293 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 303 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 312 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 329 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 346 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 358 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 364 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 370 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 376 "/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# 192 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 213 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 241 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 256 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 266 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 275 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 291 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 301 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 308 "/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# 21 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
234
235# 37 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
236
237# 50 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
238
239# 76 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
240
241# 91 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
242
243# 102 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
244
245# 115 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
246
247# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
248
249# 154 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
250
251# 165 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
252
253# 176 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
254
255# 187 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257# 198 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259# 208 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
260
261# 214 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263# 220 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264
265# 226 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
266
267# 232 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
268
269# 234 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270# 235 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
271! New line at end of file is required for FYPP
272# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
273
274# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
275
276! Caution:
277! This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI rank.
278! That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0.
279! For an example see misc/nvidia_uvm/bind.sh.
280# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
281
282# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
283
284# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
285
286# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
287
288# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
289
290# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
291
292# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
293
294# 167 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
295! New line at end of file is required for FYPP
296# 3 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp" 2
297
298!>
299!! @file
300!! @brief Contains module m_helper
301
302!> @brief Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions
303module m_helper
304
305 use m_derived_types !< definitions of the derived types
306
307 use m_global_parameters !< definitions of the global parameters
308
309 use ieee_arithmetic !< for checking nan
310
311 implicit none
312
313 private;
314 public :: s_comp_n_from_prim, &
315 s_comp_n_from_cons, &
316 s_initialize_bubbles_model, &
317 s_initialize_nonpoly, &
318 s_simpson, &
319 s_transcoeff, &
320 s_int_to_str, &
321 s_transform_vec, &
322 s_transform_triangle, &
323 s_transform_model, &
324 s_swap, &
325 f_cross, &
326 f_create_transform_matrix, &
327 f_create_bbox, &
328 s_print_2d_array, &
329 f_xor, &
330 f_logical_to_int, &
331 unassociated_legendre, &
332 associated_legendre, &
333 spherical_harmonic_func, &
334 double_factorial, &
335 factorial, &
336 f_cut_on, &
337 f_cut_off, &
338 s_downsample_data, &
339 s_upsample_data
340
341contains
342
343 !> Computes the bubble number density n from the primitive variables
344 !! @param vftmp is the void fraction
345 !! @param Rtmp is the bubble radii
346 !! @param ntmp is the output number bubble density
347 !! @param weights is the quadrature weights
348 subroutine s_comp_n_from_prim(vftmp, Rtmp, ntmp, weights)
349
350# 55 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
351#if MFC_OpenACC
352# 55 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
353!$acc routine seq
354# 55 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
355#elif MFC_OpenMP
356# 55 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
357
358# 55 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
359
360# 55 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
361!$omp declare target device_type(any)
362# 55 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
363#endif
364 real(wp), intent(in) :: vftmp
365 real(wp), dimension(nb), intent(in) :: rtmp
366 real(wp), intent(out) :: ntmp
367 real(wp), dimension(nb), intent(in) :: weights
368
369 real(wp) :: r3
370
371 r3 = dot_product(weights, rtmp**3._wp)
372 ntmp = (3._wp/(4._wp*pi))*vftmp/r3
373
374 end subroutine s_comp_n_from_prim
375
376 !> @brief Computes the bubble number density from the conservative void fraction and weighted bubble radii.
377 subroutine s_comp_n_from_cons(vftmp, nRtmp, ntmp, weights)
378
379# 70 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
380#if MFC_OpenACC
381# 70 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
382!$acc routine seq
383# 70 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
384#elif MFC_OpenMP
385# 70 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
386
387# 70 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
388
389# 70 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
390!$omp declare target device_type(any)
391# 70 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
392#endif
393 real(wp), intent(in) :: vftmp
394 real(wp), dimension(nb), intent(in) :: nrtmp
395 real(wp), intent(out) :: ntmp
396 real(wp), dimension(nb), intent(in) :: weights
397
398 real(wp) :: nr3
399
400 nr3 = dot_product(weights, nrtmp**3._wp)
401 ntmp = sqrt((4._wp*pi/3._wp)*nr3/vftmp)
402
403 end subroutine s_comp_n_from_cons
404
405 !> @brief Prints a 2D real array to standard output, optionally dividing each element by a given scalar.
406 impure subroutine s_print_2d_array(A, div)
407
408 real(wp), dimension(:, :), intent(in) :: a
409 real(wp), optional, intent(in) :: div
410
411 integer :: i, j
412 integer :: local_m, local_n
413 real(wp) :: c
414
415 local_m = size(a, 1)
416 local_n = size(a, 2)
417
418 if (present(div)) then
419 c = div
420 else
421 c = 1._wp
422 end if
423
424 print *, local_m, local_n
425
426 do i = 1, local_m
427 do j = 1, local_n
428 write (*, fmt="(F12.4)", advance="no") a(i, j)/c
429 end do
430 write (*, fmt="(A1)") " "
431 end do
432 write (*, fmt="(A1)") " "
433
434 end subroutine s_print_2d_array
435
436 !>
437 !! bubbles_euler + polytropic
438 !! bubbles_euler + non-polytropic
439 !! bubbles_lagrange + non-polytropic
440 impure subroutine s_initialize_bubbles_model()
441
442 ! Allocate memory
443 if (bubbles_euler) then
444#ifdef MFC_DEBUG
445# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
446 block
447# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
448 use iso_fortran_env, only: output_unit
449# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
450
451# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
452 print *, 'm_helper.fpp:122: ', '@:ALLOCATE(weight(nb), R0(nb))'
453# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
454
455# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
456 call flush (output_unit)
457# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
458 end block
459# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
460#endif
461# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
462 allocate (weight(nb), r0(nb))
463# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
464
465# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
466
467# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
468
469# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
470#if defined(MFC_OpenACC)
471# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
472!$acc enter data create(weight, R0)
473# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
474#elif defined(MFC_OpenMP)
475# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
476!$omp target enter data map(always,alloc:weight, R0)
477# 122 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
478#endif
479 if (.not. polytropic) then
480#ifdef MFC_DEBUG
481# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
482 block
483# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
484 use iso_fortran_env, only: output_unit
485# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
486
487# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
488 print *, 'm_helper.fpp:124: ', '@:ALLOCATE(pb0(nb), Pe_T(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))'
489# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
490
491# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
492 call flush (output_unit)
493# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
494 end block
495# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
496#endif
497# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
498 allocate (pb0(nb), pe_t(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))
499# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
500
501# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
502
503# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
504
505# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
506
507# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
508
509# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
510
511# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
512
513# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
514#if defined(MFC_OpenACC)
515# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
516!$acc enter data create(pb0, Pe_T, k_g, k_v, mass_g0, mass_v0)
517# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
518#elif defined(MFC_OpenMP)
519# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
520!$omp target enter data map(always,alloc:pb0, Pe_T, k_g, k_v, mass_g0, mass_v0)
521# 124 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
522#endif
523#ifdef MFC_DEBUG
524# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
525 block
526# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
527 use iso_fortran_env, only: output_unit
528# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
529
530# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
531 print *, 'm_helper.fpp:125: ', '@:ALLOCATE(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb))'
532# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
533
534# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
535 call flush (output_unit)
536# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
537 end block
538# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
539#endif
540# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
541 allocate (re_trans_t(nb), re_trans_c(nb), im_trans_t(nb), im_trans_c(nb))
542# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
543
544# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
545
546# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
547
548# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
549
550# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
551
552# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
553#if defined(MFC_OpenACC)
554# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
555!$acc enter data create(Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c)
556# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
557#elif defined(MFC_OpenMP)
558# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
559!$omp target enter data map(always,alloc:Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c)
560# 125 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
561#endif
562 else if (qbmm) then
563#ifdef MFC_DEBUG
564# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
565 block
566# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
567 use iso_fortran_env, only: output_unit
568# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
569
570# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
571 print *, 'm_helper.fpp:127: ', '@:ALLOCATE(pb0(nb))'
572# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
573
574# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
575 call flush (output_unit)
576# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
577 end block
578# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
579#endif
580# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
581 allocate (pb0(nb))
582# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
583
584# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
585
586# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
587#if defined(MFC_OpenACC)
588# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
589!$acc enter data create(pb0)
590# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
591#elif defined(MFC_OpenMP)
592# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
593!$omp target enter data map(always,alloc:pb0)
594# 127 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
595#endif
596 end if
597
598 ! Compute quadrature weights and nodes for polydisperse simulations
599 if (nb > 1) then
600 call s_simpson(weight, r0)
601 else if (nb == 1) then
602 r0 = 1._wp
603 weight = 1._wp
604 else
605 stop 'Invalid value of nb'
606 end if
607 r0 = r0*bub_pp%R0ref
608 end if
609
610 ! Initialize bubble variables
611 call s_initialize_bubble_vars()
612
613 end subroutine s_initialize_bubbles_model
614
615 !>
616 impure subroutine s_initialize_bubble_vars()
617
618 r0ref = bub_pp%R0ref; p0ref = bub_pp%p0ref
619 rho0ref = bub_pp%rho0ref;
620 ss = bub_pp%ss; pv = bub_pp%pv; vd = bub_pp%vd
621 mu_l = bub_pp%mu_l; mu_v = bub_pp%mu_v; mu_g = bub_pp%mu_g
622 gam_v = bub_pp%gam_v; gam_g = bub_pp%gam_g
623 if (.not. polytropic) then
624 if (bubbles_euler) then
625 m_v = bub_pp%M_v; m_g = bub_pp%M_g
626 k_v = bub_pp%k_v; k_g = bub_pp%k_g
627 end if
628 r_v = bub_pp%R_v; r_g = bub_pp%R_g
629 tw = bub_pp%T0ref
630 end if
631 if (bubbles_lagrange) then
632 cp_v = bub_pp%cp_v; cp_g = bub_pp%cp_g
633 k_vl = bub_pp%k_v; k_gl = bub_pp%k_g
634 end if
635
636 ! Input quantities
637 if (bubbles_euler .and. (.not. polytropic)) then
638 if (thermal == 2) then
639 gam_m = 1._wp
640 else
641 gam_m = gam_g
642 end if
643 end if
644
645 ! Nondimensional numbers
646 eu = p0ref
647 ca = eu - pv
648 if (.not. f_is_default(bub_pp%ss)) web = 1._wp/ss
649 if (.not. f_is_default(bub_pp%mu_l)) re_inv = mu_l
650 if (.not. polytropic) pe_c = 1._wp/vd
651
652 if (bubbles_euler) then
653 ! Initialize variables for non-polytropic (Preston) model
654 if (.not. polytropic) then
655 call s_initialize_nonpoly()
656 end if
657 ! Initialize pb based on surface tension for qbmm (polytropic)
658 if (qbmm .and. polytropic) then
659 pb0 = eu
660 if (.not. f_is_default(web)) then
661 pb0 = pb0 + 2._wp/web/r0
662 end if
663 end if
664 end if
665
666 end subroutine s_initialize_bubble_vars
667
668 !> Initializes non-polydisperse bubble modeling
669 impure subroutine s_initialize_nonpoly()
670 integer :: ir
671 real(wp), dimension(nb) :: chi_vw0, cp_m0, k_m0, rho_m0, x_vw, omegan, rhol0
672
673 real(wp), parameter :: k_poly = 1._wp !<
674 !! polytropic index used to compute isothermal natural frequency
675
676 ! phi_vg & phi_gv (phi_gg = phi_vv = 1) (Eq. 2.22 in Ando 2010)
677 phi_vg = (1._wp + sqrt(mu_v/mu_g)*(m_g/m_v)**(0.25_wp))**2 &
678 /(sqrt(8._wp)*sqrt(1._wp + m_v/m_g))
679 phi_gv = (1._wp + sqrt(mu_g/mu_v)*(m_v/m_g)**(0.25_wp))**2 &
680 /(sqrt(8._wp)*sqrt(1._wp + m_g/m_v))
681
682 ! internal bubble pressure
683 pb0 = eu + 2._wp/web/r0
684
685 ! mass fraction of vapor (Eq. 2.19 in Ando 2010)
686 chi_vw0 = 1._wp/(1._wp + r_v/r_g*(pb0/pv - 1._wp))
687
688 ! specific heat for gas/vapor mixture
689 cp_m0 = chi_vw0*r_v*gam_v/(gam_v - 1._wp) &
690 + (1._wp - chi_vw0)*r_g*gam_g/(gam_g - 1._wp)
691
692 ! mole fraction of vapor (Eq. 2.23 in Ando 2010)
693 x_vw = m_g*chi_vw0/(m_v + (m_g - m_v)*chi_vw0)
694
695 ! thermal conductivity for gas/vapor mixture (Eq. 2.21 in Ando 2010)
696 k_m0 = x_vw*k_v/(x_vw + (1._wp - x_vw)*phi_vg) &
697 + (1._wp - x_vw)*k_g/(x_vw*phi_gv + 1._wp - x_vw)
698 k_g(:) = k_g(:)/k_m0(:)
699 k_v(:) = k_v(:)/k_m0(:)
700
701 ! mixture density (Eq. 2.20 in Ando 2010)
702 rho_m0 = pv/(chi_vw0*r_v*tw)
703
704 ! mass of gas/vapor
705 mass_g0(:) = (4._wp*pi/3._wp)*(pb0(:) - pv)/(r_g*tw)*r0(:)**3
706 mass_v0(:) = (4._wp*pi/3._wp)*pv/(r_v*tw)*r0(:)**3
707
708 ! Peclet numbers
709 pe_t(:) = rho_m0*cp_m0(:)/k_m0(:)
710
711 ! natural frequencies (Eq. B.1)
712 omegan(:) = sqrt(3._wp*k_poly*ca + 2._wp*(3._wp*k_poly - 1._wp)/(web*r0))/r0/sqrt(rho0ref)
713 do ir = 1, nb
714 call s_transcoeff(omegan(ir)*r0(ir), pe_t(ir)*r0(ir), &
715 re_trans_t(ir), im_trans_t(ir))
716 call s_transcoeff(omegan(ir)*r0(ir), pe_c*r0(ir), &
717 re_trans_c(ir), im_trans_c(ir))
718 end do
719 im_trans_t = 0._wp
720
721 end subroutine s_initialize_nonpoly
722
723 !> Computes the transfer coefficient for the non-polytropic bubble compression process
724 !! @param omega natural frequencies
725 !! @param peclet Peclet number
726 !! @param Re_trans Real part of the transport coefficients
727 !! @param Im_trans Imaginary part of the transport coefficients
728 elemental subroutine s_transcoeff(omega, peclet, Re_trans, Im_trans)
729
730 real(wp), intent(in) :: omega, peclet
731 real(wp), intent(out) :: re_trans, im_trans
732
733 complex(wp) :: imag, trans, c1, c2, c3
734
735 imag = (0._wp, 1._wp)
736
737 c1 = imag*omega*peclet
738 c2 = sqrt(c1)
739 c3 = (exp(c2) - exp(-c2))/(exp(c2) + exp(-c2)) ! TANH(c2)
740 trans = ((c2/c3 - 1._wp)**(-1) - 3._wp/c1)**(-1) ! transfer function
741
742 re_trans = trans
743 im_trans = aimag(trans)
744
745 end subroutine s_transcoeff
746
747 !> @brief Converts an integer to its trimmed string representation.
748 elemental subroutine s_int_to_str(i, res)
749
750 integer, intent(in) :: i
751 character(len=*), intent(inout) :: res
752
753 write (res, '(I0)') i
754 res = trim(res)
755 end subroutine s_int_to_str
756
757 !> Computes the Simpson weights for quadrature
758 subroutine s_simpson(local_weight, local_R0)
759
760 real(wp), dimension(:), intent(inout) :: local_weight
761 real(wp), dimension(:), intent(inout) :: local_r0
762 integer :: ir
763 real(wp) :: r0mn, r0mx, dphi, tmp, sd
764 real(wp), dimension(nb) :: phi
765
766 sd = poly_sigma
767 r0mn = 0.8_wp*exp(-2.8_wp*sd)
768 r0mx = 0.2_wp*exp(9.5_wp*sd) + 1._wp
769
770 ! phi = ln( R0 ) & return R0
771 do ir = 1, nb
772 phi(ir) = log(r0mn) &
773 + (ir - 1._wp)*log(r0mx/r0mn)/(nb - 1._wp)
774 local_r0(ir) = exp(phi(ir))
775 end do
776
777# 310 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
778 dphi = phi(2) - phi(1)
779# 312 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
780
781 ! weights for quadrature using Simpson's rule
782 do ir = 2, nb - 1
783 ! Gaussian
784 tmp = exp(-0.5_wp*(phi(ir)/sd)**2)/sqrt(2._wp*pi)/sd
785 if (mod(ir, 2) == 0) then
786 local_weight(ir) = tmp*4._wp*dphi/3._wp
787 else
788 local_weight(ir) = tmp*2._wp*dphi/3._wp
789 end if
790 end do
791 tmp = exp(-0.5_wp*(phi(1)/sd)**2)/sqrt(2._wp*pi)/sd
792 local_weight(1) = tmp*dphi/3._wp
793 tmp = exp(-0.5_wp*(phi(nb)/sd)**2)/sqrt(2._wp*pi)/sd
794 local_weight(nb) = tmp*dphi/3._wp
795
796 end subroutine s_simpson
797
798 !> This procedure computes the cross product of two vectors.
799 !! @param a First vector.
800 !! @param b Second vector.
801 !! @return The cross product of the two vectors.
802 pure function f_cross(a, b) result(c)
803
804 real(wp), dimension(3), intent(in) :: a, b
805 real(wp), dimension(3) :: c
806
807 c(1) = a(2)*b(3) - a(3)*b(2)
808 c(2) = a(3)*b(1) - a(1)*b(3)
809 c(3) = a(1)*b(2) - a(2)*b(1)
810 end function f_cross
811
812 !> This procedure swaps two real numbers.
813 !! @param lhs Left-hand side.
814 !! @param rhs Right-hand side.
815 elemental subroutine s_swap(lhs, rhs)
816
817 real(wp), intent(inout) :: lhs, rhs
818 real(wp) :: ltemp
819
820 ltemp = lhs
821 lhs = rhs
822 rhs = ltemp
823 end subroutine s_swap
824
825 !> This procedure creates a transformation matrix.
826 !! @param param Parameters for the transformation.
827 !! @param center Optional center point for the transformation.
828 !! @return Transformation matrix.
829 function f_create_transform_matrix(param, center) result(out_matrix)
830
831 type(ic_model_parameters), intent(in) :: param
832 real(wp), dimension(1:3), optional, intent(in) :: center
833 real(wp), dimension(1:4, 1:4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
834
835 sc = transpose(reshape([ &
836 param%scale(1), 0._wp, 0._wp, 0._wp, &
837 0._wp, param%scale(2), 0._wp, 0._wp, &
838 0._wp, 0._wp, param%scale(3), 0._wp, &
839 0._wp, 0._wp, 0._wp, 1._wp], shape(sc)))
840
841 rz = transpose(reshape([ &
842 cos(param%rotate(3)), -sin(param%rotate(3)), 0._wp, 0._wp, &
843 sin(param%rotate(3)), cos(param%rotate(3)), 0._wp, 0._wp, &
844 0._wp, 0._wp, 1._wp, 0._wp, &
845 0._wp, 0._wp, 0._wp, 1._wp], shape(rz)))
846
847 rx = transpose(reshape([ &
848 1._wp, 0._wp, 0._wp, 0._wp, &
849 0._wp, cos(param%rotate(1)), -sin(param%rotate(1)), 0._wp, &
850 0._wp, sin(param%rotate(1)), cos(param%rotate(1)), 0._wp, &
851 0._wp, 0._wp, 0._wp, 1._wp], shape(rx)))
852
853 ry = transpose(reshape([ &
854 cos(param%rotate(2)), 0._wp, sin(param%rotate(2)), 0._wp, &
855 0._wp, 1._wp, 0._wp, 0._wp, &
856 -sin(param%rotate(2)), 0._wp, cos(param%rotate(2)), 0._wp, &
857 0._wp, 0._wp, 0._wp, 1._wp], shape(ry)))
858
859 tr = transpose(reshape([ &
860 1._wp, 0._wp, 0._wp, param%translate(1), &
861 0._wp, 1._wp, 0._wp, param%translate(2), &
862 0._wp, 0._wp, 1._wp, param%translate(3), &
863 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
864
865 if (present(center)) then
866 ! Translation matrix to move center to the origin
867 t_to_origin = transpose(reshape([ &
868 1._wp, 0._wp, 0._wp, -center(1), &
869 0._wp, 1._wp, 0._wp, -center(2), &
870 0._wp, 0._wp, 1._wp, -center(3), &
871 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
872
873 ! Translation matrix to move center back to original position
874 t_back = transpose(reshape([ &
875 1._wp, 0._wp, 0._wp, center(1), &
876 0._wp, 1._wp, 0._wp, center(2), &
877 0._wp, 0._wp, 1._wp, center(3), &
878 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
879
880 out_matrix = matmul(tr, matmul(t_back, matmul(ry, matmul(rx, matmul(rz, matmul(sc, t_to_origin))))))
881 else
882 out_matrix = matmul(ry, matmul(rx, rz))
883 end if
884
885 end function f_create_transform_matrix
886
887 !> This procedure transforms a vector by a matrix.
888 !! @param vec Vector to transform.
889 !! @param matrix Transformation matrix.
890 subroutine s_transform_vec(vec, matrix)
891
892 real(wp), dimension(1:3), intent(inout) :: vec
893 real(wp), dimension(1:4, 1:4), intent(in) :: matrix
894
895 real(wp), dimension(1:4) :: tmp
896
897 tmp = matmul(matrix, [vec(1), vec(2), vec(3), 1._wp])
898 vec = tmp(1:3)
899
900 end subroutine s_transform_vec
901
902 !> This procedure transforms a triangle by a matrix, one vertex at a time.
903 !! @param triangle Triangle to transform.
904 !! @param matrix Transformation matrix.
905 !! @param matrix_n Normal transformation matrix.
906 subroutine s_transform_triangle(triangle, matrix, matrix_n)
907
908 type(t_triangle), intent(inout) :: triangle
909 real(wp), dimension(1:4, 1:4), intent(in) :: matrix, matrix_n
910
911 integer :: i
912
913 do i = 1, 3
914 call s_transform_vec(triangle%v(i, :), matrix)
915 end do
916
917 call s_transform_vec(triangle%n(1:3), matrix_n)
918
919 end subroutine s_transform_triangle
920
921 !> This procedure transforms a model by a matrix, one triangle at a time.
922 !! @param model Model to transform.
923 !! @param matrix Transformation matrix.
924 !! @param matrix_n Normal transformation matrix.
925 subroutine s_transform_model(model, matrix, matrix_n)
926
927 type(t_model), intent(inout) :: model
928 real(wp), dimension(1:4, 1:4), intent(in) :: matrix, matrix_n
929
930 integer :: i
931
932 do i = 1, size(model%trs)
933 call s_transform_triangle(model%trs(i), matrix, matrix_n)
934 end do
935
936 end subroutine s_transform_model
937
938 !> This procedure creates a bounding box for a model.
939 !! @param model Model to create bounding box for.
940 !! @return Bounding box.
941 function f_create_bbox(model) result(bbox)
942
943 type(t_model), intent(in) :: model
944 type(t_bbox) :: bbox
945
946 integer :: i, j
947
948 if (size(model%trs) == 0) then
949 bbox%min = 0._wp
950 bbox%max = 0._wp
951 return
952 end if
953
954 bbox%min = model%trs(1)%v(1, :)
955 bbox%max = model%trs(1)%v(1, :)
956
957 do i = 1, size(model%trs)
958 do j = 1, 3
959 bbox%min = min(bbox%min, model%trs(i)%v(j, :))
960 bbox%max = max(bbox%max, model%trs(i)%v(j, :))
961 end do
962 end do
963
964 end function f_create_bbox
965
966 !> This procedure performs xor on lhs and rhs.
967 !! @param lhs logical input.
968 !! @param rhs other logical input.
969 !! @return xored result.
970 elemental function f_xor(lhs, rhs) result(res)
971
972 logical, intent(in) :: lhs, rhs
973 logical :: res
974
975 res = (lhs .and. .not. rhs) .or. (.not. lhs .and. rhs)
976 end function f_xor
977
978 !> This procedure converts logical to 1 or 0.
979 !! @param predicate A Logical argument.
980 !! @return 1 if .true., 0 if .false..
981 elemental function f_logical_to_int(predicate) result(int)
982
983 logical, intent(in) :: predicate
984 integer :: int
985
986 if (predicate) then
987 int = 1
988 else
989 int = 0
990 end if
991 end function f_logical_to_int
992
993 !> This function generates the unassociated legendre poynomials
994 !! @param x is the input value
995 !! @param l is the degree
996 !! @return P is the unassociated legendre polynomial evaluated at x
997 recursive function unassociated_legendre(x, l) result(result_P)
998
999 integer, intent(in) :: l
1000 real(wp), intent(in) :: x
1001 real(wp) :: result_p
1002
1003 if (l == 0) then
1004 result_p = 1._wp
1005 else if (l == 1) then
1006 result_p = x
1007 else
1008 result_p = ((2*l - 1)*x*unassociated_legendre(x, l - 1) - (l - 1)*unassociated_legendre(x, l - 2))/l
1009 end if
1010
1011 end function unassociated_legendre
1012
1013 !> This function calculates the spherical harmonic function evaluated at x and phi
1014 !! @param x is the x coordinate
1015 !! @param phi is the phi coordinate
1016 !! @param l is the degree
1017 !! @param m_order is the order
1018 !! @return Y is the spherical harmonic function evaluated at x and phi
1019 recursive function spherical_harmonic_func(x, phi, l, m_order) result(Y)
1020
1021 integer, intent(in) :: l, m_order
1022 real(wp), intent(in) :: x, phi
1023 real(wp) :: y, prefactor, local_pi
1024
1025 local_pi = acos(-1._wp)
1026 prefactor = sqrt((2*l + 1)/(4*local_pi)*factorial(l - m_order)/factorial(l + m_order));
1027 if (m_order == 0) then
1028 y = prefactor*associated_legendre(x, l, m_order);
1029 elseif (m_order > 0) then
1030 y = (-1._wp)**m_order*sqrt(2._wp)*prefactor*associated_legendre(x, l, m_order)*cos(m_order*phi);
1031 end if
1032
1033 end function spherical_harmonic_func
1034
1035 !> This function generates the associated legendre polynomials evaluated
1036 !! at x with inputs l and m
1037 !! @param x is the input value
1038 !! @param l is the degree
1039 !! @param m_order is the order
1040 !! @return P is the associated legendre polynomial evaluated at x
1041 recursive function associated_legendre(x, l, m_order) result(result_P)
1042
1043 integer, intent(in) :: l, m_order
1044 real(wp), intent(in) :: x
1045 real(wp) :: result_p
1046
1047 if (m_order <= 0 .and. l <= 0) then
1048 result_p = 1;
1049 elseif (l == 1 .and. m_order <= 0) then
1050 result_p = x;
1051 elseif (l == 1 .and. m_order == 1) then
1052 result_p = -(1 - x**2)**(1._wp/2._wp);
1053 elseif (m_order == l) then
1054 result_p = (-1)**l*double_factorial(2*l - 1)*(1 - x**2)**(l/2);
1055 elseif (m_order == l - 1) then
1056 result_p = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1);
1057 else
1058 result_p = ((2*l - 1)*x*associated_legendre(x, l - 1, m_order) - (l + m_order - 1)*associated_legendre(x, l - 2, m_order))/(l - m_order);
1059 end if
1060
1061 end function associated_legendre
1062
1063 !> This function calculates the double factorial value of an integer
1064 !! @param n_in is the input integer
1065 !! @return R is the double factorial value of n
1066 elemental function double_factorial(n_in) result(R_result)
1067
1068 integer, intent(in) :: n_in
1069 integer, parameter :: int64_kind = selected_int_kind(18) ! 18 bytes for 64-bit integer
1070 integer(kind=int64_kind) :: r_result
1071 integer :: i
1072
1073 r_result = product((/(i, i=n_in, 1, -2)/))
1074
1075 end function double_factorial
1076
1077 !> The following function calculates the factorial value of an integer
1078 !! @param n_in is the input integer
1079 !! @return R is the factorial value of n
1080 elemental function factorial(n_in) result(R_result)
1081
1082 integer, intent(in) :: n_in
1083 integer, parameter :: int64_kind = selected_int_kind(18) ! 18 bytes for 64-bit integer
1084 integer(kind=int64_kind) :: r_result
1085
1086 integer :: i
1087
1088 r_result = product((/(i, i=n_in, 1, -1)/))
1089
1090 end function factorial
1091
1092 !> This function calculates a smooth cut-on function that is zero for x values
1093 !! smaller than zero and goes to one. It can be used for generating smooth
1094 !! initial conditions
1095 !! @param x is the input value
1096 !! @param eps is the smoothing parameter
1097 !! @return fx is the cut-on function evaluated at x
1098 function f_cut_on(x, eps) result(fx)
1099
1100 real(wp), intent(in) :: x, eps
1101 real(wp) :: fx
1102
1103 fx = 1 - f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
1104
1105 end function f_cut_on
1106
1107 !> This function calculates a smooth cut-off function that is one for x values
1108 !! smaller than zero and goes to zero. It can be used for generating smooth
1109 !! initial conditions
1110 !! @param x is the input value
1111 !! @param eps is the smoothing parameter
1112 !! @return fx is the cut-ff function evaluated at x
1113 function f_cut_off(x, eps) result(fx)
1114
1115 real(wp), intent(in) :: x, eps
1116 real(wp) :: fx
1117
1118 fx = f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
1119
1120 end function f_cut_off
1121
1122 !> This function is a helper function for the functions f_cut_on and f_cut_off
1123 !! @param x is the input value
1124 !! @return gx is the result
1125 function f_gx(x) result(gx)
1126
1127 real(wp), intent(in) :: x
1128 real(wp) :: gx
1129
1130 if (x > 0) then
1131 gx = exp(-1._wp/x)
1132 else
1133 gx = 0._wp
1134 end if
1135
1136 end function f_gx
1137
1138 !> @brief Downsamples conservative variable fields by a factor of 3 in each direction using volume averaging.
1139 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)
1140
1141 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf, q_cons_temp
1142
1143 ! Down sampling variables
1144 integer :: i, j, k, l
1145 integer :: ix, iy, iz, x_id, y_id, z_id
1146 integer, intent(inout) :: m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds
1147
1148 m_ds = int((m + 1)/3) - 1
1149 n_ds = int((n + 1)/3) - 1
1150 p_ds = int((p + 1)/3) - 1
1151
1152 m_glb_ds = int((m_glb + 1)/3) - 1
1153 n_glb_ds = int((n_glb + 1)/3) - 1
1154 p_glb_ds = int((p_glb + 1)/3) - 1
1155
1156 do l = -1, p_ds + 1
1157 do k = -1, n_ds + 1
1158 do j = -1, m_ds + 1
1159 x_id = 3*j + 1
1160 y_id = 3*k + 1
1161 z_id = 3*l + 1
1162 do i = 1, sys_size
1163 q_cons_temp(i)%sf(j, k, l) = 0
1164
1165 do iz = -1, 1
1166 do iy = -1, 1
1167 do ix = -1, 1
1168 q_cons_temp(i)%sf(j, k, l) = q_cons_temp(i)%sf(j, k, l) &
1169 + (1._wp/27._wp)*q_cons_vf(i)%sf(x_id + ix, y_id + iy, z_id + iz)
1170 end do
1171 end do
1172 end do
1173 end do
1174 end do
1175 end do
1176 end do
1177
1178 end subroutine s_downsample_data
1179
1180 !> @brief Upsamples conservative variable fields from a coarsened grid back to the original resolution using interpolation.
1181 subroutine s_upsample_data(q_cons_vf, q_cons_temp)
1182
1183 type(scalar_field), intent(inout), dimension(sys_size) :: q_cons_vf, q_cons_temp
1184 integer :: i, j, k, l
1185 integer :: ix, iy, iz
1186 integer :: x_id, y_id, z_id
1187 real(wp), dimension(4) :: temp
1188
1189 do l = 0, p
1190 do k = 0, n
1191 do j = 0, m
1192 do i = 1, sys_size
1193
1194 ix = int(j/3._wp)
1195 iy = int(k/3._wp)
1196 iz = int(l/3._wp)
1197
1198 x_id = j - int(3*ix) - 1
1199 y_id = k - int(3*iy) - 1
1200 z_id = l - int(3*iz) - 1
1201
1202 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)
1203 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, iy + y_id, iz)
1204 temp(3) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
1205
1206 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, iy, iz + z_id)
1207 temp(2) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy + y_id, iz + z_id) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, iy + y_id, iz + z_id)
1208 temp(4) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
1209
1210 q_cons_vf(i)%sf(j, k, l) = (2._wp/3._wp)*temp(3) + (1._wp/3._wp)*temp(4)
1211
1212 end do
1213 end do
1214 end do
1215 end do
1216
1217 end subroutine s_upsample_data
1218
1219end module m_helper
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:1119
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1121
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1121
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1121
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:292
m_global_parameters
Global parameters for the post-process: domain geometry, equation of state, and output database setti...
Definition m_global_parameters.fpp.f90:19
m_global_parameters::bubbles_euler
logical bubbles_euler
Definition m_global_parameters.fpp.f90:324
m_global_parameters::weight
real(wp), dimension(:), allocatable weight
Definition m_global_parameters.fpp.f90:323
m_global_parameters::r0
real(wp), dimension(:), allocatable r0
Definition m_global_parameters.fpp.f90:323
m_global_parameters::polytropic
logical polytropic
Definition m_global_parameters.fpp.f90:326
m_global_parameters::nb
integer nb
Definition m_global_parameters.fpp.f90:321
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:303
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:349
m_helper::s_initialize_nonpoly
impure subroutine, public s_initialize_nonpoly()
Initializes non-polydisperse bubble modeling.
Definition m_helper.fpp.f90:670
m_helper::associated_legendre
recursive real(wp) function, public associated_legendre(x, l, m_order)
This function generates the associated legendre polynomials evaluated at x with inputs l and m.
Definition m_helper.fpp.f90:1042
m_helper::spherical_harmonic_func
recursive real(wp) function, public spherical_harmonic_func(x, phi, l, m_order)
This function calculates the spherical harmonic function evaluated at x and phi.
Definition m_helper.fpp.f90:1020
m_helper::s_transform_triangle
subroutine, public s_transform_triangle(triangle, matrix, matrix_n)
This procedure transforms a triangle by a matrix, one vertex at a time.
Definition m_helper.fpp.f90:907
m_helper::f_cut_on
real(wp) function, public f_cut_on(x, eps)
This function calculates a smooth cut-on function that is zero for x values smaller than zero and goe...
Definition m_helper.fpp.f90:1099
m_helper::s_transform_model
subroutine, public s_transform_model(model, matrix, matrix_n)
This procedure transforms a model by a matrix, one triangle at a time.
Definition m_helper.fpp.f90:926
m_helper::f_create_bbox
type(t_bbox) function, public f_create_bbox(model)
This procedure creates a bounding box for a model.
Definition m_helper.fpp.f90:942
m_helper::f_gx
real(wp) function f_gx(x)
This function is a helper function for the functions f_cut_on and f_cut_off.
Definition m_helper.fpp.f90:1126
m_helper::s_simpson
subroutine, public s_simpson(local_weight, local_r0)
Computes the Simpson weights for quadrature.
Definition m_helper.fpp.f90:759
m_helper::double_factorial
elemental integer(kind=int64_kind) function, public double_factorial(n_in)
This function calculates the double factorial value of an integer.
Definition m_helper.fpp.f90:1067
m_helper::s_initialize_bubble_vars
impure subroutine s_initialize_bubble_vars()
Definition m_helper.fpp.f90:617
m_helper::factorial
elemental integer(kind=int64_kind) function, public factorial(n_in)
The following function calculates the factorial value of an integer.
Definition m_helper.fpp.f90:1081
m_helper::s_upsample_data
subroutine, public s_upsample_data(q_cons_vf, q_cons_temp)
Upsamples conservative variable fields from a coarsened grid back to the original resolution using in...
Definition m_helper.fpp.f90:1182
m_helper::unassociated_legendre
recursive real(wp) function, public unassociated_legendre(x, l)
This function generates the unassociated legendre poynomials.
Definition m_helper.fpp.f90:998
m_helper::f_cut_off
real(wp) function, public f_cut_off(x, eps)
This function calculates a smooth cut-off function that is one for x values smaller than zero and goe...
Definition m_helper.fpp.f90:1114
m_helper::s_transform_vec
subroutine, public s_transform_vec(vec, matrix)
This procedure transforms a vector by a matrix.
Definition m_helper.fpp.f90:891
m_helper::s_initialize_bubbles_model
impure subroutine, public s_initialize_bubbles_model()
bubbles_euler + polytropic bubbles_euler + non-polytropic bubbles_lagrange + non-polytropic
Definition m_helper.fpp.f90:441
m_helper::s_print_2d_array
impure subroutine, public s_print_2d_array(a, div)
Prints a 2D real array to standard output, optionally dividing each element by a given scalar.
Definition m_helper.fpp.f90:407
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:729
m_helper::f_cross
pure real(wp) function, dimension(3), public f_cross(a, b)
This procedure computes the cross product of two vectors.
Definition m_helper.fpp.f90:803
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)
Downsamples conservative variable fields by a factor of 3 in each direction using volume averaging.
Definition m_helper.fpp.f90:1140
m_helper::f_xor
elemental logical function, public f_xor(lhs, rhs)
This procedure performs xor on lhs and rhs.
Definition m_helper.fpp.f90:971
m_helper::s_comp_n_from_cons
subroutine, public s_comp_n_from_cons(vftmp, nrtmp, ntmp, weights)
Computes the bubble number density from the conservative void fraction and weighted bubble radii.
Definition m_helper.fpp.f90:378
m_helper::f_create_transform_matrix
real(wp) function, dimension(1:4, 1:4), public f_create_transform_matrix(param, center)
This procedure creates a transformation matrix.
Definition m_helper.fpp.f90:830
m_helper::f_logical_to_int
elemental integer function, public f_logical_to_int(predicate)
This procedure converts logical to 1 or 0.
Definition m_helper.fpp.f90:982
m_helper::s_int_to_str
elemental subroutine, public s_int_to_str(i, res)
Converts an integer to its trimmed string representation.
Definition m_helper.fpp.f90:749
m_helper::s_swap
elemental subroutine, public s_swap(lhs, rhs)
This procedure swaps two real numbers.
Definition m_helper.fpp.f90:816