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