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